Hamilton雅可比矩阵判别系统稳定性（学了个寂寞）_雅可比矩阵判断生态系统稳定性

参考文章内容：https://github.com/Santosh-Gudlavalleti/PH3500-Classical-Physics

最近在研究端口-Hamilton系统想做个控制方面的东西所以学习了一下相关内容上述网站上的内容说它采用了哈密尔顿、拉格朗日等方法对系统的稳定性进行了一个分析，这是一个国外物理课程的实验课里面还有它的PDF介绍。虽然我觉得这个文章和Hamilton控制没什么大关系但毕竟学都学了还是记录一下。

背景介绍

目的

了解由一节非线性微分方程链接的两个状态之间的关系、稳定性和其它参数

背景

整个系统是模拟观测兔子和狼数量的场景、假设场景内只有兔子和狼两个物种，他们的数量是随时间变化的，植物的数量非常大也就是说兔子的食物是无限的。文章中说她参考了Lotka-volterra模型但是对其做了简化

系统方程

原始系统

狼和兔的系统方程如下：

a-兔子的出生速率

b-兔子被吃掉和死亡的速率

c-狼的生长速度（因为狼的食物是兔子，所以狼的生长依赖兔子）、

d-狼的死亡速率（由于年龄和相互竞争关系）

r-兔子当前的种群规模

w-狼当前的数量规模

p-植物的数量，假设数量巨大，数量变化忽略不计。

下一时刻的系统

由于上述系统对时间没有明确的依赖，所以他是一个自洽/自治系统

得出雅可比矩阵

和我们看作是系统经历了一个时间 t 后的新的状态,用原本的规模加上一个 变化率和时间的乘积，或者可以看作是一个泰勒展开,那么我们可以得到:

同样的：

这样就能得到雅可比矩阵：

计算上述雅可比矩阵的行列式

系统分析

对于整个系统的初始条件设置为： r = 100 ; w = 60 ; a = 0.18 ; b = 0.004 ; c = 0.001 ; d = 0.01

我认为在这里abcd应该可以看做是系统本身的参数是无法进行更改的

随后可以将狼和羊的数量表示出来，这个图中表示狼和羊的数量是一只发生变化并有一个对应的关系：

对于实时变化而言，我们可以看出两个钟群的数量变化关系是稳定的

如果此时我们将羊的初始值设置为500，即将r=500

那么此时

再来看看这个图可以发现这个系统的稳定性就大大下降了，兔子和狼一度都 有灭绝的风险

如果我们将兔子的初始数量设置为 r=1000，我们会发现此时狼和羊都会灭绝

通过我们最初的系统方程

我们可以计算出狼和羊数量变化同时为0的时刻所对应的点的值应该是和

当取时，文中说计算雅可比矩阵的特征值为a和-d，但是我算的并不是，我觉得原文中此时所指的雅可比矩阵应该是这个矩阵

因为只有这个矩阵我算着特征值才是和原文相同的。

当取(0,0)点时

可以算出特征值分别为通过分析Jacobi矩阵的特征值。由于符号总是相反的，我们可以断定它是鞍点。同样可以用以下逻辑推理来解释:如果狼的数量变为0，兔子的数量指数地达到无穷大，而如果兔子的数量为0，狼的数量指数地下降到0。图的性质在两个方向上是不同的，因此可以得出它是一个鞍点。

在点我们可以得到特征值为很明显，解是纯虚的(解的实部为零)和彼此的复共轭。因此系统是椭圆型的，解是周期的。系统绕这一点在椭圆上振荡。振荡的频率就是特征值的几何平均值。因此振荡的频率（这块我不懂）

我们可以观察到:当兔子的数量或狼的数量达到特定值(即r = d/c或w = a/b)时，切线斜率的符号发生变化。这些是驱动系统走向稳定的点。椭圆轨迹解释了振荡的相位差，而相位差与椭圆的偏心率有关。

随后进行推导：

上下进行除法

