数值计算方法 Chapter8. 常微分方程的数值解
Posted Espresso Macchiato
tags:
篇首语：本文由小常识网(cha138.com)小编为大家整理，主要介绍了数值计算方法 Chapter8. 常微分方程的数值解相关的知识，希望对你有一定的参考价值。
0. 问题描述
这一章节考察的问题如标题所述，即常微分方程的数值求解：
d y d x = f ( x , y ) y ( x 0 ) = y 0 \\left\\ \\beginaligned \\fracdydx &= f(x, y) \\\\ y(x_0) &= y_0 \\endaligned \\right. ⎩⎨⎧dxdyy(x0)=f(x,y)=y0
1. Euler公式
1. 向前Euler公式
Euler公式算是一个求解常微分方程数值解问题的一个比较直接的思路：
d y d x = δ y δ x = f ( x , y ) \\fracdydx = \\frac\\delta y\\delta x = f(x, y) dxdy=δxδy=f(x,y)
从而有：
y + δ y = f ( x , y ) ⋅ δ x y+\\delta y = f(x, y) \\cdot \\delta x y+δy=f(x,y)⋅δx
由此，我们不断地求解就能够近似的得到各个 x x x取值下的 y y y值。
y n + 1 = y n + h ⋅ f ( x n , y n ) y_n+1 = y_n + h\\cdot f(x_n, y_n) yn+1=yn+h⋅f(xn,yn)
给出python伪代码实现如下：
def fwd_euler_fn(f, x0, y0, step=1e-3):
def fn(x):
h = step if x >= x0 else -step
n = math.ceil((x - x0) / h)
x, y = x0, y0
for _ in range(n):
y += h * f(x, y)
x += h
return y
return fn
2. 向后Euler公式
向后Euler公式和向前Euler公式并无本质区别，不过微调公式为：
y n + 1 = y n + h ⋅ f ( x n + 1 , y n + 1 ) y_n+1 = y_n + h\\cdot f(x_n+1, y_n+1) yn+1=yn+h⋅f(xn+1,yn+1)
但是显然的在计算 y n + 1 y_n+1 yn+1时事实上无法得知 f ( x n + 1 , y n + 1 ) f(x_n+1, y_n+1) f(xn+1,yn+1)的值，因此这里需要叠加上迭代的思路：
y n + 1 ( k + 1 ) = y n ( k + 1 ) + h ⋅ f ( x n + 1 , y n + 1 ( k ) ) y_n+1^(k+1) = y_n^(k+1) + h\\cdot f(x_n+1, y_n+1^(k)) yn+1(k+1)=yn(k+1)+h⋅f(xn+1,yn+1(k))
上述迭代公式称为Picard迭代。
可以证明，当 h h h足够小时，Picard迭代收敛。
给出python伪代码如下：
def bwd_euler_fn(f, x0, y0, step=1e-3):
def fn(x, epsilon = 1e-6):
h = step if x >= x0 else -step
n = math.ceil((x - x0) / h)
xlist = [x0 + i*h for i in range(n)]
ylist = [y0 for i in range(n)]
for i in range(n-1):
ylist[i+1] = ylist[i] + h * f(xlist[i], ylist[i])
while True:
zlist = deepcopy(ylist)
for i in range(n-1):
zlist[i+1] = zlist[i] + h * f(xlist[i+1], ylist[i+1])
delta = abs(zlist[-1] - ylist[-1])
ylist = zlist
if delta < epsilon:
break
return zlist[-1]
return fn
3. 梯形公式
梯形公式本质上依然还是基于微分差商，不过不同于之前直接使用微分的形式，这里更加严格的使用了积分的表达，即：
y n + 1 = y n + ∫ x n x n + 1 f ( x , y ) d x y_n+1 = y_n + \\int_x_n^x_n+1f(x, y)dx yn+1=yn+∫xnxn+1f(x,y)dx
然后，这里使用梯形公式来近似掉其中的积分过程，有：
y n + 1 = y n + h 2 ( f ( x n , y n ) + f ( x n + 1 , y n + 1 ) ) y_n+1 = y_n + \\frach2(f(x_n, y_n) + f(x_n+1, y_n+1)) yn+1=yn+2h(f(xn,yn)+f(xn+1,yn+1))
这里同样会涉及到等号右侧的 y n + 1 y_n+1 yn+1的求解问题，不过，不同于向后Euler公式中的纯迭代思路，这里只使用一次迭代来近似，即：
y n + 1 = y n + h 2 ( f ( x n , y n ) + f ( x n + 1 , y n + h ⋅ f ( x n , y n ) ) ) y_n+1 = y_n + \\frach2(f(x_n, y_n) + f(x_n+1, y_n + h\\cdot f(x_n, y_n))) yn+1=yn+2h(f(xmatlab编程问题利用欧拉方法求常微分方程近似数值解
99插值法，函数逼近，曲线拟和，数值积分，数值微分，解线性方程组的直接方法，解线性方程组的迭代法，非线性方程求根，常微分方程的数值解法