Solución de Raíces Final

Descargar como pdf o txt
Descargar como pdf o txt
Está en la página 1de 29

Tabla de contenido

SOLUCIÓN DE RAÍCES ..................................................................................................................2


1. Métodos iterativos para resolver x = g(x) .....................................................................................2
1.1 Puntos fijos...............................................................................................................................3
1.1.1. Interpretación gráfica de la iteración de punto fijo...........................................................6
2. Los métodos de localización de raíces......................................................................................... 10
2.1. El método de bisección de Bolzano ........................................................................................ 11
3. Métodos de la regula falsi o método de la posición falsa ............................................................ 16
3.1. Covergencia del método de la regula falsi ............................................................................. 17
4. Los métodos de Newton Raphson y de la secante ....................................................................... 20
4.1. Métodos tangenciales para el cálculo de raíces; Método de Newton...................................... 20
4.1.1. Errores por la división entre cero ....................................................................................... 22
4.2. El método de la secante ........................................................................................................ 25
SOLUCIÓN DE RAÍCES

1. Métodos iterativos para resolver x = g(x)


Una técnica fundamental de cálculo científico es la iteración. Como su nombre lo dice, se trata
de repetir un proceso hasta que se obtiene un resultado. Se usan métodos iterativos para hallar
raíces de ecuaciones, soluciones de los sistemas lineales y no lineales y soluciones de
ecuaciones diferenciales.
Se necesita una regla, fórmula o función g(x), con la que calcularemos los sucesivos términos,
junto con un valor de partida p0. Lo que se produce es una sucesión {𝑝𝑘 } obtenida mediante el
proceso iterativo 𝑝𝑘+1 = 𝑔(𝑝𝑘 ). La sucesión se ajusta al siguiente patrón
𝑝0 (valor de partida)
𝑝1 = 𝑔(𝑝0 )
𝑝2 = 𝑔(𝑝1 )
.
.
.
(1.1)
𝑝𝑘 = 𝑔(𝑝𝑘−1 )
𝑝𝑘+1 = 𝑔(𝑝𝑘 )
.
.
.
¿Qué nos dice una sucesión interminable de número como ésta? Si los números tienden a un
límite, entonces se puede afirmar que se tiene algo entre manos. Pero ¿qué ocurre si los números
divergen o son periódicos? El ejemplo siguiente muestra una situación como ésta.
Ejemplo teórico 1.1. El proceso iterativo 𝑝0 =1 y 𝑝𝑘+1 = 1.001 para k = 0, 1, … produce una
sucesión divergente; sus primeros y centésimo términos son:
𝑝1 = 1.001 𝑝0 = (1.001)(1.000000) = 1.001000
𝑝2 = 1.001 𝑝1 = (1.001)(1.001000) = 1.002001
𝑝2 = 1.001 𝑝3 = (1.001)(1.002001) = 1.003003
. .
. .
. .
𝑝100 = 1.001 𝑝99 = (1.001)(1.104012) = 1.105116
Este proceso puede continuarse indefinidamente y es fácil probar que se verifica que
lim 𝑝𝑛 = +∞. Se tiene que la sucesión {𝑝𝑘 } es la solución numérica de la ecuación diferencial
𝑛→∞
𝑑𝑦 𝑑𝑦
ordinaria (EDO) y’ = 0.001 y ( 𝑑𝑡 = 0.001𝑦 → ∫ = ∫ 0.001𝑑𝑡 → 𝑙𝑛𝑦 = 0.001 𝑡 + 𝐶) . Se
𝑦

sabe que la solución de esta EDO es 𝑦(𝑥 ) = 𝑒 0.001 𝑥 . De hecho, si comparáramos el término
centésimo de la sucesión con y (100), se observa que 𝑃100 = 1.105116 ≈ 1.105171 = 𝑒 0.1 =
𝑦(100).
En esta sección estudiaremos qué tipos de funciones g(x) producen sucesiones {𝑝𝑘 } que
convergen.
1.1 Puntos fijos
Definición 1.1. (Punto fijo). Un punto fijo de una función g(x) es un número real P tal que
P=g(P).
Geométricamente hablando, los puntos fijos de una función g(x) son los puntos de intersección
de la curva y = g(x) con la recta y = x.
Definición 1.2. (Iteración de punto fijo). La iteración 𝑝𝑛+1 = 𝑔(𝑝𝑛 ) para n = 0,1, … se llama
iteración de punto fijo.
Teorema 1.1. Supongamos que g es una función continua y que {𝑝𝑛 }∞
𝑛=0 es una sucesión

generada por iteración de punto fijo. Si lim 𝑝𝑛 = 𝑃, entonces P es un punto fijo de g(x).
𝑛→∞

Demostración. Si lim 𝑝𝑛 = 𝑃, entonces lim 𝑝𝑛+1 = 𝑃. Usando esto, la continuidad de g y la


𝑛→∞ 𝑛→∞

relación 𝑝𝑛+1 = 𝑔(𝑝𝑛 ), se deduce que

𝑔(𝑃) = 𝑔 ( lim 𝑝𝑛 ) = lim 𝑔(𝑝𝑛 ) = lim 𝑝𝑛+1 = 𝑃 (1.2)


𝑛→∞ 𝑛→∞ 𝑛→∞

Por tanto, 𝑃 es un punto fijo de g(x).


Ejemplo teórico 1.2. Considerando la iteración convergente
𝑝0 = 0.5 y 𝑝𝑘+1 = 𝑒 −𝑝𝑘 para k = 0,1,…
Se puede calcular los diez primeros términos como:
𝑝1 = 𝑒 −0.500000 = 0.606531
𝑝2 = 𝑒 −0.606531 = 0.545239
𝑝3 = 𝑒 −0.545239 = 0.579703
. .
. .
. .
𝑝9 = 𝑒 −0.566409 = 0.567560
𝑝10 = 𝑒 −0.567560 = 0.566907
La sucesión converge y, calculando un poco más, se tiene:
lim 𝑝𝑛 = 0.567143
𝑛→∞

De esa manera, lo que se ha obtenido es una aproximación al punto fijo de la función 𝑔(𝑥) =
𝑒 −𝑥 .
Los dos siguientes teoremas establecen las condiciones para la existencia de un punto fijo y
para la convergencia del proceso de iteración de punto fijo.
Teorema 1.2. Supongamos que 𝑔 ∈ 𝐶 [𝑎, 𝑏]
(1.3) Si la imagen de la aplicación 𝑦 = 𝑔(𝑥) verifica que 𝑦 ∈ [𝑎, 𝑏] para cada punto 𝑥 ∈ [𝑎, 𝑏],
entonces g tiene un punto fijo en [𝑎, 𝑏].
(1.4) Suponiendo, además, que 𝑔′(𝑥) está definida en (a, b) y que |𝑔′(𝑥)| < 1 para todo 𝑥 ∈
(𝑎, 𝑏), entonces g tiene un único punto fijo P en [𝑎, 𝑏].
Demostración de ecuación 1.3. Si 𝑔(𝑎) = 𝑎 o 𝑔(𝑏) = 𝑏, entonces la conclusión es cierta.
Suponiendo, entonces, que los valores 𝑔(𝑎) y 𝑔(𝑏) verifican que 𝑔(𝑎) ∈ (𝑎, 𝑏] y que 𝑔(𝑏) ∈
[𝑎, 𝑏). La función 𝑓 (𝑥 ) ≡ 𝑥 − 𝑔(𝑥) tiene la siguiente propiedad:
𝑓 (𝑎 ) = 𝑎 − 𝑔 (𝑎 ) < 0 y 𝑓(𝑏) = 𝑏 − 𝑔(𝑏) > 0
Teorema a (Teorema del valor intermedio o de Bolzano). Supongamos que tenemos 𝑓 ∈
[𝑎, 𝑏] y que L es cualquier número entre 𝑓(𝑎) y 𝑓(𝑏). Entonces existe un número c en (𝑎, 𝑏)
tal que 𝑓 (𝑐 ) = 𝐿.
Aplicando el Teorema a, a f(x) y la constante L=0, se deduce que existe un número P en (𝑎, 𝑏)
tal que f(P)=0. Por tanto, P=g(P) y P es el punto fijo de g(x) que deberíamos encontrar.

