|
SmtC: Show me the Code
Ole Peter Smith
Instituto de Matemática e Estatística
Universidade Federal de Goiás
http://www.olesmith.com.br
|
|
E se Eva tinha optado pela cobra?
Facebook.
|
Regra dos Trapézios
Aproximamos a função com sua Polinêmio Interplador de Lagrange, [;P_1(x);].
Calculamos:
[;
\mathcal{H}_0(x)
=
\int_{x_0}^{x_1} H_0(x) dx
=
-\frac{1}{x_1-x_0} \int_{x_0}^{x_1} (x-x_1) dx
=
;]
[;
-\frac{1}{x_1-x_0}
\left\{
\frac{1}{2} \left( x_1^2-x_0^2 \right)-x_1(x_1-x_0)
\right\}
=
-\frac{1}{2}(x_1+x_0)+x_1
=
\frac{1}{2}(x_1-x_0)
;]
[;
\mathcal{H}_1(x)
=
\int_{x_0}^{x_1} H_1(x) dx
=
\frac{1}{x_1-x_0} \int_{x_0}^{x_1} (x-x_0) dx
=
;]
[;
\frac{1}{x_1-x_0}
\left\{
\frac{1}{2} \left( x_1^2-x_0^2 \right)-x_0(x_1-x_0)
\right\}
=
\frac{1}{2}(x_1+x_0)-x_0
=
\frac{1}{2}(x_1-x_0)
;]
Juntando:
[;
\int_{x_0}^{x_1} f(x) dx \simeq
\frac{1}{2}(x_1-x_0)(y_0+y_1)
;]
Python Listing: ../../Code/Integration.py.
def Trapezoid(f,x0,x1):
return 0.5*(x1-x0)*(f(x0)+f(x1))
|
Para obter uma aproximação melhor, dividimos o intervalo em [;n;] subintervalos,
[;[a=x_0,x_1,\ldots,x_{n-1},x_n=b];], do mesmo comprimento:
[;
x_i=a+i\frac{x_n-x_0}{n},
~ i=0,\ldots,n
;]
Python Listing: ../../Code/Integration.py.
def Trapezoids_1(f,x0,xn,n=100):
## Repeat Trapezoid:
## Divide in [x0,x1] em n intervals and calc
h=(xn-x0)/(1.0*n)
trapez=0.0
x=x0
for i in range(n):
trapez+=Trapezoid(f,x,x+h)
x+=h
return trapez
|
Aplicamos a Regra dos Trapézios em cada um deles:
[;
\int_a^b f(x)
=
\sum_{i=0}^{n-1}
\int_{x_i}^{x_{i+1}} f(x)
\simeq
\sum_{i=0}^{n-1}
\frac{x_{i+1}-x_i}{2} \left\{ f(x_i)+f(x_{i+1}) \right\}
;]
[;
=
\frac{1}{2} f(x_0)
+
\sum_{i=1}^{n-1}f(x_i)
+\frac{1}{2} f(x_n)
;]
Python Listing: ../../Code/Integration.py.
def Trapezoids_R(f,x0,xn,n=100):
## Repeated Trapezoids:
## Divide in [x0,x12] em n intervals and calc
h=(xn-x0)/(1.0*n)
trapez=0.5*f(x0)
x=x0
for i in range(n-1):
trapez+=f(x+h)
x+=h
trapez+=0.5*f(xn)
#print "n=",n,"trap=",h*trapez
return h*trapez
|
Python Listing: test.py.
from math import *
from time import *
from Integration import *
def f1(x):
return x
def f2(x):
return x**2
def f3(x):
return x**3
def f4(x):
return x**4
def f5(x):
return x**5
def f6(x):
return e**x
fs=[f1,f2,f3,f4,f5,f6]
Analytical=[1.0/2.0,1.0/3.0,1.0/4.0,1.0/5.0,1.0/6.0,e-1.0]
a=0.0
b=1.0
for nf in range( len(fs) ):
print "**"*20,"Trapezoid Function:",nf+1,"**"*20
text=[
"n",
"a",
"\tb",
"\tAna",
"\tTrapez 1",
"r",
"\tTrapez R",
"r",
]
print "\t".join(text)
f=fs[nf]
analytical=Analytical[nf]
for n in range(1,11):
trapez_1=Trapezoids_1(f,a,b,10*n)
trapez_r=Trapezoids_R(f,a,b,10*n)
text=[
"%d" % n,
"%.6f" % a,
"%.6f" % b,
"%.6f" % analytical,
"%.6f" % trapez_1,
"%.6e" % abs( (trapez_1-analytical)/analytical ),
"%.6f" % trapez_r,
"%.6e" % abs( (trapez_r-analytical)/analytical ),
]
print "\t".join(text)
print "\n"
|
|
Messages:
0 secs.
|