用Python實現(xiàn)四階龍格-庫塔（Runge-Kutta）方法求解高階微分方程
問題
應(yīng)用四階龍格-庫塔（Runge-Kutta）方法求解如下二階初值問題:
\begin{equation}
\left\{
\begin{aligned}
t^2x''(t)-2tx'(t)+2x(t) & = t^3\ln t, & t\in [1,5]\\
x(t) & = 1, & t=1 \\
x'(t) & = 0. & t=1
\end{aligned}
\right.
\end{equation}
要求:取步長h=0.01,給出解x(t)的圖像和在t=0,0.5,1,1.5,2,2.5,3,3.5,4,4.5,5處的近似解.
求解步驟
-
Step1. 將原問題歸結(jié)為其等價問題
引進(jìn)新的變量
y(t)=x'(t)將高階微分方程的初值問題歸結(jié)為如下一階微分方程組的初值問題來求解.\begin{equation} \left\{ \begin{aligned} x'(t) & = y, & t\in [1,5]\\ y'(t) & = \frac{t^3\ln t +2ty -2x}{t^2}, & t\in [1,5]\\ x(t) & = y, & t=1\\ y(t) & = 0. & t=1 \end{aligned} \right. \end{equation} -
Step2. 四階龍格-庫塔方法的離散格式
針對以上一階微分方程組的初值問題應(yīng)用四階龍格-庫塔方法構(gòu)造得到以下離散格式:
\begin{equation} \left\{ \begin{aligned} x_{n+1} & = x_n +\frac{h}{6}(K_1 + 2K_2 + 2K_3 + K_4),\\ y_{n+1} & = y_n +\frac{h}{6}(L_1 + 2L_2 + 2L_3 + L_4). \end{aligned} \right. \end{equation}其中
\begin{equation} \left\{ \begin{aligned} K_1 & = f(t_n,x_n,y_n),\\ K_2 & = f(t_n + \frac{h}{2},x_n + \frac{h}{2}K_1,y_n + \frac{h}{2}L_1),\\ K_3 & = f(t_n + \frac{h}{2},x_n + \frac{h}{2}K_2,y_n + \frac{h}{2}L_2),\\ K_4 & = f(t_n + h,x_n + hK_3,y_n + hL_3),\\ L_1 & = g(t_n,x_n,y_n),\\ L_2 & = g(t_n + \frac{h}{2},x_n + \frac{h}{2}K_1,y_n + \frac{h}{2}L_1),\\ L_3 & = g(t_n + \frac{h}{2},x_n + \frac{h}{2}K_2,y_n + \frac{h}{2}L_2),\\ L_4 & = f(t_n + h,x_n + hK_3,y_n + hL_3). \end{aligned} \right. \end{equation}這是一步法，利用節(jié)點
t_n上的值x_n,y_n，由 (4) 式順序計算K_1,L_1,K_2,L_2,K_3,L_3,K_4,L_4亥揖，然后代入(3)式即可求得節(jié)點t_{n+1}上的x_{n+1},y_{n+1}.- Step3. 利用龍格-庫塔法求解高階微分方程的
Python代碼如下:
# 開發(fā)者: Leo 劉 # 開發(fā)環(huán)境: macOs Big Sur # 開發(fā)時間: 2021/9/28 4:31 下午 # 郵箱 : 517093978@qq.com # @Software: PyCharm # ---------------------------------------------------------------------------------------------------------- import math import matplotlib.pyplot as plt def f(t, x, y): a = y return a def g(t, x, y): a = (t ** 3 * math.log(t) + 2 * t * y - 2 * x) / t ** 2 return a def RK4(t, x, y, h): """ :param t: t 的初始值 :param x: x 的初始值 :param y: y 的初始值 :param h: 時間步長 :return: 迭代新解 """ tarray, xarray, yarray = [], [], [] while t <= 5: tarray.append(t) xarray.append(x) yarray.append(y) t += h K_1 = f(t, x, y) L_1 = g(t, x, y) K_2 = f(t + h / 2, x + h / 2 * K_1, y + h / 2 * L_1) L_2 = g(t + h / 2, x + h / 2 * K_1, y + h / 2 * L_1) K_3 = f(t + h / 2, x + h / 2 * K_2, y + h / 2 * L_2) L_3 = g(t + h / 2, x + h / 2 * K_2, y + h / 2 * L_2) K_4 = f(t + h, x + h * K_3, y + h * L_3) L_4 = g(t + h, x + h * K_3, y + h * L_3) x = x + (K_1 + 2 * K_2 + 2 * K_3 + K_4) * h / 6 y = y + (L_1 + 2 * L_2 + 2 * L_3 + L_4) * h / 6 return tarray, xarray, yarray def main(): tarray, xarray, yarray = RK4(1, 1, 0, 0.01) print("龍格-庫塔 數(shù)值結(jié)果".center(168)) print('-' * 178) print("對象\\時刻\t", "t=0\t\t", " t=0.5\t\t", " t=1\t\t\t", " t=1.5\t\t", " t=2\t\t\t", " t=2.5\t\t\t", " t=3\t\t\t", " t=3.5\t\t", " t=4\t\t\t", " t=4.5\t\t\t", " t=5") print('-' * 178) print("x:", end='') for i in range(len(xarray)): if i % 40 == 0: print("\t\t", "%8.4f" % xarray[i], end='') print('\n', '-' * 177) print("y:", end='') for i in range(len(yarray)): if i % 40 == 0: print("\t\t", "%8.4f" % yarray[i], end='') print('\n', '-' * 177) plt.figure('龍格-庫塔 數(shù)值結(jié)果') plt.subplot(221) # plt.plot(tarray, xarray, label='x_runge_kutta') plt.scatter(tarray, xarray, label='x_scatter', s=1, c='#DC143C', alpha=0.6) # plt.xlabel('t') plt.legend() plt.subplot(222) # plt.plot(tarray, yarray, label='y_runge_kutta') plt.scatter(tarray, yarray, label='y_scatter', s=1, c='#DC143C', alpha=0.6) # plt.xlabel('t') plt.legend() plt.subplot(212) # plt.plot(tarray, xarray, label='runge_kutta') plt.scatter(tarray, xarray, label='Numerical solution scatter', s=1, c='#DC143C', alpha=0.6) plt.xlabel('t') plt.legend() plt.show() if __name__ == "__main__": main()- Step4. 代碼的運行結(jié)果如下：
龍格-庫塔 數(shù)值結(jié)果 ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- 對象\時刻 t=0 t=0.5 t=1 t=1.5 t=2 t=2.5 t=3 t=3.5 t=4 t=4.5 t=5 ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- x: 1.0000 0.8579 0.5080 0.1042 -0.1564 -0.0416 0.7105 2.3875 5.2995 9.7769 16.1685 --------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- y: 0.0000 -0.6689 -1.0161 -0.9205 -0.2855 0.9689 2.9116 5.6026 9.0952 13.4374 18.6728 --------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- - Step3. 利用龍格-庫塔法求解高階微分方程的
問題來源: 《微分方程數(shù)值解》-M.Ran
