因研究一个新项目,了解到一些数值方法,记录下来,避免忘记。
前向欧拉法 (FE)
前向欧拉法可用于计算微分方程数值解。
根据导数的定义,当 h h h 足够小时,有:
y
(
x
0
+
h
)
−
y
(
x
0
)
h
=
y
′
(
x
0
)
\frac {y(x_0+h) - y(x_0)}{h} = y'(x_0)
hy(x0+h)−y(x0)=y′(x0)
稍加变换,可得:
y
(
x
0
+
h
)
=
y
(
x
0
)
+
h
y
′
(
x
0
)
y(x_0+h) = y(x_0) + hy'(x_0)
y(x0+h)=y(x0)+hy′(x0)
设所求函数
y
y
y 的导数是关于
x
,
y
x, y
x,y 的方程, 记为
f
(
x
,
y
)
f(x, y)
f(x,y) , 则有:
y
(
x
n
+
h
)
=
y
n
+
h
f
(
x
n
,
y
n
)
y(x_n + h) = y_n + hf(x_n, y_n)
y(xn+h)=yn+hf(xn,yn)
紧凑形式:
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)
根据以上数学形式,可实现如下算法计算在
x
x
x 处原函数
y
y
y 的值
y
(
x
)
y(x)
y(x):
function FE_iterate(f, x0, y0, h, n) {
let t=x0;
let $y=y0;
for (let i=0;i<n;i++) {
t = i*h + x0;
$y += f(t, $y)*h;
t += h;
}
return [t, $y];
}
function FE(f, x0, y0, h) {
return function y(x) {
const N = Math.trunc((x-x0)/h);
let [t, $y] = N < 0 ? FE_iterate(f, x0, y0, -h, -N) : FE_iterate(f, x0, y0, h, N);
if (Math.abs(x-t)<1e-15) return $y;
return FE_iterate(f, t, $y, x - t, 1)[1];
};
}
假设有一函数 y y y 满足 y ′ = y y' = y y′=y 且 y ( 0 ) = 1 y(0) = 1 y(0)=1, 可解得 y = e x y=e^x y=ex。我们可以取 h = 0.001 h=0.001 h=0.001 用上述算法计算 y ( 1 ) y(1) y(1) :
FE((x,y)=>y, 0, 1, 0.001)(1)
得到 2.716923932235896
, 与真实值 2.718281828459045
还算接近。随着迭代次数增加误差会不可避免的增加,关于前向欧拉法的误差估计,可以参考数值方法相关的书籍。
梯形法 (Trapezoidal, TR )
梯形法从积分的计算出发,当
h
h
h 足够小时,有:
∫
x
0
x
0
+
h
y
′
(
x
)
d
x
=
1
2
[
y
′
(
x
0
)
+
y
′
(
x
0
+
h
)
]
h
\int^{x_0+h}_{x_0} y'(x) dx = \frac {1}{2}[y'(x_0)+y'(x_0+h)]h
∫x0x0+hy′(x)dx=21[y′(x0)+y′(x0+h)]h
即:
y
(
x
0
+
h
)
−
y
(
x
0
)
=
1
2
[
y
′
(
x
0
)
+
y
′
(
x
0
+
h
)
]
h
y(x_0+h) - y(x_0) = \frac {1}{2}[y'(x_0)+y'(x_0+h)]h
y(x0+h)−y(x0)=21[y′(x0)+y′(x0+h)]h
依然设所求函数
y
y
y 的导数是关于
x
,
y
x, y
x,y 的方程, 记为
f
(
x
,
y
)
f(x, y)
f(x,y) , 则有:
y
(
x
n
+
h
)
=
y
(
x
n
)
+
1
2
[
f
(
x
n
,
y
n
)
+
f
(
x
n
+
h
,
y
(
x
n
+
h
)
)
]
h
y(x_n+h) = y(x_n) + \frac {1}{2}[f(x_n, y_n)+f(x_n+h, y(x_n + h))]h
y(xn+h)=y(xn)+21[f(xn,yn)+f(xn+h,y(xn+h))]h
公式中
y
y
y 在
x
0
+
h
x_0 + h
x0+h 处的导数需要用
y
(
x
0
+
h
)
y(x_0+h)
y(x0+h) 计算,而
y
(
x
0
+
h
)
y(x_0+h)
y(x0+h) 是我们要计算的。这里可以先应用前向欧拉法估算,于是得到:
y
(
x
n
+
h
)
=
y
(
x
n
)
+
1
2
[
f
(
x
n
,
y
n
)
+
f
(
x
n
+
h
,
y
(
x
n
)
+
h
f
(
x
n
,
y
n
)
)
]
h
y(x_n+h) = y(x_n) + \frac {1}{2}[f(x_n, y_n)+f(x_n+h, y(x_n) + hf(x_n, y_n))]h
y(xn+h)=y(xn)+21[f(xn,yn)+f(xn+h,y(xn)+hf(xn,yn))]h
紧凑形式为:
y
n
+
1
=
y
n
+
1
2
[
f
(
x
n
,
y
n
)
+
f
(
x
n
+
1
,
y
n
+
h
f
(
x
n
,
y
n
)
)
]
h
y_{n+1} = y_n + \frac {1}{2}[f(x_n, y_n)+f(x_{n+1}, y_n + hf(x_n, y_n))]h
yn+1=yn+21[f(xn,yn)+f(xn+1,yn+hf(xn,yn))]h
根据以上数学形式,可以实现如下算法计算在
x
x
x 处原函数
y
y
y 的值
y
(
x
)
y(x)
y(x):
function TR_iterate(f, x0, y0, h, n) {
let t=x0;
let $y=y0;
for (let i=0;i<n;i++) {
t = i*h + x0;
const k1 = f(t, $y);
const k2 = f(t+h, $y+h*k1);
$y += 0.5*(k1+k2)*h;
t += h;
}
return [t, $y];
}
function TR(f, x0, y0, h) {
return function y(x) {
const N = Math.trunc((x-x0)/h);
let [t, $y] = N < 0 ? TR_iterate(f, x0, y0, -h, -N) : TR_iterate(f, x0, y0, h, N);
if (Math.abs(x-t)<1e-15) return $y;
return TR_iterate(f, t, $y, x - t, 1)[1];
};
}
取 h = 0.001 h=0.001 h=0.001 用上述算法计算欧拉法中的示例:
TR((x,y)=>y, 0, 1, 0.001)(1)
得到 2.7182813757517628
, 与真实值 2.718281828459045
更加接近。由此可见,梯形法比前向欧拉法更加精确。但是随着迭代次数增加,误差也会不可避免的增加,关于梯形法的误差估计,可以参考数值方法相关的书籍。
4 阶龙格库塔法 (Runger-Kutta)
龙格库塔法则具有更高精度,具体数学推导可以参考数值方法相关的书籍。这里仅给出代码实现:
function RK4_iterate(f, x0, y0, h, n) {
let t=x0;
let $y=y0;
for (let i=0;i<n;i++) {
t = i*h + x0;
const k1 = f(t, $y);
const k2 = f(t+h/2, $y+k1*h/2);
const k3 = f(t+h/2, $y+k2*h/2);
const k4 = f(t+h, $y+k3*h);
$y += h*(k1+2*k2+2*k3+k4)/6;
t += h;
}
return [t, $y];
}
function RK4(f, x0, y0, h) {
return function y(x) {
const N = Math.trunc((x-x0)/h);
let [t, $y] = N < 0 ? RK4_iterate(f, x0, y0, -h, -N) : RK4_iterate(f, x0, y0, h, N);
if (Math.abs(x-t)<1e-15) return $y;
return RK4_iterate(f, t, $y, x - t, 1)[1];
};
}
取 h = 0.001 h=0.001 h=0.001 用上述算法计算欧拉法中的示例:
Rk4((x,y)=>y, 0, 1, 0.001)(1)
得到 2.7182818284590247
, 比梯形法更加接近真实值。但同样地,它也有误差,具体可以参考数值方法相关的书籍。
牛顿迭代法 (Newton-Raphson Method)
牛顿迭代法可用于线性方程数值求解,例如求平方根。原始牛顿迭代法需要用到函数的导数,代码实现过程中可以采用差分法近似求导。代码如下:
function NR(f, x) {
const eps = 1e-15;
const h = 1e-15;
let dt = 0;
let fx = 0;
let iter = 500;
do {
x += dt;
fx = f(x);
dt = -h/(f(x+h)/fx - 1);
} while (iter-->0 && Math.abs(dt)>eps && Math.abs(fx)>eps);
return x;
}
我们用它计算 2 \sqrt2 2 :
NR(x=>x*x-2, 1)
可以得到 1.4142135623730958
, 与 Math.sqrt
求得的结果仅在最后一位有差异。调整 eps
和 iter
可以调整精度。牛顿迭代法收敛速度快,求解快速,但是能否求解与初始值的选择有很大关系。例如,
x
3
−
3
x
2
+
x
+
3
x^3-3x^2+x+3
x3−3x2+x+3 如果初始值选择为 2 则永远不收敛,迭代过程会在 1 跟 2 之间反复横跳,如下:
如果遇到反复横跳,可以在 x += dt
的计算中引入一个新参数 k
改善,如 x += k*dt
。这个改善会减慢收敛速度,但会增加收敛的概率。
牛顿迭代法可以扩展到多函数情况,用于解方程组。这时方程组的导数为雅可比矩阵,避免矩阵求逆,可以采用以下形式:
J
⋅
(
x
n
−
x
n
+
1
)
=
f
(
x
n
)
\bm {J} \cdot (\bm {x}_n - \bm {x}_{n+1}) = f(\bm {x}_n)
J⋅(xn−xn+1)=f(xn)
这时,一次迭代可以看成解一个线性方程组。
牛顿迭代法演示可以查看这里