续：6.1 - 6.11 数值积分

复化求积公式

Newton-Cotes求积公式的精度随着求积节点数的增加而增加，但求积节点数 \(\ge8\) 时，Newton-Cotes公式数值不稳定。实际应用中，将积分区间分成若干个小区间分别求积分再求和，即复化求积公式的基本思想。

在区间 \([a,b]\) 上，取等距节点 \(x_k=a+kh,\ k=0,1,\dots,n\) ，
由定积分的区间可加性得 \(\displaystyle\int_a^bf(x)dx=\sum_{k=1}^n\int_{x_{k-1}}^{x_k}f(x)dx\) 。

  • 若在每个小区间 \(x_{k-1},x_k\) 用梯形公式，则有复化梯形公式 \(T_n\) ：
    \[ \displaystyle I=\int_a^bf(x)dx\approx T_n=\frac{h}{2}\sum_{k=1}^n\left[f(x_{k-1})+f(x_k)\right]=\frac{h}{2}\left[2\sum_{k=1}^{n-1}f(x_k)+f(a)+f(b)\right] \]
    复化梯形公式的的误差为：
    \[ \begin{array}{l}\displaystyle I-T_n=-\frac{h^3}{12}[f''(\xi_1)+\dots+f''(\xi_n)]=-\frac{h^2(b-a)}{12}f''(\eta),\quad\eta\in(a,b)\\\displaystyle \left|I-T_n\right|\le\frac{(b-a)^3}{12n^2}\max_{a\le x\le b}|f''(x)|\end{array} \]
    可知复化梯形公式收敛，且要使得误差 \(\le\varepsilon\) ，只要 \(\left|I-T_n\right|\le\varepsilon\) 或 \(\displaystyle n>\sqrt{\frac{(b-a)^3\max_{a\le x\le b}|f''(x)|}{12\varepsilon}}\) 。
  • 同理，复化Simpson公式 \(S_n\) ：
    \[ \displaystyle I=\int_a^bf(x)dx\approx S_n=\frac{h}{6}\left[4\sum_{k=1}^nf(x_{k-\frac{1}{2}})+2\sum_{k=1}^{n-1}f(x_k)+f(a)+f(b)\right] \]
    复化Simpson公式的误差为：
    \[ \begin{array}{l}\displaystyle I-S_n=-\frac{h^4(b-a)}{2880}f^{4}(\eta),\quad\eta\in(a,b)\\\displaystyle \left|I-S_n\right|\le\frac{(b-a)^5}{2880n^4}\max_{a\le x\le b}|f^{(4)}(x)|\end{array} \]
    可知收敛，且要使得误差 \(\le\varepsilon\) ，只要 \(\left|I-S_n\right|\le\varepsilon\) 或 \(\displaystyle n>\sqrt[4]{\frac{(b-a)^5\max_{a\le x\le b}|f^{(4)}(x)|}{2880\varepsilon}}\) 。
  • 同理，复化Cotes公式 \(C_n\) ：
    \[ \begin{aligned}\displaystyle I&=\int_a^bf(x)dx\approx C_n\\&=\frac{h}{90}\left\{32\sum_{k=1}^n\left[f(x_{k-\frac{3}{4}})+f(x_{k-\frac{1}{4}})\right]+12\sum_{k=1}^nf(x_{k-\frac{1}{2}})+14\sum_{k=1}^{n-1}f(x_k)+7f(a)+7f(b)\right\}\end{aligned} \]
    复化Cotes公式的误差为：
    \[ \begin{array}{l}\displaystyle I-S_n=-\frac{h^6(b-a)}{1935360}f^{6}(\eta),\quad\eta\in(a,b)\\\displaystyle \left|I-S_n\right|\le\frac{(b-a)^7}{1935360n^6}\max_{a\le x\le b}|f^{(6)}(x)|\end{array} \]
    可知收敛，且要使得误差 \(\le\varepsilon\) ，只要 \(\left|I-S_n\right|\le\varepsilon\) 或 \(\displaystyle n>\sqrt[6]{\frac{(b-a)^7\max_{a\le x\le b}|f^{(6)}(x)|}{1935360\varepsilon}}\) 。

