【计算方法】数值积分与数值微分

Posted by ShawnD on November 16, 2020

引言

代数精度

如果当$f(x) = x^k \quad (k=0, 1, …, m)$时, 求积公式准确成立;而当$f(x) = x^{m+1}$时, 求积公式不准确成立, 称求积公式具有$m$次代数精度。

PS:

  • (等价定义)定义1中$f(x) = x^k(k=0, 1, …, m)$改为:$f(x)$是次数不超过$m$的多项式
  • 理论上,通过误差分析可以确定:用代数精度能够衡量求积公式的计算精度, 且$m$越大, 求积公式的绝对值误差$\mid E_n(f) \mid$也越小。

例:

$n=0$时, $\int_a^b f(x)dx \approx A_0 f(x_0) = (b - a)f(x_0)$

上面由积分中值定理推得, 矩形公式

若取$x_0 = a$, 为左矩形公式

若取$x_0 = b$, 为右矩形公式

$x_0 = \frac{a+b}{2}$, 为中矩形公式

$n=1$时, 若取$x_0 = a, x_1 = b, A_0 = A_1 = \frac{b - a}{2}$

\[\int_a^b f(x) dx \approx A_0f(x_0) + A_1f(x_1) = \frac{b-a}{2}[f(a) + f(b)]\]

上式为梯形公式。

将$f(x) = 1$代入$\int_a^b f(x)dx \approx A_0 f(x_0) = (b - a)f(x_0)$, 准确成立

将$f(x) = x$代入$\int_a^b f(x)dx \approx A_0 f(x_0) = (b - a)f(x_0)$, $\frac{1}{2}(b^2 - a^2) \neq (b - a)x_0$

因此左矩形公式、右矩形公式具有0次代数精度, 中矩形公式具有1次代数精度。

插值型求积公式及其性质

数值积分公式

一般形式:

\[I(f) \triangleq \int_a^b \rho(x)f(x)dx \approx \sum_{i=0}^n A_i f(x_i) \triangleq I_n(f)\]

插值型:

\[\int_a^b \rho(x)f(x)dx \approx \int_a^b \rho(x)L_n(x)dx = \sum_{i=0}^n A_if(x_i)\]

其中:

\[A_i = \int_a^b \rho(x)l_i(x)dx (i=0, 1, ..., n)\]

误差:

\[E_n(f) = \int_a^b \rho(x)\frac{f^{(n+1)(\xi)}}{(n+1)!} \prod_{i=0}^n(x - x_i)dx\]

牛顿——科特斯公式及余项估计

牛顿——科特斯(Newton-Cotes)公式

权函数$\rho(x) \equiv 1$, 将区间 [a, b] n等分, 步长$h = \frac{b - a}{n}$, 等距节点$x_i = a + ih(i = 0 \thicksim n)$

令$x = a + th$, 则Lagrange插值基函数为:

\[l_i(x) = \prod_{j=0, j \neq i}^n\frac{x - x_j}{x_i - x_j} = \prod_{j=0, j \neq i}^n \frac{t - j}{i - j} \qquad (i = 0 \thicksim n)\]

故求积系数为(记住吧):

\[A_i = \int_a^b l_i(x)dx = \frac{(-1)^{n-i}h}{i!(n-i)!}\int_0^n \prod_{j=0, j \neq i}^n(t - j)dt \qquad (i=0 \thicksim n)\]

\[C_i^{(n)} = \frac{A_i}{b-a} = \frac{A_i}{hn}= \frac{(-1)^{n-i}}{i!(n-i)!n}\int_0^n \prod_{j=0, j \neq i}^n (t - j)dt\]

称Cotes系数。

故:

\[A_i = (b - a)C_i^{(n)} \qquad (i=0, 1, 2, ..., n)\]

则求积公式

\[\int_a^b f(x) dx \approx (b - a)\sum_{i=0}^n C_i^{(n)}f(x_i)\]

称为n阶Newton-Cotes公式

Cotes系数与$f(x), [a, b]$无关, 仅与$n, i$有关。

几个低阶的Newton-Cotes公式

梯形公式

当$n=1$时, $2$个节点Cotes系数为$C_0^{(1)} = C_1^{(1)} = \frac{1}{2}$

求积公式:

\[\int_a^b f(x)dx \approx \frac{b - a}{2} [f(a) + f(b)] \triangleq T(f)\]

称为梯形公式, 代数精度为1.

辛普森(Simpson)公式

当$n = 2$时, $3$个节点, 步长$h = \frac{b - a}{2}$

Cotes系数为

\(C_0^{(2)} = \frac{(-1)^{2 - 0}}{2·0!·2!}\int_0^2(t-1)(t-2)dt = \frac{1}{6}\),

