GIIT_AMA_1718_T3_Tutorial 1

566 days ago by GIIT_AMA_1718

Cálculo de los ceros de una función

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.6]$.

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

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) 
       

¿El punto en el que tiene el cero está a la izquierda o a la derecha del punto medio del intervalo, $0.8$? Aunque en el dibujo parece estar a la izquierda, podemos asegurarnos viendo el valor de la función en ese punto.

c=(a+b)/2 c 
       
f(c) 
       

Como es negativo, sabemos que el cero estará a la derecha (por el Teorema de Bolzano, ya que la función es continua).

Es decir, que sabemos que la raíz está en el intervalo $[0,0.8]$. Así que podemos considerar ese intervalo y repetir el proceso tantas veces como queramos para que el intervalo sea más pequeño.

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

Ahora repetimos el proceso.

c=(a+b)/2.0 c 
       
f(c) 
       

Como es positivo, sabemos que la solución está en el intervalo $[0.4,0.8]$.

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


Practica:

1. 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.2 y 1.1.

2. Utiliza un condicional (if) para que el intervalo en el que está la raíz se escoja automáticamente.

3. Crea un bucle que realice $n$ pasos del algoritmo de la bisección (el valor de $n$ se introducirá anteriormente).

4. (avanzado) Crea una función, que reciba una función $f$, un intervalo $[a,b]$ y un número de pasos $n$ y aplique $n$ pasos de la bisección para calcular una aproximación de la raíz de $f$ en el intervalo $[a,b]$.

5. (avanzado) Crea una modificación de la función anterior para que reciba la función, el intervalo y un valor máximo de error y devuelva una aproximación de la raíz con un error menor que el valor máximo introducido. Para calcular el número de pasos se puede utilizar que $E\leq (b-a)/2^n$ y aplicando logaritmos, $n\geq log((b-a)/E)/log(2)$.


 
       

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=1.6$

f(x)=e^x-4 plot(f(x),(x,0,2)) 
       

Tomamos el punto inicial $c=1.6$.

c=1.6; 
       

Aplicamos el método.

c=c-(f(c)/f.diff(x)(c)) ; c.n() 
       

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

f(x)=e^x-4 c=1.6; for i in [1..4]: c=(c-(f(c)/f.diff(x)(c))).n() ; c 
       

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

ln(4).n() 
       

Podemos representar gráficamente el método:

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

 

 

Estimación del error mediante la diferencia entre dos valores consecutivos

Supongamos que queremos un error menor que $10^{-5}$. Una manera usual de estimarlo es ir restando a cada valor obtenido el anterior.

f(x)=e^x-4 x0=1.6; print "Punto inicial: ",x0 for i in [1..5]: xcopia=x0 x0=(x0-(f(x0)/f.diff(x)(x0))).n() ; print "Iteración ",i, ": ",x0 error=abs(x0-xcopia); print "Error:",error; print "Derivada: ",f.diff(x)(x0) 
       

A partir de los datos anteriores, podemos ver que 4 pasos serían suficientes.

Practica:

1. Aplica el método de Newton-Raphson a la función $f(x)=1-x-sen(x)$ para encontrar una raíz entre 0 y 1 con 4 cifras decimales de precisión.

2. (avanzado) Crea una función que reciba una función $f$, un valor inicial $x0$ y un número de pasos $n$ y delvuelva el resultado de aplicar $n$ pasos del método de Newton-Raphson.

3. (avanzado) Crea una función que reciba una función $f$, un valor inicial $x0$ y un error $n$ y delvuelva el resultado de aplicar $n$ pasos del método de Newton-Raphson.

4. (avanzado) Modifica el método de Newton-Raphson para obtener el método de la secante.


 
       


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) 
       

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

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 
       

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

Practica: 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$. Estima el error cometido.

 
       

Fractales (opcional)

Hemos visto que si estamos cerca de un cero simple (un cero de la función donde la derivada no se anula), entonces existe un intervalo de puntos en los que el método converge a dicho cero. Pero ¿qué forma tiene el conjunto de puntos tales que al aplicar el método obtenemos un determinado cero de la raíz? Curiosamente, tiene una forma muy compleja, de hecho, es una estructura fractal.

En vez de verlo en $\mathbb{R}$, lo vamos a dibujar en $\mathbb{C}$ (es equivalente a dibujarlo en $\mathbb{R}^2$), para poder ver mejor los fractales. 

Tomaremos como función $x^3-1$, aunque se pueden considerar otras funciones.

Calcular los fractales tiene un coste computacional bastante alto (tenemos que aplicar el método de Newton-Raphson a un montón de puntos), por lo que el código siguiente está en Cython, en lugar de Python.

%cython # cython es una versión de python que permite compilar import numpy cimport numpy # para declarar los tipos de los arrays tb tenemos que usar cimport # Nuestra función es x^3-1 (puedes probar a cambiarla) def f(double complex x): return x**3-1 # En df guardamos la derivada (hay que escribirla 'a mano') def df(double complex x): return 3*x**2 # Programamos el método de Newton-Rapson. def NR(double complex c): cdef double complex cant for i in range(15): # Permitimos como máximo 15 iteraciones cant=c if df(c)!=0: c=c-(f(c)/df(c)) if abs(c-cant) < 0.0001: #Cuando el error sea pequeño, devolveremos la parte imaginaria del número (coordenada y). return c.imag return 0.0 def NRfractal(float x0,float y0,float side,int N=200, int L=50, float R=3): '''Calcula una matriz NxN, donde el valor en cada posición será el resultado de aplicar NR al número complejo de la malla creada''' cdef double complex c, z, I cdef float delta, h cdef int j, k cdef numpy.ndarray[numpy.float64_t, ndim=2] m m = numpy.zeros((N,N), dtype=numpy.float64) I = complex(0,1) delta = side/N for j in range(N): for k in range(N): c = (x0+j*delta)+ I*(y0+k*delta) m[j,k]=NR(c) return m 
       
m=NRfractal(-3, -3, 6,1000,160) matrix_plot(matrix(m),cmap='prism',frame=false).show(figsize=(8,8))