在Matlab中使用ode45

Question

我试图模拟一个由ODEs系统控制的物理过程的时间行为。当我将输入脉冲的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在I = 0中非常简单，所以像ode45这样的自适应求解程序将采取非常大的步骤。因此，很有可能它会穿过你注射冲动的区域，而永远看不到它。如果您对为什么您的问题中的代码错过脉搏感到困惑，请参见我在这里的回答，即使您已经指定了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

另见我对类似问题的回答。