研究生課程《應用數值分析》結課了，使用python簡單記錄了下常微分方程數值解法。
y_=y_i+h_i f(x_i,y_i) \\ y_0=y(a) \end \right .
y'=x-y+1 & 0\leq x \leq 1 \\ y(0)=1 \end \right .
, yi:'
.format
(xi,yi)
) 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
⎩⎪⎪⎪⎪⎨
⎪⎪⎪⎪
⎧yi
+1=
yi+
61(
k1+
2k2
+k3
)k1
=hf(
xi,
yi)
k2=
hf(x
i+2
1h,
yi+
21k
1)k
3=h
f(xi
+h,
yi−
k1+
2k2
)y0
=y(a
),i=
0,1,
⋯,n−
1 y_=y_i+\frac(k_1+2k_2+2k_3+k_4) & \\ k_1 = hf(x_i,y_i) & \\ k_2 = hf(x_i+\frach,y_i+\frack_1) & \\ k_3 = hf(x_i+\frach,y_i+\frack_2) & \\ k_4 = hf(x_i+h,y_i+k_3) & \\ y_0=y(a),i=0,1,\cdots,n-1 \\ \end \right .
⎩⎪⎪⎪⎪⎪
⎪⎨⎪⎪
⎪⎪⎪⎪
⎧yi
+1=
yi+
61(
k1+
2k2
+2k3
+k4
)k1
=hf
(xi
,yi
)k2
=hf(
xi+
21h
,yi
+21
k1)
k3=
hf(x
i+2
1h,
yi+
21k
2)k
4=h
f(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:, yi:'
.format
(xi,yi)
) 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
⎩⎨⎧yp
=yi
+hf
(xi
,yi
)yi+
1=y
i+2
h[f
(xi
,yi
)+f(
xi,
yp)
]y0
=y(a
)i=
0,1,
⋯,n−
1求解初值問題
y'=-y &(0 \leq x \leq 1) \\ y(0)=1 \\ \end \right .
, yi:'
.format
(xi,yi)
) 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
常微分方程的數值解法系列一 常微分方程
在慣性導航以及vio等實際問題中利用imu求解位姿需要對imu測量值進行積分得到需要的位置和姿態，其中主要就是求解微分方程。但之前求解微分方程的解析方法主要是應用於一些簡單和特殊的微分方程求解中，對於一般形式的微分方程，一般很難用解析方法求出精確解，只能用數值方法求解。該系列主要介紹一些常用的常微分...
常微分方程
微分方程這裡，感覺難度明顯上來了。核心思路，消去微分 分離變數法，想方設法分離變數 齊次微分方程 對於無法直接分離變數的方程，如果是y和x的次數一樣，並且不含常數項。可以可以化為齊次，變數代換求解 一階線性微分方程 常數變易法。常數變易法我覺得關鍵是變和易，因為先當作乙個常數0，是比較容易解決的。然...
常微分方程數值解上機
二步顯式 adams 法和gear 法求解，y 0 1,步長分別為h 0.1和h 0.05 1.程式文字 二步顯式 adams法 clc y 1 1 h 0.1 y 2 y 1 2 h y 1 3 h n 1 h fori 2 n y i 1 y i 3 h y i h y i 1 3 h endt...