\[C_1^{(2)} = \frac{(-1)^{2 - 1}}{2·1!·(2-1)!}\int_0^2t(t-2)dt = \frac{4}{6}\] \[C_2^{(2)} = \frac{(-1)^{2 - 2}}{2·2!·(2-2)!}\int_0^2t(t-1)dt = \frac{1}{6}\]

求积公式:

\[\int_a^b f(x)dx \approx \frac{b-a}{6}[f(a) + 4f(\frac{a+b}{2}) + f(b)] \triangleq S(f)\]

称为Simpson公式,代数精度为3.

Simpson公式的几何意义是, S恰好是经过三点

\[(a, f(a)), (\frac{a+b}{2}, f(\frac{a+b}{2})), (b, f(b))\]

的抛物线所围成的曲边梯形的面积

科特斯(Cotes公式)

当$n=4$时,5个节点, 步长$h = \frac{b-a}{4}$

Cotes系数为$C_0^{(4)} = \frac{7}{90}, C_1^{(4)} = \frac{32}{90}, C_2^{(4)} = \frac{12}{90}, C_3^{(4)} = \frac{32}{90}, C_4^{(4)} = \frac{7}{90}$

求积公式:

\[\int_a^b f(x)dx \approx \frac{b - a}{90}[7f(a) + 32f(x_1) + 12f(x_2) + 32f(x_3) + 7f(b)] \triangleq C(f)\]

称为Cotes公式, 代数精度为5

例:

分别用梯形公式、Simpson公式、Cotes公式计算积分

\[I = \int_{0.5}^1 \sqrt{x} dx\]

解:

梯形公式:

当$n=1$时, $2$个节点Cotes系数为$C_0^{(1)} = C_1^{(1)} = \frac{1}{2}$

求积公式:

\[\int_a^b f(x)dx \approx \frac{b - a}{2} [f(a) + f(b)] = \frac{1 - 0.5}{2}(\sqrt{0.5} + \sqrt{1}) \approx 0.42677670\]

Simpson公式:

当$n = 2$时, $3$个节点, 步长$h = \frac{b - a}{2}$

Cotes系数为

\[C_0^{(2)} =\frac{(-1)^{2 - 0}}{2·0!·2!}\int_0^2(t-1)(t-2)dt = \frac{1}{6}\] \[C_1^{(2)} = \frac{(-1)^{2 - 1}}{2·1!·(2-1)!}\int_0^2t(t-2)dt = \frac{4}{6}\] \[C_2^{(2)} = \frac{(-1)^{2 - 2}}{2·2!·(2-2)!}\int_0^2t(t-1)dt = \frac{1}{6}\]

求积公式:

\[\int_a^b f(x)dx \approx \frac{b-a}{6}[f(a) + 4f(\frac{a+b}{2}) + f(b)] = \frac{1 - 0.5}{6}[\sqrt{0.5} + 4\sqrt{0.75} + \sqrt{1}] \approx 0.43093403\]

Cotes公式:

Cotes系数为$C_0^{(4)} = \frac{7}{90}, C_1^{(4)} = \frac{32}{90}, C_2^{(4)} = \frac{12}{90}, C_3^{(4)} = \frac{32}{90}, C_4^{(4)} = \frac{7}{90}$

求积公式:

\[\int_a^b f(x)dx \approx \frac{b - a}{90}[7f(a) + 32f(x_1) + 12f(x_2) + 32f(x_3) + 7f(b)] = \frac{0.5}{90}(7\sqrt{0.5} + 32\sqrt{0.625} + 12\sqrt{0.75} + 32\sqrt{0.875} + 7\sqrt{1}) \approx 0.43096407\]

偶数阶N-C公式的代数精度

因为N-C公式是插值型求积公式, 故$n$阶N-C公式至少具有$n$次代数精度。

定理1: 当$n$为偶数时, 牛顿-科特斯公式:

\[\int_a^b f(x) dx \approx (b-a)\sum_{i=0}^n C_i^{(n)}f(x_i)\]

至少具有$n+1$代数精度。

Newton-Cotes公式的积分余项

第一积分中值定理:

条件:

  • $f(x), g(x)$在闭区间$[a, b]$上可积
  • $g(x)$在$[a, b]$上不变号
  • $f(x)$连续

结论:在积分区间$[a, b]$上至少存在一点$\xi$, 使$\int_a^b f(x)g(x)dx = f(\xi)\int_a^b g(x)dx$

梯形公式的积分余项

\[E_T(f) = I(f) - T(f) = \int_a^b \frac{f''(\xi)}{2!}(x - a)(x - b)dx\]

设$f’‘(x) \in C^2[a, b]$, 由于$g(x) = (x - a)(x - b)$在区间$[a, b]$上不变号, 根据第一积分中值定理, 存在$\eta \in [a, b]$, 使:

\[E_T(f) = -\frac{(b-a)^3}{12}f''(\eta) = -\frac{h^3}{12}f''(\eta)\]

