GDIDP_AMAT_16_Practica1

hace 367 días por 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,1 
       

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= 1/2
[a,b]= [0.500000000000000, 1]
c= 1/2
[a,b]= [0.500000000000000, 1]

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.

f(x)=-x+cos(x); a,b=0,1 c=(a+b)/2 
       
for i in [1..6]: 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]
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.734375000000000, 0.750000000000000]
c= 0.742187500000000
[a,b]= [0.500000000000000, 1]
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.734375000000000, 0.750000000000000]
c= 0.742187500000000

 

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..6]: 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
Iteración 6 : c= 0.734375000000000
[a,b]= [0.734375000000000, 0.750000000000000] error= 0.0156250000000000
[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
Iteración 6 : c= 0.734375000000000
[a,b]= [0.734375000000000, 0.750000000000000] error= 0.0156250000000000

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 range(1,n+1): 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 [1..4]: 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, nos acercamos mucho 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 range(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 calcular estimaciones del error, vamos restando los sucesivos términos obtenidos. El primer error lo podemos estimar buscando un intervalo donde esté la raíz (como en bisección), y los siguientes restando las aproximanciones sucesivas que conseguimos. Tenemos que modificar ligeramente al algoritmo: 

f(x)=e^x-4 x0=2; print "Punto inicial: ",x0 error=2-1; print "Error inicial:", error for i in [1..6]: x1=(x0-(f(x0)/f.diff(x)(x0))).n() ; error=x1-x0 print "Iteración ",i, ": ",x1, "Error:",abs(error); x0=x1 
       
Punto inicial:  2
Error inicial: 1
Iteración  1 :  1.54134113294645 Error: 0.458658867053549
Iteración  2 :  1.39771625526465 Error: 0.143624877681797
Iteración  3 :  1.38635934331094 Error: 0.0113569119537131
Iteración  4 :  1.38629436323119 Error: 0.0000649800797529743
Iteración  5 :  1.38629436111989 Error: 2.11129691507494e-9
Iteración  6 :  1.38629436111989 Error: 0.000000000000000
Punto inicial:  2
Error inicial: 1
Iteración  1 :  1.54134113294645 Error: 0.458658867053549
Iteración  2 :  1.39771625526465 Error: 0.143624877681797
Iteración  3 :  1.38635934331094 Error: 0.0113569119537131
Iteración  4 :  1.38629436323119 Error: 0.0000649800797529743
Iteración  5 :  1.38629436111989 Error: 2.11129691507494e-9
Iteración  6 :  1.38629436111989 Error: 0.000000000000000

 

Practica: Aplica 4 pasos del método de Newton 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 range(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 range(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.257627640875076
0.257627640875076

Lo podemos automatizar mediante un bucle:

F(x)=e^(-x)/3 x0=1.0 for i in range(1,11): x0=F(x0).n(); x0 
       
0.122626480390481
0.294864671357105
0.248210783577719
0.260065157850550
0.257000449125435
0.257789288808731
0.257586014573722
0.257638380495925
0.257624889377776
0.257628365049042
0.122626480390481
0.294864671357105
0.248210783577719
0.260065157850550
0.257000449125435
0.257789288808731
0.257586014573722
0.257638380495925
0.257624889377776
0.257628365049042

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 range(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
f(x)=-3*x+exp(-x) 
       
plot(f,(x,0,0.5))