用Matlab求解微分方程

解微分方程有兩種解，一種是解析解，一種是數值解，這兩種分別對應不同的解法

解析解

利用dsolve函數進行求解

syms  x;
s = dsolve('eq1,eq2,...', ’cond1,cond2,...', 'v');
%eq：微分方程
%cond：條件
%v：獨立變量
%形如：方程：y'= f(t,y),初值：y(t0) = y0

1.求解析解

 dsolve('Du = 1+ u^2','t')
 ans =
 
 tan(C2 + t)
          1i
         -1i

求 的解析解

s = dsolve('D2y=3*y+2*x','x');
% D2y用以表示y的二階導數，默認是以t為自變量的，所以最好指明自變量為x.
syms y(x);
s = dsolve([diff(y,x,2) == 3*y+2*x], [y(0) == 5]) 
% diff內依次是函數、自變量、微分階數，方程用==表示相等而不是賦值

2.初值問題

求初值問題

s = dsolve('Dy = y - 2*t / y','y(0) =1');

3.邊界問題

求邊界問題

s = dsolve('x*D2y - 3*Dy =x^2','y(1)=0','y(5) = 0','x');

4.高階方程

求解方程

s=dsolve('D2y =cos(2*x) - y','y(0) =1','Dy(0) = 0','x');
simplify(s);
(eqn,cond,‘IgnoreAnalyticConstraints’,false) %設置不化簡結果

5.方程組問題

求解方程組

[f,g]= dsolve('Df = f + g','Dg = -f + g','f(0)=1','g(0) = 2','x'); 

一些例子

 dsolve('D2y+4*Dy+29*y = 0','y(0) = 0','Dy(0)= 15 ','x')
 
ans =
 
3*sin(5*x)*exp(-2*x)

[x y z] =   dsolve('Dx = 2*x-3*y+3*z','Dy = 4*x-5*y+3*z','Dz = 4*x-4*y+2*z')
 
x = 
C7*exp(2*t) + C8*exp(-t)
y =
C7*exp(2*t) + C8*exp(-t) + C9*exp(-2*t) 
z =
C7*exp(2*t) + C9*exp(-2*t)
%可以對其進行簡化操作
 x = simplify(x) 
x = C7*exp(2*t) + C8*exp(-t)
 y = simplify(y) 
y =exp(-2*t)*(C9 + C8*exp(t) + C7*exp(4*t))

數值解

%龍格庫塔法(Runge-Kutta法)
xfun=@(t,x)0.3.*x.*(1-x/8); %定義賦值函數r=0.3,k=8
[tout,xout]=ode45(fun,[0,40],0.1)  %方程數值解，四五階RK法
[tout,xout]=ode23(xfun,[t0,tfinal],x0) %二三階RK法
%%
ode系列數值求解形如  /  =   ( , )的微分方程組, 並繪圖。
xfun: 輸入參數，函數必須恰有t,x兩個變量，用函數文件定義的fun.m則用@fun或‘fun’調用。
t0：輸入參數，t的初始值。
tfinal：輸入參數，t的終值。
x0：輸入參數，x的初始值。
tout: 離散的自變量值， xout: 離散的函數值。
%%

同時也有一些其他的求解語句和輸出語句

%%
其他的求解語句
ode45        ode113        ode15s 
ode23s       ode23t        ode23tb 
其他的輸出語句
odeplot     odeprint        
odephas2    odephas3
%%

一個例子

求的數值解

首先對該方程進行換元然后建立m文件

function fyy=rhf(t,x)
    fyy=[y(1).*(1-y(2).^2)+y(2);y(1)];
end

最后計算數值解

y0=[0.25,0]’;
[t,y]=ode23(‘rhf’,[0,0.25],y0);
plot(t,y)

一些例子

%vdp1000.m

function dy = vdp1000(t,y)
dy = zeros(2,1);
dy(1) = y(2);
dy(2) = 1000*(1-y^2)*y2-y1;
end
%命令行輸入
[T,Y] = ode15s('vdp1000',[0 3000],[2 0]);%第一個參數是文件名，第二個參數是初始時間和終止時間第三個參數是y1和y2的初值
plot(T,Y(:,1),'-');
%結果是T時間 

plot(T,Y(:,1),'-k');，畫Y數組中的第一列數隨着T的變化曲線，‘-k’表示顏色黑色實線，

%定義函數
function dy=eq1(x,y)
    dy=zeros(2,1);
    dy(1)=y(2);
    dy(2)=1/5*sqrt(1+y(1)^2)/(1-x);
end

調用
x0=0;
xf=0.9999;
[x,y]=ode15s('eq1',[x0 xf],[0 0]);
plot(x,y(:,1),'-') 
hold on
y=0:0．01:2;
plot(1,y,'*')

微分方程模型

1.種群增長Logistic模型

• N(t)表示在時刻 t時刻種群數量
• r 表示種群的內稟增長率，即在沒有資源限制下的種群增長率
• K表示環境載量，反映資源環境對種群增長的制約作用

2.生物種群競爭模型

• 用 N₁(t)和N₂(t) 分別表示在時刻 甲、乙兩個種群數量。
• a₁₁ 表示種群甲自身的被抑制的情況
• a₁₂ 表示種群乙對種群甲的競爭力

參考： https://zhuanlan.zhihu.com/p/162296418