GDI_AMAT_18_Practica5

277 days ago by etlopez18

Integración Numérica

Vamos a calcular numéricamente la integral $\int_0^1 1+e^{-x}\sin(4x)\,dx$.

 Método del trapecio

El método del trapecio consiste en calcular la integral como el área del trapecio formado por el eje x y los puntos $(a,f(a))$, $(b,f(b))$. Aplicamos la fórmula:

f(x)=1+exp(-x)*sin(4*x) I1=(f(0.0)+f(1.0))/2.0;I1 
       
0.860793960474483
0.860793960474483

Vamos a dibujar la función y el trapecio anterior.

a=plot(f(x),(x,0,1),color='blue',thickness=4,ymin=0,ymax=2) v=([(0,f(0)),(1,f(1))]) b=point2d(v,size=120,color='black') c=line([(1,0),(1,f(1))],color='black',thickness=2) show(a+b+c+plot(f(0)+(x-0)*(f(1)-f(0)),(x,0,1),fill='axis', color='green',fillcolor='green',thickness=4)) #show(a+b+line([(1,0),(1,f(1))],color='black',thickness=4)+line([(0,f(0)),(1,f(1))],color='red',thickness=4)) 
       

 

Método de Simpson

Por el método de Simpson, la integral se aproxima por

$$I=\frac{h}{3}\left(f(a)+4 f(x_1)+f(b) \right)$$

Es decir:

f(x)=1+exp(-x)*sin(4*x) I2=(f(0.0)+4*f(0.5)+f(1.0))/6;I2 
       

Recordemos que consiste en interpolar la función en los extremos y el punto medio y calcular la integral de dicho polinomio. El polinomio de interpolación es:

R = PolynomialRing(RDF, 'x') px=vector([0,0.5,1]);py=vector([f(0),f(0.5),f(1)]); puntos=zip(px,py); Q(x) = R.lagrange_polynomial(puntos);expand(Q) 
       
x |--> -2.7628912307723903*x^2 + 2.4844791517213567*x + 1.0
x |--> -2.7628912307723903*x^2 + 2.4844791517213567*x + 1.0

Dibujamos ahora la función, el polinomio a integrar y el área que define la integral.

a=plot(f(x),(x,0,1),color='blue',thickness=4,ymin=0,ymax=2) v=([(0,f(0)),(0.5,f(0.5)),(1,f(1))]) b=point(v,size=120,color='black') c=line([(1,0),(1,f(1))],color='black',thickness=2) d=plot(Q(x),(x,0,1),color='green',fill='axis',fillcolor='green',thickness=4) show(a+b+c+d) 
       

Veamos que la integral obtenida por el método de Simpon según la fórmula $A=\frac{h}{3}(f(a)+4f(a+h)+f(b))$ coincide con la integral del polinomio de interpolación:

Q.integrate(x,0,1).n() 
       
1.32127583226988
1.32127583226988
f.integrate(x,0,1).n() 
       
1.30825060464267
1.30825060464267

 Trapecio compuesto

Vamos a calcular la integral mediante el método del trapecio compuesto dividiendo el intervalo en 5 subintervalos.

f(x)=1+exp(-x)*sin(4*x) I3=(f(0.0)+2*f(0.2)+2*f(0.4)+2*f(0.6)+2*f(0.8)+f(1.0))/10;I3 
       
1.29252450729011
1.29252450729011

Y representamos la función y el área calculada

a=plot(f(x),(x,0,1),color='blue',thickness=4,ymin=0,ymax=2) v=[(k,f(k)) for k in [0,0.2..1]] b=point2d(v,size=60,color='black') c=line([(1,0),(1,f(1))],color='black',thickness=2) l=sum([plot(f(k)+(x-k)*(f(k+0.2)-f(k))*5,(x,k,k+0.2),fill='axis', color='green',fillcolor='green',thickness=4) for k in [0,0.2.. 0.8]]) a+b+c+l 
       

 Integración con Sage

Podemos calcular la integral exacta utilizando el comando "integrate". Si hacemos f.integrate(x,a,b), nos devuelve el valor (exacto) de la integral de la función $f$ entre $a$ y $b$.

f.integrate(x,0,1).n() 
       

En sage tenemos el comando numerical_integral que nos permite calcular la integral numérica de una función (devuelve la integral numerica y una estimación del error cometido):

numerical_integral(f,0,1) 
       

Cuadratura adaptativa

Ejemplo: Si queremos aplicar la cuadratura adpatativa para un error E

 
       
f(x)=1+e^(-x)*sin(4*x) var('a,b') S(a,b)=(b-a)/6*(f(a)+4*f((a+b)/2)+f(b)) 
       
E1=(S(0,0.5)+S(0.5,1)-S(0,1))/15; E1.n() 
       
-0.000792744419074067
-0.000792744419074067
###si E1<E, hemos terminado I=S(0,0.5)+S(0.5,1)+E1; I.n() 
       
1.30859192156470
1.30859192156470
f.integrate(x,0,1).n()-I.n() 
       
-0.000341316922027879
-0.000341316922027879
###si E1>E, E=0.0005, SUBDIVIDIMOS EL PROBLEMA E2=(S(0,0.25)+S(0.25,0.5)-S(0,0.5))/15; E2.n() 
       
-0.000135860245704533
-0.000135860245704533
I1=S(0,0.25)+S(0.25,0.5)+E2; I1.n() 
       
0.762232054049444
0.762232054049444
E3=(S(0.5,0.75)+S(0.75,1)-S(0.5,1))/15; E3.n() 
       
0.0000649346682632604
0.0000649346682632604
I2=S(0.5,0.75)+S(0.75,1)+E3; I2.n() 
       
0.546017802695266
0.546017802695266
(I1+I2).n() 
       
1.30824985674471
1.30824985674471
(I1+I2).n()-f.integrate(x,0,1).n() 
       
-7.47897958897070e-7
-7.47897958897070e-7