4.8 微分方程
微分方程是數值計算中常見的問題，MATLAB提供了多種函數來計算微分方程的解。
4.8.1 常微分方程
眾所周知，對一些典型的常微分方程，能求解出它們的一般表達式，並用初始條件確定表達式中的任意常數。但實際中存在有這種解析解的常微分方程的範圍十分狹窄，往往只局限在線性常係數微分方程（含方程組），以及少數的線性變係數方程。對於更加廣泛的、非線性的一般的常微分方程，通常不存在初等函數解析解。由於實際問題求解的需要，求近似的數值解成為了解決問題的主要手段。常見的求數值解的方法有歐拉折線法、阿當姆斯法、龍格-庫塔法與吉爾法等。其中由於龍格-庫塔法的精度較高，計算量適中，所以使用的較廣泛。
數值解的最大優點是不受方程類型的限制，即可以求任何形式常微分方程的特解（在解存在的情況下），但是求出的解只能是數值解。
1．龍格-庫塔方法簡介
對於一階常微分方程的初值問題，在求解未知函數時，在點的值是已知的，並且根據高等數學中的中值定理，應有：
一般而言，在任意點，有：
當確定後，根據上述遞推公式能算出未知函數在點的一列數值解：
當然在遞推的過程中同樣存在著一個誤差累計的問題，實際計算中的遞推公式一般都進行過改造，龍格-庫塔公式為：
其中：
2．龍格-庫塔法的實現
基於龍格-庫塔法，MATLAB提供了ode系列函數求常微分方程的數值解。常用的有ode23 和ode45函數，其調用語法如下。
（1）[t,y]=ode23(filename,tspan,y0)：採用了二階、三階龍格-庫塔法進行計算。
（2）[t,y]=ode45(filename,tspan,y0)：採用了四階、五階龍格-庫塔法進行計算。
其中filename是定義f(t,y)的函數文件名，該函數文件必須返回一個列向量。tspan的形式為[t0,tf]，表示求解區間。y0是初始狀態列向量。t和y分別給出時間向量和相應的狀態向量。
這兩個函數分別採用了二階、三階龍格-庫塔法和四階、五階龍格-庫塔法，並採用自適應變步長的求解方法，即當解的變化較慢時採用較大的步長，從而使得計算的速度很快；當解的變化較快時步長會自動地變小，從而使得計算的精度很高。
【例4-45】 設有初值問題：
試求其數值解，並和精確解相比較，精確解為（）。
首先要建立微分方程所對應的函數文件myodefun.m，文件內容如下：
function y=myodefun(t,y)
% 建立函數文件myodefun.m
y=(y^2-t-2)/(4*(t+1));
建立myodefun函數之後，就可以調用ode23函數求解微分方程。
>> t0=0;
>> tf=10;
>> y0=2;
>> [t,y]=ode23 ('myodefun',[t0,tf],y0); % 求數值解
>> y1=sqrt(t+1)+1; % 求精確解
>> plot(t,y,'k.',t,y1,'r')
通過圖形比較，數值解用黑色圓點表示，精確解用紅色實線表示，如圖4-13所示。
【例4-46】 求下面無勁度系統微分方程組的數值解。
為了求解方程，首先要建立方程的m文件。本例中不妨建立名為rigid.m的函數文件，此文件用以描述給出的方程組，文件的內容如下：
function dy = rigid(t,y)
dy = zeros(3,1); % 一個列向量
dy(1) = y(2) * y(3);
dy(2) = -y(1) * y(3);
dy(3) = -0.51 * y(1) * y(2);
本例中，我們通過odeset函數對誤差進行控制，另外在時間[0 12]進行求解，0時刻初始條件向量為[0 1 1]。
>> options =odeset('RelTol',1e-4,'AbsTol',[1e-4 1e-4 1e-5]);% 誤差控制
>> [T,Y] = ode45(@rigid,[0 12],[0 11],options); % 求數值解
>>plot(T,Y(:,1),'-',T,Y(:,2),'-.',T,Y(:,3),'.') % 繪製結果圖
得到的結果如圖4-14所示。
圖4-13 常微分方程結果圖
圖4-14 常微分方程數值