Resolución numérica de ecuaciones

Curso de Investigación de Operadores

Uno de los problemas que más se presenta en matemáticas es el de calcular la solución de una ecuación. En algunas (pocas) ocasiones, esto puede hacerse por métodos analíticos, es decir, se puede “despejar” la incógnita para encontrar el o los valores que resuelven la ecuación. En la gran mayoría de las ocasiones con algún interés práctico esto no es posible y es necesario recurrir a un método numérico que, con la ayuda de un ordenador, nos permita calcular un valor aproximado de la solución.

Método de Bisección

Antes de comenzar a desarrollar los fundamentos de este método, repasemos un resultado matemático que nos sirve de apoyo para fundamentar la validez de nuestra solución numérica.

Teorema de Bolzano: Sea \(f\) una función real continua en un intervalo cerrado \([a,b]\) con \(f(a)\) y \(f(b)\) de signos contrarios. Entonces existe al menos un punto c del intervalo abierto \((a,b)\) con \(f(c)=0\).

El método de Bisección se apoya en el teorema anterior que asegura que (bajo las condiciones expuestas en el teorema) el método es siempre convergente aunque la velocidad de convergencia es lineal (es decir que es un método “lento”). Su simplicidad lo hace ideal para aprender a programar en R y computar numéricamente las raíces de cualquier función continua en un intervalo cerrado.

El método consiste en lo siguiente:

• Debe existir seguridad sobre la continuidad de la función \(f(x)\) en el intervalo \([a,b]\)

• A continuación se verifica que \(f(a)\cdot f(b)<0\).

• Se calcula el punto medio m del intervalo \([a,b]\) y se evalúa \(f(m)\):

  • Si ese valor es igual a cero, ya hemos encontrado la raíz buscada.

  • En caso de que no lo sea, verificamos se redefine el intervalo \([a,b]\) como: \([m,b]\) si \(\text{sign}\{f(a)\}=\text{sign}\{f(m)\}\) ó \([a,m]\) si \(\text{sign}\{f(b)\}=\text{sign}\{f(m)\}\).

• Con este nuevo intervalo se continúa sucesivamente encerrando la solución en un intervalo cada vez más pequeño, hasta alcanzar la precisión deseada.

Programando el método en R:

Vamos a implementar en R el método propuesto a través de una función que se llamará bisección y tendrá 3 argumentos: \(f\), \(a\) y \(b\):

# Halla la raíz de f (una función definida en R) en el intervalo [a,b]
biseccion <- function(f,a,b) {
plot(f,from=a,to=b) # Grafico la función y
abline(h=0,col="blue") # pinto el eje-X de color azul.

# Definimos los valores iniciales de m, n y un 'error'
# que mide la amplitud del intervalo sobre el que trabajamos.
m = (a+b)/2 ; n = 0; error<-abs(a-b)/2

while (error > 1.e-5 & n < 100) { # "Doble-tolerancia" del método
n<-n+1
if (f(m) == 0) break # Si el cero esta en el medio de [a,b] paro.
if (f(m)*f(a) < 0) {b = m} else {a = m} # La regla antes planteada.
m = (a+b)/2 # Redefino el intervalo con los valores actualizados.

# Genero output gráfico y en pantalla del método:
text(m,0,n,cex=0.8,col="red")
error = abs(a-b)/2
cat("X=",m,"\tE=",error,"\n")
                                 } # Fin del while()
}# Fin de la función biseccion().

p <- function(x){return(x-2^{-x}) }

biseccion(p,0,1)

X= 0.75     E= 0.25 
X= 0.625    E= 0.125 
X= 0.6875   E= 0.0625 
X= 0.65625  E= 0.03125 
X= 0.640625     E= 0.015625 
X= 0.6484375    E= 0.0078125 
X= 0.6445312    E= 0.00390625 
X= 0.6425781    E= 0.001953125 
X= 0.6416016    E= 0.0009765625 
X= 0.6411133    E= 0.0004882812 
X= 0.6413574    E= 0.0002441406 
X= 0.6412354    E= 0.0001220703 
X= 0.6411743    E= 6.103516e-05 
X= 0.6412048    E= 3.051758e-05 
X= 0.6411896    E= 1.525879e-05 
X= 0.6411819    E= 7.629395e-06 

Ejercicio para la clase:

1- Encuentra los ceros del polinomio \(q(x)=−x^2+0.5x+4\).

2- Encuentra los valores de \(x\) tal que \(p(x)=q(x)\) y tome \(p(x)=x2^x-1\). (Ayuda: que función se hace cero si \(p(x)=q(x)\)?)

3- Para los siguientes polinomios, se pide factorizarlos y obtener de esta forma las raíces o ceros de manera analítica:

3.1 \(p_1(x)=x^2+x−12\) (sol: -4 y 3).

3.2 \(p_2(x)=x^3−4x^2+x+6\) (sol: -1, 2, y 3).

