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?

%先对方程进行处理，为使用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%先对方程进行处理，为使用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 %先对方程进行处理，为使用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’);

Cart

Cart