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

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

求解函數(shù)采用方法適用場合

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

先看一個簡單例子，y‘=y+3x/x^2，初值y(0)=-2，求解區(qū)間為[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) 。
對于剛性問題，數(shù)值解算法必須取很小步長才能獲得滿意的結(jié)果，導致計算量會大大增加。解決剛性問題需要有專門方法。非剛性算法可以求解剛性問題，只不過需要很長的計算時間。
例、假如點燃一個火柴，火焰球迅速增大直至某個臨界體積，然后，維持這一體積不變，原因是火焰球內(nèi)部燃燒耗費的氧氣和從球表面所獲氧氣達到平衡。其簡化模型如下：
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 計算的點數(shù)' num2str(length(t))]); 時間已過 0.003704 秒。 ode45 計算的點數(shù)157 >> L=1e-5; f=@(t,y) y^2-y^3; tic;[ t,y ]=ode45(f,[0,2/L],L); toc disp (['ode45 計算的點數(shù)' num2str(length(t))]); 時間已過 1.010016 秒。 ode45 計算的點數(shù)120565 >> L=1e-5; f=@(t,y) y^2-y^3; tic;[ t,y ]=ode15s(f,[0,2/L],L); toc disp (['ode45 計算的點數(shù)' num2str(length(t))]); 時間已過 0.189977 秒。 ode45 計算的點數(shù)98

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

關(guān)于MATLAB的學習：

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

歡迎大家加入我們的MATLAB學習交流群：
953314432
掃碼關(guān)注我們
發(fā)現(xiàn)更多精彩

總結(jié)

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

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯，歡迎將生活随笔推薦給好友。