3.3 \(p_3(x)=x^4−5x^2+4\) (sol: -2, -1, 1 y 2).

4- Programa en R los polinomios del punto 3, grafícalos y utiliza el algorítmo de bisección para computar numéricamente sus raíces. Verifica que los valores numéricos que obtuviste se aproximan a los ceros de la función.

Método de Newton

Definición

El método de Newton (conocido también como el método de Newton-Raphson o el método de Newton-Fourier) es un algoritmo basado en la derivada que permite encontrar aproximaciones de los ceros o raíces de una función real derivable, o inclusive aproximar valores como \(\sqrt{3},\sqrt[4]{3},\sqrt[5]{247},\pi\) encontrando de manera aproximada las raíces de las ecuaciones \(x^2−3=0\), \(x^4−3=0\), \(x^5−247=0\), \(\cos(x)+1=0\) respectivamente.

Formalmente

Sea \(f\) una función derivable definida en un intervalo real \([a,b]\) y sea \(f(r)=0\), es decir, sea \(r\) una raiz real de \(f\). Si \(x_n\) es una aproximación a \(r\), entonces la siguiente aproximación \(x_{n+1}\) esta dada por:

\[x_{n+1}=x_n-\frac{f(x_n)}{f'(x_n)},f'(x_n)\neq 0\]

Descripción del método

Sea \(f:[a,b]→\mathbb{R}\) una función derivable

1. Se elige \(x_0\) en el eje de las \(x\), asumiendo que está cerca de la solución de \(f(x)=0\) (raíz buscada).

2. Calculamos la ecuación punto pendiente de la recta tangente a la función en \((x_0,f(x_0))\), a saber $y−f(x_0)=f′(x_0)(x−x_0) $ .\((1)\)

3. Esta recta debe intersecar al eje de las \(x\), en un punto más cercano a la raíz buscada; en el punto \((x_1,0)\).

4. El punto \((x_1,0)\) satisface la ecuación \((1)\) y sustituyendo, queda: \(0−f(x_0)=f′(x_0)(x−x_0)\). \((2)\)

5. Si \(f′(x_0)\neq0\), entonces despejando \(x_1\) en \((2)\) se obtiene \[x_1=x_0−\frac{f(x_0)}{f′(x_0)}\]

6. Repetimos el procedimiento anterior para \(x_0\), pero ahora comenzando con \(x_1\), en cuyo caso se obtiene \(\displaystyle x_2=x1−\frac{f(x_1)}{f′(x_1)}\). De forma que \(x_2\) está más cerca de la raíz buscada que \(x_1\).

7. Iterando cada vez con el número obtenido, se construye una secuencia: \(x_0,x_1,x_2,…,x_n,…\) de números cada vez más próximos a la raíz, tales que:

\[x_{n+1}=x_n−\frac{f(x_n)}{f′(x_n)}(3)\]

En resumen y de forma menos teórica el procedimiento es el siguiente:

• Observando la gráfica de la función se estima un valor adecado para la primera aproximación \(x_1\).

• Sustituyendo la primera aproximación en (3), se obtiene una segunda aproximación \(x_2\). Luego se calcula \(x_3\) sustituyendo en (3) la segunda aproximación \(x_2\); y así sucesivamente hasta que se llegue a la igualdad \(x_{n+1}=x_n\)

Código en R

El código viene dado por:

expresion <- expression (exp(x)+x-2) # escribimos el polinomio
derivada <- D(expresion, "x") # Derivada del polinomio
 
x <- 0  # Cualquier valor diferente de aprox
aprox <- 3 # valor puntoinicial

while ( x != aprox) {
  
  x <- aprox # Se le asigna el valor aproximado a x.

reemplazoexpresion <- eval(expresion) #Reemplaza el valor de x en "expresión"

reemplazoderiv <- eval(derivada) #Reemplaza el valor de x en "derivada"

#newton

aprox <- x - (reemplazoexpresion/reemplazoderiv) #Ecuación método de Newton

print(x)
}

[1] 3
[1] 2
[1] 1.119203
[1] 0.5821784
[1] 0.4488018
[1] 0.4428652
[1] 0.4428544
[1] 0.4428544

Ejercicio para la clase:

1- Para los siguientes polinomios, se pide factorizarlos y obtener de esta forma las raíces o ceros de manera analítica:

1.1 \(p_1(x)=x^2+x−12\) (sol: -4 y 3).

1.2 \(p_2(x)=x^3−4x^2+x+6\) (sol: -1, 2, y 3).

1.3 \(p_3(x)=x^4−5x^2+4\) (sol: -2, -1, 1 y 2).

2- Programa en R los polinomios del punto 1, grafícalos y utiliza el algorítmo de Newthon Rhapson para computar numéricamente sus raíces. Verifica que los valores numéricos que obtuviste se aproximan a los ceros de la función.