Trapezoidi e Simpson

1
2
3
N=100
x = linspace(0,2,N+1)
y = x**4 - 2*x + 1

1
2
3
h = abs( x[1] - x[0] )
trap = 0.5*h*(y[0] + y[-1])
trap += sum(y[1:-1]) * h

1
2
3
simp = y[0] + y[-1]
simp += 4*sum(y[1:-1:2]) + 2*sum(y[2:-1:2])
simp *= h/3 

1
2
esatto = 4.4
print "Valore esatto, Trapezoidi, 
Simpson: ", esatto, trap, simp

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
esatto = 4.4
Nmax=1000
trap_err=zeros((Nmax-1)/2.)
simp_err=zeros((Nmax-1)/2.)
i=0

for N in 
range(2,Nmax,2):
	x = linspace(0,2,N+1)
	y = x**4 - 2*x + 1
	h = abs( x[1] - x[0] )
	trap = h* ( 0.5*y[0] + 0.5*y[-1] + sum(y[1:-1]) )
	simp = h/3 *( y[0] + y[-1] + 4*sum(y[1:-1:2]) + 2*sum(y[2:-1:2]) )
	
	trap_err[i] = abs(trap - esatto)
	simp_err[i] = abs(simp - esatto)
	i=i+1

1
2
3
N_arr = range(2,Nmax,2)
loglog(N_arr, trap_err)
loglog(N_arr, simp_err)

1
2
3
4
5
6
7
8
9
deltaY = log(trap_err[350])-log(trap_err[300])
deltaX = log(N_arr[350])-log(N_arr[300])
slope_trap = deltaY / deltaX

deltaY = log(simp_err[350])-log(simp_err[300])
deltaX = log(N_arr[350])-log(N_arr[300])
slope_simp = deltaY / deltaX

print slope_simp, slope_trap