Teorema b (Teorema del valor medio o de Lagrange). Supongamos que 𝑓 ∈ [𝑎, 𝑏] y que
𝑓′(𝑥) existe para todo 𝑥 ∈ (𝑎, 𝑏). Entonces existe un número c en (𝑎, 𝑏) tal que
𝑓 (𝑏) − 𝑓(𝑎)
𝑓 ′ (𝑐 ) =
𝑏−𝑎
Geométricamente este teorema nos dice que hay al menos un número 𝑐 ∈ (𝑎, 𝑏) al que la
pendiente de la recta tangente a la curva 𝑦 = 𝑓(𝑥) en el punto (𝑐, 𝑓(𝑐)) es igual a la pendiente
de la recta secante que pasa por los puntos (𝑎, 𝑓(𝑎)) y (𝑏, 𝑓(𝑏)).
Demostración de (4). Ahora se debe probar que la solución es única. Suponiendo, por el
contrario, que existen dos puntos fijos distintos P1 y P2. Aplicando el Teorema b, el teorema
del valor medio se deduce que existe un número 𝑑 ∈ (𝑎, 𝑏) tal que
𝑔(𝑃2 ) − 𝑔(𝑃1 )
𝑔 ′ (𝑑 ) = (1.5)
𝑃2 − 𝑃1
A continuación, se usa que 𝑔(𝑃1 ) = 𝑃1 y 𝑔(𝑃2 ) = 𝑃2 para simplificar el miembro del lado
derecho de la relación
𝑃2 − 𝑃1
𝑔 ′ (𝑑 ) = =1
𝑃2 − 𝑃1
Lo que contradice la hipótesis hecha en (1.4) de que |𝑔′ (𝑥 )| < 1 en (a, b). Así que no es posible
que existan dos puntos fijos y, por tanto, bajo la condición dada en (1.4), 𝑔(𝑥 ) tiene un punto
fijo P en [a, b].
Ejemplo teórico 1.3. Se aplicará el Teorema 1.2 para probar rigurosamente que 𝑔(𝑥 ) =
cos(𝑥 ) tiene un único punto fijo en [0,1].
Claramente 𝑔 ∈ 𝐶[0,1]. Además, 𝑔(𝑥 ) = cos(𝑥 ) es una función decreciente en [0,1], por lo
que su imagen o rango 𝑔([0,1]) es [cos(1), 1)] ⊆ [0,1]. Esto prueba que se cumple la condición
(1.3) del Teorema 1.2, así que g tiene un punto fijo en [0,1]. Finalmente, si 𝑥 ∈ (0,1), entonces
|𝑔′ (𝑥 )| = |−𝑠𝑒𝑛(𝑥)| = 𝑠𝑒𝑛(𝑥) ≤ 𝑠𝑒𝑛(1) ≈ 0.8415 < 1. Por consiguiente, la condición (1.4)
del Teorema 1.2 se cumple y el punto fijo de g en [0,1] es único.
Se establecerá el teorema que se suele usar para determinar si el proceso de iteración de un
punto fijo en (1.1) produce una sucesión convergente o divergente.

Teorema 1.3. (Teorema del punto fijo). Supongamos que (i) g, g’ ∈ 𝐶 [𝑎, 𝑏], (ii) K es una
constante positiva, (iii) 𝑝0 ∈ (𝑎, 𝑏) y (iv) 𝑔(𝑥) ∈ 𝐶 [𝑎, 𝑏] para todo 𝑥 ∈ [𝑎, 𝑏]. Entonces hay un
punto fijo P de g en [𝑎, 𝑏].
(1.6) Si |𝑔′(𝑥)| ≤ 𝐾 < 1 para todo 𝑥 ∈ [𝑎, 𝑏], entonces P es el único punto fijo de g en [𝑎, 𝑏],
y la iteración 𝑝𝑛 = 𝑔(𝑝𝑛−1 ) converge a dicho punto fijo P. En este caso, se dice que P es un
punto fijo atractivo.
(1.7) Si |𝑔′(𝑃)| > 1 y 𝑝0 ≠ 𝑃 entonces la iteración 𝑝𝑛 =g (𝑝𝑛−1 ) no converge a P. En este caso,
se dice que P es un punto fijo repulsivo y la iteración presenta divergencia local.
Observación 1. Como 𝑔′ es continua en un intervalo al que pertenece P, se pueden usar hipótesis
más simples en (1.6): que |𝑔′(𝑃)| < 1 y que 𝑝0 está suficientemente cerca de P.

Demostración. El apartado (1.3) de Teorema 1.2 y las hipótesis (i) y (iv) del Teorema 1.3,
garantizan que g tiene un punto fijo P en [𝑎, 𝑏]. Para demostrar (1.6) se comienza notando que,
usando (1.3) y (1.4), las hipótesis (i) - (iv) implican, por el Teorema 1.2, que el punto fijo de g
en [𝑎, 𝑏] es único. De (iii) y (iv) se deduce, por inducción, que los puntos el punto {𝑝𝑛 }∞
𝑛=0

están todos [𝑎, 𝑏]. Seguidamente, empezando con 𝑝0 aplicando el Teorema b, el teorema del
valor medio, para deducir que existe un valor 𝑐0 ∈ (𝑎, 𝑏) tal que
|𝑃 − 𝑝1 | = |𝑔(𝑃) − 𝑔(𝑝0 )| = |𝑔′ (𝑐0 )(𝑃 − 𝑝0 )|
(1.8)
= |𝑔 ′ (𝑐0 )||(𝑃 − 𝑝0 )| ≤ 𝐾 |(𝑃 − 𝑝0 )| < |(𝑃 − 𝑝0 )|
Por tanto, 𝑝1 está más cerca de P que 𝑝0 (véase la Figura 1.1).

Figura 1.1. La relación entre 𝑃, 𝑝0 , 𝑝1 , |(𝑃 − 𝑝0 )| 𝑦 |(𝑃 − 𝑝1 )|.


Razonando de modo análogo se tiene, en general,
(1.9)
|𝑃 − 𝑝𝑛 | = |𝑔(𝑃) − 𝑔(𝑝𝑛−1 )| = |𝑔′ (𝑐𝑛−1 )(𝑃 − 𝑝𝑛−1 )|
= |𝑔′ (𝑐𝑛−1 )||(𝑃 − 𝑝𝑛−1 )| ≤ 𝐾|(𝑃 − 𝑝𝑛−1 )| < |(𝑃 − 𝑝𝑛−1 )|
Para completar la demostración de (1.6), se prueba que
(1.10)
lim |𝑃 − 𝑝𝑛 |
𝑛→∞
En primer lugar se establece por inducción la desigualdad
|𝑃 − 𝑝𝑛 | ≤ 𝐾 𝑛 |(𝑃 − 𝑝0 )| (1.11)

El caso n = 1 está ya comprobado en la relación (1.8). Usando la hipótesis de inducción


|𝑃 − 𝑝𝑛−1 | ≤ 𝐾 𝑛−1 |(𝑃 − 𝑝0 )| y la relación (1.9), se obtiene
|𝑃 − 𝑝𝑛 | ≤ 𝐾 |(𝑃 − 𝑝𝑛−1 )| ≤ 𝐾𝐾 𝑛−1 |(𝑃 − 𝑝0 )| = 𝐾 𝑛 |𝑃 − 𝑝0 |
En consecuencia, la desigualdad (1.11) se verifica, por inducción, para todo n. Puesto que
0 < K < 1, el factor 𝐾 𝑛 converge a cero cuando n tiende a infinito, por tanto
0 ≤ lim |𝑃 − 𝑝𝑛 | = lim 𝐾 𝑛 |𝑃 − 𝑝0 | = 0 (1.12)
𝑛→∞ 𝑛→∞

de aquí se deduce que lim |𝑃 − 𝑝𝑛 | = 0, con lo cual lim 𝑝𝑛 = 𝑃, lo que concluye la


𝑛→∞ 𝑛→∞

demostración del apartado (1.6) del Teorema 1.3.


Corolario 1.1. Supongamos que g verifica las hipótesis dadas en el apartado (1.6) del Teorema
1.3. Entonces las siguientes desigualdades proporcionan cotas del error que se comete cuando
se usa 𝑝𝑛 como aproximación a P:
|𝑃 − 𝑝𝑛 | ≤ 𝐾 𝑛 |𝑃 − 𝑝0 | 𝑝𝑎𝑟𝑎 𝑡𝑜𝑑𝑜 𝑛 ≥ 1, (1.13)
y
𝐾 𝑛 |𝑝1 − 𝑝0 |
|𝑃 − 𝑝𝑛 | ≤ 𝑝𝑎𝑟𝑎 𝑡𝑜𝑑𝑜 𝑛 ≥ 1, (1.14)
1−𝐾
1.1.1. Interpretación gráfica de la iteración de punto fijo
Puesto que se busca un punto fijo P de la función g(x), es necesario que la curva y = g(x) y la
recta y = x se corten en el punto (P, P). Los dos tipos simples de iteración convergente,
monótona y oscilante se muestran en la Figura 1.2 a) y 1.2 b), respectivamente.

Figura 1.2. a) Convergencia monótona cuando 0 < 𝑔 ′ (𝑃) < 1, b) Convergencia oscilante si −1 < 𝑔 ′ (𝑃) < 0.

También se tiene la divergencia monótona si 1 < 𝑔′ (𝑃) y divergencia 𝑔′ (𝑃) < −1, las cuales
no se muestran pero se puede checar en internet o en algún libro de métodos numéricos.
Ejemplo teórico 1.4. Considerando la iteración 𝑝𝑛+1 = 𝑔(𝑝𝑛 ) cuando la función viene dada
𝑥2
por 𝑔(𝑥) = 1 + 𝑥 − . Los puntos fijos pueden determinarse resolviendo la ecuación
4

𝑥 = 𝑔(𝑥 ), cuyas soluciones (los puntos fijos de g) son 𝑥 = −2 y 𝑥 = 2. La derivada de la


𝑥
función es 𝑔′(𝑥 ) = 1 − 2 y los dos casos que tenemos que considerar son:
Caso (i): 𝑃 = −2.0 Caso (ii): 𝑃 = 2.0
Inicio 𝑝0 = −2.05 Inicio 𝑝0 = 1.6
𝑝1 = −2.100625 𝑝1 = 1.96
𝑝2 = −2.20378135 𝑝2 = 1.996
𝑝3 = −2.41794441 𝑝3 = 1.99999996
⋮ ⋮
lim 𝑝𝑛 = −∞ lim 𝑝𝑛 = 2
𝑛→∞ 𝑛→∞
3 1
Como |𝑔′(𝑥)| > 2 en [−3, −1], por el Teorema 1.3, Como |𝑔′(𝑥)| < 2 en [1,3], por el Teorema 1.3, la

