Matlab ODE 求解器由于步长小而永远卡住了

Question

我对 Matlab 和这个 ODE 求解器相当陌生，下面是我的代码：

主米

format short;
tspan=[0 5];
y0=[0.30;-0.30;0;-0;0];

[t,y]=ode23s(@(t,y) pend(t,y),tspan,y0);
figure(1)
%subplot(2,1,1); 
plot(t,y(:,1),t,y(:,2),'k--')
set(gcf,'Position',[100,500,450,180]);
xlabel('time [s]');
legend('q_1','q_2')
ylabel('leg angle [rad]');
figure(2)
%subplot(2,1,2); 
plot(t,y(:,5))
set(gcf,'Position',[100,500,450,180]);
xlabel('time [s]')
ylabel('locomotion [m]')

彭德·

%the following function contains the right hand side of the
%differential equation of the form
%M(t,y)*y'=F(t,y)
%i.e. it contains F(t,y).it is also stored in a separate filenamed, pend.m.
function yp= pend(t,y)
m = 5;                  %leg masses [kg]        suggested: 5
%'shin' length [m]      suggested: 0.5
b = 0.5;                %'thigh' length [m]     suggested: 0.5
L = 2*b; 
q=y(1:2);
dq=y(3:4);
Ox=y(5);
rho=0;
k0=50; %Nm/rad
v = 0;
Hsw= L*cos(q(1));  % Height of leg1
Hst= L*cos(q(2));   % Height of leg2
H1 = L - Hsw; 
H2 = L - Hst; 
if dq(1)<0
Fid1=1;
else Fid1=0;
end
if dq(2)<0
Fid2=1;
else Fid2=0;
end    
F1 =  -15000*min(H1-0.03*L,0)*Fid1; %N
F2 =  -15000*min(H2-0.03*L,0)*Fid2;
Fc1 = F1*L*sin(q(1));
Fc2 = F2*L*sin(q(2));
Fc=[Fc1;Fc2];
M=[m*b^2 0;0 m*b^2];
Ko=k0*[1 -1; -1 1]+m*9.8*b*[1 0; 0 1];
D=M;
ddq=inv(M)*(-rho*D*dq-Ko*q+Fc);
dOx=0;
if Fid1==1
dOx=-L*dq(1)*cos(q(1))*sign(F1);
end
if Fid2==1
dOx=-L*dq(2)*cos(q(2))*sign(F2);
end
yp=[dq;ddq;dOx];

我在这里面临的问题是时间跨度 t 以极小的步骤增加，这些步长随着时间的推移而继续减少，例如 0.1 秒内 0.8 到 0.8 秒，5 分钟内 0.8 到 0.9 等等，这意味着它永远不会达到时间限制，因此卡在循环中。

我尝试过不同的求解器，如 ode45，也尝试给出不同的 RelTol 和 AbsTol 值来控制步长，但失败了。它确实对前几步有所作为，但随后又变慢了。

当我使用求解器 ode15s 时，它会发出警告

"t=9.246943e-01 时失败。 无法满足集成公差 不会将步长减小到允许的最小值以下 (1.776357e-15) 在时间 t.">

并且只绘制图形直到 0.92 秒。

欢迎任何建议或帮助解决此问题。

谢谢。

Answer 1

具有自适应步长的 ODE 求解器需要平滑的 ODE 函数。四阶求解器期望 ODE 函数至少是连续可微分的 4 倍。

您的 ODE 函数有跳跃和扭结。求解器将这些感知为非常大且混沌振荡的较高导数值。为了补偿和恢复收敛阶数，减小步长，进一步增加计算的导数值，因为您的函数实际上不可微分，直到步长增量小于最小有效浮点增量，触发您看到的误差。

使用"事件"在那些不平滑的模型更改时停止并重新启动集成过程。

Answer 2

您可以处理不同规模的数量。

例如; 将时间 (t) 计算为"毫秒"而不是"秒"。这提供了(或保留)额外的 3 位小数，以允许最小步长Δt≥ 1.776357e-15 的适当增量。

注意：如果存在另一个基本上从量t得出的物理量，如速度("米每秒")，请确保刻度/单位是等效的。在给出的例子中，计算中的速度单位必须以"米/毫秒"为单位。

额外说明; 如有必要，在返回所需的最终值时转换回适当的比例。