进行分离变量并对两侧同时进行积分：

由于 K 是常数故而：

可以得到是守恒的

通过数量模拟我们可以得到 K 关于 rw 的值的变化如下，可以看到除了边界区， 大部分区域的 K 值为稳定的在 0.64 附近

然后文章中就得出了这么一个结论：

对上述问题进行抽样，我们可以借助哈密顿动力学和其他方程在数学上预测其 行为。若雅可比矩阵行列式中 df 的系数接近于零，则称为稳定系统，(- 0.05,0.05)可为系统稳定的相当大范围。在这个范围内，与实际状态量相比不 会有脉冲变化，因此系统处于稳定区（我还是不懂）

感觉学物理专业的同学应该能明白它再讲什么，对于我这么个学通信实在是云里雾里。

一些代码如下：

% Project    : Dynamics of First order Nonlinear Differential Equations
% Instructor : Prof. S. Kasiviswanathan
% Course     : Classical Physics (PH3500)
% Project By : Santosh Gudlavalleti
% Roll Number: EE19B055
 
clc, clear
 
disp('Welcome to the Display, Please choose an option');
disp('Option 1: The Final Plots will be displayed with demo Inputs');
disp('Option 2: The Live Phase Space Plots will be displayed with demo Inputs, (Quite Slow)');
disp('Option 3: The Final Plots will be displayed with demo Inputs of population and disease factor');
disp('Option 4: The Final Plots will be displayed with Custom Inputs');
 
s = input('Option: ');
 
if s == 1
    disp('The Final Plots will be displayed with Initial Rabbit Population as 100 and Wolf Population as 60');
    disp('Please Wait while the Graph is being Made');
    disp('Use Function: "main" to re-run the program');
    figure(1)
    fun_1(0,0.18,0.004,0,0.001,0.1,200);
    fun_3(0,0.18,0.004,0,0.001,0.1,20);
    fun_5(1000,0.18,0.004,60,0.001,0.1,200);
    
end
 
if s == 2
    disp ('The Live Phase Space Plots will be displayed with Initial Rabbit Population as 100 and Wolf Population as 60');
    disp ('Please Wait while the Graph is complete(Takes Time)');
    pause(2.0);
    fun_2(100,0.18,0.004,60,0.001,0.1,200);
end
 
 
if s == 3
    disp('The Final Plots will be displayed with Initial Rabbit Population as 100 and Wolf Population as 60');
    disp('The rabbits are affected by a disease after they reach 300 i.e ~85% of the max Population');
    disp('Please Wait while the Graph is being Made');
    disp('Use Function: "main" to re-run the program')
    fun_4(100,0.18,0.004,60,0.001,0.1,200,0.008);
end
 
if s == 4
    disp('Please enter the Custom Inputs')
    r = input('Enter Initial Rabbit Population (preferred 75-150): ');
    w = input('Enter Initial Wolf Population (preferred 50-75): ');
    t = input('Enter Time Duration (preferred >50): ');
    a = input('Enter growth rate of Rabbit(a) (preferred 0.15-0.20): '); 
    b = input('Enter death rate of Rabbit(b) (preferred 0.0025-0.0075): ');    
    c = input('Enter growth rate of Wolves(c) (preferred 0.0005-0.003): ');
    d = input('Enter death rate of Wolves(d) (preferred 0.05-0.25): ');
    
    fun_1(r,a,b,w,c,d,t);
end
 
if  isempty(s) || s~=1 && s~=2 && s~=3 && s~=4
    disp ('Enter A valid Option')
 
end
 
function [out1] = fun_1(r,a,b,w,c,d,t) %Current Phase Plot
h = 0.001;
ra = (r);
wo = (w);
 
for i = 1:(t/h)
     ra(i+1) = ra(i)*(1 + h*(a - b*wo(i)));
        wo(i+1) = wo(i)*(1 + h*(c*ra(i) - d));
