研究生課程《應(yīng)用數(shù)值分析》結(jié)課了，使用python簡(jiǎn)單記錄了下常微分方程數(shù)值解法。

向前歐拉法

{ 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

需要求解的函數(shù)

a: float

求解區(qū)間起始值

b：float

求解區(qū)間終止值

ya：float

起始條件，ya=y(a)

h：float

求解步長(zhǎng)(區(qū)間[a,b]n等分)

verbose：logical，default is True

顯示迭代結(jié)果

Returns

----------

res：list like

返回向前歐拉發(fā)求解的結(jié)果

'''

# i = 0

res = []

xi = a

yi = ya

while xi<=b: # 在求解區(qū)間范圍

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

需要求解的函數(shù)

a: float

求解區(qū)間起始值

b：float

求解區(qū)間終止值

ya：float

起始條件，ya=y(a)

h：float

求解步長(zhǎng)(區(qū)間[a,b]n等分)

verbose：logical，default is True

顯示迭代結(jié)果

Returns

----------

res：list like

返回向前歐拉發(fā)求解的結(jié)果

'''

res = []

xi = a

yi = ya

while xi <= b: # 在求解區(qū)間范圍

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

改進(jìn)歐拉法

{ 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):

'''改進(jìn)歐拉法

Args

----------

f： callable function

需要求解的函數(shù)

a: float

求解區(qū)間起始值

b：float

求解區(qū)間終止值

ya：float

起始條件，ya=y(a)

h：float

求解步長(zhǎng)(區(qū)間[a,b]n等分)

verbose：logical，default is True

顯示迭代結(jié)果

Returns

----------

res：list like

返回向前歐拉發(fā)求解的結(jié)果

'''

res = []

xi = a

yi = ya

while xi <= b: # 在求解區(qū)間范圍

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

總結(jié)

以上是生活随笔為你收集整理的python 常微分方程_常微分方程数值解法——python实现的全部?jī)?nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。