几个数值分析算法

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

前向欧拉法 (FE)

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

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

\[\frac {y(x_0+h) - y(x_0)}{h} = y'(x_0) \]

稍加变换，可得：

\[y(x_0+h) = y(x_0) + hy'(x_0) \]

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

\[y(x_n + h) = y_n + hf(x_n, y_n) \]

紧凑形式：

\[y_{n+1} = y_n + hf(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+h}_{x_0} y'(x) dx = \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) + hf(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 \]

根据以上数学形式，可以实现如下算法计算在 \(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-12;
    const h = Math.sqrt(Number.EPSILON) * (1 + Math.abs(x));
    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;
}

我们用它计算 \(\sqrt2\) :

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

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

牛顿迭代法可以扩展到多函数情况，用于解方程组。这时方程组的导数为雅可比矩阵，避免矩阵求逆，可以采用以下形式：

\[\bm {J} \cdot (\bm {x}_n - \bm {x}_{n+1}) = f(\bm {x}_n) \]

这时，一次迭代可以看成解一个线性方程组。

牛顿迭代法演示

若无法看到演示可以点击这里查看。

Romberg / Richardson 法求解微分数值解

Romberg 法通过不断应用 Richardson 外推法来改善数值求解精度及加速迭代过程，可以用于求解积分与微分。以下是应用于微分的自适应算法：

// 5 点中心差分
function _Diff(f, x, h) {
  const s = 8*(f(x + h) - f(x - h)) - f(x + 2*h) + f(x - 2*h);
  return s / (12 * h);
}

function derivativeRichardson(f, x, { h0 = 0.1, maxSteps = 6, tol = 1e-12 } = {}) {
  // const R = [[centralDiff(f, x, h0)]];  // R[k][j]，k 行（更小步长），j 列（更高阶外推）
  const R1 = new Array(maxSteps).fill(0);
  const R2 = new Array(maxSteps).fill(0);
  let Rk = R1; // current Row
  let Rp = R2; // last Row
  let h = h0;
  let best = Number.POSITIVE_INFINITY;
  let bestVal = Rp[0] = _Diff(f, x, h);

  for (let k = 1; k < maxSteps; k++) {
	h /= 2;
    // 第一列：不同步长的中心差分
    Rk[0] = _Diff(f, x, h);

    // 逐列外推：误差阶次 p=4，每右移一列提高 2 阶，比例因子 2^(4j)-1
    for (let j = 1; j <= k; j++) {
      const factor = Math.pow(2, 4 * j) - 1; // 2^(4j) - 1
      Rk[j] = Rk[j - 1] + (Rk[j - 1] - Rp[j - 1]) / factor;
    }

    // 取当前最右元素作为该行最佳候选
    const candidate = Rk[k];
    const errEst = Math.abs(candidate - Rp[k - 1]);

    if (errEst < best) {
      best = errEst;
      bestVal = candidate;
    }

    // 自适应停止：相邻“最佳”变化很小则停
    if (best < tol) break;
	
	const tmp = Rk;
	Rk = Rp;
	Rp = tmp;
  }

  return bestVal;
}

posted @ 2023-03-21 22:13 1bite 阅读(207) 评论(0) 收藏举报

刷新页面返回顶部

1bite