la sucesión no converge a 𝑃 = −2. sucesión converge a 𝑃 = 2.

El Teorema 1.3 no dice qué ocurrirá si |𝑔′(𝑃)| = 1. El siguiente ejemplo se ha construido


especialmente para {𝑝𝑛 } converja cuando 𝑝0 > 𝑃 y diverja cuando 𝑝0 < 𝑃.

Ejemplo teórico 1.5. Considerando la iteración 𝑝𝑛+1 = 𝑔(𝑝𝑛 ) cuando la función viene dada
1
por 𝑔(𝑥) = 2(𝑥 − 1)2 para 𝑥 ≥ 1. Esta función solo tiene un punto fijo 𝑃 = 2 y su derivada es
1
𝑔′(𝑥) = 1 de manera que 𝑔′(2) = 1 y no se puede aplicar el Teorema 1.3. Considerando
(𝑥−1)2

los dos casos según que el punto de partida esté a la derecha o a la izquierda de 𝑃 = 2.
Caso (i): Inicio 𝑝0 = 1.5 Caso (ii): Inicio 𝑝0 = 2.5
𝑝1 = 1.41421356 𝑝1 = 2.44948974
𝑝2 = 1.28718851 𝑝2 = 2.40789513
𝑝3 = 1.07179943 𝑝3 = 2.37309514
𝑝4 = 0.53590832 𝑝4 = 2.34358284
1
𝑝5 = 2(−0.46409168) 2 ⋮
lim 𝑝𝑛 = 2
𝑛→∞

Como 𝑝4 está fuera del dominio 𝑔(𝑥), el término 𝑝5 Esta sucesión converge muy lentamente a 𝑃 = 2; de
no puede calcularse. hecho, 𝑃1000 = 2.00398714.
En el ejemplo teórico 1.5, caso (ii), la sucesión converge muy lentamente; después de 1000
iteraciones los tres siguientes términos son
𝑃1000 = 2.00398714, 𝑃1001 = 2.00398317, 𝑃1002 = 2.00397921
Esto no debería preocuparnos; después de todo, ¡podríamos calcular unos pocos miles de
términos más para encontrar una aproximación mejor! Pero ¿cuál es el criterio para detener una
iteración? Si usamos la diferencia entre dos términos consecutivos, se obtiene:
|𝑃1001 − 𝑃1002 | = |2.00398317 − 2.00397921| = 0.00000396
Sin embargo, se sabe que el error absoluto de la aproximación 𝑃1000 𝑒𝑠
|𝑃 − 𝑃1000 | = |2.00000000 − 2.00398714| = 0.00398714
Que resulta ser unas mil veces mayor que |𝑃1001 − 𝑃1002 |. Esto prueba que la cercanía entre
dos términos consecutivos no garantiza que hayamos obtenido una precisión alta; sin embargo,
éste (o bien el error relativo) es usualmente el único criterio del que se dispone y, a menudo, es
el que se usa para determinar cuándo debe detenerse un proceso iterativo.
Los coeficientes se pueden calcular con el siguiente programa en Matlab.
% Programa (Iteración de punto fijo). Aproximación a una solución de la
% ecuación x = g(x) mediante la iteración pn + a = g(pn) realizada a partir de una
% aproximación inicial p0
% LIBRO: MÉTODOS NUMÉRICOS CON MATLAB, MATHEWS Y FINK
function [k,p,err,P] = fixpt(g,p0,tol,max1)
% Valores que se van a asignar a [k,p,err,P]= fixpt(@g,1.5,0.001,10)
% Datos
% g es la función de iteración
% p0 es el punto de partida
% tol es la tolerancia
% max1 es el número máximo de iteraciones0
% Resultados
% k es el número de iteraciones realizadas
% p es la aproximación al punto fijo
% err es la diferencia entre dos términos consecutivos
% P es la sucesión {pn} completa
P(1) = p0
for k = 2:max1
P(k) = feval(g,P(k-1));
err = abs(P(k)-P(k-1));
relerr = err/(abs(P(k)) + eps)
p = P(k)
if (err<tol) | (relerr<tol),
break
end
P = P';
end

La solución para un problema es


% solución fixpt
clc
clear all
g = @ (x) 2*(x-1)^(1/2);
[k,p,err,P]= fixpt(g, 2.5, 0.001,10);

La solución para el problema antes mostrado es:


relerr =
0.006692659609693
p=
2.245162029973827
2. Los métodos de localización de raíces
Se considerará un tema de interés familiar, suponiendo que ahorramos dinero haciendo
depósitos mensuales de P euros a una tasa de interés anual I que se compone cada mes; entonces
la cantidad total A de euros acumulada después de N depósitos es
𝐼 𝐼 2 𝐼 𝑁−1
𝐴 = 𝑃 + 𝑃 (1 + 12) + 𝑃 (1 + 12) + ⋯ + 𝑃 (1 + 12) (2.1)

El primer término del miembro derecho de la expresión (2.1) es el último depósito; el segundo
es el penúltimo depósito cuya contribución, habiendo acumulado un periodo de interés, es
𝐼
𝑃 (1 + 12) euros; el antepenúltimo depósito ha acumulado dos periodos de interés y su

𝐼 2
contribución al total (1 + 12) euros, y así sucesivamente. Finalmente, el primer depósito, que

𝐼 𝑁−1
ha acumulado interés durante 𝑁 − 1 periodos, contribuye con (1 + 12) euros al total.

Recordemos la fórmula de la suma de los N términos de una progresión geométrica:


1− 𝑟 𝑁
1 + 𝑟 + 𝑟 2 + 𝑟 3 + ⋯ + 𝑟 𝑁−1 = (2.2)
1− 𝑟

Escribiendo (2.1) como


𝐼 𝐼 2 𝐼 𝑁−1
𝐴 = 𝑃 (1 + (1 + 12) + (1 + 12) + ⋯ + (1 + 12) ) (2.3)
𝐼
Y sustituyendo 𝑟 = 1 + 12 en (2.2) se obtiene

𝐼 𝑁
1 − (1 + 12)
𝐴=𝑃
𝐼
1 − (1 + 12)

Que se puede simplificar para obtener la fórmula del capital acumulado cuando el interés se
compone mensualmente:
𝑃 𝐼 𝑁
𝐴= ((1 + ) − 1)
𝐼/12 12
Se utilizará esta fórmula en el siguiente ejemplo, donde se tiene que usar varias veces para
hallar la respuesta.
Ejemplo teórico 2.1. Ahorramos 250 euros al mes durante 20 años y deseamos que el capital
acumulado al final de estos 20 años sea 250 000 euros. ¿A qué tasa de interés I se tiene que
invertir el dinero para obtener ese resultado? Si se fija N = 240, entonces A es función
únicamente de I; o sea, A = A(I). Si se comienza con dos tentativas, I0 = 0.12 e I1 = 0.13 y
realizaremos una sucesión de cálculos para ir estrechando el margen de la respuesta.
Empezando con I0 = 0.12 se obtiene
250 0.12 240
𝐴(0.12) = ((1 + ) − 1) = 247 314
0.12/12 12
Puesto que este valor es menor que el deseado, probamos con 𝐼 = 0.13
250 0.13 240
𝐴(0.13) = ((1 + ) − 1) = 282 311
0.13/12 12
Éste es un poco alto, así que probamos con la media de ambos I2 = 0.125;
250 0.125 240
𝐴(0.125) = ((1 + ) − 1) = 264 623
0.125/12 12
Éste es, de nuevo, un poco alto, por lo que se deduce que el valor deseado debe estar en el
intervalo [0.12, 0.125] y se prueba, ahora, con su punto medio I3 = 0.1225:
250 0.1225 240
𝐴(0.1225) = ((1 + ) − 1) = 255 803
0.1225/12 12
Este valor es mayor que lo que deseamos, así que ahora estrechamos el intervalo a
[0.12, 0.1225] y usamos su punto medio I4 = 0.12125 para el último cálculo
250 0.12125 240
𝐴(0.12125) = ((1 + ) − 1) = 251 518
0.12125/12 12
Se puede seguir realizando iteraciones hasta obtener tantas cifras significativas de la solución
como deseemos. El propósito de este ejemplo era encontrar el valor de I para el que la función
A(I) es igual a un cierto valor especificado L, o sea, hallar una solución de la ecuación
A(I) = L. La práctica habitual es pasar L al primer miembro y resolver la ecuación A(I) – L =
0.

