【计算方法】python求解数值积分｜梯形公式｜辛普森公式｜高斯求积公式

梯形公式

import numpy as np

def ff(x):
    return np.sqrt(x)*np.log(x)

def tixing_quad(ff,a,b,n):
    x_p = np.linspace(a,b,n+1) #linspace去得到右端点，arrange去不到
    h = (b-a)/n
    f = np.zeros(n+1)
    f[1:n+1] = ff(x_p[1:n+1])
    value = 0
    for i in range(n):
        value += (f[i]+f[i+1])*h/2
    err = abs(value - (-4/9))
    return value,err

print(tixing_quad(ff,0,1,8))

实验截图：

---

辛普森公式

import numpy as np

def ff(x):
    return np.sqrt(x)*np.log(x)

def simpson_quad(ff,a,b,n):
    x_p = np.linspace(a,b,n+1) #linspace去得到右端点，arrange去不到
    h = (b-a)/n
    f = np.zeros(n+1)
    f[1:n] = ff(x_p[1:n])
    f_m = ff(x_p+h/2)
    value = 0
    for i in range(n):
        value += (f[i]+f[i+1]+4*f_m[i])*h/6
    err = abs(value - (-4/9))
    return value,err

print(simpson_quad(ff,0,1,4))

实验截图：

---

高斯公式

import numpy as np

def ff(x):
    return np.sqrt(x)*np.log(x)

def gauss_quad(ff,a,b,m):
    
    if m == 1:
        x_p = 0
        A_p = 2
    elif m == 2:
        x_p = np.array([-1/np.sqrt(3), 1/np.sqrt(3)])
        A_p = np.array([1,1])
    elif m == 3:
        x_p = np.array([-np.sqrt(0.6), 0, np.sqrt(0.6)])
        A_p = np.array([5/9, 8/9, 5/9])
    else:
        print('gauss point error,only 1,2,3')
    
    value = np.sum(A_p*ff((a+b)/2+(b-a)/2*x_p)) * (b-a)/2
    err = abs(value-(-4/9))
    return value,err

print(gauss_quad(ff,0,1,3))

实验截图：

---

实验报告

方法	方法1	方法2	方法3
参数选择 (n)	8	4	3
数值	-0.4080900395195133	-0.42475203308844967	-0.45269478226195303
误差	0.03635440492493114	0.019692411355994754	0.008250337817508613