怎么求微分方程组的解析解

Question

主程序见图一,函数程序见图二,已经求出图像解,目标是求出此三元二阶微分方程组y与t的函数关系式,万分感谢

clear;clc;close all;
%初始参数，该参数变化，myfun程序要对应调整
wm=3890*2*pi/60;%电机输入转速rad/s
zi=4;%内转子极对数
zo=14;%外转子极对数
z1=33;%小齿轮齿数
z2=102;%大齿轮齿数
ms=4.5e-3;
r2=ms*z1/2;%小齿轮分度圆
r3=ms*z2/2;%大齿轮分度圆
I=z2/z1;%机械齿轮速比
G=(zo+zi)/zi;%磁齿轮速比
%一、计算过程

tmax=1000*2*pi/wm;% 计算时长，按输入轴的周期数计算
bc=1e-3;
y0=[0;0;0;wm;wm/G;wm/G/I];
[t,f]=odeRK4sys('myfun3',tmax,bc,y0);
figure;plot(t,f(:,1)/G-f(:,2));%磁齿轮的相对弹性角度
figure;plot(t,f(:,2)/I-f(:,3));%机械齿轮的相对弹性转角;
figure;plot(t,f(:,4),t,f(:,5),t,f(:,6))

function dydt=myfun3(t,y)
%1 几何参数
wm=3890*2*pi/60;
zi=4;%内转子极对数
zo=14;%外转子极对数
z1=33;%小齿轮齿数
z2=102;%大齿轮齿数
arf0=20/180*pi;%压力角
ms=4.5e-3;

r2=ms*z1/2;%小齿轮分度圆
r3=ms*z2/2;%大齿轮分度圆
rb2=r2*cos(arf0);%小齿轮基圆
rb3=r3*cos(arf0);%大齿轮基圆
I=z2/z1;%机械齿轮速比
G=(zo+zi)/zi;%磁齿轮速比

%2 力参
Tm=135;%电机输入扭矩
TL=Tm*I*G;%变速器负载
% T0=22;
%3 转动惯量
IM=0.081197;%电机转子转动惯量kgm2
I1=0.181197;%磁齿轮低速级转动惯量kgm2
Ig=0.013265;%磁齿轮高速级转动惯量kgm2
I2=0.0230596;%机械小齿轮转动惯量
I3=1.9825692;%机械大齿轮转动惯量
IL=116.8;%整车等效转动惯量

%4 刚度与阻尼
%齿轮的阻尼和刚度
w1=wm/G;%机械小齿轮输入转速
kp=7.9e8;
kb=3.43e8;%两个齿轮的刚度曲线为正弦曲线，
km=kp+kb*cos(w1*z1*t);
sg=0.2;%齿轮阻尼比
cm=2*sg*sqrt(kp*I2*I3/(I2*rb3^2+I3*rb2^2));

 b=1e-4;
 xx=rb2*y(2)-rb3*y(3);
 dxx=rb2*y(5)-rb3*y(6);

  if xx-b>0
      fx=xx-b;
  elseif xx+b<0
      fx=xx+b;
  else
      fx=0;
  end

F23=km*fx+cm*dxx;
T0=Tm*G*2.2;
Tc1=0.3*T0*sin(wm*t);
Tc=0.1*(zi*y(4)-(zi+zo)*y(5));
T21=T0*sin(zi*y(1)-(zi+zo)*y(2))+Tc+Tc1;

y(7)=(Tm-T21/G)/(Ig+IM);
y(8)=(T21-r2*F23)/(I2+I1);
y(9)=(-TL+r3*F23)/(I3+IL);
dydt=[y(4)
      y(5)
      y(6)
      y(7)
      y(8)
      y(9)];
end

Answer 1

🌈🌈🌈参考通义千问和郭老师的小迷弟雅思莫了-编写提供🌈🌈🌈
您可以参考如下，如果回答的不正确，及时评论区回复我，我会根据你错误描述追加回复，直到您满意为止。

您提供的代码是基于MATLAB的 odeRK4sys 函数求解一个三元二阶微分方程组的数值解，并且通过myfun3函数定义了微分方程。但是，对于解析解，一般复杂的微分方程组很难直接找到显式的解析表达式。您的这个系统包括非线性耦合项和时间相关的正弦项，因此难以得到解析解。

不过，如果您仍然希望尝试寻找解析解，可以考虑以下步骤：

1. 尝试对原微分方程组进行化简或变换，比如看看是否可以通过变量替换、分离变量、拉普拉斯变换等方法简化问题。

2. 对于部分特殊的条件或者近似情况，可能能找到解析解，但这需要深入分析具体的物理背景和数学特性。

由于该微分方程组复杂性较高，直接给出解析解的可能性较小，建议您继续采用数值方法来求解并分析结果。如果确实有特殊情况下的解析解需求，可能需要进一步具体分析该系统的特性和边界条件。