Definición 2.1 (Raíz de una ecuación, cero de una función). Suponiendo que f(x) es una
función continua. Cualquier número real r tal que f(r) = 0 se llama raíz de la ecuación
f(x) = 0; también se dice que r es un cero de la función f(x).
Por ejemplo, la ecuación 2𝑥 2 + 5𝑥 − 3 = 0 tiene dos raíces reales 𝑟1 = 0.5 y 𝑟2 = −3.0,
mientras que la función correspondiente 𝑓 (𝑥 ) = 2𝑥 2 + 5𝑥 − 3 = (2𝑥 − 1)(𝑥 + 3) tiene dos
ceros reales 𝑟1 = 0.5 y 𝑟2 = −3.0.
2.1. El método de bisección de Bolzano
En esta sección se hallarán ceros de funciones continuas, se tiene que empezar con un intervalo
de partida [𝑎, 𝑏] en el que 𝑓(𝑎) y 𝑓(𝑏) tengan distinto signo. Entonces, por el Teorema a de la
sección anterior, el teorema del valor intermedio, la gráfica 𝑦 = 𝑓(𝑥) cruzará el eje OX en un
cero x = r que está en dicho intervalo (véase Figura 2.1).
Figura 2.1. Proceso de decisión en el método de bisección. a) si f(a) y f(c) tienen signos opuestos, entonces recorta
por la derecha, b) si f(a) y f(c) tienen signos opuestos, entonces recorta por la izquierda.
El método de bisección consiste en ir acercando sistemáticamente los extremos del intervalo
hasta que se obtenga un intervalo de anchura suficientemente pequeña en el que se localiza un
cero. El proceso de decisión para subdividir el intervalo consiste en tomar el punto medio del
𝑎+𝑏
intervalo 𝑐 = y luego analizar las tres posibilidades que pueden darse:
2

Si f(a) y f(c) tienen signos opuestos, entonces hay un cero en [a,c] (2.4)
Si f(c) y f(b) tienen signos opuestos, entonces hay un cero en [c,b] (2.5)
Si f(c)=0, entonces hay un cero (2.6)
Los casos (2.4) y (2.5) se pueden observar en la Figura 2.1, el caso (2.6) no se presenta en esta
Figura, pero si se tiene el valor de c. Para continuar el proceso, renombramos el nuevo intervalo
más pequeño también como [a, b] y se repite el proceso hasta que el intervalo sea tan pequeño
como se desee. Puesto que el proceso de bisección genera una sucesión de intervalos encajados,
con sus correspondientes puntos medios, se usará la siguiente notación para tener registro de
los detalles del proceso:
𝑎0 +𝑏0
[𝑎0 , 𝑏0 ] es el intervalo de partida y 𝑐0 = es un punto medio.
2

[𝑎1 , 𝑏1 ] es el segundo intervalo en el que se localiza un cero y 𝑐1 es su punto medio; el


intervalo [𝑎1 , 𝑏1 ] es la mitad de ancho que [𝑎0 , 𝑏0 ].
(2.7)
Después de llegar al intervalo [𝑎𝑛 , 𝑏𝑛 ], en el que también se localiza un cero y cuyo
punto medio es 𝑐𝑛 , se construye el intervalo [𝑎𝑛+1 , 𝑏𝑛+1 ], en el que también se sigue
localizando un cero, y que mide la mitad que [𝑎𝑛 , 𝑏𝑛 ].
Se deja como ejercicio el demostrar que la sucesión de los extremos izquierdos de estos
intervalos es creciente y que la sucesión de los extremos derechos es decreciente, es decir,
𝑎0 ≤ 𝑎1 ≤ ⋯ ≤ 𝑎𝑛 ≤ ⋯ ≤ 𝑏𝑛 ≤ ⋯ ≤ 𝑏1 ≤ 𝑏0 (2.8)
𝑎𝑛 +𝑏𝑛
Donde 𝑐𝑛 = y si 𝑓(𝑎𝑛+1 )𝑓 (𝑏𝑛+1 ) < 0, entonces
2

[𝑎𝑛+1 , 𝑏𝑛+1 ] = [𝑎𝑛 , 𝑐𝑛 ] o bien [𝑎𝑛+1 , 𝑏𝑛+1 ] = [𝑐𝑛 , 𝑎𝑛 ] para cada n. (2.9)
Teorema 2.1 (Convergencia del método de bisección). Suponiendo que 𝑓 ∈ 𝐶[𝑎, 𝑏] y que
𝑓 (𝑎) y 𝑓(𝑏) tienen signos distintos. Sea {𝑐 }∞
𝑛=0 la sucesión de puntos medios de los intervalos

generados por el método de bisección dado en (2.8) y (2.9). Entonces existe un número 𝑟 ∈
[𝑎, 𝑏] tal que 𝑓 (𝑟) = 0 y, además,
𝑏−𝑎
|𝑟 − 𝑐𝑛 | ≤ (2.10)
2𝑛+1

En particular, la sucesión {𝑐 }∞
𝑛=0 converge al cero x = r; esto es,

lim 𝑐𝑛 = 𝑟 (2.11)
𝑛→∞

Ejemplo teórico 2.2. La función ℎ(𝑥 ) = 𝑥 𝑠𝑒𝑛(𝑥) aparece en el estudio de vibraciones


forzadas no amortiguadas. Hay que hallar el valor de 𝑥 que está dentro del intervalo [0, 2] y en
el que la función vale ℎ(𝑥 ) = 1 (el ángulo 𝑥 en la función 𝑠𝑒𝑛(𝑥) se mide en radianes). Se usa
el método de bisección para hallar un cero de la función 𝑓 (𝑥 ) = 𝑥 𝑠𝑒𝑛(𝑥 ) − 1. Empezando con
𝑎0 = 0 y 𝑏0 = 2, calculamos
𝑓 (0) = −1.000000 𝑦 𝑓(2) = 0.818595
de manera que hay una raíz de f(x) = 0 en el intervalo [0,2]. En el punto medio 𝑐0 = 1 se tiene
que f (1) = - 0.158529, luego la función cambia de signo en el intervalo [c0, b0] = [1, 2].
Para seguir, se recorta el intervalo por la izquierda y se pone 𝑎1 = 𝑐0 𝑦 𝑏1 = 𝑏0 . El nuevo punto
medio es 𝑐1 = 1.5 y se tiene 𝑓(𝑐1 ) = 0.496242. Puesto que 𝑓(1) = −0.158529 y 𝑓(1.5) =
0.496242, la raíz está en [𝑎1 , 𝑐1 ] = [1.0, 1.5] y la siguiente decisión es recortar por la derecha
y poner 𝑎2 = 𝑎1 𝑦 𝑏2 = 𝑐1 . De esta forma se obtiene una sucesión {𝑐𝑘 } que converge a
𝑟 ≈ 1.114157141. Se puede observar los cálculos de los ocho primeros pasos en la Tabla 2.1.
Tabla 2.1.
k Extremo izquierdo, 𝑎𝑘 Punto medio, 𝑐𝑘 Extremo derecho, 𝑏𝑘 Valor de la función, 𝑓(𝑐𝑘 )
0 0 1 2 -0.158529
1 1.0 1.5 2.0 0.496242
2 1.00 1.25 1.5 0.186231
3 1.000 1.125 1.125 0.015051
4 1.0000 1.0625 1.1250 -0.071827
5 1.06250 1.09375 1.12500 -0.028362
6 1.093750 1.109375 1.125000 -0.006643
7 1.1093750 1.1171857 1.1250000 0.004208
8 1.10937500 1.11328125 1.11718750 -0.001216
⋮ ⋮ ⋮ ⋮ ⋮
Una de las virtudes del método de bisección es que la fórmula proporciona una estimación
predeterminada de la precisión calculada. En el ejemplo teórico 2.2 la anchura del intervalo
era 𝑏0 − 𝑎0 = 2. Suponiendo que se hubiera continuado con la Tabla 2.1 hasta la iteración
2−0
trigésimo-primera; entonces, por (2.10), la cota del error sería |𝐸31 | ≤ ≈ 4.656613𝑥10−10 .
232

Por tanto, 𝑐31 sería una aproximación a r con nueve cifras decimales de precisión. El número
N de bisecciones sucesivas que garantizaría que el punto medio 𝑐𝑁 es una aproximación a un
cero con un error menor que el valor prefijado 𝛿 es
ln(𝑏−𝑎)−ln(𝛿)
𝑁 = 𝑒𝑛𝑡 ( ) (2.12)
ln(2)

(donde ent (x) denota la parte entera de un número x).


Los coeficientes se
% Programa (Método de bisección). Aproximación a una raíz de la
pueden calcular con % ecuación f(x)=0 en el intervalo [a,b]. Puede usarse sólo si f(x) es
% continua y f(a) y f(b) tienen distinto signo
el siguiente programa % LIBRO: MÉTODOS NUMÉRICOS CON MATLAB, MATHEWS Y FINK
en Matlab. function [c,err,yc] = bisect(f,a,b,delta)
% Valores que se van a asignar a bisect(@f,0,2,0.0001)
% Datos
% f es la función, introducida como una cadena de caractéres 'f'
% a y b son el extremo izquierdo y el extremo derecho
% delta es la tolerancia
% Resultados
% c es el cero
% y(c)=f(c)
% err es el error estimado de la aproximación a c
ya = feval(f,a);
yb = feval(f,b);
if ya*yb > 0,% El comando break siempre debe estar acompañado por un
% ciclo for o while
end
max1 = 1 + round((log(b-a) - log(delta))/log(2));
for k = 1:max1
c = (a+b)/2;
yc = feval(f,c);
if yc == 0
a = c;
b = c;
elseif yb*yc > 0
b = c;
yb = yc;
else
a = c;
ya = yc;
end
if b - a < delta,
break,
end
end
c = (a+b)/2;
err = abs(b-a);
yc = feval(f,c);
end
La solución para distintos problemas es
% solución bisección
clc
clear all
f = @(x) x*sin(x)-1;
[c,err,yc] = bisect(f,0,2,0.01)

