PDA

View Full Version : [Python]Trapezium rule


The_ouroboros
28-11-2008, 21:30
A parte l'orrida forma(che poi si puo sistemare).
Sto cercando di creare una funzione per la regola del trapezio di integrazione numerica.
http://upload.wikimedia.org/math/e/6/7/e67531409ff1359080063e7b2531cb0b.png con http://upload.wikimedia.org/math/0/c/d/0cda790da44423eae67d667a68f117ba.png k=1,..,n



#!/usr/bin/env python

def trapezium_rule(f, a, b):
"Approximate the definite integral of f from a to b by Trapezium rule."
return (b - a) * ((f(a) + f(b))/2)

"Approximate the definite integral of f from a to b(subdivided) by Trapezium rule. n is the numbers of partitions "
def trapezium_plus(f,a,b,n):
z = (b-a)/2*n
res = z
x = {0:0}

for i in range(1, n):
tmp = a + i*((b-a)/n)
x[i]= tmp

x[1] = f(x[1])
x[n] = f(x[n])

for j in range(2, n-1):
x[j] = 2 * f(x[j])

out = 0

for k in range(1, n):
out += x[k]

res = z*out
return res

from math import sin
from math import pi
print 'BEGIN'
print 'integral of sin(x) from 0 to 2*pi:\n'
print '\n'
print trapezium_plus(sin, 0, 2*pi, 30)
print 'END'


Il risultato che mi esce è però


BEGIN
integral of sin(x) from 0 to 2*pi:



Traceback (most recent call last):
File "trapezium.py", line 37, in <module>
print trapezium_plus(sin, 0, 2*pi, 30)
File "trapezium.py", line 19, in trapezium_plus
x[n] = f(x[n])
KeyError: 30



Dove sbaglio??

Grazie

The_ouroboros
28-11-2008, 21:49
in metacodice sarebbe:


FUNCTION trapez(f,a,b,n)
z = (b-a)/2*n
FOR i=1:n
x(i)=a+i*((b-a)/n)
END
res = f(x(1))+f(x(n))
FOR i=2:n-1
res = res + f(x(i))
END
res = res*z
RETURN res

cdimauro
28-11-2008, 22:20
Prova così:
def trapez(f, a, b, n):
z = (b - a) / 2 * n
for i in xrange(1, n + 1):
x[i] = a + i * ((b - a) / n)
res = f(x[1]) + f(x[n])
for i in xrange(2, n):
res += f(x[i])
return res * z

The_ouroboros
28-11-2008, 22:32
grazie... mi ero un po impantanato in una cavolata :fagiano: :fagiano:

The_ouroboros
28-11-2008, 22:41
#!/usr/bin/env python

def trapezium_rule(f, a, b):
"Approximate the definite integral of f from a to b by Trapezium rule."
return (b - a) * ((f(a) + f(b))/2)

"Approximate the definite integral of f from a to b(subdivided) by Trapezium rule. n is the numbers of partitions "
def trapezium_plus(f,a,b,n):
z = (b-a)/2*n
res = z
x = {0:0}

for i in range(1, n+1):
x[i]= a + i*((b-a)/n)

res = f(x[1]) + f(x[n])

for j in range(2, n):
x[j] = 2 * f(x[j])

for k in range(1, n):
res += x[k]

return res * z

"by cdimauro"
def trapez(f, a, b, n):
x = {0:0}
z = (b - a) / 2 * n
for i in xrange(1, n + 1):
x[i] = a + i * ((b - a) / n)
res = f(x[1]) + f(x[n])
for i in xrange(2, n):
res += f(x[i])
return res * z

from math import sin
from math import pi

n = input ('Number of iteration?')
print 'BEGIN'
print '\nWith', n, ' iteration'
print '\nIntegral of sin(x) from 0 to 2*pi(ref from Simpson):\t 2.5648942583e-16'
print '\nIntegral of sin(x) from 0 to 2*pi(trapezium_plus):\t' , trapezium_plus(sin, 0, 2*pi, n)
print '\nIntegral of sin(x) from 0 to 2*pi(trapez):\t' , trapez(sin, 0, 2*pi, n)
print '\nEND'