龙贝格加速算法及python实现

1. 龙贝格加速算法及python实现

梯形求积法简单，但精度差，收敛速度慢，如何提高收敛速度以节省计算量是这里主要探讨的话题。

1.1. 返回首页

首页地址https://kz16.top/na
本文来源https://kz16.top/na/Romberg

1.2. 梯形求积与龙贝格加速理论

1.2.1. 梯形求积公式简介

考虑函数$f(x)\in C^2(a,b)$的积分
$$\tag{1} I = \int_{a}^{b} f(x) {\rm d}x$$

把区间$(-1,1)$分成$n$等分，步长为$h = \dfrac{b-a}{n}$，

这$n+1$个等距节点$x_k = a+kh,(k=0,1,\cdots n)$对应的复化梯形求积公式为

$$\tag{2} I_n = \frac{h}{2}\sum_{k=0}^{n-1} [f(x_k) + f(x_{k+1})]$$

1.2.2. 梯形求积公式误差估计

我们知道复化梯形求积公式的本质是折线代替曲线，也就是在每个子区间上，利用子区间的两个端点做线性插值。根据拉格朗日中值定理，第$k$个子区间$(x_k, x_{k+1})$上的线性插值关系为

$$f(x) = f(x_k) + \frac{f(x_{k+1} )- f(x_k)}{h} (x - x_k) + \frac{1}{2} f''(\xi_k) (x-x_k)(x-x_{k+1})$$

因此有

$$I = I_n + \sum_{k=0}^{n-1} \int_{x_k}^{x_{k+1}} \frac{1}{2} f''(\xi_k) (x-x_k)(x-x_{k+1}) {\rm d}x$$

这里$\xi_k \in (x_k, x_{k+1})$，从而有估计

$$\begin{aligned} I-I_n &= \sum_{k=0}^{n-1} \int_{x_k}^{x_{k+1}} \frac{1}{2} f''(\xi_k) (x-x_k)(x-x_{k+1}) {\rm d}x \\ & =\sum_{k=0}^{n-1} -\frac{1}{12} f''(\xi_k) h^3 \\ & = C_n h^2. \end{aligned}$$

这里$C_n = -\frac{h}{12} \sum_{k=0}^{n-1} f''(\xi_k)$。

1.2.3. 龙贝格加速原理

由于函数$f(x) \in C^2(a,b)$，根据定积分的定义，因此有
$$\lim _{h\rightarrow0}\sum_{k=0}^{n-1} \frac{h}{2} f''(\xi_k) = \int_{a}^{b} f''(x) {\rm d}x = f'(b) -f'(a)$$

即，只要$h$足够小，就有$C_n \approx C_{2n} \approx -\frac{1}{6}[f'(b) -f'(a)]$。

于是根据等式

$$\tag{3} I - I_{n} = C_nh^2，\quad I - I_{2n} = C_{2n} \left(\frac{h}{2} \right) ^2.$$

有

$$\tag{4} \dfrac{I-I_{2n} }{ I-I_{n} } = \frac{1}{4} \frac{C_{2n}}{C_n} \approx\frac{1}{4}.$$

上式移项整理得

$$\tag{5} I \approx \frac{4}{3} I_{2n} -\frac{1}{3} I_n.$$

从上述推导可以得出结论：把梯形二分前后的两个积分值$I_n$与$I_{2n}$按照(5)式作线性组合，同样可以作为真实积分$I$的逼近，容易验证，(5)式的这种逼近恰好就是Simpson求积法的积分值$S_n$，即

$$\tag{6} S_n = \frac{4}{3} I_{2n} -\frac{1}{3} I_n.$$

同样的道理，可以把Simpson求积法作类似处理，注意到Simpson求积的截断误差大致与$h^4$成正比，因此有

$$\tag{7} \dfrac{I-S_{2n} }{ I-S_{n} } \approx\frac{1}{8}， \quad I \approx \frac{16}{15} S_{2n} -\frac{1}{15} S_n.$$

上面(6)式结果恰好是相同节点的Cotes求积法的积分值$T_n$，即

$$\tag{8} T_n = \frac{16}{15} S_{2n} -\frac{1}{15} S_n.$$

同样的手续，进一步可得龙贝格公式：

$$\tag{9} R_n = \frac{64}{63} T_{2n} -\frac{1}{63} T_n.$$

1.2.4. 龙贝格加速的条件

通过上面的推导过程，我们知道(4)式后面的约等于是有限制条件的

• 步长$h$足够小
• 函数有二阶连续的导数，即$f(x)\in C^2(a,b)$.

类似我们不难推导，后面的公式(6),(8)和(9)对函数$f(x)$的光滑性要求会更高。

1.3. 龙贝格加速数值算例

1.3.1. 问题

计算积分值
$$I = \int_{0}^{1} \dfrac{\sin x}{x}{\rm d}x$$

1.3.2. python代码

import numpy as np

class RombergQuad(object):
    def __init__(self, f, a, b):
        self.f = f
        self.a = a
        self.b = b

    def trapezoid(self, n):
        h = 1.0*(self.b - self.a)/n
        x = np.linspace(self.a+h, self.b-h, n-1)
        return 0.5*h*(f(self.a) + 2*f(x).sum() + f(self.b))
    
    def simpson(self, xn, x2n):
        return (4.0*x2n - xn)/3.0

    def cotes(self, xn, x2n):
        return (16.0*x2n - xn)/15.0
 
    def romberg(self, xn, x2n):
        return (64.0*x2n - xn)/63.0
       
def f(x):
    if isinstance(x, int) and x == 0:
        return 1
    return np.sin(x)/x

rq = RombergQuad(f, 0, 1)
t0 = rq.trapezoid(1)
t1 = rq.trapezoid(2)
t2 = rq.trapezoid(4)
t3 = rq.trapezoid(8)

s1 = rq.simpson(t0, t1)
s2 = rq.simpson(t1, t2)
s3 = rq.simpson(t2, t3)

c2 = rq.cotes(s1, s2)
c3 = rq.cotes(s2, s3)

r3 = rq.romberg(c2, c3)
print(t0, t1, t2, t3)
print(s1, s2, s3)
print(c2,c3)
print(r3)
print("I=", 0.9460831)

结果为

0.9207354924039483 0.9397932848061772 0.9445135216653896 0.9456908635827014
0.9461458822735868 0.9460869339517938 0.946083310888472
0.9460830040636742 0.9460830693509172
0.9460830703872227
I= 0.9460831

1.4. 返回首页

首页地址https://kz16.top/na
本文来源https://kz16.top/na/Romberg

---

仅供参考，部分漏洞在所难免
欢迎转载，转载请指明来源，
请勿用于商业
作者：Zhou Bingzhen
2019.11.21