GT_MM_1415_Tema_2_Practica_1

1129 days ago by trinidad1415

Precisión y error

Aritmética de coma flotante

Sage guarda en coma flotante todos los números que estén escritos en notación decimal ($1.0$ lo guardará en coma flotante, $1$ no). No guardará los números enteros o el resultado de operar con números enteros, ni los números reales $\pi$ o $e$. Para forzar a que guarde un número en coma flotante, utilizamos el método $n()$. Por defecto, utiliza 53 bits para la mantisa, pero podemos modificarlo incluyendo el número de bits entre los paréntesis.

pi.n(), pi.n(100) 
       

Vamos a ver los errores producidos por la aritmética de coma flotante. Guardamos $1/5$ en coma flotante.

a=(1/5).n(); a 
       

Comparamos ahora la representación binaria de $1/5$ en coma flotante con la representación binaria utilizando $60$ bits, para ver que el número no se guarda exacto, sino que se comete un pequeño error.

a.str(base=2) 
       
((1/5).n(60)).str(base=2) 
       

Debido a la limitación de cifras que se pueden almacenar en la mantisa, cuando realizamos una operación (por ejemplo una suma) se produce un error. Por ejemplo, si sumamos $1$ y $2^{-52}$, nos devuelve un resultado distinto de $1$, pero si sumamos $1$ y $2^{-53}$, volvemos a obtener $1$, pues la cifra la hemos sumado en una posición que no está en la mantisa.

(1.0+2.0^-52).str(base=2), (1.0+2.0^-53).str(base=2) 
       
1.0+2.0^-52,1.0+2.0^-52==1,1.0+2.0^-53,1.0+2.0^-53==1 
       
2.0^-53 
       

Propagación de errores

Vamos a determinar qué ocurre con el error cuando aplicamos la función seno a una magnitud que tienen un cierto error.

f(x)=sin(x) 
       

Suponemos que hemos medido $1,01$ cuando el valor original era $1.0$ y hemos aplicado el seno.

f(1.0),f(1.01) 
       

Vemos el error del seno de nuestra magnitud es similar al error que teníamos, multiplicado por la derivada del seno en el punto.

f(1.0)-f(1.01) 
       
f.diff(x)(1.01)*(1.0-1.01) 
       

Probamos con un error menor:

f(1.0)-f(1.001),f.diff(x)(1.001)*(1.0-1.001) 
       

Practica:

1. Calcula la menor potencia de 10 ($10^{-n}$) que sumada a $1$ no es igual a $1$. Haz lo mismo para $0.001$ y para $1000$.

2. Hemos medido un ángulo como $0.23$ radianes con un error de $\pm 0.01$ radianes. Estimar el error del coseno de dicho ángulo.

3. (avanzado) Medimos un ángulo recto $\pi/2$ con un error de $\pm 0.01$ radianes. Estimar el error del seno de dicho ángulo.

 
       

Cálculo de los ceros de una función

Métodos de dos puntos

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.0; c 
       

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

f(a),f(c),f(b) 
       

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 
       

Paso 2

Ahora repetimos el proceso.

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

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

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

Paso 3

c=(a+b)/2.0; c f(a), f(c), f(b) 
       
a,b=c,b; print[a,b]; 
       

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] 
       

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 
       

 

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

b-a 
       

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) 
       

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 
       

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 y 1. Calcula el error cometido.

2. (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]$.

3. (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étodos de un punto

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,4),thickness=3,ymin=-3,ymax=20) 
       

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 range(1,5): 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,2.1),thickness=3) 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 
       
 
       

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 
       
g(x)=f.diff(x); g 
       

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) 
       
h(x)=f.diff(x,2);h 
       

El máximo también es obvio:

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

h(1.6) 
       

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

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

Estimación del error

Entonces, para calcular estimaciones del error, haremos:

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

Y automatizando el proceso, para 5 iteraciones:

f(x)=e^x-4 x0=1.6; print "Punto inicial: ",x0 e0=1.6-1.2; 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 
       

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) Para comprobar si el error en una solución $x1$ es menor que una cota $E$, un método sencillo es comprobar si entre $x1-E$ y $x+E$ hay un cambio de signo (si la función es continua, tendrá un cero en ese intervalo, luego el error cometido es menor que $E$). Aplica este método para resolver el primer apartado.

 
       


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