La solución para el problema antes mostrado es:


c=

1.113281250000000

err =

0.007812500000000

yc =

-0.001216490418016
3. Métodos de la regula falsi o método de la posición falsa
Una de las razones de su introducción es que la velocidad de convergencia del método de
bisección es baja. Como antes, se supone que 𝑓(𝑎) 𝑦 𝑓(𝑏) tienen distinto signo. En el método
de bisección se usa el punto medio del intervalo [a, b] para llevar a cabo el siguiente paso. Suele
conseguirse una aproximación mejor usando el punto (c,0) en el que la recta secante pasa por
los puntos (𝑎, 𝑓(𝑎)) 𝑦 (𝑏, 𝑓 (𝑏)) cruza el eje OX (Figura 3.1).

Figura 3.1. El proceso de decisión para el método de la regula falsi. (a) si f(a) y f(c) tienen signos opuestos,
entonces se recorta por la derecha. (b) si f(a) y f(c) tienen signos opuestos, entonces se recorta por la izquierda.

Para hallar el punto c, se igualan dos fórmulas para la pendiente m de la recta L:


𝑓 (𝑎 ) − 𝑓 (𝑎 )
𝑚= (3.1)
𝑏−𝑎
Que resulta de usar los puntos (𝑎, 𝑓(𝑎)) y (𝑏, 𝑓(𝑏)), y
0 − 𝑓 (𝑏 ) (3.2)
𝑚=
𝑐−𝑏
Que resulta de usar los puntos (𝑐, 0) y (𝑏, 𝑓(𝑏)).
Igualando las pendientes que aparecen en (3.1) y (3.2), se tiene
De donde c puede despejarse, y se obtiene
𝑓 (𝑎 ) − 𝑓 ( 𝑎 ) 0 − 𝑓 (𝑏 )
=
𝑏−𝑎 𝑐−𝑏
De donde c puede despejarse para obtener
𝑓 (𝑏) (𝑏 − 𝑎) (3.3)
𝑐=𝑏−
𝑓 (𝑏 ) − 𝑓 (𝑎 )
Las tres posibilidades son las mismas que antes
Si f(a) y f(c) tienen signos opuestos, entonces hay un cero en [a, c] (3.4)
Si f(c) y f(b) tienen signos opuestos, entonces hay un cero en [c, b] (3.5)
Si f(c)=0, entonces hay un cero. (3.6)
3.1. Covergencia del método de la regula falsi
La fórmula (3.3) junto con el proceso de decisión descrito en (3.4) y (3.5) se usa para construir
una sucesión de intervalos {[𝑎𝑛 , 𝑏𝑛 ]} cada uno de los cuales contiene un cero. En cada paso la
aproximación al cero obtenida es:
𝑓(𝑏𝑛 ) (𝑏𝑛 − 𝑎𝑛 ) (3.7)
𝑐𝑛 = 𝑏𝑛 −
𝑓 (𝑏𝑛 ) − 𝑓(𝑎𝑛 )
Y puede probarse que la sucesión {𝑐𝑛 } converge a un cero r de la función. Sin embargo, hay
que prestar atención al siguiente hecho: aunque la anchura del intervalo 𝑏𝑛 − 𝑎𝑛 se hace más
pequeña, es posible que no tienda a cero; si la curva 𝑦 = 𝑓(𝑥) es convexa cerca de (𝑟, 0),
entonces uno de los extremos 𝑎𝑛 o 𝑏𝑛 permanece estacionario y el otro tiende a la solución
(Figura 3.2).

Figura 3.2. Un extremo estacionario en el método de la regula falsi.


Se hallará la solución de 𝑥 𝑠𝑒𝑛 (𝑥 ) − 1 = 0 usando el método de la regula falsi; observaremos
que converge más rápidamente que el de bisección y se observará también que {𝑏𝑛 − 𝑎𝑛 }∞
𝑛=0

no tiende a cero
Ejemplo 3.1. Se va a utilizar el método de regula falsi para hallar la raíz de 𝑥 𝑠𝑒𝑛(𝑥 ) − 1 = 0
que está en el intervalo [0, 2] (de nuevo con el ángulo x medido en radianes).
Empezando con 𝑎0 = 0 y 𝑏0 = 0, se tiene que 𝑓 (0) = −1.00000000 y 𝑓(0) =
−0.81859485, de manera que hay una raíz en [0, 2]. Usando la ecuación (3.7), se tiene
0.81859485 (2−0)
𝑐0 = 2 − 0.81859485 −(−1) = 1.09975017 y 𝑓(𝑐0 ) = −0.02001921

La función cambia de signo en el intervalo [𝑐0 , 𝑏0 ] = [1.09975017,2], así que se recorta por
la derecha y se pone 𝑎1 = 𝑐0 y 𝑏1 = 𝑐0 . Se utiliza otra vez la ecuación (3.7) para hallar la
siguiente aproximación:
0.81859485 (2−1.09975017)
𝑐0 = 2 − 0.81859485 −(−0.02001921) = 1.12124074 y 𝑓(𝑐1 ) = 0.00983461
Ahora 𝑓(𝑥 ) cambia de signo en [𝑎1 , 𝑐1 ] = [1.09975017, 1.12124075], así que la siguiente
decisión es recortar por la derecha y poner 𝑎2 = 𝑎1 y 𝑏2 = 𝑐1 . Estos cálculos se muestran en la
siguiente Tabla.
Tabla 3.1.
k Extremo izquierdo, Punto medio, 𝑐𝑘 Extremo derecho, Valor de la
𝑎𝑘 𝑏𝑘 función, 𝑓(𝑐𝑘 )
0 0.00000000 1.09975017 2.00000000 -0.02001921
1 1.09975017 1.12124074 2.00000000 0.00983461
2 1.09975017 1.11416120 1.12124074 0.00000563
3 1.09975017 1.11415714 1.11416120 0.00000000

El criterio de parada usado en % Programa (Método de la regula falsi o de la posición falsa). Aproxima
% ción a una raíz de la ecuación h(x)=0 en el intervalo [a,b]. Puede usarse
el método de bisección no es % sólo si f(x) es continua y f(a) y f(b) tienen distinto signo
% LIBRO:MÉTODOS NUMÉRICOS CON MATLAB, MATHEWS Y FINK
útil para el método de regula function [c,err,yc] = regula(h,a,b,delta,epsilon,max1)
% Valores que se van a asignar a [c,err,yc]= regula(@h,0,2,0.001,0.04,4)
falsi por que podría producir % Datos
% f es la función, introducida como una cadena de caractéres 'f'
iteraciones sin fin. En este % a y b son el extremo izquierdo y el extremo derecho
% delta es la tolerancia para el cero
caso, se usan tanto la cercanía % epsilon es la tolerancia para el valor de f en el cero
% max1 es el número máximo de iteraciones
entre sí de dos aproximaciones %Resultados
% c es el cero
sucesivas |𝑐𝑛 − 𝑐𝑛−1 | como el % y(c)=f(c)
% err es el error estimado de la aproximación a c
tamaño |𝑓(𝑐𝑛 )|; así se hace el
ya = feval(h,a);
programa de la regula falsi en yb = feval(h,b);
if ya*yb > 0
Matlab. disp('Note: h(a)*h(b)>0')
end
for k = 1:max1
dx = yb*(b-a)/(yb-ya);
c = b - dx;
ac = c - a;
yc = feval(h,c);
if yc == 0
break
elseif yb*yc > 0
b = c;
yb = yc;
else
a = c;
ya = yc;
end
dx = min(abs(dx),ac);
if abs(dx) < delta
break
end
if abs(yc)<epsilon
break
end
end
c;
err = abs(b-a)/2;
yc = feval(h,c);
end
La solución para el problema ℎ(𝑥 ) = 𝑥 𝑠𝑒𝑛(𝑥) − 1 es
% solución regula falsi
clc
clear all

h= @(x) x*sin(x)-1;
[c,err,yc]= regula(h,0,2,0.001,0.04,4)

La solución para el problema antes mostrado es:


c=

1.099750170294616

err =

0.450124914852692

yc =

-0.020019210242676
4. Los métodos de Newton Raphson y de la secante
A continuación se mostrarán dos métodos importantes y muy utilizados para determinar las
raíces de un polinomio.
4.1. Métodos tangenciales para el cálculo de raíces; Método de Newton
Si 𝑓(𝑥 ), 𝑓′(𝑥) y 𝑓′′(𝑥) son continuas cerca de una raíz 𝑝, esta información adicional sobre la
naturaleza de 𝑓 (𝑥) puede usarse para desarrollar algoritmos que produzcan sucesiones {𝑝𝑘 } que
converjan a 𝑝 más rápidamente que en el método de bisección o en el de regula falsi. El método
de Newton-Raphson (o, simplemente, de Newton), que descansa en la continuidad de 𝑓′(𝑥) y
𝑓′′(𝑥), es uno de los algoritmos más útiles y mejor conocidos. Se introducirá gráficamente y
luego se dará un tratamiento más riguroso basado en el teorema de Taylor.
Suponiendo que la aproximación inicial 𝑝0 está cerca de la raíz 𝑝. Entonces la curva 𝑦 = 𝑓(𝑥 )
y el eje de abscisas se cortan en el punto (𝑝, 0) y, además el punto (𝑝0 , 𝑓(𝑝0 )) está situado en
la curva cerca de (𝑝, 0)(véase 𝐅𝐢𝐠𝐮𝐫𝐚 𝟒. 𝟏). Si se define 𝑝1 como el punto de intersección del
eje de abscisas con la recta tangente a la curva en el punto (𝑝0 , 𝑓(𝑝0 )). En la Figura 4.1 muestra
que 𝑝1 estará, en este caso, más cerca de 𝑝 que 𝑝0 . Por lo que se puede encontrar la ecuación
que relacionan 𝑝1 con 𝑝0 igualando dos fórmulas distintas para la pendiente 𝑚 de la recta
tangente.

