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

Integração Numérica
E se Eva tinha optado pela cobra?
Facebook.
< Interpolação | Trapézios | Simpson >

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"
< Interpolação | Trapézios | Simpson >
Messages:
0 secs.