研究生课程《应用数值分析》结课了，使用python简单记录了下常微分方程数值解法。

向前欧拉法

{ y i + 1 = y i + h i f ( x i , y i ) y 0 = y ( a ) \left \{ \begin{array}{lr} y_{i+1}=y_i+h_i f(x_i,y_i) \\ y_0=y(a) \end{array} \right .{yi+1​=yi​+hi​f(xi​,yi​)y0​=y(a)​

from pylab import *

import warnings

warnings.filterwarnings('ignore')

求解初值问题

{ y ′ = x − y + 1 0 ≤ x ≤ 1 y ( 0 ) = 1 \left \{ \begin{array}{lr} y'=x-y+1 & 0\leq x \leq 1 \\ y(0)=1 \end{array} \right .{y′=x−y+1y(0)=1​0≤x≤1

该微分方程的精确解为：y = x + e x p ( − x ) y=x+exp(-x)y=x+exp(−x)

def f(t,y):

'''

求解的微分方程，

'''

return t-y+1

def euler_forward(f,a=0,b=1,ya=1,h=0.1,verbose=True):

'''向前欧拉法

Args

----------

f： callable function

需要求解的函数

a: float

求解区间起始值

b：float

求解区间终止值

ya：float

起始条件，ya=y(a)

h：float

求解步长(区间[a,b]n等分)

verbose：logical，default is True

显示迭代结果

Returns

----------

res：list like

返回向前欧拉发求解的结果

'''

# i = 0

res = []

xi = a

yi = ya

while xi<=b: # 在求解区间范围

y = yi + h*f(xi,yi)

if verbose:

print('xi:{:.2f}, yi:{:.6f}'.format(xi,yi))

res.append(y)

xi, yi = xi+h, y

return res

res = euler_forward(f,a=0,b=1,ya=1,h=0.1,verbose=True)

xi:0.00, yi:1.000000

xi:0.10, yi:1.000000

xi:0.20, yi:1.010000

xi:0.30, yi:1.029000

xi:0.40, yi:1.056100

xi:0.50, yi:1.090490

xi:0.60, yi:1.131441

xi:0.70, yi:1.178297

xi:0.80, yi:1.230467

xi:0.90, yi:1.287420

xi:1.00, yi:1.348678

龙格库塔(Runge-Kutta methods)

三阶龙格库塔法

{ y i + 1 = y i + 1 6 ( k 1 + 2 k 2 + k 3 ) k 1 = h f ( x i , y i ) k 2 = h f ( x i + 1 2 h , y i + 1 2 k 1 ) k 3 = h f ( x i + h , y i − k 1 + 2 k 2 ) y 0 = y ( a ) , i = 0 , 1 , ⋯ &ThinSpace; , n − 1 \left \{ \begin{array}{lr} y_{i+1}=y_i+\frac{1}{6}(k_1+2k_2+k_3) & \\ k_1 = hf(x_i,y_i) & \\ k_2 = hf(x_i+\frac{1}{2}h,y_i+\frac{1}{2}k_1) & \\ k_3 = hf(x_i+h,y_i-k_1+2k_2) & \\ y_0=y(a),i=0,1,\cdots,n-1 \\ \end{array} \right .⎩⎪⎪⎪⎪⎨⎪⎪⎪⎪⎧​yi+1​=yi​+61​(k1​+2k2​+k3​)k1​=hf(xi​,yi​)k2​=hf(xi​+21​h,yi​+21​k1​)k3​=hf(xi​+h,yi​−k1​+2k2​)y0​=y(a),i=0,1,⋯,n−1​​

四阶龙格库塔法

{ y i + 1 = y i + 1 6 ( k 1 + 2 k 2 + 2 k 3 + k 4 ) k 1 = h f ( x i , y i ) k 2 = h f ( x i + 1 2 h , y i + 1 2 k 1 ) k 3 = h f ( x i + 1 2 h , y i + 1 2 k 2 ) k 4 = h f ( x i + h , y i + k 3 ) y 0 = y ( a ) , i = 0 , 1 , ⋯ &ThinSpace; , n − 1 \left \{ \begin{array}{lr} y_{i+1}=y_i+\frac{1}{6}(k_1+2k_2+2k_3+k_4) & \\ k_1 = hf(x_i,y_i) & \\ k_2 = hf(x_i+\frac{1}{2}h,y_i+\frac{1}{2}k_1) & \\ k_3 = hf(x_i+\frac{1}{2}h,y_i+\frac{1}{2}k_2) & \\ k_4 = hf(x_i+h,y_i+k_3) & \\ y_0=y(a),i=0,1,\cdots,n-1 \\ \end{array} \right .⎩⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎧​yi+1​=yi​+61​(k1​+2k2​+2k3​+k4​)k1​=hf(xi​,yi​)k2​=hf(xi​+21​h,yi​+21​k1​)k3​=hf(xi​+21​h,yi​+21​k2​)k4​=hf(xi​+h,yi​+k3​)y0​=y(a),i=0,1,⋯,n−1​​

def runge_kutta(f,a=0,b=1,ya=1,h=0.1,verbose=True):

'''四阶龙格库塔法

Args

----------

f： callable function

需要求解的函数

a: float

求解区间起始值

b：float

求解区间终止值

ya：float

起始条件，ya=y(a)

h：float

求解步长(区间[a,b]n等分)

verbose：logical，default is True

显示迭代结果

Returns

----------

res：list like

返回向前欧拉发求解的结果

'''

res = []

xi = a

yi = ya

while xi <= b: # 在求解区间范围

k1 = h * f(xi, yi)

k2 = h * f(xi + h/2, yi + k1/2)

k3 = h * f(xi + h/2, yi + k2/2)

k4 = h * f(xi + h, yi + k3)

y = yi + 1/6 * (k1 + 2*k2 + 2*k3 + k4)

if verbose:

print('xi:{:.2f}, yi:{:.10f}'.format(xi,yi))

res.append(y)

xi, yi = xi+h, y

return res

res = runge_kutta(f,a=0,b=1,ya=1,h=0.1,verbose=True)

xi:0.00, yi:1.0000000000

xi:0.10, yi:1.0048375000

xi:0.20, yi:1.0187309014

xi:0.30, yi:1.0408184220

xi:0.40, yi:1.0703202889

xi:0.50, yi:1.1065309344

xi:0.60, yi:1.1488119344

xi:0.70, yi:1.1965856187

xi:0.80, yi:1.2493292897

xi:0.90, yi:1.3065699912

xi:1.00, yi:1.3678797744

改进欧拉法

{ y p = y i + h f ( x i , y i ) y i + 1 = y i + h 2 [ f ( x i , y i ) + f ( x i , y p ) ] i = 0 , 1 , ⋯ &ThinSpace; , n − 1 y 0 = y ( a ) \left \{ \begin{array}{lr} y_p=y_i+hf(x_i,y_i) & \\ y_{i+1} = y_i+\frac{h}{2}[f(x_i,y_i)+f(x_i,y_p)] &i=0,1,\cdots ,n-1 \\ y_0=y(a) \\ \end{array} \right .⎩⎨⎧​yp​=yi​+hf(xi​,yi​)yi+1​=yi​+2h​[f(xi​,yi​)+f(xi​,yp​)]y0​=y(a)​i=0,1,⋯,n−1​

求解初值问题

{ y ′ = − y ( 0 ≤ x ≤ 1 ) y ( 0 ) = 1 \left \{ \begin{array}{lr} y'=-y &(0 \leq x \leq 1) \\ y(0)=1 \\ \end{array} \right .{y′=−yy(0)=1​(0≤x≤1)

精确解为y = e x p ( − x ) y=exp(-x)y=exp(−x)

def f(t,y):

'''

精确解为y=exp(-x)

'''

return -y

def improved_euler(f,a=0,b=1,ya=1,h=0.1,verbose=True):

'''改进欧拉法

Args

----------

f： callable function

需要求解的函数

a: float

求解区间起始值

b：float

求解区间终止值

ya：float

起始条件，ya=y(a)

h：float

求解步长(区间[a,b]n等分)

verbose：logical，default is True

显示迭代结果

Returns

----------

res：list like

返回向前欧拉发求解的结果

'''

res = []

xi = a

yi = ya

while xi <= b: # 在求解区间范围

yp = yi + h*f(xi, yi)

y = yi + h/2 * (f(xi, yi) + f(xi, yp))

if verbose:

print('xi:{:.2f}, yi:{:.6f}'.format(xi,yi))

res.append(y)

xi, yi = xi+h, y

return res

res = improved_euler(f,a=0,b=1,ya=1,h=0.1,verbose=True)

xi:0.00, yi:1.000000

xi:0.10, yi:0.905000

xi:0.20, yi:0.819025

xi:0.30, yi:0.741218

xi:0.40, yi:0.670802

xi:0.50, yi:0.607076

xi:0.60, yi:0.549404

xi:0.70, yi:0.497210

xi:0.80, yi:0.449975

xi:0.90, yi:0.407228

xi:1.00, yi:0.368541