小跳最近在搭建一個數值模擬環境，由於需要用到python裡面的一些庫，所以不得不把simulink的模型搬過來，我們都知道在simulink裡，模擬的時候設定模擬步長和微分方程求解器是必要的步驟。但是為什麼要設定這個小跳卻早已忘記了。

一年級的時候搬磚搬多了，數分課也沒好好上，回頭一看，這麼簡單的東西，當時竟然整的稀裡糊塗的。

為什麼要用RK4

先po一張圖，直觀感受一下模擬的誤差。

對於給定線性常微分方程
\[\dot x = x\]
易得，其解是
\[x(t) = Ce^t
\]
RK4是龍格庫塔法曲線，None是一階解法\(x(t+dt) = x(t)+\dot x dt\)
可以看到，線性常微分方程誤差尚且如此之大，那麼推廣到非線性微分方程，像這種形式

\[
\dot x = f(x,t) = tx^2 - \frac{x}{t}...
\]
那肯定誤差直接起飛了。解析解求起來也挺麻煩，這裡就不再引入分析了。

接下來把定義回顧一下，貼一下程式碼，有需自取，希望對大家有所幫助。

定義回顧

數值分析中，龍格－庫塔法（Runge-Kutta methods）是用於非線性常微分方程的解的重要的一類隱式或顯式迭代法。這些技術由數學家卡爾·龍格和馬丁·威爾海姆·庫塔於1900年左右發明。該方法主要是在已知方程導數和初值資訊，利用計算機模擬時應用，省去求解微分方程的複雜過程。

令初值問題表述如下。
\[
y' = f(t,y), y(t_0) = y_0

\]
則，對於該問題的RK4由如下方程給出：
\[
y_{n+1}=y_{n}+\frac{h}{6}\left(k_{1}+2 k_{2}+2 k_{3}+k_{4}\right) \\\]
其中
\[
\begin{matrix}
k_{1}=f\left(t_{n}, y_{n}\right) \\
k_{2}=f\left(t_{n}+\frac{h}{2}, y_{n}+\frac{h}{2} k_{1}\right) \\
k_{3}=f\left(t_{n}+\frac{h}{2}, y_{n}+\frac{h}{2} k_{2}\right) \\
k_{4}=f\left(t_{n}+h, y_{n}+h k_{3}\right)\end{matrix}
\]

式中，\(h\)為模擬步長，滿足\(h<\epsilon_1 \rightarrow error<\epsilon_2\)

程式碼實現

import numpy as np
import matplotlib.pyplot as plt
import sympy as sy

# 這裡介紹一個符號運算的方法，可以用來求解方程什麼的
def diff_eq(t,x):
    return sy.diff(x(t),t,1) - x(t)
    
t = sy.symbols('t')
x = sy.Function('x')
sy.pprint(sy.dsolve(diff_eq(t,x),x(t)))

def dot_x(t,x):
    return x

def rk4(f,t,x,h):
    k1 = f(t,x);
    k2 = f(t+0.5*h,x + 0.5*h*k1)
    k3 = f(t+0.5*h,x + 0.5*h*k2)
    k4 = f(t+h,x + h*k3)
    return h/6*(k1+2*k2+2*k3+k4)

t_list = np.arange(0,5,0.1);
#print(t)
x1_list = np.exp(t_list)
x2_list = []
x3_list = []
h = 0.1
x2 = 1;
x3 = 1;

for t in t_list:
#   print(t,idx)
    x2_list.append(x2)
    x3_list.append(x3)
    x2 = x2 + rk4(dot_x,t, x2, h)
    x3 = x3 + dot_x(t,x3) * h               
error_2 = x1_list - x2_list
error_3 = x1_list - x3_list

plt.figure()
plt.subplot(2,1,1)
plt.plot(t_list,x1_list, 'b-',label='Real')
plt.plot(t_list,x2_list,'r--', label = 'RK4')
plt.plot(t_list,x3_list,'g--', label = 'None')
plt.legend()

plt.subplot(2,1,2)
plt.plot(t_list,error_2, 'r--',label='Error_RK4')
plt.plot(t_list,error_3, 'g--',label='Error_none')
plt.legend()
plt.xlabel('Time(s)')
plt.show()

閒話

這裡推薦一個提高效率的工具Matplotlib cheat sheet

對於一個經常畫圖的科研狗來說，這張圖真是太太太太有必要了，因為時常遇到以下場景，不記得colormap名字，開啟文件查一番，不記得線寬關鍵詞，開啟文件查一番，不記得marker名字，開啟文件查一番。。。。。等等等等

所以，有了這張圖，在平常畫圖的時候中遇到的95%需要查文件的問題都可以在這張圖中找到答案。

這個速查表，可以關注微信公眾號“探物及理”後臺回覆“python畫圖”領取。