文章目录
对于一阶常微分方程, 数值解法是一种用离散的点来逼近微分方程的解的方法. 常微分方程的一般形式为:
d y d x = f ( x , y ) \frac{{\rm d}y}{{\rm d}x} = f(x, y) dxdy=f(x,y)
其中 f ( x , y ) f(x, y) f(x,y)是已知的函数, y y y是未知的函数, 我们的目标是找到 y ( x ) y(x) y(x)的数值近似解.
一阶常微分方程的数值解法概括如下:
- 离散化区间 首先需要选择一个区间 [ a , b ] [a, b] [a,b], 在这个区间内求解微分方程. 通常, 这个区间会被均匀地分成若干小段, 每段的长度为 h h h.
- 初始条件 需要提供一个初始条件, 例如 y ( a ) = y 0 y(a) = y_0 y(a)=y0, 其中 y 0 y_0 y0是已知的初值.
- 选择数值方法 根据问题的性质和精确度的要求, 选择适当的数值方法, 如Euler方法, Runge-Kutta方法等.
- 迭代计算 使用选定的数值方法, 从初始条件出发, 逐步计算离散点上的近似解. 通常, 迭代的公式类似于
y i + 1 = y i + h f ( x i , y i ) y_{i+1} = y_i + hf(x_i, y_i) yi+1=yi+hf(xi,yi)
其中, x i x_i xi是离散点的横坐标, y i y_i yi是在 x i x_i xi 处的近似解, h h h是步长. - 终止条件 在每次迭代中, 可以选择一个终止条件, 如达到指定的 x x x值或达到一定的迭代次数.
- 输出结果 迭代计算得到的近似解 y ( x ) y(x) y(x)可以在所需的点上得到, 以满足问题的要求.
值得注意的是, 选择合适的步长 h h h对数值解的精确度非常重要. 太大的步长可能导致数值解的不稳定性, 而太小的步长会增加计算量. 因此, 通常需要进行步长控制, 以在保持精度的同时尽可能减少计算成本.
而对于高阶的常微分方程
y
(
n
)
+
a
1
(
x
)
y
(
n
−
1
)
+
⋯
+
a
n
−
1
(
x
)
y
′
+
a
n
(
x
)
y
=
f
(
x
)
y^{(n)}+a_1(x)y^{(n-1)}+\cdots+a_{n-1}(x)y^{\prime}+a_n(x)y=f(x)
y(n)+a1(x)y(n−1)+⋯+an−1(x)y′+an(x)y=f(x)
我们可以将其转化为等价的线性微分方程组如下:
{
d
y
1
d
x
=
y
2
d
y
2
d
x
=
y
3
⋯
⋯
d
y
n
d
x
=
−
a
n
(
x
)
y
1
−
a
n
−
1
(
x
)
y
2
−
⋯
−
a
1
(
x
)
y
n
+
f
(
x
)
\begin{cases} \dfrac{{\rm d}y_1}{{\rm d}x}=y_2\\ \dfrac{{\rm d}y_2}{{\rm d}x}=y_3\\ \cdots\cdots\\ \dfrac{{\rm d}y_n}{{\rm d}x}=-a_n(x)y_1-a_{n-1}(x)y_2-\cdots-a_1(x)y_n+f(x) \end{cases}
⎩
⎨
⎧dxdy1=y2dxdy2=y3⋯⋯dxdyn=−an(x)y1−an−1(x)y2−⋯−a1(x)yn+f(x)
关于线性微分方程组的数值解法将在下一部分介绍.
在开始之前首先要进行预处理如下:
#include <functional>
#include <math.h>
#include <vector>
#include <stdio.h>
typedef std::function<double(const double&,const double&)> func;
typedef std::vector<double> vec;
Euler方法
考虑一阶常微分方程初值问题
{
d
y
d
x
=
f
(
x
,
y
)
y
(
x
0
)
=
y
0
\begin{cases}\dfrac{{\rm d}y}{{\rm d}x}=f(x,y)\\y(x_0)=y_0\end{cases}
⎩
⎨
⎧dxdy=f(x,y)y(x0)=y0
设
f
(
x
,
y
)
∈
C
(
[
a
,
b
]
×
[
c
,
d
]
)
f(x,y)\in C([a,b]\times[c,d])
f(x,y)∈C([a,b]×[c,d]), 对
y
y
y满足Lipschitz条件, 即存在正数
L
L
L, 使得对于任意两个点
(
x
,
y
1
)
(x,y_1)
(x,y1)与
(
x
,
y
2
)
(x,y_2)
(x,y2), 有
∣
f
(
x
,
y
1
)
−
f
(
x
,
y
2
)
∣
⩽
L
∣
y
1
−
y
2
∣
|f(x,y_1)-f(x,y_2)|\leqslant L|y_1-y_2|
∣f(x,y1)−f(x,y2)∣⩽L∣y1−y2∣, 则这样初值问题的解是存在唯一的, 而且连续依赖于初始条件.
为了计算上述初值问题, 引入点列
{
x
n
}
\{x_n\}
{xn}满足
x
n
=
x
n
−
1
+
h
n
x_n=x_{n-1}+h_n
xn=xn−1+hn. 其中, 称
h
n
h_n
hn为步长. 通常考虑定长的情形, 即
∀
n
,
h
n
=
h
\forall n,h_n=h
∀n,hn=h, 此时
x
n
=
x
0
+
n
h
x_n=x_0+nh
xn=x0+nh. 根据导数的定义, 我们可以利用均差近似导数:
y
(
x
n
+
h
)
−
y
(
x
n
)
h
≈
f
(
x
n
,
y
(
x
n
)
)
\frac{y(x_n+h)-y(x_n)}h\approx f(x_n,y(x_n))
hy(xn+h)−y(xn)≈f(xn,y(xn))
令
y
n
y_n
yn为
y
(
x
n
)
y(x_n)
y(xn)的近似值, 可以将上式整理成如下递推公式:
y
n
+
1
=
y
n
+
h
f
(
x
n
,
y
n
)
y_{n+1}=y_n+hf(x_n,y_n)
yn+1=yn+hf(xn,yn)
因此我们可以逐个计算出
y
n
y_n
yn的值. 其中, 上式被称为显式Euler公式.
此外, 用均差代替后一项的导数也有
y
(
x
n
+
h
)
−
y
(
x
n
)
h
≈
f
(
x
n
+
1
,
y
(
x
n
+
1
)
)
\frac{y(x_n+h)-y(x_n)}h\approx f(x_{n+1},y(x_{n+1}))
hy(xn+h)−y(xn)≈f(xn+1,y(xn+1))
此时 y n + 1 y_{n+1} yn+1不能逐个显式计算, 因此被称为隐式Euler公式.
若将显式Euler公式和隐式Euler公式作算术平均可得梯形公式:
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 h2[f(x_n,y_n)+f(x_{n+1},y_{n+1})]
yn+1=yn+2h[f(xn,yn)+f(xn+1,yn+1)]
由于隐式公式涉及到方程求根, 不利于数值计算. 因此对于隐式公式, 通常采用估计-校正技术, 即先用显式公式计算, 得到预估值, 然后以预估值作为隐式公式的迭代初值, 用隐式公式迭代一次得到校正值, 称为预估-校正技术. 例如, 用显式Euler公式作预估, 用梯形公式作校正, 即
{
y
ˉ
n
+
1
=
y
n
+
h
f
(
x
n
,
y
n
)
y
n
+
1
=
y
n
+
h
2
[
f
(
x
n
,
y
n
)
+
f
(
x
n
+
1
,
y
ˉ
n
+
1
)
]
\begin{cases}\bar y_{n+1}=y_n+hf(x_n,y_n)\\y_{n+1}=y_n+\dfrac h2[f(x_n,y_n)+f(x_{n+1},\bar y_{n+1})]\end{cases}
⎩
⎨
⎧yˉn+1=yn+hf(xn,yn)yn+1=yn+2h[f(xn,yn)+f(xn+1,yˉn+1)]
称上式为改进的Euler公式.
由于Euler方法中
y
n
+
1
y_{n+1}
yn+1都是由
y
n
y_n
yn确定, 故称它们为单步法. 单步法可以写成如下的统一形式:
y
n
+
1
=
y
n
+
h
φ
(
x
n
,
x
n
+
1
,
y
n
,
y
n
+
1
,
h
)
y_{n+1}=y_n+h\varphi(x_n,x_{n+1},y_n,y_{n+1},h)
yn+1=yn+hφ(xn,xn+1,yn,yn+1,h)
其中
φ
\varphi
φ与
f
f
f有关. 若
φ
\varphi
φ中不含
y
n
+
1
y_{n+1}
yn+1, 则此方法是显性的, 否则为隐性的.
称
e
n
:
=
y
(
x
n
)
−
y
n
e_n:=y(x_n)-y_n
en:=y(xn)−yn为某一方法在
x
n
x_n
xn点处的整体截断误差. 由于需要使用数值方法求解的常微分方程往往不存在解析解或者解析解过于复杂, 因此整体截断误差往往不好估计. 因此可以假设
x
n
x_n
xn处的
y
n
y_n
yn没有误差, 即
y
n
=
y
(
x
n
)
y_n=y(x_n)
yn=y(xn), 考虑从
x
n
x_n
xn到
x
n
+
1
x_{n+1}
xn+1这一步的误差, 我们有局部截断误差如下:
T
n
+
1
=
y
(
x
n
+
1
)
−
y
(
x
n
)
−
h
φ
(
x
n
,
x
n
+
1
,
y
(
x
n
)
,
y
(
x
n
+
1
)
,
h
)
T_{n+1}=y(x_{n+1})-y(x_n)-h\varphi(x_n,x_{n+1},y(x_n),y(x_{n+1}),h)
Tn+1=y(xn+1)−y(xn)−hφ(xn,xn+1,y(xn),y(xn+1),h)
如果给定方法的局部截断误差
T
n
+
1
=
O
(
h
p
+
1
)
,
p
∈
Z
+
T_{n+1}=O(h^{p+1}),p\in\mathbb Z_+
Tn+1=O(hp+1),p∈Z+, 则称该方法是
p
p
p阶的, 具有
p
p
p阶精度. 进一步, 若
T
n
+
1
=
g
(
x
n
,
y
(
x
n
)
)
h
p
+
1
+
O
(
h
p
+
2
)
T_{n+1}=g(x_n,y(x_n))h^{p+1}+O(h^{p+2})
Tn+1=g(xn,y(xn))hp+1+O(hp+2)
则称第一个非零项
g
(
x
n
,
y
(
x
n
)
)
h
p
+
1
g(x_n,y(x_n))h^{p+1}
g(xn,y(xn))hp+1为该方法的局部截断误差的主项.
显式Euler法和隐式Euler法为1阶的, Euler梯形方法为2阶的.
Runge-Kutta方法
显式Euler方法是最简单的单步法, 它是一阶的, 它可以看作Taylor展开后取前两项. 因此, 得到高阶方法的一个直接想法是用Taylor展开, 如果能计算
y
(
x
)
y(x)
y(x)的高阶导数, 则可写出
p
p
p阶方法的计算公式如下:
y
n
+
1
=
y
n
+
h
y
n
′
+
h
2
2
y
n
′
′
+
⋯
+
h
p
p
!
y
n
(
p
)
y_{n+1}=y_n+hy_n^{\prime}+\frac{h^2}2y_n^{\prime\prime}+\cdots+\frac{h^p}{p!}y_n^{(p)}
yn+1=yn+hyn′+2h2yn′′+⋯+p!hpyn(p)
其中,
y
n
(
j
)
y_n^{(j)}
yn(j)为
y
(
j
)
(
x
n
)
y^{(j)}(x_n)
y(j)(xn)的近似值, 使用数学分析中的链式法则, 我们可以得到
d
2
y
d
x
2
=
∂
f
∂
x
+
∂
f
∂
y
⋅
f
\frac{{{\rm d}^2y}}{{{\rm d}x^2}} = \frac{{\partial f}}{{\partial x}} + \frac{{\partial f}}{{\partial y}} \cdot f
dx2d2y=∂x∂f+∂y∂f⋅f
d
3
y
d
x
3
=
∂
2
f
∂
x
2
+
2
∂
2
f
∂
x
∂
y
⋅
f
+
∂
f
∂
x
⋅
∂
f
∂
y
+
∂
2
f
∂
y
2
⋅
f
2
+
(
∂
f
∂
y
)
2
⋅
f
\frac{{{\rm d}^3y}}{{{\rm d}x^3}} = \frac{{\partial^2f}}{{\partial x^2}} + 2\frac{{\partial^2f}}{{\partial x\partial y}} \cdot f + \frac{{\partial f}}{{\partial x}} \cdot \frac{{\partial f}}{{\partial y}} + \frac{{\partial^2f}}{{\partial y^2}} \cdot f^2 + \left(\frac{{\partial f}}{{\partial y}}\right)^2 \cdot f
dx3d3y=∂x2∂2f+2∂x∂y∂2f⋅f+∂x∂f⋅∂y∂f+∂y2∂2f⋅f2+(∂y∂f)2⋅f
⋯
⋯
\cdots\cdots
⋯⋯
但这个方法并不实用, 因为一般情况下, 求
f
(
x
,
y
)
f(x,y)
f(x,y)的导数相当麻烦. 从计算高阶导数的公式知道, 方法的截断误差提高一阶, 需要增加的计算量很大. 因此, 我们可以用区间上若干个点的导数作线性组合得到平均斜率, 将其与解的Taylor展开式相比较, 使前面若干项吻合, 从而得到具有一定阶的方法. 这就是Runge-Kutta方法的基本思想, 其一般形式为
y
n
+
1
=
y
n
+
h
∑
i
=
1
L
λ
i
K
i
y_{n+1}=y_n+h\sum_{i=1}^L\lambda_iK_i
yn+1=yn+hi=1∑LλiKi
其中
{
K
1
=
f
(
x
n
,
y
n
)
K
i
=
f
(
x
n
+
c
i
h
,
y
n
+
c
i
h
∑
j
=
1
i
−
1
a
i
j
K
j
)
,
i
=
2
,
3
,
⋯
\begin{cases} K_1=f(x_n,y_n)\\K_i=f\left(x_n+c_ih,y_n+c_ih\sum\limits_{j=1}^{i-1}a_{ij}K_j\right),i=2,3,\cdots \end{cases}
⎩
⎨
⎧K1=f(xn,yn)Ki=f(xn+cih,yn+cihj=1∑i−1aijKj),i=2,3,⋯
c
i
⩽
1
c_i\leqslant1
ci⩽1,
∑
i
=
1
L
λ
i
=
1
\sum\limits_{i=1}^L\lambda_i=1
i=1∑Lλi=1,
∑
j
=
1
i
−
1
a
i
j
=
1
\sum_{j=1}^{i-1}\limits a_{ij}=1
j=1∑i−1aij=1. 它的局部截断误差为
T
n
+
1
=
y
(
x
n
+
1
)
−
y
(
x
n
)
−
h
∑
i
=
1
L
λ
i
K
i
∗
T_{n+1}=y(x_{n+1})-y(x_n)-h\sum_{i=1}^L\lambda_iK_i^*
Tn+1=y(xn+1)−y(xn)−hi=1∑LλiKi∗
其中, K i ∗ K_i^* Ki∗和 K i K_i Ki的区别在于用微分方程准确解 y ( x n ) y(x_n) y(xn)替代 K i K_i Ki中的 y n y_n yn就得到了 K i ∗ K_i^* Ki∗, 参数 λ i , c i , a i j \lambda_i,c_i,a_{ij} λi,ci,aij待定. 下面求出各待定参数:
- 将上式中的 y ( x n + 1 ) y(x_{n+1}) y(xn+1)在 x n x_n xn处作Taylor展开, 将 K i ∗ K_i^* Ki∗在 ( x n , y ( x n ) ) (x_n,y(x_n)) (xn,y(xn))处作二元Taylor展开
- 将展开式按 h h h的幂次进行整理
- 令 T n + 1 T_{n+1} Tn+1中 h h h的低幂次的系数为零, 使 T n + 1 T_{n+1} Tn+1首项中 h h h的幂次尽量高, 比如使 T n + 1 = O ( h p + 1 ) T_{n+1}=O(h^{p+1}) Tn+1=O(hp+1)
求出待定参数后代入上式后称为 L L L级 p p p阶Runge-Kutta方法. 下面简要介绍某些特殊的Runge-Kutta方法.
二级二阶
代入计算得方程组
{
λ
1
+
λ
2
=
1
λ
2
c
2
=
1
2
\begin{cases} \lambda_1+\lambda_2=1\\\lambda_2c_2=\dfrac12 \end{cases}
⎩
⎨
⎧λ1+λ2=1λ2c2=21
这个方程组有无穷多组解, 我们取 c 2 c_2 c2为自由未知量可以确定其余各参数的值.
取
c
2
=
1
2
c_2=\dfrac12
c2=21有中点公式:
y
n
+
1
=
y
n
+
h
f
(
x
n
+
h
2
,
y
n
+
h
2
f
(
x
n
,
y
n
)
)
y_{n+1}=y_n+hf\left(x_n+\frac h2,y_n+\frac h2f(x_n,y_n)\right)
yn+1=yn+hf(xn+2h,yn+2hf(xn,yn))
取
c
2
=
2
3
c_2=\dfrac23
c2=32有Heun公式:
y
n
+
1
=
y
n
+
h
4
[
f
(
x
n
,
y
n
)
+
3
f
(
x
n
+
2
3
h
,
y
n
+
2
3
h
f
(
x
n
,
y
n
)
)
]
y_{n+1}=y_n+\frac h4\left[f(x_n,y_n)+3f\left(x_n+\frac23h,y_n+\frac23hf(x_n,y_n)\right)\right]
yn+1=yn+4h[f(xn,yn)+3f(xn+32h,yn+32hf(xn,yn))]
取
c
2
=
1
c_2=1
c2=1得到改进的Euler公式.
三级三阶
类似于二阶方法的推导, 可以得三阶的方法, 所得系数应满足的方程组为
{
λ
1
+
λ
2
+
λ
3
=
1
λ
2
c
2
+
λ
3
c
3
=
1
2
λ
2
c
2
2
+
λ
3
c
3
2
=
1
3
λ
3
c
2
c
3
a
32
=
1
6
a
31
+
a
32
=
1
a
21
=
1
\begin{cases} \lambda_{1}+\lambda_{2}+\lambda_{3}=1\\ \lambda_{2}c_{2}+\lambda_{3}c_{3}=\dfrac12\\ \lambda_2c_2^2+\lambda_3c_3^2=\dfrac13\\ \lambda_3c_2c_3a_{32}=\dfrac16\\ a_{31}+a_{32}=1\\ a_{21}=1 \end{cases}
⎩
⎨
⎧λ1+λ2+λ3=1λ2c2+λ3c3=21λ2c22+λ3c32=31λ3c2c3a32=61a31+a32=1a21=1
该方程组的解也是不唯一的. 常见的一种三级三阶方法为
{
K
1
=
f
(
x
n
,
y
n
)
K
2
=
f
(
x
n
+
h
2
,
y
2
+
h
2
K
1
)
K
3
=
f
(
x
n
+
h
,
y
n
−
h
K
1
+
2
h
K
2
)
y
n
+
1
=
y
n
+
h
6
(
K
1
+
4
K
2
+
K
3
)
\begin{cases} K_1=f(x_n,y_n)\\ K_2=f\left(x_n+\dfrac h2,y_2+\dfrac h2K_1\right)\\ K_3=f(x_n+h,y_n-hK_1+2hK_2)\\ y_{n+1}=y_n+\dfrac h6(K_1+4K_2+K_3) \end{cases}
⎩
⎨
⎧K1=f(xn,yn)K2=f(xn+2h,y2+2hK1)K3=f(xn+h,yn−hK1+2hK2)yn+1=yn+6h(K1+4K2+K3)
四级四阶
最常用的四级四阶方法是如下的经典Runge-Kutta方法:
{
K
1
=
f
(
x
n
,
y
n
)
K
2
=
f
(
x
n
+
h
2
,
y
n
+
h
2
K
1
)
K
3
=
f
(
x
n
+
h
2
,
y
n
+
h
2
K
2
)
K
4
=
f
(
x
n
+
h
,
y
n
+
h
K
3
)
y
n
+
1
=
y
n
+
h
6
(
K
1
+
2
K
2
+
2
K
3
+
K
4
)
\begin{cases} K_1=f(x_n,y_n)\\ K_2=f\left(x_n+\dfrac h2,y_n+\dfrac h2K_1\right)\\ K_3=f\left(x_n+\dfrac h2,y_n+\dfrac h2K_2\right)\\ K_4=f(x_n+h,y_n+hK_3)\\ y_{n+1}=y_n+\dfrac h6(K_1+2K_2+2K_3+K_4) \end{cases}
⎩
⎨
⎧K1=f(xn,yn)K2=f(xn+2h,yn+2hK1)K3=f(xn+2h,yn+2hK2)K4=f(xn+h,yn+hK3)yn+1=yn+6h(K1+2K2+2K3+K4)
算法实现
显式Euler法
首先可以将算法表示为如下流程图:
向量形式
/*
* 显式欧拉法
* R:结果向量
* f:显式微分方程的右端函数
* a:区间左端点
* b:区间右端点
* y0:解在x0处的初值
* h:步长
* x0:初值点
*/
void Euler(vec &R, const func &f, const double &a, const double &b, const double &y0, const double &h, double x0)
{
if (a >= b)
throw "区间左端点小于区间右端点！";
if (x0 < a || x0 > b)
throw "初值点不在区间内！";
if (!h)
throw "步长为零！";
R = {y0};
double y(y0), x(x0);
if (h > 0)
while ((x0 += h) <= b)
{
if (isnan(y += h * f(x, y)))
throw "区间内存在奇点！";
R.push_back(y);
x = x0;
}
else
while ((x0 += h) >= a)
{
if (isnan(y += h * f(x, y)))
throw "区间内存在奇点！";
R.push_back(y);
x = x0;
}
}
/*
* 显式欧拉法
* R:结果向量
* f:显式微分方程的右端函数
* a:区间左端点
* b:区间右端点
* y0:解在x0处的初值
* h:步长
*/
inline void Euler(vec &R, const func &f, const double &a, const double &b, const double &y0, const double &h = 0.0078125)
{
return Euler(R, f, a, b, y0, h, a);
}
文件形式
/*
* 显式欧拉法
* P:结果保存路径
* f:显式微分方程的右端函数
* a:区间左端点
* b:区间右端点
* y0:解在x0处的初值
* h:步长
* x0:初值点
*/
void Euler(const char *P, const func &f, const double &a, const double &b, const double &y0, const double &h, double x0)
{
if (a >= b)
throw "区间左端点小于区间右端点！";
if (x0 < a || x0 > b)
throw "初值点不在区间内！";
if (!h)
throw "步长为零！";
FILE *file;
if (fopen_s(&file, P, "w"))
throw "文件保存失败！";
fprintf(file, "%.14f,%.14f\n", x0, y0);
double y(y0), x(x0);
if (h > 0)
while ((x0 += h) <= b)
{
if (isnan(y += h * f(x, y)))
{
fclose(file);
throw "区间内存在奇点！";
}
fprintf(file, "%.14f,%.14f\n", x = x0, y);
}
else
while ((x0 += h) >= a)
{
if (isnan(y += h * f(x, y)))
{
fclose(file);
throw "区间内存在奇点！";
}
fprintf(file, "%.14f,%.14f\n", x = x0, y);
}
fclose(file);
}
/*
* 显式欧拉法
* P:结果保存路径
* f:显式微分方程的右端函数
* a:区间左端点
* b:区间右端点
* y0:解在x0处的初值
* h:步长
* x0:初值点
*/
inline void Euler(const char *P, const func &f, const double &a, const double &b, const double &y0, const double &h = 0.0078125)
{
return Euler(P, f, a, b, y0, h, a);
}
隐式Euler法
首先将算法表示为如下流程图:
向量形式
/*
* 隐式欧拉法
* R:结果向量
* f:显式微分方程的右端函数
* a:区间左端点
* b:区间右端点
* y0:解在x0处的初值
* h:步长
* x0:初值点
*/
void imEuler(vec &R, const func &f, const double &a, const double &b, const double &y0, const double &h, double x0)
{
if (a >= b)
throw "区间左端点小于区间右端点！";
if (x0 < a || x0 > b)
throw "初值点不在区间内！";
if (!h)
throw "步长为零！";
R = {y0};
double y(y0), x(x0);
if (h > 0)
while ((x0 += h) <= b)
{
double t(y + h * f(x, y));
if (isnan(t) || isnan(y += h * f(x = x0, t)))
throw "区间内存在奇点！";
R.push_back(y);
}
else
while ((x0 += h) >= a)
{
double t(y + h * f(x, y));
if (isnan(t) || isnan(y += h * f(x = x0, t)))
throw "区间内存在奇点！";
R.push_back(y);
}
}
/*
* 隐式欧拉法
* R:结果向量
* f:显式微分方程的右端函数
* a:区间左端点
* b:区间右端点
* y0:解在x0处的初值
* h:步长
*/
inline void imEuler(vec &R, const func &f, const double &a, const double &b, const double &y0, const double &h = 0.0078125)
{
return imEuler(R, f, a, b, y0, h, a);
}
文件形式
/*
* 隐式欧拉法
* P:结果保存路径
* f:显式微分方程的右端函数
* a:区间左端点
* b:区间右端点
* y0:解在x0处的初值
* h:步长
* x0:初值点
*/
void imEuler(const char *P, const func &f, const double &a, const double &b, const double &y0, const double &h, double x0)
{
if (a >= b)
throw "区间左端点小于区间右端点！";
if (x0 < a || x0 > b)
throw "初值点不在区间内！";
if (!h)
throw "步长为零！";
FILE *file;
if (fopen_s(&file, P, "w"))
throw "文件保存失败！";
fprintf(file, "%.14f,%.14f\n", x0, y0);
double y(y0), x(x0);
if (h > 0)
while ((x0 += h) <= b)
{
double t(y + h * f(x, y));
if (isnan(t) || isnan(y += h * f(x0, t)))
{
fclose(file);
throw "区间内存在奇点！";
}
fprintf(file, "%.14f,%.14f\n", x = x0, y);
}
else
while ((x0 += h) >= a)
{
double t(y + h * f(x, y));
if (isnan(t) || isnan(y += h * f(x0, t)))
{
fclose(file);
throw "区间内存在奇点！";
}
fprintf(file, "%.14f,%.14f\n", x = x0, y);
}
fclose(file);
}
/*
* 隐式欧拉法
* P:结果保存路径
* f:显式微分方程的右端函数
* a:区间左端点
* b:区间右端点
* y0:解在x0处的初值
* h:步长
* x0:初值点
*/
inline void imEuler(const char *P, const func &f, const double &a, const double &b, const double &y0, const double &h = 0.0078125)
{
return imEuler(P, f, a, b, y0, h, a);
}
经典4阶Runge-Kutta法
算法流程图如下:
向量形式
/*
* 4阶Runge-Kutta法
* R:结果向量
* f:显式微分方程的右端函数
* a:区间左端点
* b:区间右端点
* y0:解在x0处的初值
* h:步长
* x0:初值点
*/
void Runge_Kutta4(vec &R, const func &f, const double &a, const double &b, const double &y0, const double &h, double x0)
{
if (a >= b)
throw "区间左端点小于区间右端点！";
if (x0 < a || x0 > b)
throw "初值点不在区间内！";
if (!h)
throw "步长为零！";
R = {y0};
double y(y0), x(x0);
const double h2(h / 2), h6(h2 / 3);
if (h > 0)
while ((x0 += h) <= b)
{
double k1(f(x, y));
if (isnan(k1))
throw "区间内存在奇点！";
double k2(f(x += h2, y + k1 * h2));
if (isnan(k2))
throw "区间内存在奇点！";
double k3(f(x, y + k2 * h2));
if (isnan(k3))
throw "区间内存在奇点！";
double k4(f(x0, y + h * k3));
if (isnan(k4))
throw "区间内存在奇点！";
x = x0;
R.push_back(y += h6 * (k1 + k4 + 2 * (k2 + k3)));
}
else
while ((x0 += h) >= a)
{
double k1(f(x, y));
if (isnan(k1))
throw "区间内存在奇点！";
double k2(f(x += h2, y + k1 * h2));
if (isnan(k2))
throw "区间内存在奇点！";
double k3(f(x, y + k2 * h2));
if (isnan(k3))
throw "区间内存在奇点！";
double k4(f(x0, y + h * k3));
if (isnan(k4))
throw "区间内存在奇点！";
x = x0;
R.push_back(y += h6 * (k1 + k4 + 2 * (k2 + k3)));
}
}
/*
* 4阶Runge-Kutta法
* R:结果向量
* f:显式微分方程的右端函数
* a:区间左端点
* b:区间右端点
* y0:解在x0处的初值
* h:步长
*/
inline void Runge_Kutta4(vec &R, const func &f, const double &a, const double &b, const double &y0, const double &h = 0.0078125)
{
return Runge_Kutta4(R, f, a, b, y0, h, a);
}
文件形式
/*
* 4阶Runge-Kutta法
* P:结果保存路径
* f:显式微分方程的右端函数
* a:区间左端点
* b:区间右端点
* y0:解在x0处的初值
* h:步长
* x0:初值点
*/
void Runge_Kutta4(const char *P, const func &f, const double &a, const double &b, const double &y0, const double &h, double x0)
{
if (a >= b)
throw "区间左端点小于区间右端点！";
if (x0 < a || x0 > b)
throw "初值点不在区间内！";
if (!h)
throw "步长为零！";
FILE *file;
if (fopen_s(&file, P, "w"))
throw "文件保存失败！";
fprintf(file, "%.14f,%.14f\n", x0, y0);
double y(y0), x(x0);
const double h2(h / 2), h6(h2 / 3);
if (h > 0)
while ((x0 += h) <= b)
{
double k1(f(x, y));
if (isnan(k1))
{
fclose(file);
throw "区间内存在奇点！";
}
double k2(f(x += h2, y + k1 * h2));
if (isnan(k2))
{
fclose(file);
throw "区间内存在奇点！";
}
double k3(f(x, y + k2 * h2));
if (isnan(k3))
{
fclose(file);
throw "区间内存在奇点！";
}
double k4(f(x0, y + h * k3));
if (isnan(k4))
{
fclose(file);
throw "区间内存在奇点！";
}
fprintf(file, "%.14f,%.14f\n", x = x0, y += h6 * (k1 + k4 + 2 * (k2 + k3)));
}
else
while ((x0 += h) >= a)
{
double k1(f(x, y));
if (isnan(k1))
{
fclose(file);
throw "区间内存在奇点！";
}
double k2(f(x += h2, y + k1 * h2));
if (isnan(k2))
{
fclose(file);
throw "区间内存在奇点！";
}
double k3(f(x, y + k2 * h2));
if (isnan(k3))
{
fclose(file);
throw "区间内存在奇点！";
}
double k4(f(x0, y + h * k3));
if (isnan(k4))
{
fclose(file);
throw "区间内存在奇点！";
}
fprintf(file, "%.14f,%.14f\n", x = x0, y += h6 * (k1 + k4 + 2 * (k2 + k3)));
}
fclose(file);
}
/*
* 4阶Runge-Kutta法
* P:结果保存路径
* f:显式微分方程的右端函数
* a:区间左端点
* b:区间右端点
* y0:解在x0处的初值
* h:步长
* x0:初值点
*/
inline void Runge_Kutta4(const char *P, const func &f, const double &a, const double &b, const double &y0, const double &h = 0.0078125)
{
return Runge_Kutta4(P, f, a, b, y0, h, a);
}
实例分析
例1 分别取步长
h
=
0.1
,
0.05
,
0.01
h=0.1,0.05,0.01
h=0.1,0.05,0.01, 用Euler法及经典4阶Runge-Kutta法求解初值问题
{
d
y
d
t
=
−
2
y
+
2
t
2
+
2
t
y
(
0
)
=
1
\begin{cases} \dfrac{{\rm d}y}{{\rm d}t}=-2y+2t^2+2t\\ y(0)=1 \end{cases}
⎩
⎨
⎧dtdy=−2y+2t2+2ty(0)=1
代入程序, 使用显式Euler法计算如图所示:
很容易求得方程的解析解为
y
=
e
−
2
t
+
t
2
y=e^{-2t}+t^2
y=e−2t+t2, 计算出显式Euler法的误差如图所示:
同理绘制隐式Euler法曲线图, 误差如图所示:
经典4阶Runge-Kutta法曲线图, 误差如图所示:
由图中可见, 经典4阶Runge-Kutta法误差非常小, 而显式Euler法偏小, 隐式Euler法偏大. 此外, 随着步长的减小, 误差也相应减小.
例2 用Euler法及经典4阶Runge-Kutta法取不同步长求解初值问题
{
d
y
d
x
=
−
50
y
y
(
0
)
=
1
2
,
x
∈
[
0
,
1
]
\begin{cases} \dfrac{{\rm d}y}{{\rm d}x}=-50y\\y(0)=\dfrac12,x\in[0,1] \end{cases}
⎩
⎨
⎧dxdy=−50yy(0)=21,x∈[0,1]
容易求得解析解为
y
=
e
−
50
x
2
y=\frac{e^{-50x}}2
y=2e−50x
首先使用显式Euler法求解, 取步长分别为
1
8
,
1
16
,
1
32
,
1
64
,
1
128
\dfrac18,\dfrac1{16},\dfrac1{32},\dfrac1{64},\dfrac1{128}
81,161,321,641,1281求解, 发现当
h
h
h取
1
8
\dfrac18
81和
1
16
\dfrac1{16}
161时, 算法发散, 如图所示:
当
h
=
1
32
h=\dfrac1{32}
h=321时, 算法可以收敛到零, 但误差很大, 在零点附近出现震荡; 当
h
h
h取
1
64
\dfrac1{64}
641或
1
128
\dfrac1{128}
1281时, 算法稳定, 误差较小, 如图所示:
容易知道, 当
h
h
h取
1
64
\dfrac1{64}
641或
1
128
\dfrac1{128}
1281时误差较小, 绘制误差曲线如图所示:
接着使用隐式Euler法求解, 取步长分别为
1
8
,
1
16
,
1
32
,
1
64
,
1
128
\dfrac18,\dfrac1{16},\dfrac1{32},\dfrac1{64},\dfrac1{128}
81,161,321,641,1281求解, 发现当
h
h
h取
1
8
,
1
16
\dfrac18,\dfrac1{16}
81,161和
1
32
\dfrac1{32}
321时, 算法发散, 如图所示:
当
h
h
h取
1
64
\dfrac1{64}
641或
1
128
\dfrac1{128}
1281时, 算法稳定, 如图所示:
其误差如图所示:
最后使用经典4阶Runge-Kutta法求解, 取步长分别为
1
8
,
1
16
,
1
32
,
1
64
,
1
128
\dfrac18,\dfrac1{16},\dfrac1{32},\dfrac1{64},\dfrac1{128}
81,161,321,641,1281求解, 发现当
h
h
h取
1
8
\dfrac18
81和
1
16
\dfrac1{16}
161时, 算法发散, 如图所示:
当
h
h
h取
1
32
\dfrac1{32}
321,
1
64
\dfrac1{64}
641或
1
128
\dfrac1{128}
1281时, 算法稳定, 如图所示:
其误差如图所示:
1万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
