东南大学研究生-数值分析上机题（2023）Python 5 数值积分与数值微分

最新推荐文章于 2024-09-15 22:31:42 发布

天空的蓝耀

最新推荐文章于 2024-09-15 22:31:42 发布

阅读量1k

点赞数 24

文章标签： python

本文链接：https://blog.csdn.net/qq_54716613/article/details/141815173

版权

Romberg法计算数值积分

5.1 题目

由于课本第五章给出的上机题目（26题）涉及到重积分的运算，超出课程所学范围，经与老师沟通，本单元上机题目做课件最后一页的上机题目
(1) 给定积分 $\displaystyle I(f)=\int^{b}_{a}f(x)\text{dx}$ ，取初始步长 $h$ 及精度 $\epsilon$ ，采用Romberg求积法编写计算 $I (f)$ 的通用程序，误差不超过 $\epsilon$ ；
(2) 用所编程序计算积分
$I(f)=\int^{1}_{-1}\frac{1}{1+100x^{2}}\text{dx}$
取 $\displaystyle\epsilon=\frac{1}{2}\times 10^{-7}$ .

5.2 Python 源程序

import numpy as np  
  
def f(x):  
    return 1 / (1 + 100 * x**2)  
  
def romberg(a_, b_, epsilon_):  
    h = b_ - a_  
    tn = np.array([h / 2 * (f(a_) + f(b_))])  
    sn = np.array([])  
    cn = np.array([])  
    rn = np.array([])  
    k = 0  
    result = 0  
    while 1:  
        # 计算Tn  
        h = (b_ - a_) / 2 ** k  
        x = np.array([a_ + (i + 1 / 2) * h for i in range(0, 2 ** k)])  
        tn = np.append(tn, 1 / 2 * (tn[-1] + h * sum(f(x))))  
        # 计算Sn  
        if k >= 0:  
            sn = np.append(sn, 4 / 3 * tn[k+1] - 1 / 3 * tn[k])  
        # 计算Cn  
        if k >= 1:  
            cn = np.append(cn, 16 / 15 * sn[k] - 1 / 15 * sn[k-1])  
        # 计算Rn  
        if k >= 2:  
            rn = np.append(rn, 64 / 63 * cn[k-1] - 1 / 63 * cn[k-2])  
  
        # 判断精度  
        if k == 2:  
            if 1 / 63 * abs(cn[1] - cn[0]) <= epsilon_:  
                result = rn[0]  
                break  
        elif k > 2:  
            if 1 / 255 * abs(rn[k-2] - rn[k-3]) <= epsilon_:  
                result = rn[k-2]  
                break  
        k += 1  
    return k, result, tn, sn, cn, rn  

if __name__ == '__main__':  
    a = -1  
    b = 1  
    epsilon = 0.5 * (10 ** (-7))  
  
    N, Result, Tn, Sn, Cn, Rn = romberg(a, b, epsilon)  
    print("Tn求解次数：{:d}".format(N+1+1))  
    print("在给定误差下求得积分值为：{:.7f}".format(Result))  
  
    print("Tn")  
    print(Tn)  
    print("Sn")  
    print(Sn)  
    print("Cn")  
    print(Cn)  
    print("Rn")  
    print(Rn)

5.3 程序运行结果

Tn求解次数：9
在给定误差下求得积分值为：0.2942255
Tn
[0.01980198 1.00990099 0.54341203 0.34940516 0.29832452 0.29423983
 0.29422235 0.29422474 0.29422534]
Sn
[1.33993399 0.38791571 0.2847362  0.28129764 0.29287827 0.29421652
 0.29422553 0.29422553]
Cn
[0.32444783 0.27785757 0.28106841 0.29365031 0.29430573 0.29422614
 0.29422553]
Rn
[0.27711804 0.28111937 0.29385002 0.29431614 0.29422487 0.29422553]

Romberg算法过程如下

区间等分数n	Tn(f)	Sn(f)	Cn(f)	Rn(f)
1	0.01980198	1.33993399	0.32444783	0.27711804
2	1.00990099	0.38791571	0.27785757	0.28111937
4	0.54341203	0.2847362	0.28106841	0.29385002
8	0.34940516	0.28129764	0.29365031	0.29431614
16	0.29832452	0.29287827	0.29430573	0.29422487
32	0.29423983	0.29421652	0.29422614	0.29422553
64	0.29422235	0.29422553	0.29422553
128	0.29422474	0.29422553
256	0.29422534

5.4 总结感悟

通过对Romberg算法的编程，可以将粗糙的梯形插值 $T_{n}(f)$ 逐步加工为具有较高精度的的Simpson值 $S_{n}(f)$ 、Cotes值 $C_{n}(f)$ 、和Romberg值 $R_{n}(f)$ 。Romberg算法精度相对高，收敛的速度也相对较快，计算量小，在同样的精度要求下，大大节省计算量。
Romberg算法的编程相对简单，并不涉及到非常复杂的数据结构与程序结构，只要按照算法流程进行编写即可顺利得到结果。