其中步长$h = b - a$

Simpson公式的积分余项

\[\begin{aligned} E_s(f) & = I(f) - S(f) \\ & = \int_a^b \frac{f^{(4)}(\xi)}{4!}(x - a)(x - \frac{a+b}{2})(x - b)dx \\ & = -\frac{b - a}{180} (\frac{b-a}{2})^4 f^{(4)}(\eta) \\ & = -\frac{h^5}{90}f^{(4)}(\eta), \qquad \eta \in [a, b] \end{aligned}\]

其中步长$h = \frac{b-a}{2}, f^{(4)}(x) \in C^4[a, b]$

Cotes公式的积分余项

\[\begin{aligned} E_c(f) & = I(f) - C(f) \\ & = -\frac{2(b - a)}{945}(\frac{b-a}{4})^6 f^{(6)}(\eta) \\ & = -\frac{2(b - a)}{945}h^6f^{(6)}(\eta), \qquad \eta \in [a, b] \end{aligned}\]

其中步长$h = \frac{b - a}{4}$, $f^{(6)}(x) \in C^6[a, b]$

一般的N-C公式的积分余项

定理2: 设$f(x) \in C^{n+2}[a, b]$(k是正整数)

若$n$为偶数, 则积分余项

\[R_n(f) = \frac{h^{n+3}f^{n+2}(\eta)}{(n+2)!}\int_0^n t^2(t-1)...(t - n)dt, \quad \eta \in [a, b]\]

若$n$为奇数, 则积分余项

\[R_n(f) = \frac{h^{n+2}f^{(n+1)}(\eta)}{(n+1)!}\int_0^n t(t-1)...(t-n)dt, \quad \eta \in [a, b]\]

复化(合)求积公式

常用的两种求积公式作复化求积

梯形公式$T = \frac{b-a}{2}[f(a) + f(b)]$

积分余项: $E_T = -\frac{(b-a)^3}{12}f’’(\eta), \quad \eta \in [a, b]$

代数精度: 1次代数精度

辛普森公式$S = \frac{b-a}{6}[f(a) + 4f(\frac{a+b}{2}) + f(b)]$

积分余项: $E_s = - \frac{b-a}{180}(\frac{b-a}{2})^4f^{(4)}(\eta)$

代数精度: 3次代数精度

复化梯形公式

将$[a, b]$等分成$n$个区间, $h = \frac{b-a}{n}$, $x_i = a + ih, i = 0 \thicksim n$

在每个子区间上使用梯形公式

\[\int_{x_i}^{x_{i+1}}f(x)dx = \frac{h}{2}[f(x_i) + f(x_{i+1})] - \frac{h^3}{12}f''(\xi_i), \quad x_i < \xi_i < x_{i+1}\]

则:

\[\begin{aligned} I & = \int_a^b f(x)dx = \sum_{i=0}^{n-1} \{\frac{h}{2}[f(x_i) + f(x_{i+1})] - \frac{h^3}{12}f''(\xi_i)\} \\ & = \frac{h}{2}[f(a) + 2\sum_{i=1}^{n-1}f(x_i) + f(b)] - \frac{h^3}{12}\sum_{i=0}^{n-1}f''(\xi_i) \\ & = \frac{h}{2}[f(a) + 2\sum_{i=1}^{n-1}f(a + ih) + f(b)] - \frac{b - a}{12}h^2f''(\xi) \\ & = T_n + E_{T_n} \end{aligned}\]

复化Simpson公式

将$[a, b]$分成$n$个区间, $h = \frac{b - a}{n}$, 在每个子区间$[x_k, x_{k+1}]$上使用辛普森公式, 记区间中点$x_{k+1/2}$, 得

\[\begin{aligned} I & = \int_a^b f(x)dx = \sum_{k=0}^{n-1}\int_{k=0}^{n-1}f(x)dx \\ & = \sum_{k=0}^{n-1}\{\frac{h}{6}[f(x_k) + 4f(x_{k+1/2}) + f(x_{k+1})] - \frac{h}{180}(\frac{h}{2})^4 f^{(4)}(\xi_k)\} \end{aligned}\]

\[\begin{aligned} S_n & = \frac{h}{6}\sum_{j=0}^{n-1}[f(x_k) + 4f(x_{k+1/2}) + f(x_{k+1})] \\ & = \frac{h}{6}[f(a) + 4\sum_{k=0}{n-1}f(x_{k+1/2}) + 2\sum_{k=1}^{n-1}f(x_k) + f(b)] \end{aligned}\]

为复化辛普森求积公式

例:

应用复合梯形公式计算积分:$I = \int_0^1 6e^{-x^2}dx$, 要求误差不超过$10^{-6}$, 试确定所需的步长和节点个数。

解:

