Несколько выходов решателя Matlab ODE

У меня есть следующий код Matlab ODE:

[t,y,~,~,ie] = ode23tb(@(t,y) RHSODE(t,y),[0,t_end], [i0;v0],options);

Я хочу, чтобы решатель ODE также мог дать мне результат z, который является функцией y и dy / dt, такой, что z = f (y, dy / dt).

Кто-нибудь знает, как добавить такой z в вывод решателя?


matlab ode numerical-integration
person HeyMan    schedule 28.01.2014    source источник

Ответы (1)


arrow_upward
4
arrow_downward

Есть два способа сделать это. Наиболее распространенный и, как правило, самый быстрый способ - воспользоваться функцией интеграции (RHSODE в вашем случае) и оценить функцию f после выполнения интеграции. Вы не указали много деталей в своем коде, но он может выглядеть примерно так:

ydot = RHSODE(t,y);
z = f(y,ydot);

где t и y - выходы из ode23tb. Для этого требуется, чтобы и RHSODE, и f были векторизованы (или вы можете заключить вышеуказанное в цикл for).

Другой способ требует, чтобы вы создали дополнительное уравнение (или уравнения, если z - вектор) внутри вашей функции интегрирования, RHSODE. Обычно ode23tb интегрирует что-либо в этой функции, поэтому f нужно умножить на коэффициент t, чтобы отменить это. Опять же, ваш код может выглядеть примерно так:

function ydot = RHSODE(t,y)
ydot0 = ... % Your original ODE(s)
z = f(y,ydot);
ydot = [ydot0;z*t]; % Make column vector
person horchler    schedule 28.01.2014
comment
@ user2270626: Было ли это полезно? Это ответило на ваш вопрос? Если да, примите ответ. Спасибо. - person horchler; 07.02.2014
comment
Но y в ydot = RHSODE(t,y); - это целый список y, решенных с помощью решателя од. Размер не совпадает с исходным определением y, скажем 2 by 1 column vector в RHSODE? Как справиться с этой проблемой? @horchler - person Ka-Wa Yip; 27.06.2016
comment
Если y в ydot = RHSODE(t,y); является скаляром, метод работает. Но если y в ydot = RHSODE(t,y); - это, скажем, 2 by 1 column vector. Что я должен делать? - person Ka-Wa Yip; 27.06.2016
comment
@kww: это просто означает, что ваша функция интеграции не векторизована (это не обязательно, если вы просто используете ode45 без необходимости в этом решении). Как я заявляю в своем ответе, вы можете либо векторизовать его, либо использовать цикл for для итерации ваших выходных данных. Возможно, вам также потребуется применить транспонирование в зависимости от вашего конкретного кода. - person horchler; 27.06.2016
comment
@horcher Но тут есть проблема. Что делать, если интеграция не удалась и решатель вызовет функцию для предыдущего временного шага? Нужно ли мне реализовать функцию вывода outputFcn, которая вычисляет z на начальном и всех успешных шагах одесольвера? - person Ka-Wa Yip; 30.06.2016
comment
@horcher В моем случае z = f(y,ydot,t). z также зависит от времени. В этом случае - единственный метод outputFcn, который может обнаруживать флажки и отказываться от шагов? Могу ли я использовать метод, указанный в вашем ответе, если z = f(y,ydot,t)? Спасибо за совет. - person Ka-Wa Yip; 30.06.2016
comment
@kww: Вам не нужно беспокоиться о неудачных шагах с помощью любого из этих методов. ode45 и другие решатели возвращают результат только успешных шагов. Оба эти метода должны быть более эффективными, чем использование 'OutputFcn'. - person horchler; 30.06.2016
comment
@horcher спасибо. Для второго метода в вашем ответе здесь. Вроде не получается. Скажем, в простейшем случае y как скаляр, поэтому ydot должен быть скаляром, но во втором методе выход ydot представляет собой вектор-столбец 2 на 1. Я получаю от этого ошибку в строке оды. Кроме того, должна ли запись ydot быть z/t? - person Ka-Wa Yip; 01.07.2016