MUI_IIMA_1415_Numerico_1

1127 days ago by etlopez

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

1.1 Método 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 repetición

Podemos automatizar este proceso haciendo un bucle.

for i in range(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.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 range(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
[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 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

EJERCICIO: 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. ¿Cuántos pasos son necesarios para obtener la raíz con 10 cifras decimales?

 
       

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 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 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

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 range(1,6): 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
 
       

EJERCICIO: Aplica el método de Newton-Raphson a la función $f(x)=e^{-x}-sen(x)$ para encontrar una raíz entre 0 y 1 con 4 cifras decimales de precisión.  ¿Cuántos pasos habría que realizar para tener 10 cifras significativas?

 
       

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 
       

EJERCICIO: 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 range(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 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

EJERCICIO: Aplica 4 pasos del método del punto fijo a la función $F(x)=e^{-x/4}-sen(x/4)$ para encontrar un punto fijo partiendo de $x_0=1$ con un error menor de $10^-2$

 
       
 
       

Problema: Sea dni=0.87654321.  Un objeto se mueve con ecuación de movimiento

$$(1+(t-1)^2+sin(dni\ 2\ pi  t), 1+(t-1)^3/3+cos(dni\ 2\ pi  t)), \ t\in(0,2)$$

Calcular el instante de tiempo en el que su distancia al origen es mínima (con 10 cifras de precisión). ¿Cuánto vale esa distancia (con 2 cifras de precisión)?

a) Dibujar la curva pedida usando el comando parametric_plot. Para que 't' se considere variable, puede que haga falta poner var('t') antes de definir f(t).

b) Definir una función que en cada instante de tiempo, proporcione la distancia del objeto al origen (también es válido --incluso mejor-- usar la distancia al cuadrado).

c) Dibujar dicha función y decidir en qué intervalo estará el mínimo ("a ojo").

d) Calcular la derivada de la función anterior (cuando se anule será un mínimo o un máximo).

e) Dibujar la función derivada.

f) Aplicar un método para calcular un cero de la derivada. Calcular el instante en el que la distancia es mínima

g) Calcular la mínima distancia