2017-03-27 31 views
0

我是相当新的Matlab和这个ODE求解器,下面是我的代码: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]') 

pend.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秒。

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

谢谢。

回答

0

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

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

使用"events"停止并重新启动整合过程,恰好是那些不平滑的模型更改。