如何使用 ode45 解算器为 ode 的解设置动画?

How to animate solution of ode with ode45 solver?

我正在尝试访问 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

您可以使用 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

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