Las Funciones como objetos matemáticos dentro de R

Si bien es cierto que R no fue originalmente ideado como un programa para lidiar con cálculo simbólico, es posible utilizar esta herramienta para evaluar y graficar funciones matemáticas.

Para ello necesitamos primero definir la función en cuestión. Imaginemos que nos interesa estudiar el polinomio \(p(x) = 2x^2 - 0.9x -1\) en el intervalo \([-1,2]\):

p <- function(x){return(2*x^2 -0.9*x -1) } # Declaramos el objeto.

# Podemos simplemente evaluar la función:
p(0)    # Imput un nro -> output un nro
## [1] -1
p(1.5)  # Imput un nro -> output un nro
## [1] 2.15
p(c(-1,0,1.5,2)) # Imput un vector -> output un vector
## [1]  1.90 -1.00  2.15  5.20
# Podemos graficar la función en el intervalo deseado:
plot(p, from=-1, to=2)

# Podemos graficar más de una función en mismo gráfico:
q <- function(x){return(-x^2 + 0.5*x +4) } # Declaramos el objeto.
plot(q, from=-1, to=2, col='red', add=T)

Raíces o Ceros de una función

En varios contextos aplicados necesitaremos encontrar las raíces (ceros) de una función matemática \(f\), es decir el conjunto \(R_f = \{ x\in \mathbb{R} | f(x)=0 \}\). Es bien conocido que no existen métodos generales para resolver este tipo de ecuaciones. Los métodos para aproximar raíces de ecuaciones son, en general, iterativos, es decir consisten en construir una sucesión \(\{x_n\}_{n>0}\) mediante una relación de recurrencia: dado un \(x_0\), entonces \(x_n = \varphi(x_{n-1})\) para \(n=1,2,3,\dots\).

Diremos que el método iterativo es convergente si: \[\lim_{n\rightarrow \infty} x_n = \alpha \text{ y se verifica que } f(\alpha) = 0\]

En este curso veremos y programaremos en R uno de los métodos más sencillos, conocido como Método de Bisecció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:

Un ejemplo gráfico del método (fuente: Wikipedia)

Un ejemplo gráfico del método (fuente: Wikipedia)

Nos interesa el cero de la función continua en el intervalo \([a=a_1, b = b_1]\). Chequeamos que se cumple la condición: \(sign \{f(a_1\} \neq sign\{f(b_1)\}\), luego tomamos cualquier punto intermedio, por ejemplo \(m = a_2\). Como \(sign \{f(a)\} = sign\{f(m)\}\), entonces re-definimos el intervalo de búsqueda como \([m=a_2, b]\). Volvemos a elegir algún punto intermedio ahora del intervalo \([a_2, b]\), por ejemplo \(m=b_2\). Como \(sign \{f(m)\} = sign\{f(b)\}\), redefinimos el intervalo de búsqueda a \([a_2, m=b_2]\). Volvemos a escoger un valor intermedio dentro de este último intervalo, por ejemplo \(m=a_3\) y continuamos este proceso hasta aproximarnos cada vez más a la raíz/cero de la función.

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


biseccion(p,-1,0)

## X= -0.75     E= 0.25 
## X= -0.625    E= 0.125 
## X= -0.5625   E= 0.0625 
## X= -0.53125  E= 0.03125 
## X= -0.515625     E= 0.015625 
## X= -0.5234375    E= 0.0078125 
## X= -0.5195312    E= 0.00390625 
## X= -0.5175781    E= 0.001953125 
## X= -0.5166016    E= 0.0009765625 
## X= -0.5170898    E= 0.0004882812 
## X= -0.5168457    E= 0.0002441406 
## X= -0.5169678    E= 0.0001220703 
## X= -0.5170288    E= 6.103516e-05 
## X= -0.5170593    E= 3.051758e-05 
## X= -0.5170441    E= 1.525879e-05 
## X= -0.5170364    E= 7.629395e-06

Método pre-programados en R:

Existen una batería de métodos numéricos desarrollados en matemática que nos permiten aproximar los ceros de una función y que ya se encuentran pre-programados en R. Uno de ellos es la función uniroot() que tiene 2 argumentos: una función matemática \(f\) declarada previamente en R y un intervalo \([a,b]\) en donde buscar los ceros de \(f\) bajo la condición de que \(sign \{ f(a)\} \neq sign \{ f(b) \}\)

p <- function(x){return(2*x^2 -0.9*x -1) } # Declaramos el objeto.

# Para determinar los ceros de la función hay que establecer un intervalo de búsqueda:
uniroot(p, interval=c(-1,0))
## $root
## [1] -0.5170353
## 
## $f.root
## [1] -1.716935e-05
## 
## $iter
## [1] 6
## 
## $init.it
## [1] NA
## 
## $estim.prec
## [1] 6.103516e-05

Podemos observar que el cero de \(p(x)\) en el intervalo \([-1,0]\) esta muy cerca de \(x=0.5\). Ejecuta la siguiente línea de código uniroot(p, interval=c(-1,2)), porque da errores?

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)\). (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_2(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.