Figura 4.1. Construcción geométrica de 𝑝1 y 𝑝2 en el método de Newton-Raphson.


Por un lado
0 − 𝑓(𝑝0 ) (4.1)
𝑚=
𝑝1 − 𝑝0
que es la pendiente de la línea recta que pasa por (𝑝1 ,0) y (𝑝0 , 𝑓(𝑝0 )); por otro lado,
m = 𝑓′(𝑝0 ) (4.2)
que es la pendiente de la recta tangente a la curva en el punto (𝑝0 , 𝑓 (𝑝0 )). Igualando los valores
de la pendiente m en las fórmulas (4.1) y (4.2) y despejando 𝑝1 se obtiene
𝑓(𝑝0 )
𝑝1 = 𝑝0 − (4.3)
𝑓′(𝑝0 )
Este proceso puede repetirse para obtener la sucesión {𝑝𝑘 } que converge a 𝑝. Se realizará
algunas ideas precisas.
Teorema 4.1 (Teorema de Newton-Raphson). Suponiendo que la función 𝑓 ∈ 𝐶 2 [𝑎, 𝑏] y que
existe un número 𝑝 ∈ [𝑎, 𝑏] tal que 𝑓 (𝑝) = 0. Si 𝑓′(𝑝) ≠ 0 entonces existe 𝛿 > 0 tal que la
sucesión {𝑝𝑘 }∞
𝑘=0 definida por el proceso iterativo:

𝑓 (𝑝𝑘−1 ) (4.4)
𝑝𝑘 = 𝑔(𝑝𝑘−1 ) = 𝑝𝑘−1 −
𝑓′(𝑝𝑘−1 )
converge a 𝑝 cualquiera que sea la aproximación inicial 𝑝0 ∈ [𝑝 − 𝛿, 𝑝 + 𝛿].
Observación. La función 𝑔(𝑥) definida por la relación
(4.5)
𝑓 (𝑥 )
𝑔(𝑥) = 𝑥 −
𝑓′(𝑥 )
Se llama función de iteración de Newton-Raphson. Puesto que 𝑓(𝑝) = 0, es fácil ver que
𝑔(𝑝) = 𝑝, por lo que nos dice que la iteración de Newton-Raphson para hallar una raíz de la
ecuación 𝑓 (𝑥 ) = 0 consiste en hallar un punto fijo de 𝑔(𝑥).
Corolario 4.1 (Iteración de Newton-Raphson para el cálculo de las raíces cuadradas).
Suponiendo que A > 0 es un número real y sea 𝑝0 > 0 una aproximación inicial a √𝐴. Se define
la sucesión {𝑝𝑘 }∞
𝑘=0 mediante el proceso recursivo
1 𝐴 (4.6)
𝑝𝑘 = (𝑝𝑘−1 + )
2 𝑝𝑘−1
Entonces la sucesión la sucesión {𝑝𝑘 }∞
𝑘=0 converge a √𝐴; es decir, lim 𝑝𝑘 = √𝐴.
𝑛→∞

Ejemplo ilustrativo 4.1. Se utilizará el algoritmo de Newton-Raphson dado inmediatamente


antes para calcular √5.
Empezando con 𝑝0 = 2 y usando la fórmula 4.6, además, se tiene que 𝐴 = 5 y se calcula:
1 5
𝑝1 = (2 + ) = 2.25
2 2
1 5
𝑝2 = (2.25 + ) = 2.236111111
2 2.25
1 5
𝑝3 = (2.236111111 + ) = 2.236067978
2 2.236111111
1 5
𝑝4 = (2.236067978 + ) = 2.236067978
2 2.236067978
Si se siguiera iterando se llegaría a 𝑝𝑘 ≈ 2.236067978 para 𝑘 > 4, de manera que se ha
conseguido una exactitud de 9 cifras significativas.
4.1.1. Errores por la división entre cero
Uno de los inconvenientes obvios del método de Newton-Raphson es la posibilidad de que se
divida entre cero en la fórmula 4.4, lo que ocurriría si 𝑓′(𝑝𝑘−1 ) = 0. El siguiente programa
contiene las instrucciones que permiten detectar esta situación, pero ¿de qué nos sirve en ese
caso la última aproximación calculada 𝑝𝑘−1 ? Es posible que |𝑓(𝑝𝑘−1 )| sea suficientemente
pequeño y que 𝑝𝑘−1 sea una aproximación aceptable a la raíz. Se investigará esta situación y de
paso se descubrirá un aspecto interesante: la velocidad de convergencia del método.
%Programa 2.5 (Iteración de Newton-Raphson). Aproximación a una raíz de
%t(x)=0 a partir del valor inicial p0 mediante la iteración
%pk=pk-1-f(pk-1)/f'(pk-1) para k=1,2,3,4,.....
% LIBRO: MÉTODOS NUMÉRICOS CON MATLAB, MATHEWS Y FINK
function [p0,err,k,y] = newton(f,df,p0,delta,epsilon,max1)
%Valores que se van a asignar a newton(@t,@dt,1,0.001,0.44,10)
%Datos
% f es la función, introducida como una cadena de caractéres 't'
% df es la derivada de la función f, introducida como una cadena de
% caractéres 'f'
% p0 es la aproximación inicial a un cero de f
% delta es la tolerancia para el cero
% epsilon es la tolerancia para los valores de f
% max1 es el número máximo de iteraciones
%Resultados
% p0 es la aproximación al cero, obtenida con el método de
% Newton-Raphson
% err es una estimación del error de p0
% k es el número de iteraciones realizadas
% y es el valor de la función t(p0)
for k = 1:max1
p1 = p0 - feval(f,p0)/feval(df,p0); %evalúa el método de Newton-Raphson
err = abs(p1-p0);
relerr = 2*err/(abs(p1)+delta);
p0 = p1;
y = feval(f,p0);
if (err<delta)|(relerr<delta)|(abs(y)<epsilon)
break
end
end
end

Definición 4.1 (orden de una raíz). Suponiendo que 𝑓(𝑥) y sus derivadas 𝑓 ′ (𝑥 ), … , 𝑓 𝑀 (𝑥)
están definidas y son continuas en un intervalo centrado en el punto 𝑝. Se dice que 𝑓(𝑥 ) = 0
tiene una raíz de orden 𝑀 en 𝑥 = 𝑝, si
𝑓(𝑝) = 0, 𝑓 ′(𝑝) = 0, … , 𝑓 𝑀−1 (𝑝) = 0 y 𝑓 𝑀 (𝑝) ≠ 0 (4.7)

Las raíces de orden 𝑀 = 1 se suelen llamar raíces simples, mientras que si 𝑀 > 1, entonces se
llaman raíces múltiples. En particular, las raíces de orden 𝑀 = 2 se conocen como raíces
dobles, y así sucesivamente. El siguiente resultado servirá para aclarar estos conceptos.
Lema 4.1. Si una ecuación 𝑓(𝑥 ) = 0 tiene una raíz de orden 𝑀 en 𝑥 = 𝑝, entonces existe una
función continua ℎ(𝑥) tal que 𝑓 (𝑥 ) pueda expresarse como el producto
𝑓(𝑥 ) = (𝑥 − 𝑝)𝑀 ℎ(𝑥 ), 𝑠𝑖𝑒𝑛𝑑𝑜 ℎ(𝑝) ≠ 0 (4.8)

Ejemplo teórico 4.2. La función 𝑓 (𝑥 ) = 𝑥 3 − 3𝑥 + 2 tiene una raíz simple en 𝑝 = −2 y una


raíz doble en 𝑝 = 1. Para verificar este hecho basta considerar las derivadas 𝑓′(𝑥 ) = 3𝑥 2 − 3
y 𝑓′′(𝑥 ) = 6𝑥 y evaluar. Para 𝑝 = −2, tenemos 𝑓(−2) = 0 y 𝑓 ′ (−2) = 9, de manera que 𝑀 =
1 en la definición 4.1 y, por tanto, 𝑝 = −2 es una raíz simple. Para p = 1, tenemos 𝑓(1) =
𝑓′(1) = 0 y 𝑓′′(1) = 6, de manera que 𝑀 = 2 en la definición 4.1 y, por tanto, p = 1 es una
raíz doble. Se debe notar que 𝑓(𝑥) puede factorizarse como 𝑓(𝑥) = (𝑥 + 2)(𝑥 + 1)2 .
Velocidad de convergencia
Como se demostrará, la propiedad distintiva que caracteriza cada caso es la siguiente: si 𝑝 es
una raíz simple de 𝑓 (𝑥) = 0, entonces el método de Newton-Raphson converge rápidamente,
de forma que en cada iteración doblamos (aproximadamente) el número de cifras decimales
exactas. Si, por el contrario, 𝑝 es una raíz múltiple, entonces el error en cada paso es una
fracción del error en el paso anterior. Para describir estos hechos de manera precisa, se definirá
a continuación la noción de orden de convergencia, que es una medida de la velocidad de
convergencia de una sucesión.
Definición 4.2 (orden de convergencia). Suponiendo que {𝑝𝑛 }∞
𝑛=0 converge a p y sea 𝐸𝑛 =

