10 常微分方程初值問題的數值解法
10.1 引言
包含自變量、未知函數以及未知函數導數或微分的方程稱為微分方程。在微分方程中，如果自變量的個數只有一個，就稱為常微分方程；如果自變量個數兩個及以上，就稱為偏微分方程。微分方程中出現的未知函數最高階導數的階稱為微分方程的階。如果未知函數\(y\)及其各階導數：\(y',y'',\cdots,y^{(n)}\)都是一次的，則稱它是線性的，否則是非線性的。
高等數學中給出了一些典型的常微分方程的求解方法，但實際應用中，很多微分方程往往是無法直接求出解析解的。因此，很有必要研究微分方程的數值解法。一階常微分方程的初值問題的一般形式是：
\[\left\{\begin{aligned} &y'=f(x,y),a\le x\le b\\ &y(x_0)=y_0 \end{aligned}\right.\tag{10-1-1} \]
其中函數\(f(x,y)\)的定義域\(D=\{(x,y)|a\le x\le b,c\le y\le d\}\)。
對於這一類的初值問題而言，如果\(f(x,y)\)在區域\(D\)上滿足李普希茲（Lipschitz）條件，那么初值問題在區間\([a,b]\)上具有唯一連續可微解。李普希茲條件具體表述如下：
\[\exists L>0,\forall x\in[a,b],y\in[c,d],|f(x,y_1)-f(x,y_2)|\le L|y_1-y_2| \]
所謂微分方程的數值解法，就是尋找方程(10-1-1)的解\(y(x)\)在一系列離散節點：\(x_1<x_2<\cdots<x_n\)，上的近似值\(y_1,y_2,\cdots,y_n\)。相鄰兩個節點之間的距離\(h=x_{n+1}-x_n\)稱為步長，如果所設置的步長是定值，那么就有\(x_n=x_0+nh\)。常微分方程的數值解法的基本出發點就是將連續的問題離散化，並且采用“步進式”的解法，即求解過程順着節點排序的次序一步一步向前推進。描述這一類算法，就是要用已知的\(y_i,y_{i-1},\cdots,y_0\)構建出計算\(y_{i+1}\)的遞推式，而建立遞推公式的基本方法就是在這些節點上應用數值積分、微分、泰勒展開等方法。
10.2 歐拉（Euler）方法
歐拉格式：
一種最簡單的方法便是使用差商代替微分：
\[y'(x_n)=f(x_n,y_n)\approx\frac{y_{n+1}-y_n}{h} \]
於是就可以得到歐拉格式
\[\left\{\begin{aligned} y_{n+1}&=y_n+hf(x_n,y_n)\\ y(x_0)&=y_0 \end{aligned}\right.\quad n=0,1,2,\cdots \]
隱式歐拉格式：
歐拉格式是用差商的值來近似點\((x_n,y_n)\)處的導數值，若用差商的值來近似點\((x_{n+1},y_{n+1})\)處的導數值，則可以得到隱式歐拉格式
\[\left\{\begin{aligned} y_{n+1}&=y_n+hf(x_{n+1},y_{n+1})\\ y(x_0)&=y_0 \end{aligned}\right.\quad n=0,1,2,\cdots \]
隱式歐拉格式只有在\(f(x,y)\)容易將兩個參數分離的情況，才比較容易使用。
改進的歐拉格式：
為了得到更加准確的估計，可以考慮將歐拉格式和隱式歐拉格式取平均值，得到
\[y_{n+1}=y_n+\frac{h}{2}(f(x_n,y_n)+f(x_{n+1},y_{n+1}))\tag{10-2-1} \]
上述式子稱為梯形公式，為了計算\(y_{n+1}\)常用以下迭代式
\[\left\{\begin{aligned} y_{n+1}^{(0)}&=y_n+hf(x_n,y_n)\\ y_{n+1}^{(k+1)}&=y_n+\frac{h}{2}(f(x_n,y_n)+f(x_{n+1},y_{n+1}^{(k)})) \end{aligned}\right.\quad k=0,1,2,\cdots\tag{10-2-2} \]
當\(|y_{n+1}^{(k+1)}-y_{n+1}^{(k)}|<\varepsilon\)時，取\(y_{n+1}\approx y_{n+1}^{(k+1)}\)。
定理 設函數\(f(x,y)\)對變量\(y\)滿足Lipschtz條件，\(L\)為Lipschtz常數。若步長\(h\)滿足\(0\le\frac{hL}{2}\lt1\)，由迭代式(10-2-1)產生的序列\(\{y_{n+1}^{(k)}\}\)收斂到\(y_{n+1}\)。
證明：由(10-2-1)和(10-2-2)有
\[\begin{aligned} |y_{n+1}-y_{n+1}^{(k+1)}|&=\frac{h}{2}|f(x_{n+1},y_{n+1})-f(x_{n+1},y_{n+1}^{(k)})|\\ &\le\frac{hL}{2}|y_{n+1}-y_{n+1}^{(k)}|\le\cdots\le \left(\frac{hL}{2}\right)^{k+1}|y_{n+1}-y_{n+1}^{(0)}| \end{aligned} \]
於是
\[\lim_{k\rightarrow\infty}\left(\frac{hL}{2}\right)^{k+1}=0 \Rightarrow \lim_{k\rightarrow\infty}|y_{n+1}-y_{n+1}^{(k+1)}|=0 \Rightarrow \lim_{k\rightarrow\infty}y_{n+1}^{(k+1)}=y_{n+1} \]
由於采用迭代方法計算\(y_{n+1}\)時的計算量較大，因此一般只迭代一次，構成一類預估-校正算法：
\[\left\{\begin{aligned} y^P_{n+1}&=y_n+hf(x_n,y_n)&\text{(預估)}\\ y^C_{n+1}&=y_n+\frac{h}{2}(f(x_n,y_n)+f(x_{n+1},y_{n+1}^P))&\text{(校正)} \end{aligned}\right. \]
取\(y_{n+1}=y_{n+1}^C\)，將上述式子改寫成
\[\left\{\begin{aligned} k_1&=hf(x_n,y_n)\\ k_2&=hf(x_n+h,y_n+k_1)\\ y_{n+1}&=y_n+\frac{1}{2}(k_1+k_2) \end{aligned}\right. \]
10.3 誤差、數值穩定性、收斂性
10.3.1 誤差概述
顯式單步法一般形式為
\[y_{n+1}=y_n+h\varphi(x_n,y_n,h) \]
而隱式單步法的一般形式為
\[y_{n+1}=y_n+h\varphi(x_n,y_n,x_{n+1},y_{n+1},h) \]
其中\(\varphi(\cdot)\)稱為增量函數，並且與\(f(x,y)\)相關。歐拉格式、隱式歐拉格式就是一個典型的實例。
需要注意的是，歐拉格式的基本原理就是利用了一階差商近似了一個節點上的導數值，這種近似本身就帶來了一定的誤差。也就是說，單步法本身就具有一定的誤差。
定義1 從初值\(y(x_0)=y_0\)出發，由單步法顯式、或隱式逐步計算，可以得到\(x_{n+1}\)點出的近似值\(y_{n+1}\)。此時，我們稱\(e_{n+1}=y(x_{n+1})-y_{n+1}\)為在點\(x_{n+1}\)上的整體截斷誤差。
定義2 如果第\(n\)步在\(x_n\)的計算沒有誤差，即\(y(x_n)=y_n\)。且由單步法計算出來的近似值\(\widetilde{y}_{n+1}\)，此時，我們稱\(T_{n+1}=y(x_{n+1})-\widetilde{y}_{n+1}\)為點\(x_{n+1}\)上的局部截斷誤差。
局部階段誤差的實例（以顯式單步法為例）：
\[\widetilde{y}_{n+1}=y(x_n)+h\varphi(x_n,y(x_n),h)\\ T_{n+1}=y(x_{n+1})-\widetilde{y}_{n+1} \]
定義3 若一種數值方法的局部截斷誤差是\(O(h^{p+1})\)，即\(T_{n+1}=O(h^{p+1})\)。則稱相應的數值方法是\(p\)階方法，並稱該算法具有\(p\)階精度，其中\(p\)是正整數。
一般來說，一種數值方法的局部截斷誤差階越大，那么相對精度就越高。
歐拉格式的局部截斷誤差：
對於歐拉格式，假設\(y_k=y(x_k)\)，那么
\[y_{k+1}=y(x_k)+hf(x_k,y(x_k))=y(x_k)+hy'(x_k) \]
對理論解\(y(x)\)在點\(x_k\)進行二階泰勒展開，並令\(x=x_{k+1}\)，則有
\[y(x_{k+1})=y(x_k)+hy'(x_k)+\frac{h^2}{2}y''(\xi)\quad\xi\in(x_k,x_{k+1}) \]
於是
\[T_{k+1}=y(x_{k+1})-y_{k+1}=\frac{h^2}{2}y''(\xi)\Rightarrow T_{k+1}=O(h^2) \]
即顯式歐拉格式僅具有一階精度。
梯形公式的局部截斷誤差：
假設\(y_k=y(x_k)\)，那么
\[y_{k+1}=y(x_k)+\frac{h}{2}[y'(x_k)+y'(x_{k+1})] \]
利用泰勒展開可以得到
\[y(x_{k+1})=y(x_k)+hy'(x_k)+\frac{h^2}{2}f''(x_k)+\frac{h^3}{6}f'''(x_k)+O(h^4)\\ y'(x_{k+1})=y'(x_k)+hy''(x_k)+\frac{h^2}{2}f'''(x_k)+O(h^3) \]
綜合上述式子，可以得到
\[\begin{aligned} T_{k+1}&=y(x_{k+1})-y_{k+1}\\ &=-\frac{h^3}{12}f'''(x_k)+O(h^4) \end{aligned}\Rightarrow T_{k+1}=O(h^3) \]
即梯形公式的局部截斷誤差為\(O(h^3)\)。
歐拉方法的整體截斷誤差：
下面考察使用\(n+1\)次歐拉公式所產生的累計誤差，即整體截斷誤差。
首先，記\(\bar{y}_{n+1}=y(x_n)+hf(x_n,y(x_n)),\varepsilon_{n+1}=y(x_{n+1})-y_{n+1}\)，那么
\[\begin{aligned} |\varepsilon_{n+1}|&=|y(x_{n+1})-y_{n+1}|\\ &=|y(x_{n+1})-\bar{y}_{n+1}+\bar{y}_{n+1}-y_{n+1}|\\ &\le|T_{n+1}|+|\bar{y}_{n+1}-y_{n+1}| \end{aligned} \]
其中
\[\begin{aligned} |\bar{y}_{n+1}-y_{n+1}|&=|[y(x_n)+hf(x_n,y(x_n))]-[y_n+hf(x_n,y_n)]|\\ &=|[y(x_n)-y_n]+h[f(x_n,y(x_n))-f(x_n,y_n)]|\\ &\le|y(x_n)-y_n|+h|f(x_n,y(x_n))-f(x_n,y_n)|\\ &\le|y(x_n)-y_n|+hL|y(x_n)-y_n|\\ &=(1+hL)|y(x_n)-y_n|=(1+hL)|\varepsilon_n| \end{aligned} \]
其中\(L\)是Lipschtz常數。
根據之前的局部截斷誤差的推導，可知
\[|T_{n+1}|\le\frac{Mh^2}{2},\quad M=\max_{x\in[a,b]}|y''(x)| \]
於是
\[|\varepsilon_{n+1}|\le\frac{Mh^2}{2}+(1+hL)|\varepsilon_n| \]
反復通過上述式子遞推，最終可以得到
\[\begin{aligned} |\varepsilon_{n+1}|&\le\frac{Mh^2}{2}+\frac{Mh^2}{2}(1+hL)+\cdots+\frac{Mh^2}{2}(1+hL)^n +(1+hL)^{n+1}|\varepsilon_0|\\ &=\frac{Mh^2}{2}\sum_{k=0}^n(1+hL)^k+(1+hL)^{n+1}|\varepsilon_0|\\ &=\frac{Mh}{2L}((1+hL)^{n+1}-1)+(1+hL)^{n+1}|\varepsilon_0| \end{aligned} \]
由於初值是准確的，即\(\varepsilon_0=0\)，那么
\[|\varepsilon_n|\le\frac{Mh}{2L}((1+hL)^{n}-1) \]
根據不等式\(e^x\ge1+x(x\ge0)\)，可以知道\(e^{hL}\ge1+hL\)，於是\(e^{nhL}\ge(1+hL)^n\)。
假設\(b\ge nh=x_n-x_0\)，那么\(e^{bL}\ge(1+hL)^n\)，於是
\[|\varepsilon_n|\le\frac{Mh}{2L}(e^{bL}-1) \]
所以，顯式歐拉格式的整體截斷誤差為\(\varepsilon_n=O(h)\)。可以發現此時全局截斷誤差的階恰好比局部截斷誤差的階降一級。
10.3.2 數值穩定性
用數值方法求解初值問題時，舍入誤差是不可避免地。穩定性就是研究舍入誤差傳播問題，在求解過程中，如果舍入誤差不增長，則稱算法是數值穩定的。一般地討論一個數值方法的穩定性非常復雜，因此一般是把數值方法應用於實驗方程
\[y'=\lambda y\tag{10-3-1} \]
其中\(\lambda\in\mathbb{C},Re(\lambda)<0\)。
當然，如果一個數值方法對於實驗方程是穩定的，並不能保證該方法但對於任何方程都是穩定的。但如果該數值方法連實驗方程都不穩定，那么也很難應用於其他方程的求解。
定義4 用一個數值方法求解實驗方程(10-3-1)，假設對於給定的步長\(h\)，計算\(y_n\)時引起的誤差為\(\varepsilon_n\)。若這個誤差在計算后面的\(y_{n+k}(k=1,2,\cdots)\)所引起的誤差為\(\varepsilon_{n+k}\)滿足：
\[|\varepsilon_{n+k}|\le|\varepsilon_n| \]
則稱這個數值方法對步長\(h\)和復數\(\lambda\)是絕對穩定的，使得數值方法絕對穩定的\(\bar{h}=h\lambda\)在復平面山的允許范圍稱為數值方法的絕對穩定域。
定義5 若某數值算法的絕對穩定域包含\(h\lambda\)平面上的左邊平面（即\(Re(h\lambda)<0\)），則稱該算法是A穩定的。
歐拉方法的穩定性
將歐拉方法應用到實驗方程中，可以得到
\[y_{n+1}=y_n+h(x_n,y_n)=(1+h\lambda)y_n \]
由於舍入誤差的存在，實際上應該把\(y_n\)寫成\(\widetilde{y}_n\)，因此
\[\widetilde{y}_{n+1}=(1+h\lambda)\widetilde{y}_n \]
令誤差\(\varepsilon_n=y_n-\widetilde{y}_n\)，那么
\[\begin{aligned} \varepsilon_{n+1}&=(1+h\lambda)\varepsilon_n\\ &=(1+h\lambda)^2\varepsilon_{n-1}\\ &=\cdots\\ &=(1+h\lambda)^{n+1}\varepsilon_0 \end{aligned} \]
於是
\[|\varepsilon_{n}|\le|1+h\lambda|^n|\varepsilon_0| \]
通過上述式子可以發現，只要\(|1+h\lambda|\le1\)，那么歐拉方法就是絕對穩定的。我們稱滿足\(|1+h\lambda|\le1\)的\(h\lambda\)區域為絕對穩定區域。可以發現，該絕對穩定區域實際上是復平面上的一個圓。
隱式歐拉方法的穩定性
與上述分析方法相同，可以得到
\[|\varepsilon_n|\le\left|\frac{1}{1-h\lambda}\right|^n|\varepsilon_0| \]
因為實驗方程要求\(Re(\lambda)<0\)，所以\(\left|\frac{1}{1-h\lambda}\right|<1\)是恆成立的。這也意味着隱式歐拉方法是無條件絕對穩定的，自然也是A穩定的。
梯形公式（改進歐拉格式）的穩定性
同樣可以得到
\[|\varepsilon_n|\le\left|\frac{1+\frac{h\lambda}{2}}{1-\frac{h\lambda}{2}}\right|^n|\varepsilon_0| \]
令\(\frac{h\lambda}{2}=a+bi\)，那么
\[r=\left|\frac{1+\frac{h\lambda}{2}}{1-\frac{h\lambda}{2}}\right|= \left|\frac{1+a+bi}{1-a-bi}\right|= \sqrt{\frac{(1+a)^2+b^2}{(1-a)^2+b^2}}= \sqrt{1+\frac{4a}{(1-a)^2+b^2}} \]
因為\(Re(\lambda)<0\)，所以\(a=Re(h\lambda)<0,h>0\)。因此\(r<1\)，這對於實驗方程是恆成立的。這意味着，梯形公式是絕對穩定的，並且還是A穩定的。
10.3.3 收斂性
常微分方程初值問題的求解，是將微分方程轉換成差分方程來求解，並用計算值\(y_n\)來近似代替\(y(x_n)\)。這種近似是否合理，需要看當分割區間\([x_n,x_{n+1}]\)的長度\(h\)越來越小時，即\(h=x_{n+1}-x_n\rightarrow0\)時，\(y_n\rightarrow y(x_n)\)是否成立。如果成立，則稱該方法是收斂的，否則是不收斂的。
以歐拉方法為例，結合之前在全局截斷誤差中的推導可知
\[\lim_{h\rightarrow0}|\varepsilon_n|=\lim_{h\rightarrow0}\frac{M}{2L}((1+hL)^{n}-1)h=0 \]
這說明歐拉方法是收斂的，並且具有一階收斂速度。
10.4 龍格-庫塔（Runge-Kutta）方法
歐拉公式：
\[\left\{\begin{aligned} K_1=f(x_n,y_n)\\ y_{n+1}=y_n+hK_1 \end{aligned}\right. \]
改進的歐拉公式：
\[\left\{\begin{aligned} K_1&=f(x_n,y_n)\\ K_2&=f(x_n+h,y_n+K_1)\\ y_{n+1}&=y_n+\frac{h}{2}(K_1+K_2) \end{aligned}\right. \]
可以發現上述兩組公式都使用了函數\(f(x,y)\)在某些點上的函數值的線性組合來計算\(y(x_{n+1})\)的近似值\(y_{n+1}\)。歐拉公式每一步需要計算一次\(f(x,y)\)的值，稱為一階方法。改進的歐拉公式則需要計算兩次\(f(x,y)\)的值，稱為二階方法。
於是可以考慮使用函數\(f(x,y)\)在若干個點上的函數值的線性組合來構造近似公式，構造是要求近似公式在\((x_n,y_n)\)處的泰勒展開與理論解\(y(x)\)在\(x_n\)處的泰勒展開前面幾項重合，從而提高近似公式的階數。這樣一來，既能夠避免求偏導，又能夠提高計算方法的精度階數。
簡言之，這種方法就是在區間\([x_i,x_{i+1}]\)上多預報幾個點的斜率值，然后將其加權平均作為平均斜率，以此構造更高精度的計算格式，這就是龍格-庫塔方法的基本思想。
一般地，在RK方法中，將近似公式設為
\[\left\{\begin{aligned} y_{n+1}&=y_n+h\sum_{i=1}^pc_iK_i\\ K_1&=f(x_n,y_n)\\ K_i&=f(x_n+a_ih,y_n+h\sum_{j=1}^{i-1}b_{ij}K_j)\;(i=2,3,\cdots,p)\\ \end{aligned}\right. \]
其中\(a_i,b_{ij},c_i\)都是待定參數，確定它們的原則是使近似公式在\((x_n,y_n)\)處的泰勒展開式與解\(y(x)\)在\(x_n\)處的泰勒展開式的前幾項盡可能多地重合。
二階方法： 當\(p=2\)時，近似公式設為
\[\left\{\begin{aligned} y_{n+1}&=y_n+h(c_1K_1+c_2K_2)\\ K_1&=f(x_n,y_n)\\ K_2&=f(x_n+a_2h,y_n+hb_{21}K_1) \end{aligned}\right. \]
在點\((x_n,y_n)\)處對\(f(x,y)\)應用泰勒展開
\[\begin{aligned} y_{n+1}&=y_n+h(c_1f(x_n,y_n)+c_2f(x_n+a_2h,y_n+hb_{21}K_1))\\ &=y_n+h\left(\color{red}{c_1f(x_n,y_n)}+\color{blue}{c_2\left[f(x_n,y_n)+a_2hf'_x(x_n,y_n)+hb_{21}f'_y(x_n,y_n)f(x_n,y_n)+O(h^2)\right]}\right)\\ &=y_n+(c_1+c_2)f(x_n,y_n)h+c_2[a_2f'_x(x_n,y_n)+b_{21}f'_y(x_n,y_n)f(x_n,y_n)]h^2+O(h^3)\\ \end{aligned} \]
\(y(x_{n+1})\)在點\((x_n,y_n)\)處的泰勒展開為（這里需要注意\(y'=f(x,y)\)）
\[\begin{aligned} y(x_{n+1})&=y(x_n)+hy'(x_n)+\frac{h^2}{2}y''(x_n)+O(h^3)\\ &=y_n+hf(x_n,y_n)+\frac{h^2}{2}(f'_x(x_n,y_n)+f'_y(x_n,y_n)f(x_n,y_n))+O(h^3)\\ \end{aligned} \]
對比上下兩式的泰勒展開的系數，可得
\[c_1+c_2=1\\ c_2a_2=\frac12\\ c_2b_{21}=\frac12 \]
顯然三個方程有四個未知數，具有無窮多組解。
如果取\(c_1=c_2=1/2,a_2=b_{21}=1\)，即可得到改進的歐拉公式。
如果取\(c_1=0,c_2=1,a_2=b_{21}=1/2\)，就得到了常用的二階RK公式（也稱為中點公式）：
\[\left\{\begin{aligned} y_{n+1}&=y_n+hK_2\\ K_1&=f(x_n,y_n)\\ K_2&=f(x_n+h/2,y_n+hK_1/2) \end{aligned}\right. \]
三階方法： 當\(p=3\)時，需要計算三個K值，此時近似公式為
\[\left\{\begin{aligned} y_{n+1}&=y_n+h(c_1K_1+c_2K_2+c_3K_3)\\ K_1&=f(x_n,y_n)\\ K_2&=f(x_n+a_2h,y_n+hb_{21}K_1)\\ K_3&=f(x_n+a_3h,y_n+hb_{31}K_1+hb_{32}K_2) \end{aligned}\right. \]
仍然按照上述分析方法，此時要求泰勒展開式的前四項相同，得到
\[\left\{\begin{aligned} &c_1+c_2+c_3=1\\ &a_2=b_{21}\\ &a_3=b_{31}+b_{32}\\ &c_2a_2+c_3a_3=\frac12\\ &c_2a_2^2+c_3a_3^2=\frac13\\ &c_3a_2b_{32}=\frac16 \end{aligned}\right. \]
令
\[c_1=\frac16,c_2=\frac46,c_3=\frac16\\ a_2=\frac12,a_3=1\\ b_{21}=\frac12,b_{31}=-1,b_{32}=2 \]
就得到了常用的三階RK公式：
\[\left\{\begin{aligned} y_{n+1}&=y_n+\frac{h}6(K_1+4K_2+K_3)\\ K_1&=f(x_n,y_n)\\ K_2&=f(x_n+\frac{h}2,y_n+\frac{h}2K_1)\\ K_3&=f(x_n+h,y_n-hK_1+2hK_2) \end{aligned}\right. \]
四階方法： 在實際應用中，最常用的時四階龍格-庫塔公式。但是該方法的推導十分繁瑣，需要求解含有13個參數的11個方程組成的方程組。因此，這里就不加以推導地給出計算公式。
\[\left\{\begin{aligned} y_{n+1}&=y_n+\frac{h}6(K_1+2K_2+2K_3+K_4)\\ K_1&=f(x_n,y_n)\\ K_2&=f(x_n+\frac{h}2,y_n+\frac{h}2K_1)\\ K_3&=f(x_n+\frac{h}2,y_n+\frac{h}2K_2)\\ K_4&=f(x_n+h,y_n+hK_3) \end{aligned}\right. \]
其局部截斷誤差為\(O(h^5)\)。
四階RK方法的穩定性：
將四階RK方法應用到實驗方程中可以得到
\[y_{n+1}=[1+(h\lambda)+\frac{(h\lambda)^2}{2!}+\frac{(h\lambda)^3}{3!}+\frac{(h\lambda)^4}{4!}]y_n \]
於是穩定區域由下式確定
\[\left|1+(h\lambda)+\frac{(h\lambda)^2}{2!}+\frac{(h\lambda)^3}{3!}+\frac{(h\lambda)^4}{4!}\right|\le1 \]
當\(\lambda<0\)時，可以解得：\(-2.78<h\lambda<0\)。因此，當\(|\lambda|\)的值比較大時，只有當\(h\)比較小的時候才能夠保證算法的穩定性。
一階方程組：
四階RK方法可以拓展到多個函數的情況，假設有
\[\left\{\begin{aligned} y'_1&=f_1(x,y_1,y_2,\cdots,y_n)\\ y'_2&=f_2(x,y_1,y_2,\cdots,y_n)\\ &\vdots\\ y'_n&=f_n(x,y_1,y_2,\cdots,y_n) \end{aligned}\right. \]
並且有初值條件
\[\left\{\begin{aligned} y_1(x_0)&=y_{10}\\ y_2(x_0)&=y_{20}\\ &\vdots\\ y_n(x_0)&=y_{n0} \end{aligned}\right. \]
為了簡化記號，可以先將其上面的條件改寫成向量的形式：
\[\left\{\begin{aligned} &y'=f(x,y)\\ &y(x_0)=y_0 \end{aligned}\right. \]
其中
\[y'=\begin{bmatrix} y_1'\\y_2'\\\vdots\\y_n' \end{bmatrix},y=\begin{bmatrix} y_1\\y_2\\\vdots\\y_n \end{bmatrix}\\ f(x,y)=\begin{bmatrix} f_1\\f_2\\\vdots\\f_n \end{bmatrix}(x,y)=\begin{bmatrix} f_1(x,y)\\f_2(x,y)\\\vdots\\f_n(x,y) \end{bmatrix}=\begin{bmatrix} f_1(x,y_1,\cdots,y_n)\\f_2(x,y_1,\cdots,y_n)\\\vdots\\f_n(x,y_1,\cdots,y_n) \end{bmatrix}\\ y_0=y(x_0)=\begin{bmatrix} y_1\\y_2\\\vdots\\y_n \end{bmatrix}(x_0)=\begin{bmatrix} y_1(x_0)\\y_2(x_0)\\\vdots\\y_n(x_0) \end{bmatrix} \]
有了上述向量形式的記號，於是
\[K_1=f(x_k,y_k)=\begin{bmatrix} f_1\\f_2\\\vdots\\f_n \end{bmatrix}(x_k,y_k)\\ K_2=f(x_k+\frac{h}2,y_k+\frac{h}2K_1)\\ K_3=f(x_k+\frac{h}2,y_k+\frac{h}2K_2)\\ K_4=f(x_k+h,y_k+hK_3) \]
於是形式上就和四階方法相同
\[\left\{\begin{aligned} y_{k+1}&=y_k+\frac{h}6(K_1+2K_2+2K_3+K_4)\\ K_1&=f(x_k,y_k)\\ K_2&=f(x_k+\frac{h}2,y_k+\frac{h}2K_1)\\ K_3&=f(x_k+\frac{h}2,y_k+\frac{h}2K_2)\\ K_4&=f(x_k+h,y_k+hK_3) \end{aligned}\right. \]
但需要注意的是，需要將上述的變量當成向量來對待。
高階方程化為一階方程組：
對於高階微分方程的初值問題，原則上總可以歸結成一階方程組來求解。例如，對於如下m階微分方程
\[y^{(m)}=f(x,y,y',\cdots,y^{(m-1)}) \]
初始條件為
\[y(x_0)=y_0,y'(x_0)=y'_0,\cdots,y^{(m-1)}(x_0)=y_0^{(m-1)} \]
引入記號
\[y_1=y,y_2=y',\cdots,y_m=y^{(m-1)} \]
可以得到
\[\left\{\begin{aligned} &y_1'=y_2\\&y_2'=y_3\\&\cdots\\&y_{m-1}'=y_m\\&y'_m=f(x,y_1,y_2,\cdots,y_m) \end{aligned}\right. \]
然后使用一階方程組的四階RK方法即可。