Romberg求积公式

复化求积公式对步长有较高要求。

由复化梯形公式推导 由梯形公式 \(\begin{array}{l}\displaystyle I-T_n=-\frac{(b-a)^3}{12n^2}f''(\eta)\quad\eta\in(a,b)\\\displaystyle I-T_{2n}=-\frac{(b-a)^3}{36n^2}f''(\tilde{\eta})\quad\tilde{\eta}\in(a,b)\end{array}\) ，视 \(f''(\eta)\approx f''(\tilde{\eta})\) ，
得 \(\displaystyle I\approx\frac{4T_{2n}-T_n}{3}\) 或 \(\displaystyle I-T_{2n}\approx\frac{T_{2n}-T_n}{3}\) 。

该二式表述精度更高，代入有 \(\displaystyle\frac{4T_{2n}-T_n}{3}=\frac{h}{6}\left[4\sum_{k=1}^nf(x_{k-\frac{1}{2}})+2\sum_{k=1}^{n-1}f(x_k)+f(a)+f(b)\right]=S_n\) ，即为Simpson公式。且 \(\displaystyle T_{2n}=\frac{T_n}{2}+\frac{h}{2}\sum_{k=1}^nf(x_{k-\frac{1}{2}})\) 。

  • 由此得逐次分半的复化梯形公式的递推公式： \[ \left\{\begin{array}{l} \displaystyle T_{2^0}=T_1=\frac{b-a}{2}[f(a)+f(b)]\\ \displaystyle T_{2^k}=\frac{T_{2^{k-1}}}{2}+\frac{b-a}{2^k}\sum_{i=1}^{2^{k-1}}f(a+\frac{(2i-1)(b-a)}{2^k}),\quad k=1,2,\dots \end{array}\right. \] 且要使得 \(|I-T_{2^k}|<\varepsilon\) ，只要 \(|T_{2^k}-T_{2^{k-1}}|<3\varepsilon\) 。
  • 同理得逐次分半的复化Simpson公式的递推公式： \[ \left\{\begin{array}{l} \displaystyle T_{2^0}=T_1=\frac{b-a}{2}[f(a)+f(b)]\\ \displaystyle T_{2^k}=\frac{T_{2^{k-1}}}{2}+\frac{b-a}{2^k}\sum_{i=1}^{2^{k-1}}f(a+\frac{(2i-1)(b-a)}{2^k}),\quad k=1,2,\dots\\ \displaystyle S_{2^{k-1}}=\frac{4T_{2^k}-T_{2^{k-1}}}{3},\quad k=1,2,\dots \end{array}\right. \]

由复化Simpson公式推导 类似的推法得到 \(\displaystyle I\approx\frac{16S_{2n}-S_n}{15}\) 或 \(\displaystyle I-S_{2n}\approx\frac{S_{2n}-S_n}{15}\) 。且有 \(\displaystyle\frac{16S_{2n}-S_n}{15}=C_n\) 。

由复化Cotes公式推导

同理得 \(\displaystyle I\approx\frac{64C_{2n}-C_n}{63}\) 或 \(\displaystyle I-C_{2n}\approx\frac{C_{2n}-C_n}{63}\) 。记 \(\displaystyle R_n=\frac{64C_{2n}-C_n}{63}\) ，即为Romberg求积公式。

一般化推论

令 \(T_{2^k}^{(1)}=T_{2^k},\ T_{2^k}^{(1)}=S_{2^k},\ T_{2^k}^{(2)}=C_{2^k},\ T_{2^k}^{(3)}=R_{2^k}\) ，则有： \[ \left\{\begin{array}{l} \displaystyle T_{2^0}=T_1=\frac{b-a}{2}[f(a)+f(b)]\\ \displaystyle T_{2^k}^{(0)}=\frac{T_{2^{k-1}}}{2}+\frac{b-a}{2^k}\sum_{i=1}^{2^{k-1}}f(a+\frac{(2i-1)(b-a)}{2^k}),\quad k=1,2,\dots\\ \displaystyle T_{2^k}^{(m)}=\frac{4^mT_{2^{k+1}}^{(m-1)}-T_{2^k}^{(m-1)}}{4^m-1},\quad k=1,2,\dots,\ m=1,2,\dots \end{array}\right. \] 且要使得 \(|I-T_{2^k}^{(m)}|<\varepsilon\) ，只要 \(|T_{2^k}^{(m)}-T_{2^{k-1}}^{(m)}|<(4^{m+1}-1)\varepsilon,\ m=0,1,\dots\) 。

正交多项式

函数内积：若 \(f(x),g(x)\in C[a,b]\) ，则称 \(\displaystyle\int_a^bf(x)g(x)dx\) 为 \(f(x)\) 和 \(g(x)\) 的内积，记为： \((f,g)\) ，其满足： + \((f,g)=(g,f)\) + \((cf,g)=c(f,g)\) + \((f_1+f_2,g)=(f_1,g)+(f_2,g)\) 若 \((f,g)=0\) ，则称 \(f(x)\) 与 \(g(x)\) 正交，记为 \(f\perp g\) 。

利用内积可定义函数的平方模 \(\displaystyle\|f\|_2=\sqrt{(f,f)}=\sqrt{\int_a^bf^2(x)dx}\) ，其满足： + \(\|f\|_2\ge0,\ \|f\|_2=0\Leftrightarrow f(x)=0\) + \(\|cf\|_2=|c|\|f\|_2\) + \(\|f,g\|_2\le\|f\|_2+\|g\|_2\) + \(\|(f,g)\|_2\le\|f\|_2\|g\|_2\)

常引进加权形式定义： \(\displaystyle(f,g)=\int_a^b\rho(x)f(x)g(x)dx\) ，则 \(\displaystyle\|f\|_2=\sqrt{\int_a^b\rho(x)f^2(x)dx}\) 。

定理1：若 \(f_0(x),f_1(x),\dots,f_n(x)\) 为 \(C[a,b]\) 上的一组线性无关函数，则由 \(f_k(x)\) 线性组合可得到 \(C[a,b]\) 上的一组两两正交的函数组 \(g_0(x),g_1(x),\dots,g_n(x)\) ，单位化（平方模变为 \(1\) ）为规范正交组 \(e_0(x),r_1(x),\dots,r_n(x)\) 。（Schemite正交化）

将 \(P_n\) 上由线性无关函数 \(1,x,x^2,\dots,x^n\) 经过Schemite正交化得到的多项式 \(p_0(x),p_1(x),\dots,p_n(x)\) 称为 \([a,b]\) 上的正交多项式。

若 \(p_0(x),p_1(x),\dots,p_n(x)\) 为 \([a,b]\) 上权函数为 \(\rho(x)\) 的正交多项式，则满足： + \(p_k(x)\) 时首项系数不为零的 \(k\) 次多项式 + \(p_0(x),p_1(x),\dots,p_n(x)\) 构成 \(P_n\) 上的一组正交基 + \(p_n(x)\) 与不高于 \(n-1\) 次的多项式正交, \(p_n(x)\perp P_{n-1}\) + 方程 \(p_n(x)=0\) 在 \([a,b]\) 上有 \(n\) 个单根 + 方程 \(p_{n-1}(x)=0\) 的根与方程 \(p_n(x)=0\) 在 \([a,b]\) 上交错分布

常用正交多项式系

  1. Legendre多项式 \(\displaystyle L_n(x)=\frac{1}{2^nn!}\frac{d^n}{dx^n}(x^2-1)^n\quad x\in[-1,1],\ n=0,1,\dots\) ，权函数 \(\displaystyle\rho(x)=1\)
    • \(\displaystyle (L_m,L_n)=\left\{\begin{aligned}&0&m\neq n\\&\frac{2}{2n+1}&m=n\end{aligned}\right.\)
    • 有三项递推关系 \(\left\{\begin{array}{l}\displaystyle (n+1)L_{n+1}(x)=(2n+1)xL_n(x)-nL_{n-1}(x),\quad n\ge1\\L_0(x)=1,\ L_1(x)=x\end{array}\right.\)
  2. Chebyshev多项式 \(T_n(x)=\cos(n\arccos x)\quad x\in[-1,1],\ n=0,1,\dots\) ，权函数 \(\displaystyle\rho(x)=\frac{1}{\sqrt{1-x^2}}\)
    • \(\displaystyle (T_m,T_n)=\left\{\begin{array}{l}0\quad\quad m\neq n\\\pi\quad\quad m=n=0\\\pi/2\quad m=n\neq0\end{array}\right.\)
    • 有三项递推关系 \(\left\{\begin{array}{l}\displaystyle T_{n+1}(x)=2xT_n(x)-nT_{n-1}(x),\quad n=1,2,\dots\\T_0(x)=1,\ T_1(x)=x\end{array}\right.\)
    • \(T_n(x)\) 在 \([-1,1]\) 上的 \(n\) 个零点为 \(\displaystyle x_k^{(n)}=\cos\frac{2k-1}{2n}\pi,\quad k=1,2,\dots,n\)
  3. Laguere多项式 \(\displaystyle L_n(x)=e^x\frac{d^n}{dx^n}(x^ne^{-x}),\quad0<x<+\infty,\quad n=0,1,2,\dots\) ，权函数 \(\rho(x)=e^{-x}\)
    • \(\displaystyle (L_m,L_n)=\left\{\begin{aligned}&0&m\neq n\\&(n!)^2&m=n\end{aligned}\right.\)
    • 有三项递推关系 \(\left\{\begin{array}{l}\displaystyle L_{n+1}(x)=(2n+1-x)L_n(x)-n^2L_{n-1}(x),\quad n=1,2,\dots\\L_0(x)=1,\ L_1(x)=1-x\end{array}\right.\)
  4. Hermite多项式 \(\displaystyle H_n(x)=(-1)^ne^{x^2}\frac{d^n}{dx^n}(e^{-x^2}),\quad-\infty<x<+\infty,\quad n=0,1,2,\dots\) ，权函数 \(\rho(x)=e^{-x^2}\)
    • \(\displaystyle (L_m,L_n)=\left\{\begin{aligned}&0&m\neq n\\&2^nn!\pi&m=n\end{aligned}\right.\)
    • 有三项递推关系 \(\left\{\begin{array}{l}\displaystyle H_{n+1}(x)=2xH_n(x)-nH_{n-1}(x),\quad n\ge1\\H_0(x)=1,\ H_1(x)=2x\end{array}\right.\)

Gauss型求积公式

对一个求积公式而言，若不固定节点位置，节点数不变的情况下，代数精度如何提高。

定理1：区间 \([a,b]\) 上权函数为 \(\rho(x)\) 的具有 \(n\) 个节点的数值积分公式代数精度不超过 \(2n-1\) 次。

Gauss型求积公式：使求积公式具有 \(2n-1\) 次代数精度的节点 \(x_1,x_2,\dots,x_n\) 称为Gauss点，此时的插值型求积公式称为Gauss型求积公式： \(I\approx\displaystyle\sum_{i=1}^nA_if(x_i)\) 。

定理2：取区间 \([a,b]\) 上权函数为 \(\rho(x)\) 的正交多项式 \(p_n(x)\) 的 \(n\) 个零点 \(x_1,x_2,\dots,x_n\) 恰为Gauss点。

因此构造Gauss型求积公式的方法为：

  • 求出区间 \([a,b]\) 上权函数为 \(\rho(x)\) 的正交多项式 \(p_n(x)\) ；
  • 求出 \(p_n(x)\) 的 \(n\) 个零点；
  • 计算积分系数 \(A_i=\displaystyle\int_a^bl_i(x)\rho(x)dx\) 。

定理3：设 \(f(x)\in C^{2n}[a,b]\) ，则Gauss公式的误差为： \(\displaystyle R[f]=\int_a^bf(x)\rho(x)dx-\sum_{i=1}^nA_if(x_i)=\frac{f^{(2n)}(\eta)}{(2n)!}\int_a^b\rho(x)\omega^2(x)dx,\quad \eta\in(a,b)\)

常见Gauss型求积公式

  • Gauss-Legendre求积公式：区间 \([-1,1]\) 上权函数 \(\rho(x)=1\) 的Guass型求积公式 \(\displaystyle\int_{-1}^1f(x)dx\approx\displaystyle\sum_{i=1}^nA_if(x_i)\) ，其Gauss点为Legendre多项式的零点。可通过数学用表查询对应Gauss点和求积系数。余项为 \(\displaystyle R[f]=\frac{2^{2n+1}(n!)^4}{[(2n)!]^3(2n+1)}f^{(2n)}(\eta),\quad\eta\in(-1,1)\) 。

    利用积分变换 \(\displaystyle\int_a^bf(x)dx=\frac{b-a}{2}\int_{-1}^1(\frac{a+b}{2}+\frac{b-a}{2}t)dt,\quad(x=\frac{(a+b)+(b-a)t}{2})\)
    可用Gauss-Legendre求积公式求任意区间的数值积分，其在 \([a,b]\) 上权函数 \(\rho(x)=1\) 的求积公式为 \(\displaystyle\int_a^bf(x)dx\approx\frac{b-a}{2}\sum_{i=1}^{n}A_if(\frac{a+b}{2}+\frac{b-a}{2}x_i)\) 。余项为 \(\displaystyle R[f]=\frac{(b-a)^{2n+1}(n!)^4}{[(2n)!]^3(2n+1)}f^{(2n)}(\eta),\quad\eta\in(a,b)\) 。

  • Gauss-Laguerre求积公式：区间 \([0,+\infty)\) 上权函数 \(\rho(x)=e^{-x}\) 的Guass型求积公式 \(\displaystyle\int_0^\infty e^{-x}f(x)dx\approx\sum_{i=1}^nA_if(x_i)\) ，其Gauss点为Laguerre多项式的零点。可通过数学用表查询对应Gauss点和求积系数。余项为 \(\displaystyle R[f]=\frac{(n!)^2}{(2n)!}f^{(2n)}(\eta),\quad\eta\in(0,+\infty)\) 。

    对 \([0,+\infty)\) 上权函数 \(\rho(x)=1\) 的积分，也可构造Gauss-Laguerre求积公式（再乘一个 \(e^x\) ）： \(\displaystyle\int_0^\infty f(x)dx\approx\sum_{i=1}^nA_ie^{x_i}f(x_i)\) 。

  • Gauss-Hermite求积公式：区间 \((-\infty,+\infty)\) 上权函数 \(\rho(x)=e^{-x^2}\) 的Guass型求积公式 \(\displaystyle\int_{-\infty}^\infty e^{-x^2}f(x)dx\approx\sum_{i=1}^nA_if(x_i)\) ，其Gauss点为Herimite多项式的零点。可通过数学用表查询对应Gauss点和求积系数。余项为 \(\displaystyle R[f]=\frac{n!\sqrt\pi}{2^n(2n)!}f^{(2n)}(\eta),\quad\eta\in(-\infty,+\infty)\) 。

    同理，对 \((-\infty,+\infty)\) 上权函数 \(\rho(x)=1\) 的积分，也可构造Gauss-Hermite求积公式（再乘一个 \(e^{x^2}\) ）： \(\displaystyle\int_{-\infty}^\infty f(x)dx\approx\sum_{i=1}^nA_ie^{x_i^2}f(x_i)\) 。

6.12 - 6.13 数值微分

数值微分是指用函数值的线性组合近似函数在某点的导数值。

差商型数值微分

  • 向前差商数值微分公式： \(\displaystyle f'(x_0)\approx\frac{f(x_0+h)-f(x_0)}{h}\)

    Taylor展开得 \(\displaystyle f(x_0+h)=f(x_0)+f'(x_0)h+\frac{h^2}{2}f''(x_0+\theta h)\quad0\le\theta\le1\)
    可得误差 \(\displaystyle f'(x_0)-\frac{f(x_0+h)-f(x_0)}{h}=-\frac{h}{2}f''(x_0+\theta h)\quad0\le\theta\le1\)

  • 向后差商数值微分公式： \(\displaystyle f'(x_0)\approx\frac{f(x_0)-f(x_0-h)}{h}\)

    Taylor展开得 \(\displaystyle f(x_0-h)=f(x_0)-f'(x_0)h+\frac{h^2}{2}f''(x_0-\theta h)\quad0\le\theta\le1\)
    可得误差 \(\displaystyle f'(x_0)-\frac{f(x_0)-f(x_0-h)}{h}=\frac{h}{2}f''(x_0-\theta h)\quad0\le\theta\le1\)

  • 中心差商数值微分公式： \(\displaystyle f'(x_0)\approx\frac{f(x_0+h)-f(x_0-h)}{2h}\)

    可得误差 \(\begin{aligned}\displaystyle f'(x_0)-\frac{f(x_0+h)-f(x_0-h)}{2h}&=-\frac{h^2}{12}[f'''(x_0+\theta_1 h)+f'''(x_0-\theta_2 h)]\\&=-\frac{h^2}{6}f'''(x_0+\theta h)\quad\quad-1\le\theta\le1\end{aligned}\)

  • 二阶中心差商数值微分公式： \(\displaystyle f''(x_0)\approx\frac{f(x_0+h)-2f(x_0)+f(x_0-h)}{h^2}\)

    Taylor展开得 \(\begin{array}{l}\displaystyle f(x_0+h)=f(x_0)+hf'(x_0)+\frac{h^2}{2}f''(x_0)+\frac{h^3}{6}f'''(x_0)+\frac{h^4}{24}f^{(4)}(x_0+\theta_1 h)\\\displaystyle f(x_0-h)=f(x_0)-hf'(x_0)+\frac{h^2}{2}f''(x_0)-\frac{h^3}{6}f'''(x_0)+\frac{h^4}{24}f^{(4)}(x_0-\theta_2 h)\end{array}\)
    两式相加可得误差 \(\displaystyle f'(x_0)-\frac{f(x_0+h)-2f(x_0)+f(x_0-h)}{h^2}=-\frac{h^2}{12}f^{(4)}(x_0+\theta h)\quad\theta\in(-1,1)\)

从截断误差的角度看，步长 \(h\) 越小，计算越精确；但是 \(h\) 过小时，计算过程中有相近的数相减，会严重损失数值精度。实际应用中，可采用步长逐次减半的方法确定最终补偿。记 \(G(h),\ G(h/2)\) 分别为对应步长取值时的差商公式，对给定精度 \(\varepsilon>0\) ，若 \(|G(h)-G(h/2)|<\varepsilon\) 就取步长为 \(h/2\) ，反之取 \(h\) 。

插值型数值微分

建立插值多项式 \(L_n(x)\) ，取 \(L_n'(x)\) 作为 \(f'(x)\) 的近似。误差余项 \(\displaystyle f'(x_k)-L_n'(x_k)=\frac{f^{(n+1)}(\xi)}{(n+1)!}\omega_{n+1}(x_k)\) （ \(x\neq x_k\) 时难以分析）

仅考察节点处的导数值，假定所给节点等距。

  • 两点公式：线性插值函数 \(\displaystyle L_1(x)=\frac{x-x_1}{x_0-x_1}f(x_0)+\frac{x-x_0}{x_1-x_0}f(x_1)\) ，对公式两端求导，记 \(x_1-x_0=h\) ，有：
    \(\displaystyle L_1'(x)=\frac{1}{h}[f(x_1)-f(x_0)]\) ，
    得到两个两点公式： \(\displaystyle L_1'(x_0)=L_1'(x_1)=\frac{1}{h}[f(x_1)-f(x_0)]\)

  • 三点公式：二次插值函数 \(L_2(x)\) ，对其求导，记 \(x=x_0+th\) ，有：
    \(\displaystyle L_2'(x)=\frac{1}{2h}[(2t-3)f(x_0)-(4t-4)f(x_1)+(2t-1)f(x_2)]\)
    得到三个三点公式： \(\begin{array}{l}\displaystyle L_2'(x_0)=\frac{1}{2h}[-3f(x_0)+4f(x_1)-f(x_2)]\\L_2'(x_1)=\frac{1}{2h}[-f(x_0)+f(x_2)]\\L_2'(x_2)=\frac{1}{2h}[f(x_0)-4f(x_1)+3f(x_2)]\end{array}\)

    再次求导，可以得到更高阶的数值微分公式 \(\displaystyle L_2''(x)=\frac{1}{h^2}[f(x_0)-2f(x_1)+f(x_2)]\)