GDIDP_AMAT_1415_Practica1

1093 days ago by etlopez

MÉTODOS PARA EL CÁLCULO DE LOS CEROS DE UNA FUNCIÓN

1. Métodos de la bisección

Comencemos definiendo una función para aplicar el método y un intervalo donde sea continua y tenga un cero. Vamos a aplicar el algoritmo para calcular la raíz de f(x)=-x+cos(x) en el intervalo [0,1].

f(x)=-x+cos(x); a,b=0.0,1.0 
       

Representamos la función para ver que tiene un cero en ese intervalo. Nótese que no haría falta representarla sino comprobar que cambia de signo en los extremos del intervalo.

plot(f(x),(x,a,b),thickness=3) 
       

Cálculo paso a paso

Paso 1

Calculamos el punto medio del intervalo

c=(a+b)/2; c 
       
0.500000000000000
0.500000000000000

Vemos cuál es el signo de f en c y comparamos con el signo en los extremos

f(a),f(c),f(b) 
       
(1.00000000000000, 0.377582561890373, -0.459697694131860)
(1.00000000000000, 0.377582561890373, -0.459697694131860)

Podemos observar que f(c) tiene signo distinto a f(b), así que la raíz estará entre c y b. Tomamos ese intervalo.

a,b=c,b; a,b 
       
(0.500000000000000, 1.00000000000000)
(0.500000000000000, 1.00000000000000)

Paso 2

Ahora repetimos el proceso.

c=(a+b)/2; c 
       
0.750000000000000
0.750000000000000
f(a),f(c),f(b) 
       
(0.377582561890373, -0.0183111311261791, -0.459697694131860)
(0.377582561890373, -0.0183111311261791, -0.459697694131860)

Ahora hemos de coger el intervalo [a,c] y volver a repetir el proceso.

a,b=a,c; print[a,b]; 
       
[0.500000000000000, 0.750000000000000]
[0.500000000000000, 0.750000000000000]

Paso 3

c=(a+b)/2; c f(a), f(c), f(b) 
       
(0.377582561890373, 0.185963119505218, -0.0183111311261791)
(0.377582561890373, 0.185963119505218, -0.0183111311261791)
a,b=c,b; print[a,b]; 
       
[0.625000000000000, 0.750000000000000]
[0.625000000000000, 0.750000000000000]

Automatizando la selección del intervalo

En el proceso anterior, todas las instrucciones se repiten salvo la elección del intervalo siguiente. Vamos a intentar automatizar este paso.

Volvamos a los datos iniciales:

f(x)=-x+cos(x); a,b=0.0,1.0 
       

Ahora hay que calcular $c$ y elegir entre los intervalos $[a,c]$ y $[c,b]$. Podemos automatizar la selección:

c=(a+b)/2; print "c=", c if f(c)*f(a)<0: b=c.n() else: a=c.n() print "[a,b]=",[a,b] 
       
c= 0.500000000000000
[a,b]= [0.500000000000000, 1.00000000000000]
c= 0.500000000000000
[a,b]= [0.500000000000000, 1.00000000000000]

Para hacer un nuevo paso basta volver a ejecutar la celda anterior. Y así sucesivamente hasta realizar tantos pasos como nos indiquen.

Automatizando la repetición

Ejecutando la celda anterior, de manera sucesiva, vamos obteniendo los intervalos del método de bisección. Podemos automatizar este proceso haciendo un bucle.

for i in [1..5]: if f(c)*f(a)<0: b=c.n() else: a=c.n() print "[a,b]=",[a,b] c=(a+b)/2; print "c=", c 
       
[a,b]= [0.500000000000000, 1.00000000000000]
c= 0.750000000000000
[a,b]= [0.500000000000000, 0.750000000000000]
c= 0.625000000000000
[a,b]= [0.625000000000000, 0.750000000000000]
c= 0.687500000000000
[a,b]= [0.687500000000000, 0.750000000000000]
c= 0.718750000000000
[a,b]= [0.718750000000000, 0.750000000000000]
c= 0.734375000000000
[a,b]= [0.500000000000000, 1.00000000000000]
c= 0.750000000000000
[a,b]= [0.500000000000000, 0.750000000000000]
c= 0.625000000000000
[a,b]= [0.625000000000000, 0.750000000000000]
c= 0.687500000000000
[a,b]= [0.687500000000000, 0.750000000000000]
c= 0.718750000000000
[a,b]= [0.718750000000000, 0.750000000000000]
c= 0.734375000000000

 

