上篇博客介紹了Matlab求解常微分方程組解析解的方法：博客地址
微分方程組復雜時，無法求出解析解時，就需要求其數值解，這里來介紹。
以下內容按照Matlab官方文檔提供的方程來展開（提議多看官方文檔）

介紹一下核心函數ode45()

一般形式：[t,y] = ode45(odefun,tspan,y0) 其中 tspan = [t0 tf]
功能介紹：求微分方程組 y′=f(t,y) 從 t0 到 tf 的積分，初始條件為 y0。解數組 y 中的每一行都與列向量 t 中返回的值相對應。

1. 一階微分方程求解(簡單調用即可)

方程：y’=2*t
代碼：

tspan=[1 6]; %定義自變量x的取值空間為1-6 y0=0;%定義因變量的初值，當x=1（x取值空間的第一個數）時，y0=0 [t,y]=ode45(@(t,y) 2*t,tspan,y0); %定義函數y'=2*t,使用ode45求解 plot(t,y,'-o'); %繪制求得的數值曲線

說明：簡單的odefun參數就是這個形式，@(x,y) fun 中間是空格，不能用逗號
結果：

2. 二階微分方程求解(引入函數文件)

方程：范德波爾方程 y1’’-u(1-y12)*y1’+y1=0;（這里設u=1）
代碼：
定義輸入的方程，以函數形式定義

function dydt=odefun(t,y) %二階方程為y1''-(1-y12)*y1'+y1=0; %降階為兩個方程：y1'=y2; % y2'=(1-y12)*y2-y1; %t雖然沒有使用，但必須要作為參數寫入 dydt=[y(2);(1-y(1)^2)*y(2)-y(1)];

求解作圖

tspan=[0 20]; %定義自變量x的取值空間為0-20 y0=[2;0];%定義因變量的初值，當x=0時，y1=2,y2=y1'=0; [t,y]=ode45(@odefun,tspan,y0); %使用ode45求解 %%下面為作圖過程，不解釋 plot(t,y(:,1),'-o',t,y(:,2),'-o') title('Solution of van der Pol Equation (\mu = 1) with ODE45'); xlabel('Time t'); ylabel('Solution y'); legend('y_1','y_2')

說明：先進行降階，使y2=y1’，這時候y2’=y1’’，然后把兩個方程寫成列向量形式就行了
結果：

3. 求解微分方程組（和2類似）

這里就和求解二階方程類似的，只不過不需要降階，仍舊需要一個函數來定義方程組。我們這里不用官方文檔的例子，用同學的循壞擺問題來進行演示。
方程：
給定的初值（w接近0，但實際上不能設置為0）：

代碼：
定義輸入的方程

function dRvw=func(t,Rvw) %% 函數功能：為ode45提供微分方程 %輸入：t:時間序列，就是θ；Rvw:因變量,Rvw(1)代表R,Rvw(2)代表v,Rvw(3)代表w %輸出：dRvw:因變量的一階微分,dRvw(1)代表dR,dRvw(2)代表dv,dRvw(3)代表dw %% 初始化因變量的一階微分，3×1的向量 dRvw=zeros(3,1); %% 參數初始化 r=0.01;u=0.1;g=9.8;M=10;m=1; %% 輸入微分方程式 dRvw(1)=-Rvw(2)/Rvw(3)-r; dRvw(2)=(M*g-(m*Rvw(3)^2*Rvw(1)+m*g*sin(t)) *exp(u*(t+pi/2)))/(Rvw(3)*(M+m)); dRvw(3)=Rvw(3)*r/Rvw(1)+2*Rvw(2)/Rvw(1)+g*cos(t)/(Rvw(3)*Rvw(1));

求解作圖

%% 參數初始化 start_Theta是θ的初始值 end_Theta是θ的結束值 %R是半徑初值；v是線速度初值；w是角速度初值 start_Theta=0;end_Theta=2*pi;R=1;v=0;w=1e-5; %% 使用ode45方法計算微分方程組func的數值解 %func是帶有方程組的函數 %[start_Theta end_Theta]是自變量范圍 %[R;v;w]是方程初值 %T是自變量的數組，Rvw是對應的因變量的數值。Rvw(:,1)=R;Rvw(:,2)=v;Rvw(:,3)=w; [T,Rvw]=ode45(@func,[start_Theta end_Theta],[R;v;w]); %% 組合自變量和因變量。TRvw(:,1)=θ;TRvw(:,2)=R;TRvw(:,3)=v;Trvw(:,4); TRvw=[T,Rvw]; %% 繪制自變量-因變量圖像.figure1是R-θ,figure2是v-θ,figure3是w-θ plot(T,Rvw(:,1),'-o',T,Rvw(:,2),'-o',T,Rvw(:,3),'-o') title('自變量-因變量圖像'); xlabel('Angle θ'); ylabel('Solution'); legend('R','v','w')

說明：注釋的應該是比較清楚的，把三個方程寫成列向量的形式就行

PS:有些人和我說不能運行，然后我看了他們出錯的原因，有點兒哭笑不得。出錯的基本上都是運行上面的dRvw=func（t，Rvw）這個函數的。說明一下，這是有參數的函數，不給參數不能直接運行的。下面的求解作圖腳本才是需要運行的哈，它調用了函數，才得到的結果。如果大家還發現什么問題，歡迎私戳或評論。

PS+ 有了PS之后，還是很多人問我參數的問題，我在這里直接把文件給大家：cupt.zip
提取碼：6k8n。
解壓后，有兩個文件，大家直接運行cupt.m即可。

結果：

4. 更多形式

講到這里，大部分我們用到的微分方程形式都可以求解了，Matlab還支持帶有時變項和額外參數的微分方程求解，這里不再贅述，大家可以自行參閱官方文檔。

總結

以上是生活随笔為你收集整理的Matlab求常微分方程组的数值解的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。