怎么画出微分方程的相图

wwzqbx

比说说如何画出下面这个微分方程的相图。dx/dt=x+y;dy/dt=4x-2y。求源码！最好带注释的那种，谢谢各位了！！！

博大广阔

运用数值计算方（比如R—K）法求解出x和Y，在画x与Y的图

伤痕累累

那怎么画x，y分别与t的时域图呢？我看了，LZ的微分方程还是耦合的。求高人指点。

博大广阔

%二阶常微分方程组的求解——画相位图
function Runge_Kutta()
clc
clear all

%算法参数预设
h=0.1;Tmax=0.6;                                      %步长为h,计算时间为20s,n用于计数
xold=0;yold=2;                                     %设置初始值
%四阶龙格算法
for t=0:h:Tmax
t=0;
  m1=f(t,xold,yold)*h;k1=g(t,xold,yold)*h;         %将h放在这里减少计算舍去误差
  m2=f(t+0.5*h,xold+0.5*m1,yold+0.5*k1)*h;k2=g(t+0.5*h,xold+0.5*m1,yold+0.5*k1)*h;
  m3=f(t+0.5*h,xold+0.5*m2,yold+0.5*k2)*h;k3=g(t+0.5*h,xold+0.5*m2,yold+0.5*k2)*h;
   m4=f(t+h,xold+m3,yold+k3)*h;k4=g(t+h,xold+m3,yold+k3)*h;

%保存
if t==0; n=1;end
TT(n,:)=t;XX(n,:)=xold;YY(n,:)=yold;n=n+1;
%更新
xold=xold+(m1+2*m2+2*m3+m4)/6;
yold=yold+(k1+2*k2+2*k3+k4)/6;
t=t+h;
end
%计算结果输出
%plot(TT,XX,'b');grid on;
plot(XX,YY,'r');grid on
end


%要求的位移的导数表达式dx=f(t,x,y);
function dx=f(t,x,y)
dx=2*x+4*y;                 %不同的系统，需要修改
end
%要求的速度的导数表达式dy=g(t,x,y);
function dy=g(t,x,y)
%t是时间，x是位移，或者速度的导数，y是速度，或者速度的导数
dy=-x+6*y;                     %不同的系统，需要修改
end

伤痕累累

回复 4 # 博大广阔 的帖子

程序有错误。

博大广阔

t=0; 这行得删了。。具体你自己修改下。。我平常只是用该程序的循环部分。。还不行话给你个简洁的

伤痕累累

回复 6 # 博大广阔 的帖子

你好，请问你用两个function, x和y怎么耦合的呢？这样表示的是耦合吗？

博大广阔

给你个其他类型的吧。。function a=R_K()
%%%%%%%%%%%%%%%%%%%%%%%%%%%%计算初始值及计算时间长%%%%%%%%%%%%%%%%%%%%%%
u=0; %控制规律
t=0; tmax=10; h=0.01;
x10=[0;0];       %初始条件
%%%%%%%%%%%%%%%%%%%%%四阶龙格计算%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

n=1;        %循环程序需要的初始值
%%%%%%%%%%%%%%%%%%%%%%%主程序%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
while t<tmax   
[u,dX]=derives(t,x10,u);
       m1=h*dX;             t=t+0.5*h;X=x10+0.5*m1; %更新导数     
     [u,dX]=derives(t,X,u);
      m2=h*dX;            X=x10+0.5*m2;      %更新导数
    [u,dX]=derives(t,X,u);
     m3=h*dX;          t=t+0.5*h; X=x10+m3;     %更新导数
    [u,dX]=derives(t,X,u);
    m4=h*dX;        
     x10=x10+(1/6)*(m1+2*m2+2*m3+m4);
       y=x10(1);                                %输出
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%        
     XX(n)=y;  Time(n)=t;                 %用于保存
     n=n+1;  
   end
plot(Time,XX,'b');grid on;
end

%求导数
function [u,dX]=derives(t,X,u)
fi=10*exp(-0.5*t).*cos(3*t);          %振动系统对象
A=[-14.6254 0;1 0]; B=[1;0]; f=u+fi;
dX=A*X+B*f;   
end

伤痕累累

回复 8 # 博大广阔 的帖子

真是好人啊。厚道！！！

wzfly830613

厚道。呀。

WF1987

没看懂  惭愧

ChaChing

好奇, 不是有现成的ode45,ode23,...
为何不直接使用? 对求解有差异吗?

伤痕累累

回复 4 # 博大广阔 的帖子

谢谢你给的程序，我觉得用你的方法去画相图，有点麻烦了。如果用ode45，在时间项控制一下步长，效果不是一样的么。而且速度 也不用再去用个 g 函数另外求。

WF1987

楼上和楼上的楼上的方法可以试试

博大广阔

回复 13 # 伤痕累累 的帖子

效果是一样的，可能直接用matlab的效果还会好一些。因为他考虑的情况较多，比如变步长，数值溢出等。但是自己编也有其优点。