在 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 的参数。这会导致非常 stiff system and less accurate, or even completely wrong, results. In this case, because the your ODE is so simple when I = 0, an adaptive solver like ode45 will take very large steps. Thus, there's a high probability that it will step right over the region where you inject the impulse and never see it. See 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。

在 Matlab 中使用 ode45

Using ode45 in Matlab

matlab

numerical-methods

numerical-integration