如何使用 ode45 求解器对 ode 的解进行动画处理？

Question

我正在尝试访问 ode45 提供的当前解决方案。我见过的所有例子，他们解决了颂歌，然后存储解决方案以供进一步操作。虽然这是一个替代解决方案，但我需要在每次迭代时访问该解决方案，以便我可以为动画目的绘制它们。这是简单摆的代码。

clear 
clc

theta0 = [pi/2; 0]; 
tspan = [0 10];
[t, theta] = ode45(@odeFun, tspan, theta0);
plot(t, theta(:,1));

function dydx = odeFun(t, theta)
    g = 9.8;
    l = 1;
    dydx = zeros(2, 1);
    dydx(1) = theta(2);
    dydx(2) = -g/l*theta(1);
end

我的目标是实现以下目标

clear 
clc

theta0 = [pi/2; 0]; 
t=0;
while t < 10
    [t, theta] = ode45(@odeFun, ..., ...);
    plot(t,theta)
    drawnow
    t = t + 0.1;
end

[t, theta] = ode45(@odeFun, tspan, theta0);
plot(t, theta(:,1));

function dydx = odeFun(t, theta)
    g = 9.8;
    l = 1;
    dydx = zeros(2, 1);
    dydx(1) = theta(2);
    dydx(2) = -g/l*theta(1);
end

Answer 1

您可以使用ode output function 选项。 ODE 求解器在每个成功的时间步后调用输出函数。下面是一个使用输出函数的非常简单的例子：

clear 
clc
close all
theta0 = [pi/2; 0]; 
tspan = [0 10];
options = odeset('OutputFcn',@plotfun,'MaxStep',0.0001);
axes; hold on
[t, theta] = ode45(@odeFun, tspan, theta0, options);

function status = plotfun(t,y,flag)
    plot(t, y(1,:),'k-',t,y(2,:),'m-');
    drawnow
    status = 0;
end

function dydx = odeFun(t, theta)
    g = 9.8;
    l = 1;
    dydx = zeros(2, 1);
    dydx(1) = theta(2);
    dydx(2) = -g/l*theta(1);
end

最大积分步长设置为 0.0001，以便不会立即绘制解图。

或者，您可以通过将时间间隔分成几部分并使用前一部分的终点作为下一部分的起始点来计算解：

clear 
clc
close all
axes; hold on
theta0 = [pi/2; 0]; 
tstart=0; dt= 0.1;
options= odeset('MaxStep',0.0001);
while tstart < 10
    [t, theta] = ode45(@odeFun,[tstart tstart+dt],theta0,options);
    plot(t,theta(:,1),'k-',t,theta(:,2),'m-');
    drawnow
    tstart = tstart + dt;
    theta0= theta(end,:);
end

function dydx = odeFun(t, theta)
    g = 9.8;
    l = 1;
    dydx = zeros(2, 1);
    dydx(1) = theta(2);
    dydx(2) = -g/l*theta(1);
end

当然，您可以将这两种方法结合起来。