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?
1 Ansicht (letzte 30 Tage)
Ältere Kommentare anzeigen
fayu
am 23 Mai 2024
Kommentiert: Star Strider
am 24 Mai 2024
%先对方程进行处理,为使用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
2 Kommentare
Torsten
am 23 Mai 2024
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.
Akzeptierte Antwort
Star Strider
am 23 Mai 2024
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);
toc
%绘图
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
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
.
4 Kommentare
Weitere Antworten (0)
Siehe auch
Kategorien
Mehr zu Ordinary Differential Equations finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!