首页/文章/ 详情

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

过冷水

2年前浏览2467

我们在前面的多摆模型中，利用多摆的微分方程模型，求解出了多摆每时每刻的位置随时间的变化。当然那是一个高度复杂的非线性模型，难以上手分析。

这篇文章，我们首先利用一个二阶的线性模型进行求解，并引入微分方程定性分析中常用的工具——相图。

首先考虑下面这个经典的二阶阻尼振动方程：

将它整理为线性系统，如dx=Ax形式的样子：

矩阵A是一个二阶矩阵。我们取k=0.925，c=0.3。则矩阵A的特征根可以求得为：λ1=-0.15 0.95i，λ2=-0.15-0.95i。特征值的实部都是负数，代表系统收敛。特征值的虚部不为零，存在震荡现象。

我们取三组不同的位置和速度，作为初始状态点进行计算，可以得到下面的时间-位移曲线：

从时域曲线上看，这些曲线之间似乎没有什么关系。

但是如果我们绘制出 位移-速度 图，可以得到：

三条曲线以一种相同的方式逐渐向内汇集。我们把这种曲线叫做轨迹。这些轨迹如同在流场中的水流一般，在一个场空间内运动。我们把这个二维场叫做相空间。这个场每一个点的向量为[dx,ddx]。这个向量场组成的相空间，能够完全的描述出整个系统，沿着向量场的箭头逐步连线，便可以预测出任意一点未来的运动状态。把这个向量场绘制出来，如下图：

我们把这种螺旋点叫做稳定的焦点，特点是特征根的实根为负数，虚部不为零。这种相空间在时域上，对应着震荡收敛的解。

同理，当两个特征根都为正，虚部不为零时，则会出现发散的螺旋点，称为不稳定的焦点。

当特征根不存在虚部，也对应着4种情况：特征根同号，特征根异号，特征根一个为0，特征根两个为0。

特征根同号的话，还是以前面的阻尼举例子。继续增大阻尼，使得系统达到过阻尼状态。我们取k=0.96，c=2，此时特征根等于λ1=-0.8，λ2=-1.2。这种稳定点叫做稳定的结点。同样取三个不同的初始值，计算出相轨迹，与计算得到的相空间叠加，如下图：

如果把前面单自由度震荡系统中的弹簧去掉，变成只有阻尼c的滑块。此时取k=0，c=2，此时特征根为-2和0。代表只有一个吸引分量，另一个维度保持既不吸引也不排斥的状态。这时候表现为滑块速度逐渐降低，最终停到哪里算哪里，整个dx=0这根线都是它的稳定位置。

同样，如果不要阻尼项，只有弹簧，则变成了永远震荡下去的弹簧振子。此时取k=1，c=0，特征根为-1i和1i两个虚数。代表只有震荡分量，没有吸引或排斥分量。此时系统为一个标准的正圆。各个轨线在各自的圆上不停的运动互不影响。

还剩一个特征根是一正一负的，用弹簧振子举例不是直观，所以就不拿弹簧说事了。此时的中心点叫做鞍点，由于特征根的特性，导致了一部分吸引一部分排斥的特点。很多常见的不稳定平衡点，比如单摆的180°最高点。表现为“你以为我要在这里停下，其实我只是路过”。

绘制二维线性系统的相平面代码如下：

这篇文章，我们首先利用一个二阶的线性模型进行求解，并引入微分方程定性分析中常用的工具——相图。

首先考虑下面这个经典的二阶阻尼振动方程:

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

将它整理为线性系统，如dx=Ax形式的样子：
矩阵A是一个二阶矩阵。我们取k=0.925，c=0.3。则矩阵A的特征根可以求得为：λ1=-0.15 0.95i，λ2=-0.15-0.95i。特征值的实部都是负数，代表系统收敛。特征值的虚部不为零，存在震荡现象。
我们取三组不同的位置和速度，作为初始状态点进行计算，可以得到下面的时间-位移曲线：
从时域曲线上看，这些曲线之间似乎没有什么关系。
但是如果我们绘制出 位移-速度 图，可以得到：
三条曲线以一种相同的方式逐渐向内汇集。我们把这种曲线叫做轨迹。这些轨迹如同在流场中的水流一般，在一个场空间内运动。我们把这个二维场叫做相空间。这个场每一个点的向量为[dx,ddx]。这个向量场组成的相空间，能够完全的描述出整个系统，沿着向量场的箭头逐步连线，便可以预测出任意一点未来的运动状态。把这个向量场绘制出来，如下图：
我们把这种螺旋点叫做稳定的焦点，特点是特征根的实根为负数，虚部不为零。这种相空间在时域上，对应着震荡收敛的解。
同理，当两个特征根都为正，虚部不为零时，则会出现发散的螺旋点，称为不稳定的焦点。
当特征根不存在虚部，也对应着4种情况：特征根同号，特征根异号，特征根一个为0，特征根两个为0。
特征根同号的话，还是以前面的阻尼举例子。继续增大阻尼，使得系统达到过阻尼状态。我们取k=0.96，c=2，此时特征根等于λ1=-0.8，λ2=-1.2。这种稳定点叫做稳定的结点。同样取三个不同的初始值，计算出相轨迹，与计算得到的相空间叠加，如下图：
如果把前面单自由度震荡系统中的弹簧去掉，变成只有阻尼c的滑块。此时取k=0，c=2，此时特征根为-2和0。代表只有一个吸引分量，另一个维度保持既不吸引也不排斥的状态。这时候表现为滑块速度逐渐降低，最终停到哪里算哪里，整个dx=0这根线都是它的稳定位置。
同样，如果不要阻尼项，只有弹簧，则变成了永远震荡下去的弹簧振子。此时取k=1，c=0，特征根为-1i和1i两个虚数。代表只有震荡分量，没有吸引或排斥分量。此时系统为一个标准的正圆。各个轨线在各自的圆上不停的运动互不影响。
还剩一个特征根是一正一负的，用弹簧振子举例不是直观，所以就不拿弹簧说事了。此时的中心点叫做鞍点，由于特征根的特性，导致了一部分吸引一部分排斥的特点。很多常见的不稳定平衡点，比如单摆的180°最高点。表现为“你以为我要在这里停下，其实我只是路过”。
绘制二维线性系统的相平面代码如下：
参考资料：
[1]Morris W. Hirsch. 微分方程、动力系统与混沌导论[M]. 人民邮电出版社, 2008.
[2]刘秉正. 非线性动力学与混沌基础[M]. 东北师范大学出版社, 1994.

参考资料：

[1]Morris W. Hirsch. 微分方程、动力系统与混沌导论[M]. 人民邮电出版社, 2008.

[2]刘秉正. 非线性动力学与混沌基础[M]. 东北师范大学出版社, 1994.

理论科普

著作权归作者所有，欢迎分享，未经许可，不得转载

首次发布时间：2021-12-08