几个数值分析算法

因研究一个新项目，了解到一些数值方法，记录下来，避免忘记。

前向欧拉法 (FE)

前向欧拉法可用于计算微分方程数值解。

根据导数的定义，当 $h$ 足够小时，有：

\frac{y (x_{0} + h) - y (x_{0})}{h} = y^{'} (x_{0})

稍加变换，可得：

y (x_{0} + h) = y (x_{0}) + h y^{'} (x_{0})

设所求函数 $y$ 的导数是关于 $x, y$ 的方程, 记为 $f (x, y)$ , 则有：

y (x_{n} + h) = y_{n} + h f (x_{n}, y_{n})

紧凑形式：

y_{n + 1} = y_{n} + h f (x_{n}, y_{n})

根据以上数学形式，可实现如下算法计算在 $x$ 处原函数 $y$ 的值 $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 (0) = 1$ , 可解得 $y = e^{x}$ 。我们可以取 $h = 0.001$ 用上述算法计算 $y (1)$ :

FE((x,y)=>y, 0, 1, 0.001)(1)

得到 2.716923932235896, 与真实值 2.718281828459045 还算接近。随着迭代次数增加误差会不可避免的增加，关于前向欧拉法的误差估计，可以参考数值方法相关的书籍。

梯形法（Trapezoidal, TR ）

梯形法从积分的计算出发，当 $h$ 足够小时，有：

\int_{x_{0}}^{x_{0} + h} y^{'} (x) d x = \frac{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$ 的导数是关于 $x, y$ 的方程, 记为 $f (x, y)$ , 则有：

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$ 在 $x_{0} + h$ 处的导数需要用 $y (x_{0} + h)$ 计算，而 $y (x_{0} + h)$ 是我们要计算的。这里可以先应用前向欧拉法估算，于是得到：

y (x_{n} + h) = y (x_{n}) + \frac{1}{2} [f (x_{n}, y_{n}) + f (x_{n} + h, y (x_{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} + h f (x_{n}, y_{n}))] h

根据以上数学形式，可以实现如下算法计算在 $x$ 处原函数 $y$ 的值 $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$ 用上述算法计算欧拉法中的示例:

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$ 用上述算法计算欧拉法中的示例:

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;
}

我们用它计算 $\sqrt{2}$ :

NR(x=>x*x-2, 1)

可以得到 1.4142135623730958, 与 Math.sqrt 求得的结果仅在最后一位有差异。调整 eps 和 iter 可以调整精度。牛顿迭代法收敛速度快，求解快速，但是能否求解与初始值的选择有很大关系。例如， $x^{3} - 3 x^{2} + x + 3$ 如果初始值选择为 2 则永远不收敛，迭代过程会在 1 跟 2 之间反复横跳，如下：
如果遇到反复横跳，可以在 x += dt 的计算中引入一个新参数 k 改善，如 x += k*dt。这个改善会减慢收敛速度，但会增加收敛的概率。