在Matlab中使用ode45

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的参数。这导致了非常不准确的结果，甚至完全错误的结果。在这种情况下，因为您的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。