Análisis del error: podemos estimar el error cometido mediante el cálculo de la longitud de los intervalos.

b-a 
       
0.0312500000000000
0.0312500000000000

Algoritmo completo

Juntándolo todo:

f(x)=-x+cos(x); a,b=0.0,1.0 print "[a,b]=",[a,b], "error=",b-a for i in [1..5]: c=(a+b)/2; print "Iteración",i,": c=", c if f(c)*f(a)<0: b=c else: a=c print "[a,b]=",[a,b],"error=",(b-a) 
       
[a,b]= [0.000000000000000, 1.00000000000000] error= 1.00000000000000
Iteración 1 : c= 0.500000000000000
[a,b]= [0.500000000000000, 1.00000000000000] error= 0.500000000000000
Iteración 2 : c= 0.750000000000000
[a,b]= [0.500000000000000, 0.750000000000000] error= 0.250000000000000
Iteración 3 : c= 0.625000000000000
[a,b]= [0.625000000000000, 0.750000000000000] error= 0.125000000000000
Iteración 4 : c= 0.687500000000000
[a,b]= [0.687500000000000, 0.750000000000000] error= 0.0625000000000000
Iteración 5 : c= 0.718750000000000
[a,b]= [0.718750000000000, 0.750000000000000] error= 0.0312500000000000
[a,b]= [0.000000000000000, 1.00000000000000] error= 1.00000000000000
Iteración 1 : c= 0.500000000000000
[a,b]= [0.500000000000000, 1.00000000000000] error= 0.500000000000000
Iteración 2 : c= 0.750000000000000
[a,b]= [0.500000000000000, 0.750000000000000] error= 0.250000000000000
Iteración 3 : c= 0.625000000000000
[a,b]= [0.625000000000000, 0.750000000000000] error= 0.125000000000000
Iteración 4 : c= 0.687500000000000
[a,b]= [0.687500000000000, 0.750000000000000] error= 0.0625000000000000
Iteración 5 : c= 0.718750000000000
[a,b]= [0.718750000000000, 0.750000000000000] error= 0.0312500000000000

Respresentación gráfica

Vamos a representar los pasos anteriores:

f(x)=-x+cos(x); a,b=0.0,1.0 P = plot(f, a, b, color='red',thickness=3) print "[a,b]=",[a,b], "error=", b-a n=5 for i in [1..n]: c=(a+b)/2; print "Iteración",i,": c=", c if f(c)*f(a)<0: b=c.n() else: a=c.n() print "[a,b]=",[a,b],"error=",b-a h = (P.ymax() - P.ymin())/ (1.5*n) P += line([(a,h*i), (b,h*i)],thickness=3) P += line([(a,h*i-h/4), (a,h*i+h/4)],thickness=3) P += line([(b,h*i-h/4), (b,h*i+h/4)],thickness=3) P 
       
[a,b]= [0.000000000000000, 1.00000000000000] error= 1.00000000000000
Iteración 1 : c= 0.500000000000000
[a,b]= [0.500000000000000, 1.00000000000000] error= 0.500000000000000
Iteración 2 : c= 0.750000000000000
[a,b]= [0.500000000000000, 0.750000000000000] error= 0.250000000000000
Iteración 3 : c= 0.625000000000000
[a,b]= [0.625000000000000, 0.750000000000000] error= 0.125000000000000
Iteración 4 : c= 0.687500000000000
[a,b]= [0.687500000000000, 0.750000000000000] error= 0.0625000000000000
Iteración 5 : c= 0.718750000000000
[a,b]= [0.718750000000000, 0.750000000000000] error= 0.0312500000000000
[a,b]= [0.000000000000000, 1.00000000000000] error= 1.00000000000000
Iteración 1 : c= 0.500000000000000
[a,b]= [0.500000000000000, 1.00000000000000] error= 0.500000000000000
Iteración 2 : c= 0.750000000000000
[a,b]= [0.500000000000000, 0.750000000000000] error= 0.250000000000000
Iteración 3 : c= 0.625000000000000
[a,b]= [0.625000000000000, 0.750000000000000] error= 0.125000000000000
Iteración 4 : c= 0.687500000000000
[a,b]= [0.687500000000000, 0.750000000000000] error= 0.0625000000000000
Iteración 5 : c= 0.718750000000000
[a,b]= [0.718750000000000, 0.750000000000000] error= 0.0312500000000000
 
       

