非线性可视化（2）非线性相图

%1二维相空间
%非线性clear
clc
close all
%1多平衡点的非线性系统
%参考 非线性系统（中文翻译第三版） Khalil P32
[y,dy]=meshgrid(-0.5:0.02:1.5,-0.5:0.02:1.5);%初始化网格
u=zeros(size(y));v=u;
for k=1:numel(y)    
    %计算网格上每一个点的上的方向
    
    F=Fdydx(0,[y(k);dy(k)],1);
    
    u(k)=F(1);
    
    v(k)=F(2);
    
end
figure()
streamslice(y,dy,u,v,4)
xlabel('y')
ylabel('dy')
box on%2极限环
[y,dy]=meshgrid(-4:0.1:4,-4:0.1:4);%初始化网格
u=zeros(size(y))
;v=u;
for k=1:numel(y)
    %计算网格上每一个点的上的方向
    
    F=Fdydx(0,[y(k);dy(k)],2);
    
    u(k)=F(1);
    
    v(k)=F(2);
    
end
figure()
streamslice(y,dy,u,v,3)
xlabel('y'
)ylabel('dy')
box on
xlim([-4,4]);ylim([-4,4])figure()
streamslice(y,dy,u,v,3)
xlabel('y')
ylabel('dy'
)box on
hold on
t2=0:0.1:50;
[y2,~]=ODE_RK4_hyh(t2,t2(2)-t2(1),[-4;4],2);
plot(y2(1,:),y2(2,:),'r','linewidth',2)
[y3,~]=ODE_RK4_hyh(t2,t2(2)-t2(1),[-0.1;0.1],2);
plot(y3(1,:),y3(2,:),'g','linewidth',2)
hold off
xlim([-4,4]);ylim([-4,4])
figure()
plot(t2,y2(1,:),'r','linewidth',2);ylim([-4,4])
figure()
plot(t2,y3(1,:),'g','linewidth',2);ylim([-4,4])%3不同时域信号在相图上的样子
figure()
t1=0:0.01:1;
y1=sin(10*pi*t1);
dy1=diffhyh(y1,2);
subplot(2,3,1)plot(t1,y1);ylim([-1.5,1.5])
xlabel('t');ylabel('y')
subplot(2,3,4)
plot(y1,dy1);xlim([-1.5,1.5])
xlabel('y');ylabel('dy')
%非线性极限环t2=0:0.01:18;
[y,Output]=ODE_RK4_hyh(t2,t2(2)-t2(1),[-1;1.161],2);
y2=y(1,:);
dy2=y(2,:);
subplot(2,3,2)
plot(t2,y2);ylim([-2.5,2.5]);xlim([0,18])
xlabel('t');ylabel('y')
subplot(2,3,5)
plot(y2,dy2);xlim([-2.5,2.5])
xlabel('y');ylabel('dy')%混沌（模拟的信号，并没有给出实际的系统）
%二维系统下，相平面是可以完全描述系统的，不可能出现交叉。
%如果出现轨线交叉，则说明系统的维度大于2，相平面展示的图形只是高维系统在2维上的投影。
t3=0:0.01:2;
Gas=@(t,a,b,c) c*exp(-b*(t-a).^2);
y3=Gas(t3,0,120,1) Gas(t3,1,120,1) Gas(t3,2,120,1);
y3=y3 Gas(t3,0.3,90,0.5) Gas(t3,1.3,90,0.5);
y3=y3-Gas(t3,0.6,300,0.3)-Gas(t3,1.6,300,0.3);
dy3=diffhyh(y3,2);
subplot(2,3,3)
plot(t3,y3);ylim([-0.5,1.5])
xlabel('t');ylabel('y')
subplot(2,3,6)
plot(y3,dy3);xlim([-0.5,1.5]);ylim([-0.15,0.15]
)xlabel('y');ylabel('dy')function [F,Output]=Fdydx(x,y,Input)
%形式为Y'=F(x,Y)的方程，参见数值分析求解常系数微分方程相关知识
%高次用列向量表示，F=[dy(1);dy(2)];y(1)为函数，y(2)为函数导数
switch Input
    %示例1
    
    case 1
    
        p=[83.72,-226.31,229.62,-103.79,17.76,0];
        
        dy(1)=0.5*(-polyval(p,y(1)) y(2));
        
        dy(2)=0.2*(-y(1)-1.5*y(2) 1.2);
        
        F=[dy(1);dy(2)];
        
    %示例2
    
    case 2
    
        dy(1)=y(2);
        
        dy(2)=-y(1) 0.3*(1-y(1)^2)*y(2);%0.3
        
        F=[dy(1);dy(2)];endOutput=[];
endfunction [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=[];
endfunction Fdiff=diffhyh(F,dim)
%采用2阶中心差分，边缘采用一阶向前或向后差分（边缘没有搞迎风差分，精度要求不高）
%dim，差分的维度方向，dim=1对应着矩阵向下方向的差分，dim=2对应的向右
diffF1=diff(F,1,dim);if dim==1%向下
    Fdiff=([zeros(1,size(F,2));diffF1] [diffF1;zeros(1,size(F,2))])/2;
    
    Fdiff(1,:)=diffF1(1,:);
    
    Fdiff(end,:)=diffF1(end,:);elseif dim==2%向右
    Fdiff=([zeros(size(F,1),1),diffF1] [diffF1,zeros(size(F,1),1)])/2;
    
    Fdiff(:,1)=diffF1(:,1);
    
    Fdiff(:,end)=diffF1(:,end);
    
end
end