MATLAB非线性可视化之线性系统相图

clear

clc

close all

%线性系统，用来展示相空间的用途
%小阻尼震荡（负实，共轭虚部）

A=[0,1;

    -0.925,-0.3];
    
%大阻尼震荡（负实，无虚部）

 A=[0,1;
     -0.96,-2];
     
%只有阻尼（一个为0一个为负）

% A=[0,1;

%     0,-2];

%只有弹簧（实部都为0）

% A=[0,1;

%     -1,0];

%鞍点

% A=[0.15,0;

%     0,-0.15];
[y,dy,u,v]=PhaseSpace_Linear(-5:0.05:5,-4:0.05:4,A);
h=1e-2;

x0=0:h:20;N=length(x0);

%相同的微分方程，不同的初始状态

[y1,~]=ODE_RK4_hyh(x0,h,[4;-2],A);

[y2,~]=ODE_RK4_hyh(x0,h,[0;-4],A);

[y3,~]=ODE_RK4_hyh(x0,h,[-3;3],A);


figure()

hold on

plot(x0,y1(1,:),'color',[1,0,0])

plot(x0,y2(1,:),'color',[0,1,0])

plot(x0,y3(1,:),'color',[0,0,1])

%绘制箭头

xy_lim=axis;

xy_ratio=get(gca, 'DataAspectRatio');xy_ratio=xy_ratio(1:2);

N_Arrow=5;

for k=1:N_Arrow

    N_k=round(N/N_Arrow*(k-0.83));
    
    plot_arrow( xy_lim,xy_ratio,[x0(N_k),y1(1,N_k)]...
    
        ,[1,0,0],1, [x0(N_k)-x0(N_k-10),y1(1,N_k)-y1(1,N_k-10)] )
        
   plot_arrow( xy_lim,xy_ratio,[x0(N_k),y2(1,N_k)]...
   
        ,[0,1,0],1, [x0(N_k)-x0(N_k-10),y2(1,N_k)-y2(1,N_k-10)] )
        
    plot_arrow( xy_lim,xy_ratio,[x0(N_k),y3(1,N_k)]...
    
        ,[0,0,1],1, [x0(N_k)-x0(N_k-10),y3(1,N_k)-y3(1,N_k-10)] )
        
end

hold off

legend("条件1","条件2","条件3")

box on

xlabel('t');ylabel('x');

figure()

hold on

plot(y1(1,:),y1(2,:),'color',[0.8,0,0],'linewidth',1)

plot(y2(1,:),y2(2,:),'color',[0,0.8,0],'linewidth',1)

plot(y3(1,:),y3(2,:),'color',[0,0,0.8],'linewidth',1)

xlim([-5,5]);ylim([-4,4])

%绘制箭头xy_lim=axis;

xy_ratio=get(gca, 'DataAspectRatio');
xy_ratio=xy_ratio(1:2);

N_Arrow=3;

for k=1:N_Arrow

    N_k=round(N/N_Arrow*(k-0.9));
    
    plot_arrow( xy_lim,xy_ratio,[y1(1,N_k),y1(2,N_k)]...
    
        ,[0.8,0,0],2, [y1(1,N_k)-y1(1,N_k-10),y1(2,N_k)-y1(2,N_k-10)] )
        
    plot_arrow( xy_lim,xy_ratio,[y2(1,N_k),y2(2,N_k)]...
    
        ,[0,0.8,0],2, [y2(1,N_k)-y2(1,N_k-10),y2(2,N_k)-y2(2,N_k-10)] )
        
    plot_arrow( xy_lim,xy_ratio,[y3(1,N_k),y3(2,N_k)]...
    
        ,[0,0,0.8],2, [y3(1,N_k)-y3(1,N_k-10),y3(2,N_k)-y3(2,N_k-10)] )
        
end

hold off

xlim([-5,5]);ylim([-4,4])

box on

xlabel('x');ylabel('dx');
% figure()hold onh=streamslice(y,dy,u,v);xlim([-5,5]);ylim([-4,4])box onxlabel('x');ylabel('dx');


function [F,Output]=Fdydx(x,y,Input)

%形式为Y'=F(x,Y)的方程，参见数值分析求解常系数微分方程相关知识%高次用列向量表示，F=[dy(1);dy(2)];y(1)为函数，y(2)为函数导数A=Input;F=A*y;Output=[];end
function [y,dy,u,v]=PhaseSpace_Linear(y1,y2,A)[y,dy]=meshgrid(y1,y2);%初始化网格u=zeros(size(y));v=u;for k=1:numel(y)    %计算网格上每一个点的上的方向    F=Fdydx(0,[y(k);dy(k)],A);    u(k)=F(1);    v(k)=F(2);endend

function [y,Output]=ODE_RK4_hyh(x,h,y0,Input)%4阶RK方法%h间隔为常数的算法y=zeros(size(y0,1),size(x,2));y(:,1)=y0;for ii=1:length(x)-1    yn=y(:,ii);    xn=x(ii);    [K1,~]=Fdydx(xn    ,yn       ,Input);    [K2,~]=Fdydx(xn h/2,yn h/2*K1,Input);    [K3,~]=Fdydx(xn h/2,yn h/2*K2,Input);    [K4,~]=Fdydx(xn h  ,yn h*K3  ,Input);    y(:,ii 1)=yn h/6*(K1 2*K2 2*K3 K4);endOutput=[];end
function plot_arrow(xy_lim,xy_ratio,xy_arrow,arrow_color,arrow_width,arrow_direction)%初始化箭头形状（归一化的形状）arrow_0=[0,0;-1,0.5;-1,-0.5];%对方向进行归一化a_dn=arrow_direction(:)./xy_ratio(:);a_dn=a_dn/sqrt(sum(a_dn.^2));d=(xy_lim(4)-xy_lim(3) xy_lim(2)-xy_lim(1))/2;%箭头对窗口缩放arrow_1=arrow_0*arrow_width*0.03*d;%箭头旋转arrow_2=arrow_1*[a_dn(1),a_dn(2);-a_dn(2),a_dn(1)];%箭头变形xy_ratio_n=xy_ratio/sqrt(sum(xy_ratio.^2));%对比例尺归一化arrow_3=arrow_2.*xy_ratio_n xy_arrow;fill(arrow_3(:,1),arrow_3(:,2),arrow_color,'EdgeColor','none')end