Есть ли способ вывести промежуточные вычисления решателя оды в MATLAB?

Question

У меня есть система жестких ODE, которые я решаю с помощью ode15s. Производная по времени имеет тип:

yt = A*y - u.*y;

Где yt, y и u - векторы столбцов, а A - разреженная матрица. Вычисление u должно выполняться каждый раз, когда вычисляется, и оно требует значительных вычислительных ресурсов (u нелинейно, и это является причиной жесткости). Я не хочу предварительно вычислять u, потому что я хотел бы сохранить свободу решателя, чтобы выбрать любое желаемое время позже. Однако меня интересует не решение y (t), а решение u (t). * Y (t), которое уже рассчитано для вычисления yt. Если я получу решение в виде структуры sol, которую затем смогу оценить с помощью deval, я легко вычислю y (t), но u (t) придется пересчитать. Мне было интересно, есть ли какой-нибудь способ сохранить результат u.*y при использовании решателя, чтобы я мог вызывать его так же легко, как и deval(y,t). Самая близкая альтернатива, которую я нашел, состоит в том, чтобы полагаться на deval, также предоставляя производную по времени, поэтому сделайте:

[y,yt] = deval(sol,t); %t is a row vector
uy     = A*y - yt;     %uy = u.*y

Что не очень неэффективно, но тратит некоторые вычисления в операциях, которые уже были вычислены.

Есть идеи, как это сделать?

Answer 1

Поскольку я исследовал эту точно такую ​​же проблему, у меня было несколько месяцев go, это невозможно из-за встроенных функций. Есть некоторые обходные пути для сохранения промежуточных вычислений (например, использование постоянных или глобальных переменных), но они создают больше проблем.

На форуме MATLAB много дискуссий. Этот ответ может быть особенно полезен. Другие комментарии можно найти здесь или здесь .

Проблема сохранения промежуточных вычислений с помощью встроенных в MATLAB функций ode заключается в том, что они являются адаптивными размер шага и интегрируемая функция вызывается несколько раз за шаг. Если вы передадите t как вектор с более чем 2 элементами , MATLAB вернет решение в указанные временные шаги. Это, однако, не означает, что он вычисляет решение по этим временным шагам - он просто интерполирует их встроенное решение.

Лучший вариант - это то, что вы упомянули, для расчета желаемого результата после решения оды. , что означает пересчет u. Если вы правильно векторизуете свои функции, этот шаг не должен быть проблемным временем c.

Answer 2

Вы можете "взломать" решатели ODE и включить управляющие переменные u в вектор состояния. Тогда ваша система ODE становится (тривиальной, явной) системой DAE

y' = f(t,y,u)
0  = u - g(t,x)

Единственный «трюк» теперь заключается в том, что вам также необходимо передать матрицу масс в опциях, которая является диагональной матрицей с записями 1 для дифференциальных уравнений и 0 для алгебраических уравнений c. Проверьте документацию, если используемый вами решатель поддерживает DAE, это должны обеспечивать все неявные методы.

Может быть некоторая потеря эффективности, поскольку новые уравнения также будут влиять на управление размером шага, но я не ожидаю большой разницы .

В результате значения элемента управления будут сохранены в последних компонентах вектора состояния, поэтому их извлечение является простым разрезанием этого вектора.