Использование ode45 в Matlab

Question

Я пытаюсь моделировать поведение времени для физического процесса, управляемого системой ODE. Когда я переключаю width входного импульса от 20 до 19, нет никакого истощения состояния y(1), что физически не имеет смысла. Что я делаю не так? Использую ли я ode45 неправильно?

function test 

width = 20; 
center = 100; 

tspan = 0:0.1:center+50*(width/2); 

[t,y] = ode45(@ODEsystem,tspan,[1 0 0 0]); 

plot(t,y(:,1),'k*',t,y(:,2),'k:',t,y(:,3),'k--',t,y(:,4),'k'); 
hold on; 
axis([center-3*(width/2) center+50*(width/2) -0.1 1.1]) 
xlabel('Time') 
ylabel('Relative values') 
legend({'y1','y2','y3','y4'}); 

    function dy = ODEsystem(t,y) 
     k1 = 0.1; 
     k2 = 0.000333; 
     k3 = 0.1; 

     dy = zeros(size(y)); 

     % rectangular pulse 
     I = rectpuls(t-center,width); 

     % ODE system 
     dy(1) = -k1*I*y(1); 
     dy(2) = k1*I*y(1) - k2*y(2); 
     dy(3) = k2*y(2) - k3*I*y(3); 
     dy(4) = k3*I*y(3); 
    end 
end 

Answer 1

Вы меняете параметры своих ODE прерывисто вовремя. Это приводит к очень stiff system и менее точным или даже совершенно неправильным результатам. В этом случае, поскольку ваш ODE настолько прост, когда I = 0, адаптивный решатель, такой как ode45, предпримет очень большие шаги. Таким образом, существует высокая вероятность того, что он будет проходить прямо над регионом, где вы вводите импульс и никогда не увидите его. См. my answer here, если вы смущены тем, почему код в вашем вопросе пропускает пульс, даже если вы указали tspan, чтобы иметь (выходные) шаги всего 0.1.

В целом это плохая идея, чтобы иметь какие-либо разрывы (if заявления, abs, min/max, функции, такие как rectpuls и т.д.) в вашей функции интеграции. Вместо этого вам нужно разбить интеграцию и вычислить ваши результаты кусочно во времени. Вот модифицированная версия кода, который реализует этот:

function test_fixed 

width = 19; 
center = 100; 

t = 0:0.1:center+50*(width/2); 
I = rectpuls(t-center,width); % Removed from ODE function, kept if wanted for plotting 

% Before pulse 
tspan = t(t<=center-width/2); 
y0 = [1 0 0 0]; 
[~,out] = ode45(@(t,y)ODEsystem(t,y,0),tspan,y0); % t pre-calculated, no need to return 
y = out; 

% Pulse 
tspan = t(t>=center-width/2&t<=center+width/2); 
y0 = out(end,:); % Initial conditions same as last stage from previous integration 
[~,out] = ode45(@(t,y)ODEsystem(t,y,1),tspan,y0); 
y = [y;out(2:end,:)]; % Append new data removing identical initial condition 

% After pulse 
tspan = t(t>=center+width/2); 
y0 = out(end,:); 
[~,out] = ode45(@(t,y)ODEsystem(t,y,0),tspan,y0); 
y = [y;out(2:end,:)]; 

plot(t,y(:,1),'k*',t,y(:,2),'k:',t,y(:,3),'k--',t,y(:,4),'k'); 
hold on; 
axis([center-3*(width/2) center+50*(width/2) -0.1 1.1]) 
xlabel('Time') 
ylabel('Relative values') 
legend({'y1','y2','y3','y4'}); 

    function dy = ODEsystem(t,y,I) 
     k1 = 0.1; 
     k2 = 0.000333; 
     k3 = 0.1; 

     dy = zeros(size(y)); 

     % ODE system 
     dy(1) = -k1*I*y(1); 
     dy(2) = k1*I*y(1) - k2*y(2); 
     dy(3) = k2*y(2) - k3*I*y(3); 
     dy(4) = k3*I*y(3); 
    end 
end 

Смотрите также my answer to a similar question.