清空筛选条件
清空筛选条件

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?

3 次查看（过去 30 天）
fayu
fayu 2024-5-23
评论： Star Strider 2024-5-24
%先对方程进行处理，为使用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 个评论
Torsten
Torsten 2024-5-23
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.

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Star Strider
Star Strider 2024-5-23
Ran in:
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
.
  4 个评论
显示 2更早的评论

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by