本節將介紹常微分方程初值問題的數值求解，主要內容分為三個部分：常微分方程數值求解的概念、求解函數及剛性問題。
一、常微分方程數值求解的一般概念
首先，凡含有參數，未知函數和未知函數導數 (或微分) 的方程，稱為微分方程，有時簡稱為方程，未知函數是一元函數的微分方程稱作常微分方程，未知函數是多元函數的微分方程稱作偏微分方程。微分方程中出現的未知函數最高階導數的階數，稱為微分方程的階。
常微分方程數值求解，指研究求解初值問題各類數值方法的構造、理論分析與數值實現問題。研究的主要對象為一階方程組初值問題：

其中，其中y=y(x)是未知函數,y(x0)=y0是初值條件,而f(x,y)是給定的二元函數。
所謂其數值解法，就是求y(x)在離散結點xn處的函數近似值yn 的方法 ，yn≈y(xn)。這些近似值稱為常微分方程初值問題的數值解。相鄰兩個結點之間的距離稱為步長。
這里主要介紹兩種：單步法和多步法
單步法：在計算y(n+1)時只用到前一步的y(n)，因此在有了初值之后就可以逐步往下計算，其代表是龍格- - 庫塔（ Runge- - Kutta ） 法
多步法：在計算y(n+1)時，除了用到前一步的值y(n)之外, , 還要用到y(n-p)（ p=1，2 ， … ，k，k＞0）的值, , 即前面的k步。其代表就是亞當斯 (Adams) 法
更多介紹可參見這個鏈接：
https://wenku.baidu.com/view/cafe161b9a6648d7c1c708a1284ac850ad0204fc.html
二、常微分方程求解函數
MATLAB 提供了多個求常微分方程初值問題數值解的函數，一般調用格式為：
[t,y]=solver(filename,tspan,y0,option)
其中，t和y分別給出時間向量和相應的數值解。solver為求常微分方程數值解的函數。filename 是定義 f(t ，y) 的函數名，該函數必須返回一個列向量。
tspan 形式為 [t0，tf]，表示求解區間。 y0 是初始狀態向量。Option 是可選參數，用于設置求解屬性，常用的屬性包括相對誤差值 RelTol(默認值是10的-3次方)和絕對誤差值 AbsTol( ( 默認值是10的-6次方)。
常微分方程數值求解函數的統一命名格式：
odennx
ode是Ordinary Differential Equation 的縮寫，是常微分方程的意思。
nn 是數字，代表所用方法的階數。例如，ode23采用2階龍格- - 庫塔（ Runge- - Kutta ）算法 ，用3階公式做誤差估計來調節步長，具有低等精度。ode45 采用4階龍格- - 庫塔算 法 ，用5階公式做誤差估計來調節步長，具有中等精度。
x是字母，用于標注方法的專門特征。例如， ode15s 、ode23s 中的“s”代表（ Stiff ），表示函數適用于剛性方程。
下表列出了求常微分方程數值解的函數：

求解函數采用方法適用場合

ode23	2 階或3階 Runge- - Kutta 算法，低精度	非剛性
ode45	4 階或5階 Runge- - Kutta 算法，中精度	非剛性
ode113	Adams 算法，精度可到10的-3次方至10的-6次方	非剛性，計算時間比 ode45
ode23t	梯形算法	適度剛性
ode15s	Gear’s 反向數值微分算法，中精度	剛性
ode23s	2階 Rosebrock 算法，低精度	剛性，當精度較低時，計算時間比 ode15s
ode23tb	梯形算法，低精度	剛性，當精度較低時，計算時間比 ode15s

先看一個簡單例子，y‘=y+3x/x^2，初值y(0)=-2，求解區間為[1,4]。

>> odefun=@(x,y) (y+3*x)/x^2; tspan=[1 4]; y0=-2; [x y]=ode45(odefun,tspan,y0) plot(x,y)

二、剛性問題
有一類常微分方程，其解的分量有的變化很快，有的變化很慢，且相差懸殊，這就是所謂的剛性問題 (Stiff) 。
對于剛性問題，數值解算法必須取很小步長才能獲得滿意的結果，導致計算量會大大增加。解決剛性問題需要有專門方法。非剛性算法可以求解剛性問題，只不過需要很長的計算時間。
例、假如點燃一個火柴，火焰球迅速增大直至某個臨界體積，然后，維持這一體積不變，原因是火焰球內部燃燒耗費的氧氣和從球表面所獲氧氣達到平衡。其簡化模型如下：
y’=y^2-y3，y(0)=L，0<=t<=2/L
其中， y(t) 代表火焰球半徑，初始半徑是λ ，它很小。分析 λ 的大小對方程求解過程的影響。

>> L=0.01; f=@(t,y) y^2-y^3; tic;[ t,y ]=ode45(f,[0,2/L],L); toc disp (['ode45 計算的點數' num2str(length(t))]); 時間已過 0.003704 秒。 ode45 計算的點數157 >> L=1e-5; f=@(t,y) y^2-y^3; tic;[ t,y ]=ode45(f,[0,2/L],L); toc disp (['ode45 計算的點數' num2str(length(t))]); 時間已過 1.010016 秒。 ode45 計算的點數120565 >> L=1e-5; f=@(t,y) y^2-y^3; tic;[ t,y ]=ode15s(f,[0,2/L],L); toc disp (['ode45 計算的點數' num2str(length(t))]); 時間已過 0.189977 秒。 ode45 計算的點數98

注：tic 和 toc 函數 用來記錄 微分方程求解 命令執行的時間 ，使用 tic 函數啟動計時器，使用 toc 函數顯示從計時器啟動到當前所經歷的時間。
上述采用了三種不同方法，可以發現，第一種的運行結果表明這時常微分方程不算很剛性；第二種這時計算時間明顯加長，計算的點數劇增，表明這時常微分方程表現為剛性；因此對于剛性問題，我們需要改變求解算法，我們選擇以“s”結尾的函數，例如第三種方法他們專門用于求解剛性方程。計算時間大大縮短，計算的點數大大減少。
常微分方程數值解法的研究已發展得相當成熟，理論上也頗為完善，小編由于能力有限，只能了解個大概，本文講解也就寫的是比較基礎的一些方面，如果大家有更多需求，可以自己查閱相關資料。

關于MATLAB的學習：

大家可以關注我們的知乎專欄——數據可視化和數據分析中matlab的使用：
https://zhuanlan.zhihu.com/c_1131568134137692160

歡迎大家加入我們的MATLAB學習交流群：
953314432
掃碼關注我們
發現更多精彩

總結

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

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