高斯型求积公式与数值积分的龙贝格方法——Python代码

高斯型求积公式的介绍

• Gauss型求积公式

高斯-勒让德求积公式


#高斯-勒让德求积公式
from sympy import *
from scipy.special import perm,comb  #排列,组合
x,t = symbols('x,t')
#积分区间
a = 1
b = 2
#需要求积的目标函数
def f(x):
    f = 1/x
    return f

n = 2  #n次多项式正交，n越大精度越高(n=0,1,2,...)
#勒让德多项式
def L(n):
    df = diff((x ** 2 - 1) ** (n + 1), x, n + 1)
    # Python内置阶乘函数factorial
    L = 1 /2**(n+1)/factorial(n+1) * df
    return L

#高斯点x求取
def Gauss_point(n):
    x_k_list = solve(L(n))   #求得零点解集
    return x_k_list

#求积系数A
def Quadrature_coefficient(x_k_list):
    A_list = []
    for x_k in x_k_list:
        A = 2/(1-x_k**2)/(diff(L(n),x,1).subs(x,x_k))**2
        A_list.append(A)
    return A_list

result = 0
x_k_list = Gauss_point(n)
A_list = Quadrature_coefficient(Gauss_point(n))
sum = len(A_list)
#区间变换
if a == -1 and b == 1:
    for i in range(sum):
        result += (A_list[i] * f(x_k_list[i])).evalf()
    print(result)
#将求求粉公式中的区间(a,b)转换为[-1,1]
else:
    for i in range(sum):
        X = (b-a)/2 * x_k_list[i] + (a+b)/2  #区间变换
        result += (b-a)/2 * (A_list[i] * f(X)).evalf()
    print(result)



数值积分的龙贝格方法

## 导入相关包
import numpy as np
import math # 绝对值函数
from scipy import integrate # 求积分

## 求梯形值(返回用k阶复化梯形公式估计的积分)
def Trap(f,a,b,Iold,k):
    if k == 1:
        Inew = (f(a)+f(b))*(b-a)/2
        print ('二分'+str(k-1)+'次后的梯形值为'+'%.6f'%Inew)
    else:
        n = pow(2,k-2)
        h = (b-a)/n # 步长
        x = a+(h/2) #第一步的中心点
        sum_k = 0
        for i in range(n):
            sum_k = sum_k + f(x) # 求和
            x = x + h # 下一个点
        Inew = (Iold+h*sum_k)/2 # 递推公式
        print ('二分'+str(k-1)+'次后的梯形值为'+'%.6f'%Inew)
    return Inew

## 求加速值(运用理查森外推加速算法)
def Richardson(R,k):
    for i in range(k-1,0,-1):
        c = pow(2,2*(k-i))
        R[i] = (c*R[i+1]-R [i])/(c-1) # 龙贝格求积算法
    for a in sorted(R.keys(),reverse=True)[1:]: # 逆序输出
        print ('第'+str(k-1)+'次二分的第'+str(k-a)+'次加速值为'+'%.6f'%R[a])
    return R

## 龙贝格求积分
def romberg(f, a, b, eps):
    T = {} # 定义空字典
    k = 1
    print ('区间[a,b]的二分次数为：'+str(k-1))
    T[1] = Trap(f,a,b,0.0,1)
    former_R = T[1]
    while True:
        k += 1
        print ('\n区间[a,b]的二分次数为：'+str(k-1))
        # 求梯形值
        T[k] = Trap(f,a,b,T[k-1],k)
        # 求加速值
        T = Richardson(T,k)
        # 判断是否满足终止条件
        if abs(T[1] - former_R) < eps:
            return T[1]
        former_R = T[1] #最后一个值置为初始值
        
## 定义函数
def f(x):
    return x**(3/2) 

## 给定参数
a = 0 # 积分上限
b = 1 # 积分下限
eps = 10**-5  # 给定精度

## 龙贝格求积分值
I = romberg(f,a,b,eps) 
print('\n龙贝格求积分结果为: {:.6f}'.format(I))

## 计算机参考值
I_exact, Error = integrate.quad(f, a, b)
print('计算机参考值: {:.6f}'.format(I_exact))

print('相对误差(与计算机参考值相比): {:.6f}%'.format((I-I_exact)/I_exact*100))