龙贝格求积公式

利用龙贝格法求f(x)=x^3-2x*x+7x-5的积分值

2019.4.17

# -*- coding: utf-8 -*-
"""
Created on Wed Apr 17 19:29:54 2019

@author: xh216
"""
import numpy as np
a=1
b=3
N=800

R=[0.038271604938272,
   0.180599647266314,
   0.062081128747795,
   0.180599647266314,
   0.076895943562610,
   0.180599647266314,
   0.062081128747795,
   0.180599647266314,
   0.038271604938272]   #龙贝格系数
R=np.transpose(R)
K=N/8   #分8个小区间
Nk=N/K
aa=[]
tt=[]

#tt=[(i/Nk) for i in range(Nk+1)] 步长必须为整数
for i in range(8+1):
    tt.append(i/Nk)
for i in range(100+1):
    aa.append(a+i*(b-a)/K) 
x=[]
w=[]
x0=[]
tt=np.array(tt)
for k in range(100):
    x0=aa[k]+(aa[k+1]-aa[k])/2*(tt+1)
    x.extend(x0)
    w0=(aa[k+1]-aa[k])*R
    #w=[w,w0]
    w.extend(w0)

#print(x)
def f(x):
    N=len(x)
    y=np.zeros(900)
    for i in range(N):
        y[i]=x[i]*x[i]*x[i]-2*x[i]*x[i]+7*x[i]-5
    return y
y=f(x)
w=np.transpose(w)
I=np.matmul(y,w)
error=abs(I-20-2/3)
print(error)