跳转至

Chap 4: Numerical Differentiation and Integration⚓︎

约 2156 个字 预计阅读时间 11 分钟

Numerical Differentiation⚓︎

目标：对于给定的 \(x_0\)，近似计算 \(f'(x_0)\)（即数值微分(numerical differentiation)）

微分计算公式：\(f'(x_0) = \lim\limits_{h \rightarrow 0} \dfrac{f(x_0 + h) - f(x_0)}{h}\)

  • 向后：\(f'(x_0) \approx \dfrac{f(x_0 + h) - f(x_0)}{h}\)

  • 向前：\(f'(x_0) \approx \dfrac{f(x_0) - f(x_0 - h)}{h}\)

现在我们用 \(f(x)\) 的带有插值点 \(x_0, x_0 + h\) 的拉格朗日多项式来近似表示 \(f(x)\)：

\[ \begin{align} f(x) & = \dfrac{f(x_0)(x - x_0 - h)}{x_0 - x_0 - h} + \dfrac{f(x_0 + h)(x - x_0)}{x_0 + h - x_0} \notag \\ & + \dfrac{(x - x_0)(x - x_0 - h)}{2} f''(\xi_x) \notag \end{align} \]
\[ \begin{align} f'(x) & = \dfrac{f(x_0 + h) - f(x_0)}{h} + \dfrac{2(x - x_0) - h}{2} f''(\xi_x) \notag \\ & + \dfrac{(x - x_0)(x - x_0 - h)}{2} \cdot \dfrac{d}{dx} [f''(\xi_x)] \notag \end{align} \]

因此 \(f'(x_0) = \dfrac{f(x_0 + h) - f(x_0)}{h} - \dfrac{h}{2}f''(\xi)\)

接下来用插值点为 \(\{x_0, x_1, \dots, x_n\}\) 的拉格朗日多项式来近似表示 \(f(x)\)

\[ \begin{align} f(x) & = \sum\limits_{k=0}^n f(x_k) L_k(x) + \dfrac{(x - x_0) \dots (x - x_n)}{(n+1)!} f^{(n+1)}(\xi_x) \notag \\ f'(x_j) & = \sum\limits_{k=0}^n f(x_k)L_k'(x_j) + \dfrac{f^{(n+1)}(\xi_j)}{(n+1)!} \prod\limits_{\substack{k = 0 \\ k \ne j}}^n (x_j - x_k) \notag \end{align} \]

注

  • 一般来说，更多的评估点（即这里的插值点）会带来更大的近似精度
  • 但另一方面，随着评估点的增加，舍入误差也在变大，因此数值微分是不稳定的！
例子

给定三个点 \(x_0, x_0 + h, x_0 + 2h\)，请求得关于每个点的三点公式，然后选出对于 \(f'(x)\) 而言最佳的三点公式。

\[ \begin{align} f'(x_0) & = \dfrac{1}{h}\Big[-\dfrac{3}{2} f(x_0) + 2f(x_0 + h) - \dfrac{1}{2}f(x_0 + 2h) \Big] + \dfrac{h^2}{3} f^{(3)}(\xi_0) \notag \\ f'(x_0 + h) & = \dfrac{1}{h} \Big[ -\dfrac{1}{2} f(x_0) + \dfrac{1}{2} f(x_0 + 2h) \Big] - \dfrac{h^2}{6} f^{(3)}(\xi_1) \notag \\ \Rightarrow f'(x_0) & = \dfrac{1}{2h} \Big[f(x_0 + h) - f(x_0 - h) \Big] - \dfrac{h^2}{6} f^{(3)} (\xi_1) \notag \end{align} \]

寻找近似计算 \(f''(x_0)\) 的方式。

考虑在 \(x_0\) 处的 \(f(x_0 + h), f(x_0 - h)\) 的泰勒展开式：

\[ \begin{align} f(x_0 + h) & = f(x_0) + f'(x_0) h + \dfrac{1}{2}f''(x_0)h^2 + \dfrac{1}{6} f'''(x_0) h^3 + \dfrac{1}{24} f^{(4)} (\xi_1) h^4 \notag \\ f(x_0 - h) & = f(x_0) - f'(x_0) h + \dfrac{1}{2}f''(x_0)h^2 - \dfrac{1}{6} f'''(x_0) h^3 + \dfrac{1}{24} f^{(4)} (\xi_{-1}) h^4 \notag \end{align} \]

因此：\(f''(x_0) = \dfrac{1}{h^2} [f(x_0 - h) - 2f(x_0) + f(x_0 + h)] - \dfrac{h^2}{12} f^{(4)} (\xi)\)

Elements of Numerical Integration⚓︎

目标：近似计算 \(I = \int_a^b f(x) dx\)（即数值求积(numerical quadratrue)）

思路：使用 \(f(x)\) 的拉格朗日插值多项式——从区间 \([a, b]\) 上选择一组不同的点 \(a \le x_0 < x_1 \dots < x_n \le b\)。拉格朗日多项式为 \(P_n(x) = \sum\limits_{k=0}^n f(x_k) L_k(x)\)，因此：

\[ \int_a^b f(x) dx \approx \sum\limits_{k=0}^n f(x_k) \int_a^b L_k(x) dx \]

令 \(A_k = \int_a^b L_k(x) dx = \int_a^b \prod\limits_{j \ne k} \dfrac{x - x_j}{x_k - x_j} dx\)

误差 \(R[f]\) 为：

\[ \begin{align} & R[f] \notag \\ = & \int_a^b f(x) dx - \sum\limits_{k=0}^n A_kf(x_k) \notag \\ = & \int_a^b [f(x) - P_n(x)] dx = \int_a^b R_n(x) dx \notag \\ = & \int_a^b \dfrac{f^{(n+1)}(\xi_x)}{(n+1)!} \prod\limits_{k=0}^n(x - x_k)dx \notag \end{align} \]

定义

求积公式的精度(degree of accuracy/precision) 为最大的正整数 \(n\)，使得公式对于每个 \(x^k(k = 0, 1, \dots, n)\) 都是精确的。

例子

考虑在 \([a, b]\) 上的线性插值，我们有 \(P_1(x) = \dfrac{x - b}{a - b}f(a) + \dfrac{x - a}{b - a}f(b)\)。可以得到：

  • \(A_1 = A_2 = \dfrac{b - a}{2}\)
  • \(\int_a^b f(x) dx \approx \dfrac{b - a}{2}[f(a) + f(b)]\)

请计算上述公式的精度。

考虑 \(x^k(k = 0, 1, \dots)\)：

  • \(x^0\)：\(\int_a^b 1dx = b - a = \dfrac{b - a}{2}[1 + 1]\)
  • \(x^1\)：\(\int_a^b xdx = \dfrac{b^2 - a^2}{2} = \dfrac{b - a}{2}[a + b]\)
  • \(x^2\)：\(\int_a^b x^2dx = \dfrac{b^3 - a^3}{3} = \dfrac{b - a}{2}[a^2 + b^2]\)

因此精度阶数 = 1

对于等间距的节点 \(x_i = a + ih, h = \dfrac{b - a}{n}, i = 0, 1, \dots, n\)

\[ \begin{align} A_i & = \int_{x_0}^{x_n} \prod\limits_{j \ne i} \dfrac{x - x_j}{x_i - x_j} dx \notag \\ & = \int_0^n \prod\limits_{i \ne j} \dfrac{(t - j)h}{(i - j)h} \times h dt = \dfrac{(b - a)(-1)^{n-i}}{ni!(n - i)!} \int_0^n \prod\limits_{i \ne j} (t - j) dt \notag \end{align} \]

注

科茨系数不取决于 \(f(x)\) 或 \([a, b]\)，而仅由 \(n, i\) 决定。因此我们可以从一张表中找出这些系数。上述公式称为牛顿 - 科茨公式(Newton-Cotes formula)

  • \(n = 1\)
    • \(C_0^{(1)} = C_1^{(1)} = \dfrac{1}{2}\)
    • \(\int_a^b f(x) dx \approx \dfrac{b - a}{2}[f(a) + f(b)]\)（称为梯形法则(trapezoidal rule)）
    • \(R[f] = \int_a^b \dfrac{f''(\xi_x)}{2!}(x-a)(x-b) dx = -\dfrac{1}{12}h^3f''(\xi)\)
    • \(\xi \in [a, b], h = \dfrac{b - a}{1}\)
    • 精度 = 1
  • \(n = 2\)
    • \(C_0^{(2)} = \dfrac{1}{6}, C_1^{(2)} = \dfrac{2}{3}, C_2^{(2)} = \dfrac{1}{6}\)
    • \(\int_a^b f(x) dx \approx \dfrac{b - a}{6}[f(a) + 4f(\dfrac{a+b}{2}) + f(b)]\)（称为辛普森法则(Simpson's rule)）
    • \(R[f] = -\dfrac{1}{90} h^5 f^{(4)}(\xi)\)
    • \(\xi \in (a, b), h = \dfrac{b - a}{2}\)
    • 精度 = 3
  • \(n = 3\)：辛普森 3/8 法则，精度 = 3，\(R[f] = -\dfrac{3}{80}h^5 f^{(5)}(\xi)\)
  • \(n = 4\)：科茨法则，精度 = 5，3，\(R[f] = -\dfrac{8}{945}h^7 f^{(6)}(\xi)\)

定理

对于使用 \(n+1\) 个点的牛顿 - 科茨公式，\(\exists \xi \in (a, b)\)，使得：

\[ \int_a^b f(x) dx = \sum\limits_{k=0}^n A_k f(x_k) + \dfrac{h^{n+3}f^{(n+2)}(\xi)}{(n+2)!} \int_0^n t^2(t - 1) \dots (t - n) dt \]
  • 如果 \(n\) 为偶数，那么 \(f \in C^{n+2}[a, b]\) 且 \(\int_a^b f(x) dx = \sum\limits_{k=0}^n A_k f(x_k) + \dfrac{h^{n+2}f^{(n+1)}(\xi)}{(n+1)!} \int_0^n t(t - 1) \dots (t - n) dt\)
  • 如果 \(n\) 为奇数，那么 \(f \in C^{n+1}[a, b]\)

Gaussian Quadrature⚓︎

目标：构造一个公式 \(\int_a^b w(x) f(x) dx \approx \sum\limits_{k=0}^n A_k f(x_k)\)，对于 \(n+1\) 个点而言精度为 \(2n+1\)。

思路：确定 \(2n+2\) 个未知量 \(x_0, \dots, x_n;\ A_0, \dots, A_n\)，使得公式在 \(f(x) = 1, x, x^2, \dots, x^{2n+1}\) 上都是精确的。点 \(x_0, \dots, x_n\) 被称为高斯点(Gaussian points)，这个方法被称为高斯求积(Gaussian quadrature)

例子

TBD

假设 \(\int_0^1 \sqrt{x} f(x) dx \approx A_0 f(x_0) + A_1 f(x_1)\)。

公式必须在 \(f(x) = 1, x, x^2, x^3\) 上精确表示。计算 \(\begin{cases}\frac{2}{3} = A_0 + A_1 \\ \frac{2}{5} = A_0 x_0 + A_1 x_1 \\ \frac{2}{7} = A_0 x_0^2 + A_1 x_1^2 \\ \frac{2}{9} = A_0x_0^3 + A_1 x_1^3\end{cases}\)，解得：\(\begin{cases}x_0 \approx 0.8212 \\ x_1 \approx 0.2899 \\ A_0 \approx 0.3891 \\ A_1 \approx 0.2776\end{cases}\)

定理

当且仅当 \(W(x) = \prod\limits_{k=0}^n(x - x_k)\) 与所有阶数不超过 \(n\) 的多项式正交时，\(x_0, \dots, x_n\) 是高斯点。

证明

  • 若 \(x_0, \dots, x_n\) 是高斯点，则公式 \(\int_a^b w(x) f(x) dx \approx \sum\limits_{k=0}^nA_kf(x_k)\) 的精度至少为 \(2n+1\)。那么对于任意多项式 \(P_m(x)\ (m \le n)\)，\(P_m(x) W(x)\) 的阶数不超过 \(2n+1\)。因此上述公式对于 \(P_m(x) W(x)\) 而言是精确的，也就是说：

    \[ \int_a^b w(x) P_m(x) W(x) dx = \sum\limits_{k=0}^n A_k P_m(x_k) W(x_k) = 0 \]

  • 要证明 \(x_0, \dots, x_n\) 是高斯点，我们需要证明公式对任意多项式 \(P_m(x)\ (m \le 2n + 1)\) 是精确的。令 \(P_m(x) = W(x) q(x) + r(x)\)，那么 $\int_a^n w(x) P_m(x) dx = $\int_a^n w(x) W(x) q(x) dx + \(\int_a^n w(x) r(x) dx = \sum\limits_{k=0}^n A_k r(x_k) = \sum\limits_{k=0}^n A_k P_m(x_k)\)

正交多项式的集合 \(\{\varphi_0, \varphi_1, \dots, \varphi_n, \dots\}\) 是线性独立的，且 \(\varphi_{n+1}\) 和任何多项式 \(P_m(x)\ (m \le n)\) 正交。所以，如果我们拿 \(\varphi_{n+1}\) 作为 \(W(x)\)，那么 \(\varphi_{n+1}\) 的根就是高斯点了。

例子

使用高斯求积来近似计算 \(\int_0^1 \sqrt{x} f(x) dx\)，其中 \(n = 1\)

假设 \(\int_0^1 \sqrt{x} f(x) dx \approx A_0 f(x_0) + A_1 f(x_1)\)。

  1. 构造正交多项式 \(\varphi_2\)。令 \(\varphi_0(x) = 1, \varphi_1(x) = x + a, \varphi_2(x) = x^2 + bx + c\)

    • \((\varphi_0, \varphi_1) = 0 \Rightarrow \int_0^1 \sqrt{x} (x + a) dx = 0 \Rightarrow a = - \dfrac{3}{5}\)
    • \((\varphi_0, \varphi_2) = 0 \Rightarrow \int_0^1 \sqrt{x} (x^2 + bx + c) dx = 0 \Rightarrow a = - \dfrac{10}{9}\)
    • \((\varphi_1, \varphi_2) = 0 \Rightarrow \int_0^1 \sqrt{x} (x - \dfrac{3}{5})(x + bx + c) dx = 0 \Rightarrow a = \dfrac{5}{21}\)
    • \(\therefore \varphi_2(x) = x^2 - \dfrac{10}{9}x + \dfrac{5}{21}\)

  2. 找到 \(\varphi_2\) 的两个根，作为高斯点 \(x_0, x_1\)：\(x_{0;1} = \dfrac{\frac{10}{9} \pm \sqrt{(\frac{10}{9})^2 - \frac{20}{21}}}{2}\)

  3. 因为这个公式必须在 \(f(x) = 1, x\) 上是精确的，所以我们能比较容易地求解 \(A_0, A_1\) 的线性方程组

最终结果为：\(x_0 \approc 0.8212, x_1 \approx 0.2899, A_0 \approx 0.3891, A_1 \approx 0.2776\)

现在使用上述结果来近似计算 \(\int_0^1 \sqrt{x} e^x dx\)

\[ $\int_0^1 \sqrt{x} f(x) dx$ \approc A_0 e^{x_0} + A_1 e^{x_1} = 0.3891 \times e^{0.8212} + 0.2776 \times e^{0.2899} \approx 1.2555 \]

而 \(\int_0^1 \sqrt{x} (2x-1) dx = \dfrac{2}{15}\) 是精确的

一些特殊的正交多项式：

  • 勒让德多项式(Legendre polynomials)：定义在 \([-1, 1]\) 上且 \(w(x) \equiv 1\)

    \[ P_k(x) = \dfrac{1}{2^k k!} \dfrac{d^k}{dx^k}(x^2 - 1)^k \quad \quad (P_k, P_l) = \begin{cases}0 & k \ne l \\ \dfrac{2}{2k+1} & k = l\end{cases} \]

    \(P_0 = 1, P_1 = x, (k + 1)P_{k+1} = (2k + 1)xP_k - kP_{k-1}\)

    使用 \(P_{n+1}\) 的根的公式称为高斯 - 勒让德求积公式。

  • 切比雪夫多项式(Chebyshev polynomials)：定义在 \([-1, 1]\) 上且 \(w(x) \equiv \dfrac{1}{\sqrt{1 - x^2}}\)

    \[ T_k(x) = \cos (k \times \arccos x) \]

    \(T_{n+1}\) 的根为 \(x_k = \cos \Big(\dfrac{2k + 1}{2n + 2} \pi \Big) \quad (k = 0, \dots, n)\)

    公式 \(\int_{-1}^1 \dfrac{1}{\sqrt{1 - x^2}} f(x) dx = \sum\limits_{k=0}^n A_k f(x_k)\) 称为高斯 - 切比雪夫求积公式。

评论区

如果大家有什么问题或想法，欢迎在下方留言~