令$f(x) = 6 e^{-x^2}$, 则$f’(x) = -12xe^{-x^2}$, $f’‘(x) = 12e^{-x^2}(2x^2 - 1)$, $f’’’^{x} = 24xe^{-x^3}(3 - 2x^2) \neq 0, x \in (0,1)$

由于$f’‘(x)$在$[0,1]$上为单调函数, 所以:

\[\max_{x \in [0, 1]}\mid f''(x) \mid = \max\{\mid f''(0) \mid, \mid f''(1) \mid\} = \max\{\mid -12 \mid, \mid 12e^{-1} \mid\} = 12\]

由复化梯形的误差公式:

\[\begin{aligned} E_n(f) & = -\frac{h^2(1 - 0)}{12} f''(\xi), 0 < \xi < 1 \\ & \Rightarrow \mid E_n(f) \mid \leq \frac{h^2}{12} \max_{x \in [0, 1]}\mid f''(x) \mid \end{aligned}\]

要使$\mid E_n(f) \mid \leq 10^{-6} \Rightarrow h \leq 10^{-3}$

取$h = 10^{-3}$, 由$h = \frac{b - a}{n} = \frac{1}{n} \Rightarrow n = 10^3$

故可取节点数为1001

龙贝格(Romberg)积分法

梯形法的逐次半分递推公式

将区间$[a, b]$ n等分, 节点$a = x_0 < x_1 < … < x_{n-1} < x_n = b$, 记$x_i = a + ih (i = 0 \thicksim n)$, 步长$h = \frac{b - a}{n}$, 则复化梯形公式:

\[T_n = \sum_{i=0}^{n-1} \frac{h}{2} [f(x_i) + f(x_{i+1})]\]

在每个子区间$[x_i, x_{i+1}]$上取重点$x_{i + \frac{1}{2}} = \frac{1}{2}(x_i + x_{i+1}), i=0 \thicksim n-1$, 将区间$[a, b]$半分成了$2n$等分。

每个子区间$x_i, x_{i+1}$上的积分应用复化梯形公式, 计算的

\[\frac{h}{4}[f(x_i) + 2f(x_{i + 1/2}) + f(x_{i+1})]\]

$h$是二分前的步长

各子区间上复化梯形值相加, 求得半分后的复化梯形值为

半分后:

\[T_{2n} = \frac{h}{4} \sum_{i=0}^{n-1}[f(x_i) + f(x_{i+1})] + \frac{h}{2} \sum_{i=0}^{n-1}f(x_{i+1/2})\]

半分前:

\[T_n = \sum_{i=0}^{n-1}\frac{h}{2}[f(x_i) + f(x_{i+1})]\]

\[T_{2n} = \frac{1}{2}T_n + \frac{h}{2}\sum_{i=0}^{n-1}f(x_{i + 1/2})\]

梯形法的半分递推公式, 也称为变步长的梯形公式

龙贝格(Romberg)算法

首先 根据复合梯形法求积公式的余项

\[I - T_n = -\frac{b-a}{12}h^2f''(\xi_n), \quad \xi_n \in (a, b)\] \[I - T_{2n} = -\frac{b-a}{12} (\frac{h}{2})^2 f''(\xi_{2n}), \quad \xi_{2n} \in (a, b)\]

当$f’‘(x) \in C[a, b]$, 且$n$充分大时$f’’(\xi_n) \approx f’’(\xi_{2n})$, 有:

\[I - T_{2n} \approx -\frac{b - a}{12} (\frac{h}{2})^2 f''(\xi_)\]

由逐次二分公式

\[\begin{cases} T_1 = \frac{b - a}{2}[f(a) + f(b)] \\ T_{2^k} = \frac{1}{2}T_{2^{k-1}} + \frac{b-a}{2^k} \end{cases}\]

高斯型求积公式

$[-1, 1]$上两点的高斯公式:

\[\int_{-1}^1 f(x)dx \approx f(-\frac{1}{\sqrt{3}}) + f(\frac{1}{\sqrt{3}})\]

若$t \in [-1, 1]$, 那么$x \in [a, b]$

\[x = \frac{b-a}{2}t + \frac{a+b}{2}\]

$[a, b]$上的两点高斯公式:

\[\begin{align} \int_a^b f(x)dx &= \frac{b - a}{2} \int_{-1}^1 f(\frac{b-a}{2}t + \frac{a+b}{2})dt \\ &\approx \frac{b-a}{2}[f(-\frac{b-a}{2\sqrt{3}} + \frac{a+b}{2}) + f(\frac{b-a}{2\sqrt{3}} + \frac{a+b}{2})] \end{align}\]

因此

\[A = B = \frac{b-a}{2}, x_0 = -\frac{b-a}{2\sqrt{3}} + \frac{a + b}{2}, x_1 = \frac{b - a}{2\sqrt{3}} + \frac{a+b}{2}\]