Решатели Matlab ода: изменение состояния и заданного времени

Question

Я решаю набор ODE (dy/dt) при t = 0, все начальные условия t=0 y_0=(0,0,0). Могу ли я добавить некоторое число к значениям y в разное время (например, при t=10 к этому числу следует добавить y1; при t=20 к этому числу следует добавить y2 и т. Д.) И решить уравнения?

Answer 1

Вставка больших разрывов в ваш ODE способом, который вы предлагаете (и способ, показанный @macduff), может привести к снижению точности и увеличению времени вычислений (особенно с ode45 - ode15s может быть лучшим вариантом или, по крайней мере, убедитесь, что ваши абсолютные и относительные допуски подходят). Вы эффективно создали очень жесткую систему. Если вы хотите добавить некоторое число в ODE, начиная с определенного времени, имейте в виду, что решатель оценивает эти уравнения только в определенные моменты времени по своему выбору. (Не вводите в заблуждение тот факт, что вы можете получить фиксированный выходной размер, указав tspan как более двух элементов - все решатели Matlab являются решателями с переменным размером шага и выбирают свои истинные шаги на основе критериев ошибки.)

Лучшим вариантом является частичная интеграция системы и добавление итоговых результатов каждого прогона вместе:

% t = 0 to t = 10, pass parameter a = 0 to add to ODEs
a = 0;
tspan = [0 10];
[T,Y] = ode45(@(t,y)myfun(t,y,a),tspan,y_0);

% t = 10 to t = 20, pass parameter a = 10 to add to ODEs
a = 10;
[t,y] = ode45(@(t,y)myfun(t,y,a),tspan+T(end),Y(end,:));
T = [T;t(2:end)];
Y = [Y;y(2:end,:)];

% t = 20 to t = 30, pass parameter a = 20 to add to ODEs
a = 20;
[t,y] = ode45(@(t,y)myfun(t,y,a),tspan+T(end),Y(end,:));
T = [T;t(2:end)];
Y = [Y;y(2:end,:)];

Редактор Matlab может пожаловаться на массив T а также Y не распределяются заранее и / или не растут, но в этом случае это нормально, так как они растут большими кусками всего несколько раз. В качестве альтернативы, если вам нужны фиксированные выходные размеры шагов, вы можете сделать это:

dt = 0.01;
T = 0:dt:30;
Y = zeros(length(T),length(y_0));

% t = 0 to t = 10, pass parameter a = 0 to add to ODEs
a = 0;
[~,Y(1:10/dt+1,:)] = ode45(@(t,y)myfun(t,y,a),T(1:10/dt+1),y_0);

% t = 10 to t = 20, pass parameter a = 10 to add to ODEs
a = 10;
[~,Y(10/dt+1:20/dt+1,:)] = ode45(@(t,y)myfun(t,y,a),T(10/dt+1:20/dt+1),Y(10/dt+1,:));

% t = 20 to t = 30, pass parameter a = 20 to add to ODEs
a = 20;
[~,Y(20/dt+1:end,:)] = ode45(@(t,y)myfun(t,y,a),T(20/dt+1:end),Y(20/dt+1,:));

Можно легко преобразовать оба вышеупомянутых блока кода в более компактный for петли при желании.

В обоих случаях ваша функция ODE myfun включает в себя параметр a сюда:

function ydot = myfun(t,y,a)
    y(1) = ... % add a however you like
    ...

Answer 2

Хорошо, как говорит Саймон МакКензи, нам действительно нужно больше информации по вашей неотложной проблеме, но я думаю, что смогу помочь. Из того, что вы нам дали, я предполагаю, что у вас есть функция myfun что вы переходите к чему-то вроде ode45

y_0 = [0,0,0];
% here Tfinal is the time in seconds that you want to simulate to
% and specify the tspan so that you will get a solution at each whole
% number of seconds, I believe
[t,y] = ode45(@myfun,[0:0.1:Tfinal],y_0);

Где-то вы определили свою функцию, здесь я назову ее myfun

function dy = myfun(t,y)
  % here let's check to see if it's time to add your offsets
  % still, you may want to put a little fudge factor for the time, but
  % you'll have to experiment, I'll set it up though
  EPS = 1e-12;
  if( t < 10 + EPS || t > 10 - EPS )
    y(1) = y(1) + 10;
  .
  .
  .
  .
  % this only makes sense if you then use y(1) in the compuatation

В противном случае просто добавьте смещение к возвращенному вектору решения, т.е.

idx10 = find( t == 10 ); % should be there since it's in the tspan
y(idx10:end) = y(idx10:end) + 10; % I guess you add it from that point on?