在 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 和不太准确，甚至完全错误的结果。在这种情况下，由于您的 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。