当前位置:网站首页>数值计算方法 Chapter8. 常微分方程的数值解
数值计算方法 Chapter8. 常微分方程的数值解
2022-07-05 17:45:00 【Espresso Macchiato】
0. 问题描述
这一章节考察的问题如标题所述,即常微分方程的数值求解:
{ d y d x = f ( x , y ) y ( x 0 ) = y 0 \left\{ \begin{aligned} \frac{dy}{dx} &= f(x, y) \\ y(x_0) &= y_0 \end{aligned} \right. ⎩⎨⎧dxdyy(x0)=f(x,y)=y0
1. Euler公式
1. 向前Euler公式
Euler公式算是一个求解常微分方程数值解问题的一个比较直接的思路:
d y d x = δ y δ x = f ( x , y ) \frac{dy}{dx} = \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+1}}f(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 + \frac{h}{2}(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 + \frac{h}{2}(f(x_n, y_n) + f(x_{n+1}, y_{n} + h\cdot f(x_n, y_n))) yn+1=yn+2h(f(xn,yn)+f(xn+1,yn+h⋅f(xn,yn)))
同样,我们给出python伪代码实现如下:
def tapezoid_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/2 * (f(x, y) + f(x+h, y+h*f(x, y)))
x += h
return y
return fn
2. Runge-Kutta方法
1. 二阶Runge-Kutta方法
Runge-Kutta方法较之之前的Euler公式是一个相对而言精度更高的方法。
Euler公式来源于微分的定义,而Runge-Kutta方法则考察了Taylor展开,有:
y n + 1 = y n + ∑ i h i i ! y ( i ) ( x n ) = y n + h ⋅ y ′ ( x n ) + h 2 2 ! ⋅ y ′ ′ ( x n ) + . . . \begin{aligned} y_{n+1} &= y_n + \sum_{i} \frac{h^i}{i!} y^{(i)}(x_n) \\ &= y_n + h \cdot y'(x_n) + \frac{h^2}{2!} \cdot y''(x_n) + ... \end{aligned} yn+1=yn+i∑i!hiy(i)(xn)=yn+h⋅y′(xn)+2!h2⋅y′′(xn)+...
我们保留二阶项即有:
y n + 1 = y n + h ⋅ y ′ ( x n ) + h 2 2 y ′ ′ ( x n ) = y n + h ⋅ f ( x n , y n ) + h 2 2 ( ∂ f ∂ x ( x n , y n ) + ∂ f ∂ y ( x n , y n ) ⋅ f ( x n , y n ) ) \begin{aligned} y_{n+1} &= y_n + h \cdot y'(x_n) + \frac{h^2}{2} y''(x_n) \\ &= y_n + h \cdot f(x_n, y_n) + \frac{h^2}{2}(\frac{\partial f}{\partial x}(x_n, y_n) + \frac{\partial f}{\partial y}(x_n, y_n) \cdot f(x_n, y_n)) \end{aligned} yn+1=yn+h⋅y′(xn)+2h2y′′(xn)=yn+h⋅f(xn,yn)+2h2(∂x∂f(xn,yn)+∂y∂f(xn,yn)⋅f(xn,yn))
但是,这里的一个问题在于说两个偏导函数事实上也不一定好解,因此上述表达式没有办法直接调用。
而Runge-Kutta方法则是使用一个近似的位移公式来对其进行估计,使得二者在二阶导范围内没有误差。
具体而言:
Δ = y n + − y n = h ⋅ y ′ ( x n ) + h 2 2 y ′ ′ ( x n ) = h ⋅ f ( x n , y n ) + h 2 2 ( ∂ f ∂ x ( x n , y n ) + ∂ f ∂ y ( x n , y n ) ⋅ f ( x n , y n ) ) = c 1 ⋅ h ⋅ f ( x n , y n ) + c 2 ⋅ h ⋅ f ( x n + a h , y n + b h f ( x n , y n ) ) = c 1 ⋅ h ⋅ f ( x n , y n ) + c 2 ⋅ h ⋅ [ f ( x n , y n ) + ∂ f ∂ x ( x n , y n ) ⋅ a h + ∂ f ∂ y ( x n , y n ) ⋅ b h f ( x n , y n ) ] \begin{aligned} \Delta &= y_{n+} - y_n \\ &= h \cdot y'(x_n) + \frac{h^2}{2} y''(x_n) \\ &= h \cdot f(x_n, y_n) + \frac{h^2}{2}(\frac{\partial f}{\partial x}(x_n, y_n) + \frac{\partial f}{\partial y}(x_n, y_n) \cdot f(x_n, y_n)) \\ &= c_1\cdot h\cdot f(x_n, y_n) + c_2 \cdot h \cdot f(x_n + ah, y_n + bhf(x_n, y_n)) \\ &= c_1 \cdot h \cdot f(x_n, y_n) + c_2 \cdot h \cdot [f(x_n, y_n) + \frac{\partial f}{\partial x}(x_n, y_n) \cdot ah + \frac{\partial f}{\partial y}(x_n, y_n) \cdot bhf(x_n, y_n)] \end{aligned} Δ=yn+−yn=h⋅y′(xn)+2h2y′′(xn)=h⋅f(xn,yn)+2h2(∂x∂f(xn,yn)+∂y∂f(xn,yn)⋅f(xn,yn))=c1⋅h⋅f(xn,yn)+c2⋅h⋅f(xn+ah,yn+bhf(xn,yn))=c1⋅h⋅f(xn,yn)+c2⋅h⋅[f(xn,yn)+∂x∂f(xn,yn)⋅ah+∂y∂f(xn,yn)⋅bhf(xn,yn)]
或者说,给出参数 c 1 , c 2 , a , b c_1, c_2, a, b c1,c2,a,b使得满足如下条件:
{ y n + 1 = y n + h ( c 1 ⋅ k 1 + c 2 ⋅ k 2 ) k 1 = f ( x n , y n ) k 2 = f ( x n + a h , y n + b h k 1 ) \left\{ \begin{aligned} y_{n+1} &= y_n + h(c_1 \cdot k_1 + c_2 \cdot k_2) \\ k_1 &= f(x_n, y_n) \\ k_2 &= f(x_n + ah, y_n + bhk_1) \end{aligned} \right. ⎩⎪⎨⎪⎧yn+1k1k2=yn+h(c1⋅k1+c2⋅k2)=f(xn,yn)=f(xn+ah,yn+bhk1)
使得 y n + 1 y_{n+1} yn+1在二阶导数范围内没有误差。
亦即:
{ ( c 1 + c 2 ) h ⋅ f ( x n , y n ) = h c 2 a h 2 ⋅ ∂ f ∂ x ( x n , y n ) = h 2 2 ∂ f ∂ x ( x n , y n ) c 2 b h 2 ⋅ ∂ f ∂ y ( x n , y n ) = h 2 2 ∂ f ∂ y ( x n , y n ) \left\{ \begin{aligned} (c_1 + c_2)h \cdot f(x_n, y_n) &= h \\ c_2 a h^2 \cdot \frac{\partial f}{\partial x}(x_n, y_n) &= \frac{h^2}{2} \frac{\partial f}{\partial x}(x_n, y_n) \\ c_2 b h^2 \cdot \frac{\partial f}{\partial y}(x_n, y_n) &= \frac{h^2}{2} \frac{\partial f}{\partial y}(x_n, y_n) \end{aligned} \right. ⎩⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎧(c1+c2)h⋅f(xn,yn)c2ah2⋅∂x∂f(xn,yn)c2bh2⋅∂y∂f(xn,yn)=h=2h2∂x∂f(xn,yn)=2h2∂y∂f(xn,yn)
简化上式即可得到:
{ c 1 + c 2 = 1 2 c 2 a = 1 2 c 2 b = 1 \left\{ \begin{aligned} c_1 + c_2 &= 1 \\ 2c_2 a &= 1 \\ 2c_2 b &= 1 \end{aligned} \right. ⎩⎪⎨⎪⎧c1+c22c2a2c2b=1=1=1
只要满足上述方程组,对应的参数 a , b , c 1 , c 2 a,b,c_1,c_2 a,b,c1,c2均可以令上两式在二阶导范围内没有误差。
给出两组常见的二阶Runge-Kutta公式如下:
c 1 = 1 2 , c 2 = 1 2 , a = 1 , b = 1 c_1=\frac{1}{2}, c_2=\frac{1}{2}, a=1, b=1 c1=21,c2=21,a=1,b=1
{ k 1 = f ( x n , y n ) k 2 = f ( x n + h , y n + h k 1 ) y n + 1 = y n + h 2 ( k 1 + k 2 ) \left\{ \begin{aligned} k_1 &= f(x_n, y_n) \\ k_2 &= f(x_n + h, y_n + hk_1) \\ y_{n+1} &= y_n + \frac{h}{2}(k_1 + k_2) \end{aligned} \right. ⎩⎪⎪⎪⎨⎪⎪⎪⎧k1k2yn+1=f(xn,yn)=f(xn+h,yn+hk1)=yn+2h(k1+k2)
c 1 = 0 , c 2 = 1 , a = 1 2 , b = 1 2 c_1=0, c_2=1, a=\frac{1}{2}, b=\frac{1}{2} c1=0,c2=1,a=21,b=21
{ k 1 = f ( x n , y n ) k 2 = f ( x n + h 2 , y n + h 2 k 1 ) y n + 1 = y n + h k 2 \left\{ \begin{aligned} k_1 &= f(x_n, y_n) \\ k_2 &= f(x_n + \frac{h}{2}, y_n + \frac{h}{2}k_1) \\ y_{n+1} &= y_n + hk_2 \end{aligned} \right. ⎩⎪⎪⎪⎨⎪⎪⎪⎧k1k2yn+1=f(xn,yn)=f(xn+2h,yn+2hk1)=yn+hk2
2. 高阶Runge-Kutta方法
同样的,我们仿照上述的思路,给出一般情况下的高阶Runge-Kutta方法的表达式如下:
{ y n + 1 = y n + h ( c 1 ⋅ k 1 + c 2 ⋅ k 2 + . . . + c m ⋅ k m ) k 1 = f ( x n , y n ) k 2 = f ( x n + a 1 h + y n + b 1 , 1 h ⋅ k 1 ) . . . k m = f ( x n + a m − 1 h + y n + ∑ i = 1 m − 1 b m − 1 , i h ⋅ k i ) \left\{ \begin{aligned} y_{n+1} &= y_n + h(c_1 \cdot k_1 + c_2 \cdot k_2 + ... + c_m \cdot k_m) \\ k_1 &= f(x_n, y_n) \\ k_2 &= f(x_n + a_{1} h + y_n + b_{1,1} h \cdot k_1) \\ &... \\ k_m &= f(x_n + a_{m-1} h + y_n + \sum_{i=1}^{m-1}b_{m-1,i} h \cdot k_{i}) \end{aligned} \right. ⎩⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎧yn+1k1k2km=yn+h(c1⋅k1+c2⋅k2+...+cm⋅km)=f(xn,yn)=f(xn+a1h+yn+b1,1h⋅k1)...=f(xn+am−1h+yn+i=1∑m−1bm−1,ih⋅ki)
参数 a , b , c a,b,c a,b,c可以使得 y n y_n yn具有 O ( h m + 1 ) O(h^{m+1}) O(hm+1)阶精度。
亦即,Taylor展开直到 h m h^{m} hm都有两侧的展开系数相同。
1. 三阶Runge-Kutta方法
我们给出三阶Runge-Kutta方法的两组典型的系数如下:
系数组合(一)
{ y n + 1 = y n + h 6 ( k 1 + 4 k 2 + k 3 ) k 1 = f ( x n , y n ) k 2 = f ( x n + 1 2 h , y n + 1 2 h k 1 ) k 3 = f ( x n + h , y n − h k 1 + 2 h k 2 ) \left\{ \begin{aligned} y_{n+1} &= y_n + \frac{h}{6}(k_1 + 4k_2 + k_3) \\ k_1 &= f(x_n, y_n) \\ k_2 &= f(x_n + \frac{1}{2}h, y_n + \frac{1}{2}hk_1) \\ k_3 &= f(x_n + h, y_n - hk_1 + 2hk_2) \end{aligned} \right. ⎩⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎧yn+1k1k2k3=yn+6h(k1+4k2+k3)=f(xn,yn)=f(xn+21h,yn+21hk1)=f(xn+h,yn−hk1+2hk2)
系数组合(二)
{ y n + 1 = y n + h 4 ( k 1 + 3 k 3 ) k 1 = f ( x n , y n ) k 2 = f ( x n + 1 3 h , y n + 1 3 h k 1 ) k 3 = f ( x n + 2 3 h , y n + 2 3 h k 2 ) \left\{ \begin{aligned} y_{n+1} &= y_n + \frac{h}{4}(k_1 + 3k_3) \\ k_1 &= f(x_n, y_n) \\ k_2 &= f(x_n + \frac{1}{3}h, y_n + \frac{1}{3}hk_1) \\ k_3 &= f(x_n + \frac{2}{3}h, y_n + \frac{2}{3}hk_2) \end{aligned} \right. ⎩⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎧yn+1k1k2k3=yn+4h(k1+3k3)=f(xn,yn)=f(xn+31h,yn+31hk1)=f(xn+32h,yn+32hk2)
系数组合(三)
{ y n + 1 = y n + h 9 ( 2 k 1 + 3 k 2 + 4 k 3 ) k 1 = f ( x n , y n ) k 2 = f ( x n + 1 2 h , y n + 1 2 h k 1 ) k 3 = f ( x n + 3 4 h , y n + 3 4 h k 2 ) \left\{ \begin{aligned} y_{n+1} &= y_n + \frac{h}{9}(2k_1 + 3k_2 + 4k_3) \\ k_1 &= f(x_n, y_n) \\ k_2 &= f(x_n + \frac{1}{2}h, y_n + \frac{1}{2}hk_1) \\ k_3 &= f(x_n + \frac{3}{4}h, y_n + \frac{3}{4}hk_2) \end{aligned} \right. ⎩⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎧yn+1k1k2k3=yn+9h(2k1+3k2+4k3)=f(xn,yn)=f(xn+21h,yn+21hk1)=f(xn+43h,yn+43hk2)
2. 四阶Runge-Kutta方法
同样的,我们可以给出两组典型的四阶Runge-Kutta公式如下:
系数组合(一)
{ y n + 1 = y n + h 6 ( k 1 + 2 k 2 + 2 k 3 + k 4 ) k 1 = f ( x n , y n ) k 2 = f ( x n + 1 2 h , y n + 1 2 h k 1 ) k 3 = f ( x n + 1 2 h , y n + 1 2 h k 2 ) k 4 = f ( x n + h , y n + h k 3 ) \left\{ \begin{aligned} y_{n+1} &= y_n + \frac{h}{6}(k_1 + 2k_2 + 2k_3 + k_4) \\ k_1 &= f(x_n, y_n) \\ k_2 &= f(x_n + \frac{1}{2}h, y_n + \frac{1}{2}hk_1) \\ k_3 &= f(x_n + \frac{1}{2}h, y_n + \frac{1}{2}hk_2) \\ k_4 &= f(x_n + h, y_n + hk_3) \end{aligned} \right. ⎩⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎧yn+1k1k2k3k4=yn+6h(k1+2k2+2k3+k4)=f(xn,yn)=f(xn+21h,yn+21hk1)=f(xn+21h,yn+21hk2)=f(xn+h,yn+hk3)
系数组合(二)
{ y n + 1 = y n + h 8 ( k 1 + 3 k 2 + 3 k 3 + k 4 ) k 1 = f ( x n , y n ) k 2 = f ( x n + 1 3 h , y n + 1 3 h k 1 ) k 3 = f ( x n + 2 3 h , y n + 1 3 h k 1 + h k 2 ) k 4 = f ( x n + h , y n + h k 1 − h k 2 + h k 3 ) \left\{ \begin{aligned} y_{n+1} &= y_n + \frac{h}{8}(k_1 + 3k_2 + 3k_3 + k_4) \\ k_1 &= f(x_n, y_n) \\ k_2 &= f(x_n + \frac{1}{3}h, y_n + \frac{1}{3}hk_1) \\ k_3 &= f(x_n + \frac{2}{3}h, y_n + \frac{1}{3}hk_1 + hk_2) \\ k_4 &= f(x_n + h, y_n + hk_1 - hk_2 + hk_3) \end{aligned} \right. ⎩⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎧yn+1k1k2k3k4=yn+8h(k1+3k2+3k3+k4)=f(xn,yn)=f(xn+31h,yn+31hk1)=f(xn+32h,yn+31hk1+hk2)=f(xn+h,yn+hk1−hk2+hk3)
3. python伪代码实现
下面,我们来给出Runge-Kutta方法的python伪代码实现如下:
def runge_kutta_fn(f, x0, y0, a, b, c, step=1e-3):
def fn(x):
h = step if x >= x0 else -step
n = math.ceil((x - x0) / h)
m = len(c)
assert(len(a) == m-1 and all(len(b[i]) == i+1 for i in range(m-1)))
k = [0 for _ in range(m)]
x, y = x0, y0
for _ in range(n):
k[0] = f(x, y)
for i in range(m-1):
xi = x + a[i] * h
yi = y
for j in range(i+1):
yi += h * b[i][j] * k[j]
k[i+1] = f(xi, yi)
x += h
y += sum(ci * ki * h for ci, ki in zip(c, k))
return y
return fn
3. 线性多步法
1. 基本思路
线性多步法的思路来源依然还是积分公式:
y ( x ) = y ( x 0 ) + ∫ x 0 x y ′ ( t ) d t y(x) = y(x_0) + \int_{x_0}^{x}y'(t)dt y(x)=y(x0)+∫x0xy′(t)dt
之前,Euler公式和Runge-Kutta方法都是直接对积分的值进行估计。
其中,Euler公式是直接令 ∫ x x + h y ′ ( t ) d t = y ′ ( x ) h \int_{x}^{x+h}y'(t)dt = y'(x)h ∫xx+hy′(t)dt=y′(x)h,而Runge-Kutta方法则是通过引入一系列的偏移量使得 y ( x ) y(x) y(x)在 n n n阶Taylor展开当中没有误差。
而线性多步法的近似思路则是用采用之前的插值公式的思路,来对 y ′ ( x ) y'(x) y′(x)来进行拟合,然后用这个拟合函数来计算后面这个积分值。
给出书中的描述如下:
若用积分节点 x n , x n − 1 , . . . , x n − q x_n, x_{n-1}, ..., x_{n-q} xn,xn−1,...,xn−q构造插值多项式近似 y ′ ( x ) y'(x) y′(x),在区间 [ x n − p , x n + 1 ] [x_{n-p}, x_{n+1}] [xn−p,xn+1]上计算数值积分 ∫ x n − p x n + 1 y ′ ( x ) d x \int_{x_{n-p}}^{x_{n+1}}y'(x)dx ∫xn−pxn+1y′(x)dx,则称构造计算 y n + 1 y_{n+1} yn+1的方法为线性多步法。
特别的:
- 如果使用 x n , x n − 1 , . . . , x n − q x_{n}, x_{n-1}, ..., x_{n-q} xn,xn−1,...,xn−q构造插值多项式,则称拟合函数为显式Adams公式;
- 如果使用 x n + 1 , x n , . . . , x n + 1 − q x_{n+1}, x_{n}, ..., x_{n+1-q} xn+1,xn,...,xn+1−q构造插值多项式,则称拟合函数为隐式Admas公式;
2. Adams公式
这里,我们舍去推导,直接给出一些常见的Adams公式如下:
二阶显式Adams公式
y n + 1 = y n + h 2 ( 3 f ( x n , y n ) − f ( x n − 1 , y n − 1 ) ) y_{n+1} = y_{n} + \frac{h}{2}(3f(x_{n}, y_{n}) - f(x_{n-1}, y_{n-1})) yn+1=yn+2h(3f(xn,yn)−f(xn−1,yn−1))
二阶隐式Adams公式
y n + 1 = y n + h 2 ( f ( x n , y n ) + f ( x n + 1 , y n + 1 ) ) y_{n+1} = y_{n} + \frac{h}{2}(f(x_{n}, y_{n}) + f(x_{n+1}, y_{n+1})) yn+1=yn+2h(f(xn,yn)+f(xn+1,yn+1))
三阶显示Adams公式
y n + 1 = y n + h 12 ( 23 f ( x n , y n ) − 16 f ( x n − 1 , y n − 1 ) + 5 f ( x n − 2 , y n − 2 ) ) y_{n+1} = y_{n} + \frac{h}{12}(23f(x_{n}, y_{n}) - 16f(x_{n-1}, y_{n-1}) + 5f(x_{n-2}, y_{n-2})) yn+1=yn+12h(23f(xn,yn)−16f(xn−1,yn−1)+5f(xn−2,yn−2))
三阶隐示Adams公式
y n + 1 = y n + h 12 ( 5 f ( x n + 1 , y n + 1 + 8 f ( x n , y n ) − f ( x n − 1 , y n − 1 ) ) y_{n+1} = y_{n} + \frac{h}{12}(5f(x_{n+1}, y_{n+1} + 8f(x_{n}, y_{n}) - f(x_{n-1}, y_{n-1})) yn+1=yn+12h(5f(xn+1,yn+1+8f(xn,yn)−f(xn−1,yn−1))
四阶显式Adams公式
y n + 1 = y n + h 24 ( 55 f ( x n , y n ) − 59 f ( x n − 1 , y n − 1 ) + 37 f ( x n − 2 , y n − 2 ) − 9 f ( x n − 3 , y n − 3 ) ) y_{n+1} = y_{n} + \frac{h}{24}(55f(x_{n}, y_{n}) - 59f(x_{n-1}, y_{n-1}) + 37f(x_{n-2}, y_{n-2}) - 9f(x_{n-3}, y_{n-3})) yn+1=yn+24h(55f(xn,yn)−59f(xn−1,yn−1)+37f(xn−2,yn−2)−9f(xn−3,yn−3))
四阶隐式Adams公式
y n + 1 = y n + h 24 ( 9 f ( x n + 1 , y n + 1 + 19 f ( x n , y n ) − 5 f ( x n − 1 , y n − 1 ) + f ( x n − 2 , y n − 2 ) ) y_{n+1} = y_{n} + \frac{h}{24}(9f(x_{n+1}, y_{n+1} + 19f(x_{n}, y_{n}) - 5f(x_{n-1}, y_{n-1}) + f(x_{n-2}, y_{n-2})) yn+1=yn+24h(9f(xn+1,yn+1+19f(xn,yn)−5f(xn−1,yn−1)+f(xn−2,yn−2))
4. 常微分方程组的数值解法
1. 一阶常微分方程组的数值解法
我们给出一阶常微分方程的初值问题表达如下:
{ d y 1 d x = f 1 ( x , y 1 , y 2 , . . . , y m ) d y 2 d x = f 2 ( x , y 1 , y 2 , . . . , y m ) . . . d y m d x = f m ( x , y 1 , y 2 , . . . , y m ) y 1 ( x 0 ) = η 1 y 2 ( x 0 ) = η 2 . . . y m ( x 0 ) = η m \left\{ \begin{aligned} \frac{dy_{1}}{dx} &= f_1(x, y_{1}, y_{2}, ..., y_{m}) \\ \frac{dy_{2}}{dx} &= f_2(x, y_{1}, y_{2}, ..., y_{m}) \\ &... \\ \frac{dy_{m}}{dx} &= f_m(x, y_{1}, y_{2}, ..., y_{m}) \\ y_1(x_0) &= \eta_1 \\ y_2(x_0) &= \eta_2 \\ &... \\ y_m(x_0) &= \eta_m \\ \end{aligned} \right. ⎩⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎧dxdy1dxdy2dxdymy1(x0)y2(x0)ym(x0)=f1(x,y1,y2,...,ym)=f2(x,y1,y2,...,ym)...=fm(x,y1,y2,...,ym)=η1=η2...=ηm
这类方程组的解法问题其实完全可以原模原样照搬常微分方程的数值解法,倒是也没啥必要详细展开就是了。
我们直接以二元方程组为例,给出一些常见的解:
{ d y d x = f ( x , y , z ) d z d x = g ( x , y , z ) y ( x 0 ) = y 0 z ( x 0 ) = z 0 \left\{ \begin{aligned} \frac{dy}{dx} &= f(x, y, z) \\ \frac{dz}{dx} &= g(x, y, z) \\ y(x_0) &= y_0 \\ z(x_0) &= z_0 \end{aligned} \right. ⎩⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎧dxdydxdzy(x0)z(x0)=f(x,y,z)=g(x,y,z)=y0=z0
前向Euler公式
{ y n + 1 = y n + h f ( x n , y n , z n ) z n + 1 = z n + h g ( x n , y n , z n ) \left\{ \begin{aligned} y_{n+1} &= y_n + hf(x_n, y_n, z_n) \\ z_{n+1} &= z_n + hg(x_n, y_n, z_n) \end{aligned} \right. { yn+1zn+1=yn+hf(xn,yn,zn)=zn+hg(xn,yn,zn)
梯形公式
{ y n + 1 ^ = y n + h f ( x n , y n , z n ) z n + 1 ^ = z n + h g ( x n , y n , z n ) y n + 1 = y n + h 2 ( f ( x n , y n , z n ) + f ( x n + 1 , y n + 1 ^ , z n + 1 ^ ) ) z n + 1 = z n + h 2 ( g ( x n , y n , z n ) + g ( x n + 1 , y n + 1 ^ , z n + 1 ^ ) ) \left\{ \begin{aligned} \hat{y_{n+1}} &= y_n + hf(x_n, y_n, z_n) \\ \hat{z_{n+1}} &= z_n + hg(x_n, y_n, z_n) \\ y_{n+1} &= y_n + \frac{h}{2}(f(x_n, y_n, z_n) + f(x_{n+1}, \hat{y_{n+1}}, \hat{z_{n+1}})) \\ z_{n+1} &= z_n + \frac{h}{2}(g(x_n, y_n, z_n) + g(x_{n+1}, \hat{y_{n+1}}, \hat{z_{n+1}})) \end{aligned} \right. ⎩⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎧yn+1^zn+1^yn+1zn+1=yn+hf(xn,yn,zn)=zn+hg(xn,yn,zn)=yn+2h(f(xn,yn,zn)+f(xn+1,yn+1^,zn+1^))=zn+2h(g(xn,yn,zn)+g(xn+1,yn+1^,zn+1^))
四阶Runge-Kutta方法
{ y n + 1 = y n + h 6 ( k 11 + 2 k 12 + 2 k 13 + k 14 ) z n + 1 = z n + h 6 ( k 21 + 2 k 22 + 2 k 23 + k 24 ) k 11 = f ( x n , y n , z n ) k 21 = g ( x n , y n , z n ) k 12 = f ( x n + h 2 , y n + h 2 k 11 , z n + h 2 k 21 ) k 22 = g ( x n + h 2 , y n + h 2 k 11 , z n + h 2 k 21 ) k 13 = f ( x n + h 2 , y n + h 2 k 12 , z n + h 2 k 22 ) k 23 = g ( x n + h 2 , y n + h 2 k 12 , z n + h 2 k 22 ) k 14 = f ( x n + h , y n + h k 13 , z n + h k 23 ) k 24 = g ( x n + h , y n + h k 13 , z n + h k 23 ) \left\{ \begin{aligned} y_{n+1} &= y_n + \frac{h}{6}(k_{11} + 2k_{12} + 2k_{13} + k_{14}) \\ z_{n+1} &= z_n + \frac{h}{6}(k_{21} + 2k_{22} + 2k_{23} + k_{24}) \\ k_{11} &= f(x_n, y_n, z_n) \\ k_{21} &= g(x_n, y_n, z_n) \\ k_{12} &= f(x_n + \frac{h}{2}, y_n + \frac{h}{2}k_{11}, z_n + \frac{h}{2}k_{21}) \\ k_{22} &= g(x_n + \frac{h}{2}, y_n + \frac{h}{2}k_{11}, z_n + \frac{h}{2}k_{21}) \\ k_{13} &= f(x_n + \frac{h}{2}, y_n + \frac{h}{2}k_{12}, z_n + \frac{h}{2}k_{22}) \\ k_{23} &= g(x_n + \frac{h}{2}, y_n + \frac{h}{2}k_{12}, z_n + \frac{h}{2}k_{22}) \\ k_{14} &= f(x_n + h, y_n + hk_{13}, z_n + hk_{23}) \\ k_{24} &= g(x_n + h, y_n + hk_{13}, z_n + hk_{23}) \end{aligned} \right. ⎩⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎧yn+1zn+1k11k21k12k22k13k23k14k24=yn+6h(k11+2k12+2k13+k14)=zn+6h(k21+2k22+2k23+k24)=f(xn,yn,zn)=g(xn,yn,zn)=f(xn+2h,yn+2hk11,zn+2hk21)=g(xn+2h,yn+2hk11,zn+2hk21)=f(xn+2h,yn+2hk12,zn+2hk22)=g(xn+2h,yn+2hk12,zn+2hk22)=f(xn+h,yn+hk13,zn+hk23)=g(xn+h,yn+hk13,zn+hk23)
2. 高阶微分方程数值方法
这里,我们再来考察一下一元高阶微分方程的数值解法。
{ d m y ( x ) d x m = f ( x , y , y ′ , . . . , y ( m − 1 ) ) y ( x 0 ) = η 0 y ′ ( x 0 ) = η 1 . . . y ( m − 1 ) ( x 0 ) = η m − 1 \left\{ \begin{aligned} \frac{d^{m}y(x)}{dx^{m}} &= f(x, y, y', ..., y^{(m-1)}) \\ y(x_0) &= \eta_0 \\ y'(x_0) &= \eta_1 \\ & ... \\ y^{(m-1)}(x_0) &= \eta_{m-1} \end{aligned} \right. ⎩⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎧dxmdmy(x)y(x0)y′(x0)y(m−1)(x0)=f(x,y,y′,...,y(m−1))=η0=η1...=ηm−1
这一类问题事实上可以作为上述一阶常微分方程组的一个应用实例,我们只需要做如下变换就可以将问题完全转换为一个一阶常微分方程组,然后就可以运用之前的一阶常微分方程组的数值解法进行求解了。
{ d y 0 d x = y 1 ( x ) d y 1 d x = y 2 ( x ) . . . d y m − 2 d x = y m − 1 ( x ) d y m − 1 d x = f ( x , y 0 , y 1 , . . . , y m − 1 ) y 0 ( x 0 ) = η 0 y 1 ( x 0 ) = η 1 . . . y m − 1 ( x 0 ) = η m − 1 \left\{ \begin{aligned} \frac{dy_0}{dx} &= y_1(x) \\ \frac{dy_1}{dx} &= y_2(x) \\ &... \\ \frac{dy_{m-2}}{dx} &= y_{m-1}(x) \\ \frac{dy_{m-1}}{dx} &= f(x, y_0, y_1, ..., y_{m-1}) \\ y_{0}(x_0) &= \eta_0 \\ y_{1}(x_0) &= \eta_1 \\ & ... \\ y_{m-1}(x_0) &= \eta_{m-1} \end{aligned} \right. ⎩⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎧dxdy0dxdy1dxdym−2dxdym−1y0(x0)y1(x0)ym−1(x0)=y1(x)=y2(x)...=ym−1(x)=f(x,y0,y1,...,ym−1)=η0=η1...=ηm−1
边栏推荐
- Sentinel flow guard
- 外盘黄金哪个平台正规安全,怎么辨别?
- Knowledge points of MySQL (6)
- This 17-year-old hacker genius cracked the first generation iPhone!
- Ordinary programmers look at the code, and top programmers look at the trend
- 7 pratiques devops pour améliorer la performance des applications
- Abnormal recovery of virtual machine Oracle -- Xi Fenfei
- ISPRS2022/云检测:Cloud detection with boundary nets基于边界网的云检测
- mybash
- Disabling and enabling inspections pycharm
猜你喜欢
2022新版PMP考试有哪些变化?
"Xiaodeng in operation and maintenance" is a single sign on solution for cloud applications
flask接口响应中的中文乱码(unicode)处理
网络威胁分析师应该具备的十种能力
Zabbix
Ten top automation and orchestration tools
mybash
What are the changes in the 2022 PMP Exam?
nacos -分布式事务-Seata** linux安装jdk ,mysql5.7启动nacos配置ideal 调用接口配合 (保姆级细节教程)
Mask wearing detection based on yolov3
随机推荐
证券网上开户安全吗?证券融资利率一般是多少?
Customize the theme of matrix (I) night mode
Tkinter window preload
Tencent music launched its new product "quyimai", which provides music commercial copyright authorization
2022 information system management engineer examination outline
Webapp development - Google official tutorial
如何修改mysql字段为自增长字段
Server configuration jupyter environment
[BeanShell] there are many ways to write data locally
Disabling and enabling inspections pycharm
Short the command line via jar manifest or via a classpath file and rerun
Simple query cost estimation
ELK日志分析系统
Data access - entityframework integration
使用QT设计师界面类创建2个界面,通过按键从界面1切换到界面2
PMP认证需具备哪些条件啊?费用多少啊?
较文心损失一点点性能提升很多
Compter le temps d'exécution du programme PHP et définir le temps d'exécution maximum de PHP
[JMeter] advanced writing method of JMeter script: all variables, parameters (parameters can be configured by Jenkins), functions, etc. in the interface automation script realize the complete business
c#图文混合,以二进制方式写入数据库