Solución de Raíces Final
Solución de Raíces Final
Solución de Raíces Final
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).
𝑛→∞
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.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
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
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.
𝐼 𝑁
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 ] = [𝑎𝑛 , 𝑐𝑛 ] 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)
𝑛→∞
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)
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.
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)
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.
𝑓 (𝑝𝑘−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 𝑝𝑘 = √𝐴.
𝑛→∞
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)
𝑝 − 𝑝𝑛 para cada 𝑛 ≥ 0. Si existen dos constantes positivas 𝐴 > 0 y 𝑅 > 0 tales que
|𝑝 − 𝑝𝑛+1 | |𝐸𝑛+1 | (4.9)
lim 𝑅
= lim =𝐴
𝑛→∞ |𝑝 − 𝑝𝑛 | 𝑛→∞ |𝐸𝑛 |𝑅
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
𝑓 ′ (𝑝𝑘 )