%     if ra(i) > 1 && wo(i) > 1
%         ra(i+1) = ra(i)*(1 + h*(a - b*wo(i)));
%         wo(i+1) = wo(i)*(1 + h*(c*ra(i) - d));
%     else
%         ra(i+1) = 0;
%         wo(i+1) = 0;
%     end
end
   out1 = [ra wo];
   
   plot(ra, wo)
   title('Plot of population of Wolves vs population of Rabbits')
   xlabel('Population of Rabits')
   ylabel('Population of Wolves')
   grid on;
   legend({'Direction of flow: Anticlockwise'},'Location','northeast')
   hold on
 
function [out2] = fun_2(r,a,b,w,c,d,t) %Live Phase plot 
h = 0.001;
ra = (r);
wo = (w);
curve = animatedline
for i = 1:(t/h)
    if ra(i) > 0 && wo(i) > 0
        ra(i+1) = ra(i)*(1 + h*(a - b*wo(i)));
        wo(i+1) = wo(i)*(1 + h*(c*ra(i) - d));
    else
        ra(i+1) = 0;
        wo(i+1) = 0;
    end
 
    title('Plot of population of Wolves vs population of Rabbits')
    xlabel('Population of Rabits')
    ylabel('Population of Wolves')
    set(gca,'XLim',[50 150],'YLim',[25 75]);
    grid on;
    addpoints(curve, ra(i), wo(i));
    drawnow limitrate
 
end
 
function [out3] = fun_3(r,a,b,w,c,d,t)  %Variation with Time
 
h = 0.001;
ra = (r);
wo = (w);
 
for i = 1:(t/h)
    if ra(i) > 0 && wo(i) > 0
        ra(i+1) = ra(i)*(1 + h*(a - b*wo(i)));
        wo(i+1) = wo(i)*(1 + h*(c*ra(i) - d));
    else
        ra(i+1) = 0;
        wo(i+1) = 0;
    end
end
   figure(3)
   plot(0:t/h,ra)
   hold on
   plot (0:t/h,wo)
   title('Plots of population of Wolves vs Time and Population of Rabbits vs Time')
   xlabel('Time (in days)')
   ylabel('Population')
   legend({'y = Rabbits Count', 'y = Wolves Count'},'Location','northeast')
   set(gca,'XTickLabel',(0:(t/10):t));
   
   hold off 
 
    
function [out4] = fun_4(r,a,b,w,c,d,t,s) %Diseased, Variation with time
h = 0.001;
ra = (r);
wo = (w);
 
for i = 1:(t/h)
    if ra(i) > 0 && wo(i) > 0
        if ra(i) > 130
            ra(i+1) = ra(i)*(1 + h*(a - s*wo(i)));
            wo(i+1) = wo(i)*(1 + h*(c*ra(i) - d));
            
        else 
            ra(i+1) = ra(i)*(1 + h*(a - b*wo(i)));
            wo(i+1) = wo(i)*(1 + h*(c*ra(i) - d));
        end
    else
        ra(i+1) = 0;
        wo(i+1) = 0;
    end
end
   out4 = [ra wo];
    
   figure(4)
   plot(0:t/h,ra)
   hold on
   plot (0:t/h,wo)
   title('Plots of population of Wolves vs Time and Population of Rabbits vs Time')
   xlabel('Time (in days)')
   ylabel('Population')
   legend({'y = Rabbits Count', 'y = Wolves Count'},'Location','northeast')
   set(gca,'XTickLabel',(0:(t/10):t));
   
   hold off 
 
function [out5] = fun_5(r,a,b,w,c,d,t) %Demo Phase plots 
 
   figure(5)
   
   fun_1(50,0.18,0.004,30,0.001,0.1,200);
   hold on;
   fun_1(100,0.18,0.004,60,0.001,0.1,200);
   hold on;
   fun_1(150,0.18,0.004,90,0.001,0.1,200);
   title('Phase-Space plot of population of Wolves vs population of Rabbits')
   hold off;