matlab練習程序（常微分方程組向量場）

過去有畫過常微分方程的向量場，通過向量場能夠很形象的看出方程解的狀態。

最近過節在家刷視頻刷到了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)];

結果：