非线性可视化（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=[];

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

end
function 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