𝑝 − 𝑝𝑛 para cada 𝑛 ≥ 0. Si existen dos constantes positivas 𝐴 > 0 y 𝑅 > 0 tales que
|𝑝 − 𝑝𝑛+1 | |𝐸𝑛+1 | (4.9)
lim 𝑅
= lim =𝐴
𝑛→∞ |𝑝 − 𝑝𝑛 | 𝑛→∞ |𝐸𝑛 |𝑅

entonces se dice que la sucesión converge a p con orden de convergencia R y el número A se


llama constante asintótica del error. Los casos R = 1, 2 merecen una consideración especial:
Si R = 1, la convergencia de {𝑝𝑛 }∞
𝑛=0 se llama lineal
(4.10)
Si R = 2, la convergencia de {𝑝𝑛 }∞
𝑛=0 se llama cuadrática
(4.11)
Si R es grande, entonces la sucesión {𝑝𝑛 } converge rápidamente a p; esto es, la relación (4.9)
implica que para valores grandes de n se tiene que la aproximación |𝐸𝑛+1 | ≈ 𝐴|𝐸𝑛 |𝑅 . Por
ejemplo, suponiendo que R = 2 y que |𝐸𝑛 | ≈ 10−2 , entonces cabe esperar que |𝐸𝑛+1 | ≈
𝐴𝑥10−4 .
Algunas sucesiones convergen con un orden que no es un número natural; se verá, por ejemplo,
(1+√5)
que el orden de convergencia del método de la secante es 𝑅 = ≈ 1.618033989.
2
Ejemplo teórico 4.3 (Convergencia cuadrática cuando la raíz es simple). Partiendo de 𝑝0 =
−2.4 y usando la iteración de Newton-Raphson, se aproximará a la raíz simple 𝑝 = −2 del
polinomio 𝑓 (𝑥 ) = 𝑥 3 − 3𝑥 + 2.
Utilizando el método de Newton-Raphson, con la ecuación:
𝑓 (𝑥 ) 𝑥 3 − 3𝑥 + 2 𝑥(3𝑥 2 − 3) − (𝑥 3 − 3𝑥 + 2)
𝑔 (𝑥 ) = 𝑥 − = 𝑥 − =
𝑓 ′(𝑥) 3𝑥 2 − 3 3𝑥 2 − 3
3𝑥 3 − 3𝑥 − 𝑥 3 + 3𝑥 − 2 2𝑥 3 − 2
= = 2
3𝑥 2 − 3 3𝑥 − 3
La fórmula de iteración en términos de 𝑝, para calcular {𝑝𝑘 } es
2 𝑝𝑘3 − 2 (4.12)
𝑝𝑘+1 = 𝑔(𝑝𝑘 ) =
3 𝑝𝑘2 − 3
Usando la ecuación 4.9 con R = 2 para comprobar si la convergencia es cuadrática se obtienen
los valores que se muestran en la Tabla 4.1.
Tabla 4.1. El método de Newton-Raphson converge cuadráticamente en una raíz simple.
k 𝑝𝑘 𝑝𝑘+1 − 𝑝𝑘 𝐸𝑘 = |𝑝 − 𝑝𝑘 | |𝐸𝑘+1 |
|𝐸𝑘 |
0 -2.400000000 0.323809524 0.400000000 0.476190475
1 -2.076190476 0.072594465 0.076190476 0.619469086
2 -2.003596011 0.003587422 0.003596011 0.664202613
3 -2.000008589 0.000008589 0.000008589
4 -2.000000000 0.000000000 0.000000000

Analizando con atención la velocidad de convergencia en el ejemplo teórico 4.3 se puede ver
que el error en una iteración es proporcional al cuadrado del error en la iteración previa;
concretamente,
|𝑝 − 𝑝𝑘+1 | ≈ 𝐴|𝑝 − 𝑝𝑘 |2
Donde 𝐴 ≈ 2/3. Para comprobar este hecho, se utiliza, por ejemplo, que
|𝑝 − 𝑝3 | = 0.000008589 y |𝑝 − 𝑝2 |2 = |0.003596011|2 = 0.000012931
Con lo cual
2
|𝑝 − 𝑝3 | = 0.000008589 ≈ 0.000008621= |𝑝 − 𝑝2 |2
3

Ejemplo teórico 4.4 (Convergencia lineal cuando la raíz es doble). Partiendo de 𝑝0 = 1.2 y
usando la iteración de Newton-Raphson, se aproximará a la raíz doble 𝑝 = 1 del polinomio
𝑓 (𝑥 ) = 𝑥 3 − 3𝑥 + 2.
Utilizando la ecuación (4.19) con R = 1 para comprobar si la convergencia es lineal, se obtienen
los resultados que se muestran en la Tabla 4.2.
Tabla 4.2. El método de Newton-Raphson converge linealmente en una raíz doble.
k 𝑝𝑘 𝑝𝑘+1 − 𝑝𝑘 𝐸𝑘 = |𝑝 − 𝑝𝑘 | |𝐸𝑘+1 |
|𝐸𝑘 |
0 1.200000000 -0.096969697 -0.200000000 0.515151515
1 1.103030303 -0.050673883 -0.103030303 0.508165253
2 1.052356420 -0.025955609 -0.052356420 0.496751115
3 1.026400811 -0.013143081 -0.026400811 0.509753688
4 1.013257730 -0.006614311 -0.013257730 0.501097775
5 1.006643419 -0.003318055 -0.006643419 0.500550093
⋮ ⋮ ⋮ ⋮ ⋮

Se debe notar que el método de Newton-Raphson converge a la raíz doble, pero con una
velocidad bastante baja. Los valores de 𝑓(𝑝𝑘 ) en el ejemplo teórico 4.4 convergen a cero más
𝑓(𝑝𝑘 )
rápidamente que los valores de 𝑓′(𝑝𝑘 ), de manera que el cociente que aparece en la
𝑓 ′ (𝑝𝑘 )

ecuación 4.4 puede calcularse cuando 𝑝𝑘 ≠ 𝑝. La sucesión converge linealmente y el error


decrece en cada paso hasta, aproximadamente, la mitad de su tamaño, en términos de su orden
de convergencia, para raíces simples y dobles.
Teorema 4.2 (Orden de convergencia del método de Newton-Raphson). Suponiendo que el
método de Newton-Raphson genera una sucesión {𝑝𝑛 }∞
𝑛=0 que converge a un cero 𝑝 de la

función 𝑓(𝑥). Si 𝑝 es una raíz simple, entonces la convergencia es cuadrática:


|𝑓′′(𝑝)|
|𝐸𝑛+1 | ≈ |𝐸𝑛 |2 para n suficientemente grande (4.13)
2|𝑓′(𝑝)|

Si 𝑝 es una raíz múltiple de orden M > 1, entonces la convergencia es lineal:


𝑀−1
|𝐸𝑛+1 | ≈ |𝐸𝑛 |2 para n suficientemente grande (4.14)
𝑀

4.2. El método de la secante


En el algoritmo de Newton-Raphson hay que evaluar dos funciones en cada iteración,
𝑓 (𝑝𝑘−1 ) 𝑦 𝑓′(𝑝𝑘−1 ). Tradicionalmente, el cálculo de la derivada de una función elemental
puede llegar a suponer un esfuerzo considerable. Sin embargo, con los modernos paquetes de
cálculo simbólico, esto no es un problema serio hoy en día. El método de la secante necesita
solo una evaluación de 𝑓 (𝑥 ) por paso y en una raíz simple tiene un orden de convergencia
𝑅 ≈ 1.618033989, es casi tan veloz como el de Newton, cuyo orden de convergencia es de 2.
La ecuación de iteración del método de la secante es la misma que la que aparece en el método
de la regula falsi, la diferencia entre ambos estriba en la estructura lógica de la forma de decidir
cómo se elige el siguiente término. Se parte de dos puntos iniciales (𝑝0 , 𝑓 (𝑝0 )) y (𝑝1 , 𝑓 (𝑝1 ))
cercanos al punto (𝑝, 0) como se muestra en la Figura 4.2, y se define 𝑝2 como la abscisa del
punto de intersección de la recta que pasa por estos dos puntos con el eje OX. En la Figura 4.2,
se muestra que 𝑝2 estará más cerca de 𝑝 que 𝑝0 y que 𝑝1 .

Figura 4.2. Construcción geométrica de 𝑝2 en el método de la secante.

La ecuación que relaciona 𝑝2 , 𝑝1 𝑦 𝑝0 se halla escribiendo la pendiente de la recta en cuestión:


