matlab可以将微分方程的解显示为相图形式吗

Question

如下图所示,解微分方程组的程序见图二和图三,目标是将微分方程组求解出的时域响应图转换为右图所示的相图,请问有修改的建议吗?非常感谢

clear;clc;close all;
%初始参数，该参数变化，myfun程序要对应调整
syms t;
wm_func = @(t)((10+2.67*t*3.6)*9550/(3600*0.4064/(0.93*14)))*pi*2/60;
zi=4;%内转子极对数
zo=14;%外转子极对数k
z1=32;%小齿轮齿数
z2=99;%大齿轮齿数
ms=4.62e-3;
I=z2/z1;%机械齿轮速比
G=(zo+zi)/zi;%磁齿轮速比
%一、计算过程
tmax=4.16;
tspan = 0:1e-4:tmax;
% tmax=900*2*pi/wm;% 计算时长，按输入轴的周期数计算
bc=1e-4;
y0=[0;0;0;wm_func(0);wm_func(0)/G;wm_func(0)/(G*I)];
% y0=[0;0;0;wm;wm/G;wm/(G*I)];
[t,y] = ode23(@(t,y)myfun3(t,y,wm_func), tspan, y0);
figure;plot(t,y(:,1)/G-y(:,2));
signal = y(:,1)/G - y(:,2);
set(gca,'FontSize',16);
title('时域响应图','FontSize',16);
xlabel('时间 (s)','FontSize',16);
ylabel('扭转振幅(rad)','FontSize',16);%磁齿轮的相对弹性角度
figure;plot(t,y(:,2)/I-y(:,3));
signal1=y(:,2)/I-y(:,3);
set(gca,'FontSize',16);
title('时域响应图','FontSize',16);
xlabel('时间 (s)','FontSize',16);
ylabel('扭转振幅(rad)','FontSize',16);%机械齿轮的相对弹性转角;
figure;plot(t,wm_func(t));

function dydt=myfun3(t,y,wm_func)
wm=wm_func(t);
% wm = 3891*2*pi/60;
zi=4;%内转子极对数
zo=14;%外转子极对数
z1=32;%小齿轮齿数
z2=99;%大齿轮齿数
arf0=20/180*pi;%压力角
ms=4.62e-3;
r2=ms*z1/2;%小齿轮分度圆
r3=ms*z2/2;%大齿轮分度圆
rb2=r2*cos(arf0);%小齿轮基圆
rb3=r3*cos(arf0);%大齿轮基圆
I=z2/z1;%机械齿轮速比
G=(zo+zi)/zi;%磁齿轮速比
Tm =(1460*9.8*0.012+1.04*1460*2.67+0.32*2.25*(10+2.67*t*3.6)^2/21.15)*0.4064/(0.93*14);%起步加速
TL=Tm*I*G;%变速器负载
%3 转动惯量
IM=0.081197;%电机转子转动惯量kgm2
I0=0.181197;%磁齿轮低速级转动惯量kgm2
I1=0.013663;%磁齿轮高速级转动惯量kgm2
I2=0.0230596;%机械小齿轮转动惯量
I3=1.9825692;%机械大齿轮转动惯量
IL=21;%整车等效转动惯量
%4 刚度与阻尼
%齿轮的阻尼和刚度
w1=wm*G;%机械小齿轮输入转速
kp=10.51e8;
kb=0.25e8;%两个齿轮的刚度曲线为正弦曲线，
km=kp+kb*sin(w1*t*z1);
sg=0.16;%齿轮阻尼比
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=1.7*Tm*G;
Tc1=0.1*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)/(I0+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

该回答结合ChatGPT4o及杨同学*共同作答, 如有帮助，还请采纳。
要在MATLAB中将微分方程的解显示为相图形式，可以使用plot函数绘制相位图。相图是一种表示一个系统状态变量与其导数之间关系的图形。下面是如何修改你的代码，以生成类似于右图所示的相图：

首先，需要对原始代码稍作修改，以便从时域响应图生成相图。你已经有了解微分方程组的时域响应，现在我们需要绘制DTE相图。

以下是修改后的代码示例：

clear;clc;close all;

% 初始参数
syms t;
wm_func = @(t)((10+2.67*t*3.6)*9550/(3600*0.4064/(0.93*14)))*pi*2/60;
zi=4; % 内转子极对数
zo=14; % 外转子极对数
z1=32; % 小齿轮齿数
z2=99; % 大齿轮齿数
ms=4.62e-3;
I=z2/z1; % 机械齿轮速比
G=(zo+zi)/zi; % 磁齿轮速比

% 计算过程
tmax=4.16;
tspan = 0:1e-4:tmax;
y0=[0;0;0;wm_func(0);wm_func(0)/G;wm_func(0)/(G*I)];
[t,y] = ode23(@(t,y)myfun3(t,y,wm_func), tspan, y0);

% 绘制时域响应图
figure;
plot(t,y(:,1)/G-y(:,2));
title('时域响应图');
xlabel('时间 (s)');
ylabel('扭转振幅 (rad)');

% 计算相图
DTE = y(:,1)/G - y(:,2); % DTE (电磁弹性角度)
dDTE_dt = gradient(DTE, t); % DTE的一阶导数

% 绘制相图
figure;
plot(DTE, dDTE_dt);
title('DTE相图');
xlabel('DTE (rad)');
ylabel('d(DTE)/dt (rad/s)');

function dydt = myfun3(t, y, wm_func)
    wm = wm_func(t);
    zi = 4; % 内转子极对数
    zo = 14; % 外转子极对数
    z1 = 32; % 小齿轮齿数
    z2 = 99; % 大齿轮齿数
    arf0 = 20/180*pi; % 压力角
    ms = 4.62e-3;
    r2 = ms*z1/2; % 小齿轮分度圆
    r3 = ms*z2/2; % 大齿轮分度圆
    rb2 = r2*cos(arf0); % 小齿轮基圆
    rb3 = r3*cos(arf0); % 大齿轮基圆
    I = z2/z1; % 机械齿轮速比
    G = (zo+zi)/zi; % 磁齿轮速比
    Tm = (1460*9.8*0.012+1.04*1460*2.67+0.32*2.25*(10+2.67*t*3.6)^2/21.15)*0.4064/(0.93*14); % 起步加速
    TL = Tm*I*G; % 变速器负载
    IM = 0.081197; % 电机转子转动惯量kgm2
    I0 = 0.181197; % 磁齿轮低速级转动惯量kgm2
    I1 = 0.013663; % 磁齿轮高速级转动惯量kgm2
    I2 = 0.0230596; % 机械小齿轮转动惯量
    I3 = 1.9825692; % 机械大齿轮转动惯量
    IL = 21; % 整车等效转动惯量
    w1 = wm*G; % 机械小齿轮输入转速
    kp = 10.51e8;
    kb = 0.25e8; % 两个齿轮的刚度曲线为正弦曲线
    km = kp + kb*sin(w1*t*z1);
    sg = 0.16; % 齿轮阻尼比
    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 = 1.7*Tm*G;
    Tc1 = 0.1*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)/(I0 + 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

修改后的代码包含以下几个主要部分：

1. 计算DTE相图所需的信号（DTE）和其导数（dDTE/dt）。
2. 使用MATLAB的gradient函数计算DTE的一阶导数。
3. 使用plot函数绘制DTE相图。

通过这段代码，你可以得到DTE的相图，类似于你提供的右图。希望这可以满足你的需求。如果有任何其他问题或进一步的修改建议，请告诉我。