I want to solve the equation in the picture and use ode45 to solve it, but after running for 5 hours there is still no result, I would like to ask if the program is written wrong or is there any optimization?

2 views (last 30 days)
fayu on 23 May 2024
Commented: Star Strider on 24 May 2024 at 11:15
%先对方程进行处理，为使用ode45准备
% clc
% clear
% syms r dx2 rm l varphi dvarphi2 dvarphi1 theta1 theta C1 dv1 Lb C2 dvarphi
% syms g Cp dVR omega1 omega2 v V dv2 z %z=(3*pi-8)/(2*pi)避免其中转化为小数
% eqn1=dv2+r*dx2+rm*l*(dvarphi2-varphi*dvarphi1^2)+omega1^2*v+theta1*V+C1*dv1*Lb*z==0;
% eqn2=C2*dvarphi1+dvarphi2+dv2/l+dx2/l+omega2^2*varphi==0;
% eqns=[eqn1 eqn2];
% S=solve(eqns,[dv2,dvarphi2]);
% S.dv2;
% S.dvarphi2;
%变量赋值 其中还没有确定对应的变量的值
clc
clear
m=0.013;
z=(3*pi - 8)/(2*pi);
Lb=0.19;
Ab=0.008*0.0006;
Lp=0.02;
Ap=0.008*0.0002;
l=0.0275;
rho_b=7850;
rho_p=7500;
Eb=109e9;
Ep=106e9;
theta=190e-12;
Cp=12.7e-9;
Ib=1/12*0.008*0.0006^3;
Ip=1/12*0.008*0.0002^3;
R=10000;
Cb=0.0117; %正常摆动时候的阻尼系数
Ch=0.00117; %旋转时候的阻尼系数 论文中没有出现
g=9.81;
dx2=0.2*9.81;
syms y
psi=@(y) 1-cos(pi*y/(2*Lb));
M1=double(rho_b*Ab*int((psi(y))^2,y,[0,Lb])+rho_p*Ap*int((psi(y))^2,y,[0,Lp]));
M2=double(rho_b*Ab*int((psi(y)),y,[0,Lb])+rho_p*Ap*int((psi(y)),y,[0,Lp]));
K=double(Eb*Ib*int(diff(psi(y),y,2)^2,y,[0,Lb])+Ep*Ip*int(diff(psi(y),y,2)^2,y,[0,Lp]));
C1=Cb/(M1+m);
C2=Ch/(m*l^2);
theta1=theta/(M1+m);
rm=m/(M2+m);
r=(M2+m)/(M1+m);
omega1=sqrt(K/M1);
omega2=sqrt(g/l);
%传递参数准备
params = struct('z',z,'Lb',Lb,'dx2',dx2,'l',l,'theta',theta,'Cp', Cp, ...
'R',R,'C1',C1,'C2',C2,'theta1',theta1,'rm',rm,'r',r,'omega1',omega1,'omega2',omega2);
%采用ode求解
y0=[0 0 0 0 0];
tspan=0:0.001:2;
[t,f]=ode45(@(t,f) odesys(t,f,params),tspan,y0);
%绘图
figure;
plot(t,f(:,1));
title('v.time');
xlabel('time');
ylabel('v');
figure;
plot(t, f(:,3));
title('varphi.time');
xlabel('time');
ylabel('varphi');
figure;
plot(t, f(:,5));
title('V.time');
xlabel('time');
ylabel('V');
function df=odesys(t,f,params)
% 从参数结构体中提取参数
z=params.z;
Lb=params.Lb;
dx2=params.dx2;
l=params.l;
Cp=params.Cp;
R=params.R;
C1=params.C1;
C2=params.C2;
theta=params.theta;
theta1=params.theta1;
rm=params.rm;
r=params.r;
omega1=params.omega1;
omega2=params.omega2;
% 状态变量
v= f(1);
dv1=f(2);
varphi=f(3);
dvarphi1=f(4);
V=f(5);
%微分方程
df=zeros(5,1);
df(1)=dv1; % df(1) 代表 dv/dt
df(3)=dvarphi1; % df(3) 代表 dvarphi/dt
df(2)=(-l*rm*df(3)^2-C2*l*rm*df(3)+v*omega1^2-...
l*rm*varphi*omega2^2+V*theta1+dx2*r-dx2*rm+C1*Lb*df(1)*z)/(rm-1); % df(2) 代表 dv1/dt
df(4)=(l*rm*varphi*df(3)^2+C2*l*df(3)-v*omega1^2+l*varphi*omega2^2+...
dx2-V*theta1-dx2*r-C1*Lb*df(1)*z)/(l*(rm - 1)); % df(4) 代表 dvarphi1/dt
df(5)=(V/R-theta*df(1))/Cp; % df(5) 代表 dV/dt
end
Torsten on 23 May 2024 at 17:51
Your model contains very small and very big parameters. To be honest, I would have been surprised if it worked without problems.
Check again your equations and the values of the model parameters (units, physical sensefulness). Maybe switching from ode45 to ode15s can help.
fayu on 24 May 2024 at 2:40
Thanks for your reply, I will double check the parameters in it.

