作者：a313291036

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 on
h=streamslice(y,dy,u,v);
xlim([-5,5]);ylim([-4,4])
box on
xlabel('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);
end
end


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);
end
Output=[];
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

本文出自：https://mp.weixin.qq.com/s/Bzm14vTWzQcCwoPY7oSO4Q