Practica: Aplica 4 pasos del método de la bisección a la función $f(x)=e^{-x}-sen(x)$ para encontrar una raíz entre 0 y 1. Calcula el error cometido.

 
       

 

 

2. Métodos de un punto

2.1 Método de Newton-Raphson

Apliquemos el método para calcular la raíz de la función $f(x)=e^x-4$, partiendo de $x_0=2$

f(x)=e^x-4 plot(f(x),(x,0,4),thickness=3,ymin=-3,ymax=20) 
       
c=2; 
       
c=c-(f(c)/f.diff(x)(c)) ; c.n() 
       
1.54134113294645
1.54134113294645

Evaluando la celda de manera iterada, vamos obteniendo las sucesivas aproximaciones. 

f(x)=e^x-4 c=2; for i in range(1,5): c=(c-(f(c)/f.diff(x)(c))).n() ; c 
       
1.54134113294645
1.39771625526465
1.38635934331094
1.38629436323119
1.54134113294645
1.39771625526465
1.38635934331094
1.38629436323119

Como podemos comprobar, llegaremos al valor exacto de la raíz, que no es otra que ln4

ln(4).n() 
       
1.38629436111989
1.38629436111989

Podemos representar gráficamente el método:

f(x)=e^x-4 P=plot(f(x),(x,1,2.1),thickness=3) c=2; for i in [1..5]: P+=line([(c,f(c)),(c-(f(c)/f.diff(x)(c)),0)], color='black',thickness=2)+line([(c,0),(c,f(c))], color='red',thickness=2) c=(c-(f(c)/f.diff(x)(c))).n(); P 
       
 
       

