美文网首页
12-04:Gauss求积公式

12-04:Gauss求积公式

作者: zengw20 | 来源:发表于2020-12-04 11:10 被阅读0次
    import numpy as np
    import matplotlib.pyplot as plt
    plt.rcParams['font.sans-serif'] = ['SimHei']  # 用来正常显示中文标签
    plt.rcParams['axes.unicode_minus'] = False  # 用来正常显示负号
    

    Part 1:Gauss-Legendre求积公式

    def f(x):
        return ((10/x)**2)*np.sin(10/x)
    
    # 真实值
    from scipy import integrate
    def f(x):
        return ((10/x)*8*2)*np.sin(10/x)
    v, err = integrate.quad(f, 1, 3)
    print (v,err)
    
    -1.426024756346262 1.674207356451877e-08
    

    两点

    a,b = 1,3
    A0 = (b-a)/2
    
    t = np.array([-0.5773502692,0.5773502692])
    A = np.array([1,1])
    
    X = (a+b)/2+((b-a)/2)*t
    
    F = np.sum(A*A0*f(X))
    print(F)
    
    23.399780052805777
    
    def Two(a,b):
        t = np.array([-0.5773502692,0.5773502692])
        A = np.array([1,1])
        X = (a+b)/2+((b-a)/2)*t
        F = (b-a)*np.sum(A*f(X))/2
        return F
    a,b = 1,3
    Two(a,b)
    
    23.399780052805777
    

    三点

    t = np.array([-0.7745966692,0,0.7745966692])
    A = np.array([5/9,8/9,5/9])
    
    X = (a+b)/2+((b-a)/2)*t
    
    F = np.sum(A*A0*f(X))
    print(F)
    
    10.74235460649307
    
    def Three(a,b):
        t = np.array([-0.7745966692,0,0.7745966692])
        A = np.array([5/9,8/9,5/9])
        X = (a+b)/2+((b-a)/2)*t
        F = (b-a)*np.sum(A*f(X))/2
        return F
    a,b = 1,3
    Three(a,b)
    
    10.74235460649307
    

    四点

    t = np.array([-0.8611363116,-0.3399810436,0.3399810436,0.8611363116])
    A = np.array([0.3478548451,0.6521451549,0.6521451549,0.3478548451])
    
    X = (a+b)/2+((b-a)/2)*t
    
    F = np.sum(A*A0*f(X))
    print(F)
    
    -2.212858309779213
    

    五点

    t = np.array([-0.9061798459,-0.5384693101,0,0.5384693101,0.9061798459])
    A = np.array([0.2369268851,0.4786286705,0.568888889,0.4786286705,0.2369268851])
    
    X = (a+b)/2+((b-a)/2)*t
    
    F = np.sum(A*A0*f(X))
    print(F)
    
    -2.379448453092436
    

    七点

    t = np.array([-0.9602898565,-0.79666647774,-0.5255324099,-0.1834346425,0.1834346425,0.5255324099,0.79666647774,0.9602898565])
    A = np.array([0.1012285361,0.2223810345,0.3137066459,0.3626837834,0.3626837834,0.3137066459,0.2223810345,0.1012285361])
    
    X = (a+b)/2+((b-a)/2)*t
    
    F = np.sum(A*A0*f(X))
    print(F)
    
    -1.4206298545026215
    
    def Seven(a,b):
        t = np.array([-0.9602898565,-0.79666647774,-0.5255324099,-0.1834346425,0.1834346425,0.5255324099,0.79666647774,0.9602898565])
        A = np.array([0.1012285361,0.2223810345,0.3137066459,0.3626837834,0.3626837834,0.3137066459,0.2223810345,0.1012285361])
        X = (a+b)/2+((b-a)/2)*t
        F = (b-a)*np.sum(A*f(X))/2
        return F
    a,b = 1,3
    Seven(a,b)
    
    -1.4206298545026215
    

    Part 2 :复合Gauss-Legendre求积公式

    分段二阶

    a,b,n = 1,3,100
    x = np.linspace(a,b,n)
    re = []
    for i in range(len(x)-1):
        xi = x[i]
        xii = x[i+1]
        
        t = np.array([-0.5773502692,0.5773502692])
        A = np.array([1,1])
        
        X = (xi+xii)/2+((xii-xi)/2)*t
        F = (xii-xi)*np.sum(A*f(X))/2
        re.append(F)
    
    print(np.sum(re))
    
    -1.426029291763135
    
    def TwoPlus(a,b,n):
        x = np.linspace(a,b,n)
        re = []
        for i in range(len(x)-1):
            xi = x[i]
            xii = x[i+1]
            F = Two(xi,xii)
            re.append(F)
        return np.sum(re)
    a,b,n = 1,3,100
    TwoPlus(a,b,n)
    
    -1.426029291763135
    

    分段三阶

    a,b,n = 1,3,10
    x = np.linspace(a,b,n)
    re = []
    for i in range(len(x)-1):
        xi = x[i]
        xii = x[i+1]
        
        t = np.array([-0.7745966692,0,0.7745966692])
        A = np.array([5/9,8/9,5/9])
        
        X = (xi+xii)/2+((xii-xi)/2)*t
        F = (xii-xi)*np.sum(A*f(X))/2
        re.append(F)
    
    print(np.sum(re))
    
    -1.4268535580329276
    
    def ThreePlus(a,b,n):
        x = np.linspace(a,b,n)
        re = []
        for i in range(len(x)-1):
            xi = x[i]
            xii = x[i+1]
            F = Three(xi,xii)
            re.append(F)
        return np.sum(re)
    a,b,n = 1,3,20
    ThreePlus(a,b,n)
    
    -1.426030077865446
    

    分段七阶

    def SevenPlus(a,b,n):
        x = np.linspace(a,b,n)
        re = []
        for i in range(len(x)-1):
            xi = x[i]
            xii = x[i+1]
            F = Seven(xi,xii)
            re.append(F)
        return np.sum(re)
    a,b,n = 1,3,5
    SevenPlus(a,b,n)
    
    -1.4260247638153247
    

    复合三阶随分段数变化的积分值变化

    n = np.arange(5,55,2)
    a,b = 1,3
    result = []
    for i in n:
        res = ThreePlus(a,b,i)
        result.append(res)
    result = np.array(result)
    
    plt.figure(figsize=(8,6))
    plt.plot(n,result-v)
    plt.scatter(n,result-v)
    plt.xlabel('分段数')
    plt.ylabel('误差')
    plt.title('复合三阶Gauss-Legendre求积误差与分段数关系')
    plt.grid()
    
    image.png

    相关文章

      网友评论

          本文标题:12-04:Gauss求积公式

          本文链接:https://www.haomeiwen.com/subject/upibwktx.html