目录
1、概述
2、代码
import numpy as np
def funcval(x):
return np.e**(1/x)
def funcval1(x):
a,b= 8755,6810
return (a**2*(np.sin(x)**2)+b**2*(np.cos(x)**2))**(1/2.)
#梯形公式
def trapezoid(a,b,n,f):
sum = 0
sum = (f(a)+f(b))*(b-a)/2
return sum
#中矩形公式
def rectangle(a,b,n,f):
sum = 0
sum = f((a+b)/2)*(b-a)
return sum
#simpson
def simpson(a,b,n,f):
sum = 0
sum = 4*f((a+b)/2)+f(a)+f(b)
return sum*(b-a)/6
#复化梯形
def ftrapezoid(a,b,n,f):
h = (b-a)/(n-1)
x=np.zeros((n,1),dtype=np.float32)
for i in range(n):
x[i] = a+i*h
sum =0
sum = f(x[0])+f(x[-1])
sum = sum+2*np.sum(f(x[1:-1]))
sum = sum*h/2
return sum
#复化simpson
def fsimpson(a,b,n,f):
if n%2 !=1:
n = n+1
h = (b-a)/(n-1)
x=np.zeros((n,1),dtype=np.float32)
for i in range(n):
x[i] = a+i*h
h = h*2
sum = f(x[0])+f(x[-1])
if x.shape[0] >= 2:
#grid
fg = f(x[2:-1:2])
#half grid
fh =f(x[1:-1:2])
sum = sum + 2*np.sum(fg)+4*np.sum(fh)
return sum*h/6
def main():
a,b =1,2
n = 20000
tr = trapezoid(a,b,n,funcval)
sip =simpson(a,b,n,funcval)
fs = fsimpson(a,b,n,funcval)
ft = ftrapezoid(a,b,n,funcval)
print(tr,ft,sip,fs)
a,b = 0,np.pi/2
fs = fsimpson(a,b,n,funcval1)
print(fs*4)
if __name__ == '__main__':
main()