𝑓(𝑝1 ) − 𝑓(𝑝0 ) 0 − 𝑓(𝑝1 ) (4.15)
𝑚= 𝑦 𝑚=
𝑝1 − 𝑝0 𝑝2 − 𝑝1
El primer valor de m en 4.15 es la pendiente de la recta secante que pasa por los dos puntos
iniciales y el segundo valor es la pendiente de la recta que pasa por (𝑝1 , 𝑓(𝑝1 )) y (𝑝2 , 0).
Igualando los miembros derechos de las dos ecuaciones de 4.15 y despejando 𝑝1 = 𝑔(𝑝1 , 𝑝0 ),
se obtiene
𝑓 (𝑝1 )(𝑝1 − 𝑝0 )
𝑝2 = 𝑔(𝑝1 , 𝑝0 ) = 𝑝1 −
𝑓(𝑝1 ) − 𝑓(𝑝0 ) (4.16)
El término general de la sucesión generada por este método viene dado por la ecuación general
de iteración de dos puntos
𝑓(𝑝𝑘 )(𝑝𝑘 − 𝑝𝑘−1 ) (4.17)
𝑝𝑘+1 = 𝑔(𝑝𝑘 , 𝑝𝑘−1 ) = 𝑝𝑘 −
𝑓(𝑝𝑘 ) − 𝑓(𝑝𝑘−1 )
Ejemplo 2.16 (Método de la secante de una raíz simple). Empezando con 𝑝0 = −2.6 y 𝑝1 =
−2.4, se utilizará el método de la secante para aproximarse a la raíz 𝑝 = −2 de la función
𝑓 (𝑥 ) = 𝑥 3 − 3𝑥 + 2.
En este caso, la ecuación de iteración (4.17) es
(𝑝𝑘3 − 3𝑝𝑘 + 2)(𝑝𝑘 − 𝑝𝑘−1 )
𝑝𝑘+1 = 𝑔(𝑝𝑘 , 𝑝𝑘−1 ) = 𝑝𝑘 −
(𝑝𝑘3 − 3𝑝𝑘 + 2) − (𝑝𝑘−1
3
− 3𝑝𝑘−1 + 2)
(𝑝𝑘3 − 3𝑝𝑘 + 2)(𝑝𝑘 − 𝑝𝑘−1 ) (4.18)
= 𝑝𝑘 −
𝑝𝑘3 − 𝑝𝑘−1
3
− 3𝑝𝑘 + 3𝑝𝑘−1
esta ecuación se puede manipular algebraicamente
𝑝𝑘+1 = 𝑔(𝑝𝑘 , 𝑝𝑘−1 )
𝑝𝑘 [𝑝𝑘3 − 𝑝𝑘−1
3
− 3𝑝𝑘 + 3𝑝𝑘−1 ] − [𝑝𝑘4 − 3𝑝𝑘2 + 2𝑝𝑘 − {𝑝𝑘3 𝑝𝑘−1 − 3𝑝𝑘 𝑝𝑘−1 + 2𝑝𝑘−1 }]
=
𝑝𝑘3 − 𝑝𝑘−1
3
− 3𝑝𝑘 + 3𝑝𝑘−1
Reduciendo la ecuación anterior, se tiene:
𝑝𝑘+1 = 𝑔(𝑝𝑘 , 𝑝𝑘−1 )
[𝑝𝑘4 − 𝑝𝑘 𝑝𝑘−1
3
− 3𝑝𝑘2 + 3𝑝𝑘 𝑝𝑘−1 ] − [𝑝𝑘4 − 3𝑝𝑘2 + 2𝑝𝑘 − 𝑝𝑘3 𝑝𝑘−1 + 3𝑝𝑘 𝑝𝑘−1 − 2𝑝𝑘−1 ]
=
𝑝𝑘3 − 𝑝𝑘−1
3
− 3𝑝𝑘 + 3𝑝𝑘−1
𝑝𝑘+1 = 𝑔(𝑝𝑘 , 𝑝𝑘−1 )
𝑝𝑘4 − 𝑝𝑘 𝑝𝑘−1
3
− 3𝑝𝑘2 + 3𝑝𝑘 𝑝𝑘−1 − 𝑝𝑘4 + 3𝑝𝑘2 − 2𝑝𝑘 + 𝑝𝑘3 𝑝𝑘−1 − 3𝑝𝑘 𝑝𝑘−1 + 2𝑝𝑘−1
=
𝑝𝑘3 − 𝑝𝑘−1
3
− 3𝑝𝑘 + 3𝑝𝑘−1
Ya la forma simplificada es
3
−𝑝𝑘 𝑝𝑘−1 − 2𝑝𝑘 + 𝑝𝑘3 𝑝𝑘−1 + 2𝑝𝑘−1
𝑝𝑘+1 = 𝑔(𝑝𝑘 , 𝑝𝑘−1 ) =
𝑝𝑘3 − 𝑝𝑘−1
3
− 3𝑝𝑘 + 3𝑝𝑘−1
𝑝𝑘3 𝑝𝑘−1 − 𝑝𝑘 𝑝𝑘−1
3
− 2𝑝𝑘 + 2𝑝𝑘−1
= 3 3
𝑝𝑘 − 𝑝𝑘−1 − 3𝑝𝑘 + 3𝑝𝑘−1
Debe dar
𝑝𝑘2 𝑝𝑘−1 − 𝑝𝑘 𝑝𝑘−1
2
−2
𝑝𝑘+1 = 𝑔(𝑝𝑘 , 𝑝𝑘−1 ) = 2 2
(4.19)
𝑝𝑘 + 𝑝𝑘 𝑝𝑘−1 − 𝑝𝑘−1 − 3
La sucesión generada aparece en la Tabla 4.3.
Tabla 4.3. Convergencia del método de la secante en una raíz simple.
k 𝑝𝑘 𝑝𝑘+1 − 𝑝𝑘 𝐸𝑘 = |𝑝 − 𝑝𝑘 | |𝐸𝑘+1 |
|𝐸𝑘 |
0 -2.600000000 0.200000000 0.600000000 0.914152831
1 -2.400000000 0.293401015 0.400000000 0.469497765
2 -2.106598985 0.083957573 0.106598985 0.847290012
3 -2.022641412 0.021130314 0.022641412 0.693608922
4 -2.001511098 0.001488561 0.001511098 0.825841116
5 -2.000022537 0.000022515 0.000022515 0.727100987
6 -2.000000022 0.000000022 0.000000022
7 -2.000000000 0.000000000 0.000000000
El método de la secante y el método de Newton-Raphson están relacionados de la siguiente
manera: para una función polinomial 𝑓 (𝑥 ), la ecuación de dos puntos del método de la secante
𝑝𝑘+1 = 𝑔(𝑝𝑘 , 𝑝𝑘−1 ) se reduce, después de simplificar a la ecuación de un solo punto de
Newton-Raphson 𝑝𝑘+1 = 𝑔(𝑝𝑘 ) cuando 𝑝𝑘 sustituye a 𝑝𝑘−1 por 𝑝𝑘 en 4.19, entonces el
miembro derecho de esta ecuación se convierte en el miembro derecho de la ecuación 4.20 dada
en el ejemplo 4.3.
Los términos de la sucesión de los errores verifican, cuando la raíz es simple, la relación
0.618
|𝑓′′(𝑝)| (4.20)
|𝐸𝑘+1 | ≈ |𝐸𝑛 |1.618 | |
2|𝑓′(𝑝)|
1+√5
Siendo el orden de convergencia 𝑅 = ≈ 1.618. Pueden encontrarse demostraciones de
2

este hecho en los libros de análisis numérico avanzado.


Aceleración de la convergencia
Cabe esperar que existan métodos de cálculo de las raíces que tengan un orden de convergencia
mejor que el lineal cuando p sea una raíz múltiple de orden M. El último resultado que se
presentó en esta sección muestra una modificación del método de Newton-Raphson cuya
convergencia pasa a ser cuadrática en una raíz múltiple.
Teorema 2.7 (Iteración de Newton-Raphson acelerada). Suponiendo que el algoritmo de
Newton-Raphson produce una sucesión que converge linealmente a una raíz x = p de orden
M > 1. Entonces la ecuación de iteración de Newton-Raphson acelerada
𝑀 𝑓 (𝑝𝑘−1 )
𝑝𝑘 = 𝑝𝑘−1 −
𝑓′(𝑝𝑘−1 ) (4.21)
Genera una sucesión {𝑝𝑘 }∞
𝑘=0 que converge cuadráticamente a 𝑝.

A continuación se muestra el programa de la secante


% Programa (Método de la secante). Aproximación a una raíz de
%d(x)=0 a partir de unos valores iniciales p0 y p1 mediante la iteración
%pk+1=pk-d(pk)-(pk-pk-1)/(f(pk)-f(pk-1)) para k=1,2,3,4,.....
%LIBRO: MÉTODOS NUMÉRICOS CON MATLAB, MATHEWS Y FINK
function [p1,err,k,y] = secant(d,p0,p1,delta,epsilon,max1)
%Valores que se van a asignar a secant(@d,0,1,0.001,0.44,10)
%Datos
% d es la función, introducida como una cadena de caractéres 'd'
% p0 y p1 son las aproximaciones iniciales a un cero de d
% delta es la tolerancia para p1
% epsilon es la tolerancia para los valores de la función d
% max1 es el número máximo de iteraciones
%Resultados
% p1 es la aproximación al cero, obtenida con el método de la secante
% err es una estimación del error de p1
% k es el número de iteraciones realizadas
% y es el valor de la función d(p1)
for k = 1:max1
p2 = p1-feval(d,p1)*(p1-p0)/feval(d,p0)-feval(d,p0);
err = abs(p2-p1);
relerr = 2*err/(abs(p2)+delta);
p0 = p1;
p1 = p2;
y = feval(d,p1);
if (err<delta)|(relerr<delta)|(abs(y)<epsilon),
break,
end
end
end

También podría gustarte