數值積分原理與應用

1 理論

計算定積分時，可以通過微元法逼近，如下：

 若取值均勻，公式如下：

 通過上述公式可以看到，定積分的值可以通過取樣點函數值的線性組合得到。當取樣點已經确定時，定積分的精度取決于每個取樣點所配置設定的權值。梯形公式、Simpson公式、Cotes公式、Romberg公式正是通過調整權值來逐漸逼近積分值。

（1）複合梯形公式

（2） Simpson公式

（3）Cotes公式

（4）Romberg公式

2 實驗

借助如下公式，使用 Romberg 求積公式計算 

import numpy as np

def fun(x):
    return 4/(1+x*x)

def romberg(a,b,m): # m>=3
    p=np.zeros(m+1,dtype='int32') #p[i]=2**i
    p[0]=1;
    for i in range(1,m+1):
        p[i]=p[i-1]*2
    n=p[m]
    
    d=b-a
    f=np.random.uniform(a,b,n+1) #儲存函數值
    for i in range(n+1):
        x=a+d/n*i
        f[i]=fun(x)
    print("f=\n",f)
    
    T=np.random.uniform(a,b,m) #儲存梯度值
    T[0]=(f[0]+f[n])/2
    for i in range(1,m):
        s=0.0
        for k in np.arange(1,p[i-1]+1):
            s+=f[(2*k-1)*p[m-i]]
        T[i]=T[i-1]/2+d/p[i]*s
    print("T=\n",T)
    
    S=np.random.uniform(a,b,m-1) 
    for i in range(m-1):
        S[i]=(4*T[i+1]-T[i])/3 #Simpson公式
    print("S=\n",S)
    
    C=np.random.uniform(a,b,m-2) 
    for i in range(m-2):
        C[i]=(16*S[i+1]-S[i])/15 #Cotes公式
    print("C=\n",C)
    
    R=np.random.uniform(a,b,m-3) 
    for i in range(m-3):
        R[i]=(64*C[i+1]-C[i])/63 #Romberg公式
    print("R=\n",R)
    return R[-1]
    
s=romberg(0,1,8)
           

T=
 [3.         3.1        3.13117647 3.13898849 3.14094161 3.14142989
 3.14155196 3.14158248]
S=
 [3.13333333 3.14156863 3.1415925  3.14159265 3.14159265 3.14159265
 3.14159265]
C=
 [3.14211765 3.14159409 3.14159266 3.14159265 3.14159265 3.14159265]
R=
 [3.14158578 3.14159264 3.14159265 3.14159265 3.14159265]
           

 從運作結果可以看出，在取樣點相同的情況下，梯形公式、Simpson公式、Cotes公式、Romberg公式的精度依次增加。