复化梯形公式、复化Simpon公式、Romberg算法（python）

1、代码

author:hewang
qq:207962168


import numpy as np
PI25DT=3.141592653589793238462643

def funEval(x):
    fx = 4.0/(1+x**2)
    return fx

def fhtan(x0,f):
    a=x0[0]
    b=x0[-1]
    fx = f(x0)
    y=2*np.sum(fx)-fx[0]-fx[-1]
    tn=((b-a)/(fx.shape[0]-1))*y/2
    return tn

def fhsimpson(x0,f):
    a=x0[0]
    b=x0[-1]
    fx = f(x0)
    f1=fx[1::2].copy()
   
    f2=fx[0::2].copy()
  
   
    if fx.shape[0] % 2== 1:
        y=4*np.sum(f1)+2*np.sum(f2)-fx[0]-fx[-1]
    else:
        y=4*np.sum(f1)+2*np.sum(f2)-fx[0]-2*fx[-1]
    sn=2*((b-a)/(fx.shape[0]-1))*y/6
    
    return sn

def romberg(x0,f):
    k=0
    n = x0.shape[0]
    fx = f(x0)
    xlb=np.zeros((4,n+3))
    for i in range(n+3):
        xlb[0,i]=fhtan(x0,f)
    for i in range(n+2):
        xlb[1,i+1]=4*xlb[0,i+1]/3-xlb[0,i]/3
    for i in range(n+1):
        xlb[2,i+2]=16*xlb[1,i+2]/15-xlb[1,i+1]/15
    while k<n:
        xlb[3,k+3]=64*xlb[2,k+3]/63-xlb[2,k+2]/63
        k+=1
    k=np.arange(n+3)
    xl=np.vstack([k,xlb])
    # print(xl.T)
    # print(xlb[3][3:])
    return xlb[3][3:]


x=np.array([0,1])
x0=np.linspace(x[0],x[1],2**10+1)

print(fhtan(x0,funEval))

print(fhsimpson(x0,funEval))

print(romberg(x0,funEval))


      

2、结果

3.141592494644074
3.141592653589793
[3.14159249 3.14159249 3.14159249 ... 3.14159249 3.14159249 3.14159249]

Process finished with exit code 0