过去有画过常微分方程的向量场，通过向量场能够很形象的看出方程解的状态。

最近过节在家刷视频刷到了3Blue1Brown介绍微分方程的视频。

视频中对钟摆建立的微分方程组通过向量场的形式也很形象的表达了系统状态。

这里用matlab也实现一下，同时对三维情况也做了一个实现。

绘制的方法就是计算方程在二维或三维某个点的方向，然后把方向归一化，画出归一化的向量即可。

二维微分方程组如下：

三维微分方程组如下：

matlab代码如下：

二维情况：

clear all;close all;clc;

mu = 0.1;  
gdl = 1;  

x = -2:0.4:16;
y = -5:0.4:5;
[x,y] = meshgrid(x,y);

dx = y;
dy= -mu*y-gdl*sin(x);

d = sqrt(dx.^2+dy.^2);
dx = dx./d;
dy = dy./d;

quiver(x,y,dx,dy);
hold on;
[t,h] = ode45(@test,[0 100],[0.01 3]);    
plot(h(:,1),h(:,2),'r')

[t,h] = ode45(@test,[0 100],[0.01 2]);   
plot(h(:,1),h(:,2),'g')
grid on;
axis equal;

test.m：

function dy=test(t,x)
mu = 0.1;  
gdl = 1;  
dy=[x(2);
    -mu*x(2)-gdl*sin(x(1))];

结果：

红色和绿色的线为方程组的两个特解。

三维情况：

clear all;close all;clc;
a = 16;
b = 4;
c = 45;

x = -40:6:40;
y =-40:6:40;
z = 0:6:80;
[x,y,z] = meshgrid(x,y,z);

dx = a*(y-x);
dy = c*x - x.*z-y;
dz = x.*y-b*z;

d = sqrt(dx.^2+dy.^2+dz.^2);
dx = dx./d;
dy = dy./d;
dz = dz./d;
quiver3(x,y,z,dx,dy,dz);

hold on;
[t,h] = ode45(@test2,[0 30],[12 4 0]);
plot3(h(:,1),h(:,2),h(:,3));
grid on;
axis equal;

test2.m：

function dy=test2(t,y)
a = 16;
b = 4;
c = 45;

dy=[a*(y(2)-y(1));
    c*y(1)-y(1)*y(3)-y(2);
    y(1)*y(2)-b*y(3)];

结果：