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

[Torsten]

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]

Thanks for your reply, I will double check the parameters in it.

Answer 1

[Star Strider]

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');
set(gca, 'YScale','log')                                            % <— ADDED

figure;
plot(t, f(:,3));
title('varphi.time');
xlabel('time');
ylabel('varphi');
set(gca, 'YScale','log')                                            % <— ADDED

figure;
plot(t, f(:,5));
title('V.time');
xlabel('time');
ylabel('V');
set(gca, 'YScale','log')                                            % <— ADDED

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]

Thank you for your reply, the parameters are indeed a bit unreasonable, I'm trying to change them, in the meantime I'd like to ask apart from the parameter problem is there any error in solving this second order equation using ode45 in this way?

[Star Strider]

As always, my pleasure!

There is no error with respect to using ode45. Indeed, ode45 will probably be able to integrate your system if given enough time (perhaps a week or so, considering that it has not managed to get to 0.09 in several hours). A stiff solver is designed to deal with these problems more efficiently.

[fayu]

Thanks, I will try it with ode15s.

[Star Strider]

As always, my pleasure!

The ode23s function is also an option.