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

Question

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

main.m

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]')

待定.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 随着极小的步长而增加，随着时间的推移继续减少，例如20秒0.1到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 函数。 4 阶求解器期望 ODE 函数至少 4 次连续可微。

您的 ODE 函数有跳跃和扭结。求解器将这些感知为非常大且混乱地振荡的较高导数值。为了补偿和恢复收敛顺序，减小了步长，进一步增加了计算的导数值，因为您的函数并不是真正可微的，直到步长增量变得小于最小有效浮点增量，从而触发您看到的错误。

使用"events" 来停止和重新启动集成过程，恰好是那些不顺利的模型更改。

Answer 2

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

例如； 将时间 (t) 计算为“毫秒”而不是“秒”。这给出（或保留）额外的 3 位小数，以允许最小步长的适当增量，Δt ≥ 1.776357e-15。

注意：如果有另一个物理量从根本上派生自量 t，例如速度（“米每秒”），请确保比例/单位是等价的.在给出的示例中，计算中的速度单位必须以“米/毫秒”为单位。

补充说明； 返回所需的最终值时，如有必要，转换回适当的比例。