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

前文我们介绍了线性系统的相图绘制。

这篇文章里，我们用几个例子，来介绍非线性系统的相图的绘制方法。

这里举的例子都是自治系统方程的例子，也就是方程结果与t0的初始取值无关（时不变系统），不含外部周期性驱动力之类的与t相关的量。因为描述自治系统，只需要知道系统的空间上的各个变量的导数，然后组成相空间即可。而时变系统各个状态都会随时间变化，无法用静态的相图去定性分析。

所以对于自治二阶系统，二阶的相平面已经可以完全的描述出系统的运动状态，无论线性还是非线性。通常的步骤可以分为两步：（1）计算出每一个点的dy和ddy导数，（2）根据每个点得到的向量，绘制出向量场对应的流线图。

以《非线性系统》这本书中给出的一个例子作为展示。其中二阶非线性方程的公式如下：

绘制出空间中每一个点的系统导数，绘制出流线，即可得到这个非线性系统的相图。

可以看到，非线性系统的相平面，可能拥有不止一个平衡点。但是这些平衡点附近的邻域，也都可以线性化，用上一节中的线性系统的相图来解释其行为。比如上图，就拥有两个稳定结点，和一个鞍点。

接下来再介绍一种只有在非线性条件下，才会出现的一种经典相平面图案：极限环。

以经典的Van der Pol方程为例，这个方程的形式如下：

后面的ε为一个常数，ε越大方程的非线性越大。取ε=0.3，绘制出相空间图：

整个空间只有一个平衡点，即极限环内部发散的螺旋点。为了直观的了解极限环的特性，选取不同的初始条件，对方程进行时域的求解，得到下面的结果：

在极限环外的点，刚开始呈现出收敛的波形，但是随后振幅保持恒定。在极限环里面的点，刚开始呈现出发散的波形，但是随后振幅也保持稳定。两者都被空间中的极限环所吸引，进行周而复始的画圈圈运动中。

如果方程中某些系数发生变化（比如上面方程的ε由负变为正），会导致整个相空间发生改变。这种改变会使得平衡点性质发生变化，导致整个系统的性质发生变化，我们称系统在这个位置发生了分岔。比如Van der Pol方程的ε由负变为正，平衡点由稳定变为发散，导致空间中稳定位置由一个点变为一个运动的极限环。

如果实验中能够观测到振动信号，也可以用绘制相平面的方法，观测信号的特性。下面用三幅图给出了测量过程中常见周期信号的相平面图：

第一幅图为典型的线性振动图，波形为正弦曲线，相平面为一个圆。

第二幅图为典型的极限环振动。这类振动有个特点，就是即使受到扰动，振幅也能始终保持在一个稳定的值。在实际实验中如果发现这个现象，而且没有外部周期驱动力（方程为自治系统），则基本能确定其为极限环运动。

第三幅图为典型的高维非线性。因为相平面内的流线不会交叉。这种交叉曲线是高维空间在二维平面上的投影。图中展示的是高维非线性中的倍周期现象的模拟。这个在后面文章中会介绍到。

当然，实际试验中由于噪声，测量的结果不会这么好看。比如最终的结果如果杂乱无章没有规律，理论上肯定是混沌现象出现了，但实际上很有可能是噪声太大而已。

后面附上本章绘图用到的matlab代码：

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

参考资料：

[1] Khalil H K . 非线性系统(第3版)[M].电子工业出版社,2005.

[2]Morris W. Hirsch. Differential equations, dynamical systems, and an introduction to chaos [M] Academic Press