Análisis del error: para este método, podemos usar la fórmula $E_{n+1}\leq C E_n^2$, donde $$C=\frac{\max |f''(x)|}{2 \min |f'(x)|}$$

Vamos a tomar como intervalo inicial [1.2,1.6]. Entonces, $E_0=1.6-1.2$ pues es la longitud del intervalo inicial donde vamos a buscar la raíz.

e0=1.6-1.2;e0 
       
0.400000000000000
0.400000000000000
g(x)=f.diff(x); g 
       
x |--> e^x
x |--> e^x

Cálculo exacto del máximo y el mínimo

En este caso, el mínimo de la función $e^x$ en el intervalo $[1.2,1.6]$ es fácil de calcular, pues la función $e^x$ es creciente:

$\min |f'(x)|=e^{1.2}$

 

g(1.2) 
       
3.32011692273655
3.32011692273655
h(x)=f.diff(x,2);h 
       
x |--> e^x
x |--> e^x

El máximo también es obvio:

$\max |f''(x)| =e^{1.6}$

h(1.6) 
       
4.95303242439511
4.95303242439511

Por tanto, $C=\frac{e^{1.6}}{2e^{1.2}}$

C=h(1.6)/(2*g(1.2)); C.n() 
       
0.745912348820635
0.745912348820635

Aproximando numéricamente el máximo y el mínimo

También podríamos haber calculado el máximo o el mínimo usando los métodos de Sage para ello. El primer valor de gmin es el mínimo alcanzado por la función y el segundo valor es en qué punto lo alcanza.

gmin=g.find_local_minimum(1.2,1.6);gmin 
       
Traceback (click to the left of this block for traceback)
...
AttributeError: 'sage.symbolic.expression.Expression' object has no
attribute 'find_local_minimum'
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "_sage_input_68.py", line 10, in <module>
    exec compile(u'open("___code___.py","w").write("# -*- coding: utf-8 -*-\\n" + _support_.preparse_worksheet_cell(base64.b64decode("Z21pbj1nLmZpbmRfbG9jYWxfbWluaW11bSgxLjIsMS42KTtnbWlu"),globals())+"\\n"); execfile(os.path.abspath("___code___.py"))
  File "", line 1, in <module>
    
  File "/tmp/tmpX7UL9a/___code___.py", line 3, in <module>
    exec compile(u'gmin=g.find_local_minimum(_sage_const_1p2 ,_sage_const_1p6 );gmin
  File "", line 1, in <module>
    
  File "element.pyx", line 331, in sage.structure.element.Element.__getattr__ (sage/structure/element.c:2892)
AttributeError: 'sage.symbolic.expression.Expression' object has no attribute 'find_local_minimum'
hmax=h.find_local_maximum(1.2,1.6);hmax 
       
Traceback (click to the left of this block for traceback)
...
AttributeError: 'sage.symbolic.expression.Expression' object has no
attribute 'find_local_maximum'
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "_sage_input_38.py", line 10, in <module>
    exec compile(u'open("___code___.py","w").write("# -*- coding: utf-8 -*-\\n" + _support_.preparse_worksheet_cell(base64.b64decode("aG1heD1oLmZpbmRfbG9jYWxfbWF4aW11bSgxLjIsMS42KTtobWF4"),globals())+"\\n"); execfile(os.path.abspath("___code___.py"))
  File "", line 1, in <module>
    
  File "/tmp/tmpJOC4jN/___code___.py", line 3, in <module>
    exec compile(u'hmax=h.find_local_maximum(_sage_const_1p2 ,_sage_const_1p6 );hmax
  File "", line 1, in <module>
    
  File "element.pyx", line 331, in sage.structure.element.Element.__getattr__ (sage/structure/element.c:2892)
AttributeError: 'sage.symbolic.expression.Expression' object has no attribute 'find_local_maximum'
hmax[0] 
       
Traceback (click to the left of this block for traceback)
...
NameError: name 'hmax' is not defined
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "_sage_input_39.py", line 10, in <module>
    exec compile(u'open("___code___.py","w").write("# -*- coding: utf-8 -*-\\n" + _support_.preparse_worksheet_cell(base64.b64decode("aG1heFswXQ=="),globals())+"\\n"); execfile(os.path.abspath("___code___.py"))
  File "", line 1, in <module>
    
  File "/tmp/tmpobfAJG/___code___.py", line 3, in <module>
    exec compile(u'hmax[_sage_const_0 ]
  File "", line 1, in <module>
    
NameError: name 'hmax' is not defined
C=hmax[0]/(2*gmin[0]); C 
       
Traceback (click to the left of this block for traceback)
...
NameError: name 'hmax' is not defined
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "_sage_input_40.py", line 10, in <module>
    exec compile(u'open("___code___.py","w").write("# -*- coding: utf-8 -*-\\n" + _support_.preparse_worksheet_cell(base64.b64decode("Qz1obWF4WzBdLygyKmdtaW5bMF0pOyBD"),globals())+"\\n"); execfile(os.path.abspath("___code___.py"))
  File "", line 1, in <module>
    
  File "/tmp/tmpLhrh46/___code___.py", line 3, in <module>
    exec compile(u'C=hmax[_sage_const_0 ]/(_sage_const_2 *gmin[_sage_const_0 ]); C
  File "", line 1, in <module>
    
NameError: name 'hmax' is not defined

Estimación del error

Entonces, para calcular estimaciones del error, haremos:

e1=C*e0^2; e1.n() 
       
0.119345975811302
0.119345975811302
e2=C*e1^2; e2.n() 
       
0.0106243741527570
0.0106243741527570

Y automatizando el proceso, para 5 iteraciones:

f(x)=e^x-4 x0=2; print "Punto inicial: ",x0 e0=2-1; print "Error inicial:", e0 for i in [1..5]: x0=(x0-(f(x0)/f.diff(x)(x0))).n() ; print "Iteración ",i, ": ",x0 e0=C*(e0^2); print "Error:",e0 
       
Punto inicial:  2
Error inicial: 1
Iteración  1 :  1.54134113294645
Error: 0.745912348820635
Iteración  2 :  1.39771625526465
Error: 0.415014615342069
Iteración  3 :  1.38635934331094
Error: 0.128473802899196
Iteración  4 :  1.38629436323119
Error: 0.0123116697232891
Iteración  5 :  1.38629436111989
Error: 0.000113063313764671
Punto inicial:  2
Error inicial: 1
Iteración  1 :  1.54134113294645
Error: 0.745912348820635
Iteración  2 :  1.39771625526465
Error: 0.415014615342069
Iteración  3 :  1.38635934331094
Error: 0.128473802899196
Iteración  4 :  1.38629436323119
Error: 0.0123116697232891
Iteración  5 :  1.38629436111989
Error: 0.000113063313764671

Practica: Aplica 4 pasos del método de la bisección a la función $f(x)=e^{-x}-sen(x)$ para encontrar una raíz entre 0 y 1. Calcula el error cometido.

 
       
 
       


2.2 Método de la secante

En este caso, partimos de la función $f(x)=e^x-4$ y dos valores iniciales, x0 y x1.

f(x)=e^x-4 x0=2.3 x1=2.0 
       

De nuevo obtenemos los valores iterando el método:

c=(x0*f(x1)-x1*f(x0))/(f(x1)-f(x0)); x0=x1; x1=c; x1.n() 
       
1.60670517038170
1.60670517038170

Y podemos automatizar el proceso usando un bucle:

f(x)=e^x-4 x0=2.3 x1=2.0 for i in [1..5]: c=(x0*f(x1)-x1*f(x0))/(f(x1)-f(x0)); x0=x1; x1=c; x1.n() 
       
1.60670517038170
1.44525009388054
1.39249605344681
1.38647519413630
1.60670517038170
1.44525009388054
1.39249605344681
1.38647519413630

Y también representar los valores obtenidos gráficamente:

f(x)=e^x-4 x0=2.3 x1=2.0 P=plot(f(x),(x,1,2.4),thickness=3) for i in [1..5]: c=((x0*f(x1)-x1*f(x0))/(f(x1)-f(x0))).n(); P+=line([(x0,f(x0)),(c,0)], color='red',thickness=2) x0=x1.n(); x1=c.n(); x1.n() P 
       

Practica: Aplica 5 pasos del método de la secante a la función $f(x)=e^{-x}-sen(x)$ para encontrar una raíz entre 0 y 1.

 
       

 

3. Método del punto fijo

Consideramos el problema f(x)=-3x+e^(-x)=0, partiendo de x0=1. La función para el método de punto fijo será F(x)=e^(-x)/3.

F(x)=e^(-x)/3 
       
x0=1.0; F(x0) 
       
0.122626480390481
0.122626480390481

El método consiste en aplicar la función F de forma sucesiva. Por ello, bastará con evaluar sucesivamente la siguiente celda:

x0=F(x0); x0.n() 
       
0.122626480390481
0.122626480390481

Lo podemos automatizar mediante un bucle:

F(x)=e^(-x)/3 x0=1.0 for i in [1..5]: x0=F(x0).n(); x0 
       
0.122626480390481
0.294864671357105
0.248210783577719
0.260065157850550
0.122626480390481
0.294864671357105
0.248210783577719
0.260065157850550

Y gráficamente:

F(x)=e^(-x)/3 g1=plot(F(x),(x,0,1),thickness=3)+plot(x,(x,0,1),color='red',thickness=2) g1=g1+line([(0,0),(0,1)], color='yellow',linestyle='--',thickness=3)+line([(0,1),(1,1)], color='yellow',linestyle='--',thickness=3)+line([(1,1),(1,0)], color='yellow',linestyle='--',thickness=3)+line([(1,0),(0,0)], color='yellow',linestyle='--',thickness=3) x0=1.0 for i in [1..5]: g1+=line([(x0,x0),(x0,F(x0))], color='green',thickness=2)+line([(x0,F(x0)),(F(x0),F(x0))], color='green',thickness=2) x0=F(x0).n(); x0 g1 
       

Podemos observar que $F(x)\in [0,1]$ para todo $x\in[0,1]$ pues la gráfica (azul) está contenida en el cuadrado amarillo. Por otra parte, vamos a ver que la derivada en valor absoluto es menor que $1$ por lo que tenemos asegurada la convergencia para todo punto inicial en $[0,1]$:

plot(F.derivative(x),(x,0,1)) 
       

Se puede observar que el máximo del valor absoluto se alcanza en $0$. Luego $K=F'(0)=1/3$ y el error cometido en tres pasos es como máximo:

K=F.diff(x)(0) (K^2*abs(F(1)-1)/(1-K)).n() 
       
0.0731144599674599
0.0731144599674599