Star Strider on 23 May 2024 at 17:47
The system is ‘stiff’ probably because of the large variation in the magnitudes of the parameters.
Using a stiff solver (ode15s here, there are others), it solves in seconds. (There are still problems, because the function encounters a singularity at about 0.09 secoinds, and the integration then stops.)
Try something like this —
m=0.013;
z=(3*pi - 8)/(2*pi);
Lb=0.19;
Ab=0.008*0.0006;
Lp=0.02;
Ap=0.008*0.0002;
l=0.0275;
rho_b=7850;
rho_p=7500;
Eb=109e9;
Ep=106e9;
theta=190e-12;
Cp=12.7e-9;
Ib=1/12*0.008*0.0006^3;
Ip=1/12*0.008*0.0002^3;
R=10000;
Cb=0.0117; %正常摆动时候的阻尼系数
Ch=0.00117; %旋转时候的阻尼系数 论文中没有出现
g=9.81;
dx2=0.2*9.81;
syms y
psi=@(y) 1-cos(pi*y/(2*Lb));
M1=double(rho_b*Ab*int((psi(y))^2,y,[0,Lb])+rho_p*Ap*int((psi(y))^2,y,[0,Lp]));
M2=double(rho_b*Ab*int((psi(y)),y,[0,Lb])+rho_p*Ap*int((psi(y)),y,[0,Lp]));
K=double(Eb*Ib*int(diff(psi(y),y,2)^2,y,[0,Lb])+Ep*Ip*int(diff(psi(y),y,2)^2,y,[0,Lp]));
C1=Cb/(M1+m);
C2=Ch/(m*l^2);
theta1=theta/(M1+m);
rm=m/(M2+m);
r=(M2+m)/(M1+m);
omega1=sqrt(K/M1);
omega2=sqrt(g/l);
%传递参数准备
params = struct('z',z,'Lb',Lb,'dx2',dx2,'l',l,'theta',theta,'Cp', Cp, ...
'R',R,'C1',C1,'C2',C2,'theta1',theta1,'rm',rm,'r',r,'omega1',omega1,'omega2',omega2);
%采用ode求解
y0=[0 0 0 0 0];
tspan=0:0.001:2;
tic
[t,f]=ode15s(@(t,f) odesys(t,f,params),tspan,y0);
Warning: Failure at t=9.147712e-02. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (2.220446e-16) at time t.
toc
Elapsed time is 0.356515 seconds.
%绘图
figure;
plot(t,f(:,1));
title('v.time');
xlabel('time');
ylabel('v');
figure;
plot(t, f(:,3));
title('varphi.time');
xlabel('time');
ylabel('varphi');
figure;
plot(t, f(:,5));
title('V.time');
xlabel('time');
ylabel('V');
Warning: Negative data ignored
Warning: Negative data ignored
function df=odesys(t,f,params)
% 从参数结构体中提取参数
z=params.z;
Lb=params.Lb;
dx2=params.dx2;
l=params.l;
Cp=params.Cp;
R=params.R;
C1=params.C1;
C2=params.C2;
theta=params.theta;
theta1=params.theta1;
rm=params.rm;
r=params.r;
omega1=params.omega1;
omega2=params.omega2;
% 状态变量
v= f(1);
dv1=f(2);
varphi=f(3);
dvarphi1=f(4);
V=f(5);
%微分方程
df=zeros(5,1);
df(1)=dv1; % df(1) 代表 dv/dt
df(3)=dvarphi1; % df(3) 代表 dvarphi/dt
df(2)=(-l*rm*df(3)^2-C2*l*rm*df(3)+v*omega1^2-...
l*rm*varphi*omega2^2+V*theta1+dx2*r-dx2*rm+C1*Lb*df(1)*z)/(rm-1); % df(2) 代表 dv1/dt
df(4)=(l*rm*varphi*df(3)^2+C2*l*df(3)-v*omega1^2+l*varphi*omega2^2+...
dx2-V*theta1-dx2*r-C1*Lb*df(1)*z)/(l*(rm - 1)); % df(4) 代表 dvarphi1/dt
df(5)=(V/R-theta*df(1))/Cp; % df(5) 代表 dV/dt
end
.
fayu on 24 May 2024 at 6:40
Thanks, I will try it with ode15s.
Star Strider on 24 May 2024 at 11:15
As always, my pleasure!
The ode23s function is also an option.