Métodos

 Numéricos  
Capítulo  4.  La  resolución  de  ecuaciones  y  de  sistemas  de  
ecuaciones  

Carlos  Beltrán  Álvarez  

Departamento  de  Matemá.cas,  Estadís.ca  y  
Computación  

Este  tema  se  publica  bajo  Licencia:  

Crea.ve  Commons  BY-­‐NC-­‐SA  4.0  
Capı́tulo 4

La resolución de ecuaciones y de
sistemas de ecuaciones

Cuando en el instituto aprendemos a resolver ecuaciones solemos enfrentarnos a ecuaciones cuadráticas de la

forma ax2 + bx + c = 0, para cuyas soluciones existe una conocida fórmula general. Sin embargo, en el ámbito
práctico de la fı́sica la mayorı́a de ecuaciones no son tan sencillas. Podemos encontrarnos con multitud de
situaciones más complicadas como las siguientes:

i) La ecuación a resolver contiene polinomios de grado mayor que 2, o funciones trascendentes (trigonométri-
cas, exponenciales, etc.) con lo que no disponemos de una fórmula para su resolución. Por ejemplo,
podemos tratar de comprender un muelle con rozamiento cuya ecuación (simplificada) viene dada por
ẍ = −cẋ − kx/m donde c es una constante de rozamiento. La solución de esta ecuación depende de
los valores de c, k, m. Si por ejemplo se tiene c2 = 4k/m entonces la solución es x(t) = e−ct/2 (A + Bt)
para algunas (cualesquiera) constantes A, B. Si suponemos que partimos del reposo con elongación de 1
metro, esto es, x(0) = 1, ẋ(0) = 0 la solución es x(t) = e−ct/2 (1 + ct/2). Pongamos que hemos medido
c = 0,01 ¿Podrı́amos decir en qué momento o momentos se alcanza la elongación 1/2? Nos enfrentamos a
la equación
e−0,005t (1 + 0,005t) = 0,5,
que no es un polinomio cuadrático ni se le parece. Siempre podemos intentar dibujar la función con un
paquete de software numérico y tratar de localizar visualmente el punto de corte de la función con el eje de
abscisas pero ¿es eficaz este método cuando es posible que tengamos que resolver el problema un número
grande de veces para distintos parámetros? Por ejemplo, pongamos que quiero comprender, en función
del rozamiento c, el valor de t para el que la elongación es 1/2. Necesito resolver una ecuación parecida
a la anterior para muchos valores de c (pongamos para 100 valores de c entre 0 y 1). ¿Estás dispuesto a
dibujar 100 gráficas y mirar a ojo el punto en el que pasan por x = 1/2 o prefieres aprender a resolver
esta clase de problemas numéricamente con un método eficaz?

ii) Otro ejemplo nos viene dado por una partı́cula que parte del reposo en un plano inclinado uniforme sin
rozamiento, cuyo ángulo θ cambia con el tiempo a velocidad constante ω. Se puede demostrar que después
de t segundos la partı́cula ha recorrido
g
x(t) = − (sinh ωt − sin ωt) metros,
2ω 2

69
70 CAPÍTULO 4. LA RESOLUCIÓN DE ECUACIONES Y DE SISTEMAS DE ECUACIONES

donde aparece una función trigonométrica y otra hiperbólica. Si medimos que en 1 segundo ha recorrido
0,5 metros, ¿cuál es la velocidad de cambio de ángulo ω? De nuevo nos vemos obligados a resolver una
ecuación para la que no hay fórmulas sencillas.

iii) No tenemos una fórmula de la ecuación a resolver, sino que sabemos que tal función existe y solo podemos
calcularla en algunos pocos puntos de nuestra elección. Por ejemplo, nos preguntamos por el peso que
ha de añadirse a un globo de helio en un laboratorio para hacerlo caer al suelo desde 1 metro de altura
de forma que la caı́da dure 2 segundos. En un contexto real en que tenemos una fuerza de rozamiento
(desconocida), no existe una fórmula que en función del peso nos devuelva el tiempo que tarda el globo
en caer, pues el rozamiento depende de la forma concreta del globo. ¿Hay una opción mejor que “ir
probando”? Realizar experimentos es costoso en tiempo (y por lo tanto en dinero) y en materiales (o sea,
en dinero), de modo que mejor realizamos el menor número posible de experimentos. ¿Cómo lo hacemos
para no perder nuestros recursos innecesariamente?
Ilustramos este punto con otro ejemplo: en el diseño de los coches modernos se utiliza la simulación
computacional para tratar de comprender el comportamiento del habitáculo frente a un choque lateral,
frontal, etc. Posteriormente, se fabrican varios prototipos y se realizan pruebas “reales” en las que, con
maniquı́es en lugar de seres humanos, se evaluan los efectos que unas u otras decisiones sobre el diseño
hayan tenido en la cantidad de daño recibido por los ocupantes del vehı́culo. Pongamos que todos los
componentes de un vehı́culo se han decidido ya y que nos preguntamos a qué velocidad máxima puede
producirse un choque frontal sin que pase algo concreto, pongamos a qué velocidad máxima la cabeza
del maniquı́ no se separa de su cuerpo. Realizaremos una serie de pruebas incluyendo algunas velocidades
bajas y otras altas, observando cómo en las bajas todo va bien y en las altas no. Cada prueba lleva
obviamente un gasto muy elevado, ¿qué podemos hacer para tratar de encontrar la velocidad máxima sin
hacer más pruebas de las estrictamente necesarias? Nótese que aquı́ tenemos que encontar el mayor x en
un intervalo tal que f (x) = 0, donde f (x) vale 0 hasta un momento y 1 de ahı́ en adelante. Se trata pues
de encontrar la mayor solución de f (x) = 0. Este problema (ası́ como otros muchos) se resuelve con los
métodos estudiados en este capı́tulo.

iv) La función que define nuestra ecuación puede calcularse en cualquier punto de nuestra elección pero solo
a través de un proceso numérico. Este es un caso mezclado de los dos anteriores, y es difı́cil entenderlo
sin estudiar antes la resolución numérica de ecuaciones diferenciales ordinarias que se verá en este curso.
Por ejemplo, en el estudio del Sistema Solar comprendemos muy bien el movimiento de los planetas, que
se pueden tabular con gran precisión con años de anterioridad, pero los asteroides pequeños son más
difı́ciles de comprender. En la estupenda novela “El martillo de Dios” de Arthur C. Clarke, un equipo
de astronautas trata de desviar un asteroide cuyo curso parece ser directo hacia nuestro querido planeta.
Dado un “vector” aplicado al asteroide, los programas pueden calcular qué efectos tendrá el cambio
en esa velocidad en su posición meses más tarde (cuando se encuentre próximo a la Tierra), esto es,
dado un vector aplicado al asteroide se puede predecir con precisión a qué distancia pasará de la Tierra.
Considerando 40,000 kilómetros una distancia de seguridad, ¿cómo podrı́amos calcular la magnitud del
vector que debemos aplicar? Nuevamente nos encontramos ante el mismo problema que antes, solo que en
lugar de experimentación en laboratorio tenemos resultados computacionales, que (podemos suponer) son
más fáciles y baratos hoy en dı́a que los que nos exigen experimentación en laboratorio. El impacto de un
asteroide potencialmente devastador es una certeza, no una posibilidad, a la que la humanidad habrá de
enfrentarse en el futuro si no es precedida por otra catástrofe comparable. ¡Afortunadamente muchas de
las matemáticas necesarias para saber lo que hay que hacer (que no para poder hacerlo, harina de otro
costal) son relativamente elementales!
4.1. ECUACIONES DE UNA VARIABLE REAL 71

4.1. Ecuaciones de una variable real

En esta sección suponemos que nos es dada una función f : [a, b] → R de la que hemos de encontrar una raı́z,
esto es un ζ tal que f (ζ) = 0. Comenzamos con el método más sencillo.

4.1.1. El método de bisección

Sea f ∈ C[a, b] tal que α = f (a) y β = f (b) tienen signos diferentes. El Corolario 1.1.8 nos garantiza que existe
una solución de f (x) = 0 en el intervalo [a, b], y es nuestra voluntad encontrarla. Para ello, calculemos el punto
medio c y el valor de la función en dicho punto:
a+b
c= , γ = f (c).
2
Si γ y α tienen distinto signo, entonces de nuevo por el Corolario 1.1.8 hay una solución de f (x) = 0 en [a, c]. En
caso contrario, por el mismo argumento, sabemos que hay una solución en [c, b]. Puede parecer que no hemos
avanzado, pero si prestamos atención veremos que sı́: antes buscábamos una solución en un intervalo de longitud
b − a y ahora la buscamos en un intervalo de la mitad de longitud. El método de bisección consiste en repetir
este proceso hasta que se cumple una condición (condición de parada) suficientemente satisfactoria.
La principal propiedad del método de bisección es la siguiente. Llamemos c1 , c2 , . . . a los puntos medios calculados
por el algoritmo. Entonces:
Teorema 4.1.1 Si f ∈ C[a, b] y f (a) y f (b) tienen signo diferente, entonces el método de bisección descrito
arriba genera una sucesión cn , n ≥ 1, que aproxima una raı́z ζ de f en el sentido de que
b−a
|cn − ζ| ≤ , n ≥ 1.
2n
Demostración. La longitud del intervalo de búsqueda se divide por dos en cada paso, luego tras n pasos la
hemos dividido por 2n . Vemos una ilustración gráfica de este método en la Figura 4.1.
Ha de notarse que la cota dada al error es siempre cierta, pero puede ser demasiado conservadora (esto es, puede
suceder que las aproximaciones sean mucho mejores), y también puede suceder que cn sea mejor aproximación
que cn+1 para algún valor de n. Es más, bisección solo calcula la aproximación de una solución pero no nos
dice nada sobre la existencia de otras soluciones, ni tenemos a priori control sobre a qué solución convergerá el
método, por ejemplo, no podemos saber a priori si, de haber varias soluciones, bisección convergerá a la mayor,
la menor, la más cercana al punto medio... no tenemos ninguna información sobre ese punto.
√
Ejemplo 4.1.2 Calculemos una aproximación de 2 usando el método de bisección para la ecuación f (x) =
x2 − 2 en el intervalo [0, 3]. Escribimos en una tabla cada intervalo de búsqueda, el punto medio del intervalo,
y el valor de la función en el punto medio.
Describimos ahora, en código Matlab, el proceso descrito. Lógicamente debemos incluir una condición de parada
razonable, en este caso, dada por el usuario con la precisión .

function c=biseccion(f,a,b,epsilon,M)

% Si los signos de f(a) y f(b) son diferentes se utiliza el metodo de biseccion

% para encontrar, con precision epsilon, una solucion de f(x)=0.
% Tambien incluimos un contador para hacer a lo mAs M iteraciones.
72 CAPÍTULO 4. LA RESOLUCIÓN DE ECUACIONES Y DE SISTEMAS DE ECUACIONES

Figura 4.1: Ilustración de la primera iteración del método de bisección. Comenzamos en el intervalo de búsqueda
[a, b] y en una iteración reducimos el intervalo de búsqueda a [a, c] que tiene la mitad de longitud.

[a, b] c = a+b
2 f (c)
[0, 3] 1,5 0,25
[0, 1,5] 0,75 −1,4375
[0,75, 1,5] 1,125 −0,7344
[1,125, 1,5] 1,3125 −0,2773
[1,3125, 1,5] 1,4062 −0,0225
[1,4062, 1,5] 1,4531 0,1116
[1,4062, 1,4531] 1,4297 0,044

Cuadro 4.1: Método de bisección aplicado a f (x) = x2 − 2 con a = 0, b = 3. La convergencia del punto medio
al valor exacto (aprox. 1,4142) es algo lenta pero sucede de manera segura.

alpha=f(a);
beta=f(b);
if sign(alpha)==sign(beta)
error('El signo de la funciOn debe cambiar en los extremos del intervalo');
end
mierror=1;
contador=1;
while (mierror>epsilon) && (contador<=M)
c=a+(b-a)/2;
gamma=f(c);
if sign(gamma)~=sign(alpha)
b=c;
else
a=c;
alpha=gamma;
end
4.1. ECUACIONES DE UNA VARIABLE REAL 73

mierror=max(abs(gamma),b-a);
contador=contador+1;
end
if contador==(M+1)
error('Se ha alcanzado el numero maximo de iteraciones');
end

Si bien el método de bisección es el de elección en el caso de que la función sea fácil de calcular y conozcamos
dos puntos en que los valores tienen distinto signo, no siempre poseemos esta información. Pongamos que en el
problema f (x) = 0 que estamos tratando de resolver, f (x) tiene la forma aproximada de la función g(x) = x2 −.
Por mucho que busquemos un punto en que f sea negativa, lo tenemos muy difı́cil si  << 1. En el caso extremo
 = 0, es directamente imposible.
Además, el método de bisección depende fuertemente de la propiedad de que los números reales están bien
ordenados. Esto no se cumple ni en el caso de los complejos ni en el caso de que estemos tratando de resolver
sistemas de ecuaciones. Por ello necesitamos otras alternativas.

4.1.2. Métodos de punto fijo

Un punto fijo de una función g es un punto tal que g(p) = p. Observamos que todo problema de “encontrar un
punto fijo” puede escribirse como “resolver la ecuación g(p) − p = 0”, con lo que econtrar puntos fijos deberı́a
ser “más sencillo que” resolver ecuaciones. No obstante, si quiero resolver f (x) = 0 puedo en su lugar buscar un
punto fijo de la función g(p) = p − f (p), con lo que descubrimos que buscar puntos fijos y resolver ecuaciones son
en realidad problemas equivalentes: un método para uno de esos problemas nos proporciona inmediatamente
un método para el otro.
Resulta que los métodos de punto fijo son frecuentemente más efectivos y fáciles de estudiar, por lo que son
una opción muy eficaz para resolver ecuaciones. Comenzamos con un resultado de existencia y, bajo ciertas
condiciones, unicidad:
Teorema 4.1.3 Sea g ∈ C[a, b], y supongamos que g(x) ∈ [a, b] para todo x ∈ [a, b]. Entonces, g tiene un punto
fijo en [a, b]. Adicionalmente, si g es derivable en (a, b) y |g 0 (x)| < 1 para x ∈ (a, b), entonces el punto fijo es
único.
Demostración. Para la primera parte, si g(a) = a o g(b) = b, ya tenemos lo que queremos. En otro caso,
necesariamente g(a) > a y g(b) < b con lo que la función h(x) = g(x) − x es continua, positiva en x = a y
negativa en x = b, luego por el Corolario 1.1.8 tiene un cero, esto es, g tiene un punto fijo, en (a, b).
Para la segunda parte, si hubiese dos puntos fijos p, q por el Teorema 1.1.9 habrı́a un número c ∈ (a, b) tal que

g(p) − g(q)
g 0 (c) = = 1,
p−q
lo que contradice las hipótesis, luego no pueden existir dos puntos fijos diferentes.

El método (o iteración) de punto fijo consiste en tomar un candidato inicial p0 ∈ [a, b] y generar la secuencia

pn = g(pn−1 ), n ≥ 1,

esto es calcular p0 , g(p0 ), g(g(p0 )), . . . Si la secuencia converge a algún punto p ∈ [a, b] y asumiendo que g es
continua tenemos por el Teorema 1.1.4

p = lı́m pn = lı́m g(pn−1 ) = g(p),

n→∞ n→∞
74 CAPÍTULO 4. LA RESOLUCIÓN DE ECUACIONES Y DE SISTEMAS DE ECUACIONES

luego p es un punto fijo de g.

El siguiente resultado (y su corolario) nos da una herramienta para garantizar la convergencia de este método
en algunos casos.

Teorema 4.1.4 Sea g ∈ C[a, b] derivable en (a, b), y supongamos que g(x) ∈ [a, b] para todo x ∈ [a, b]. Si existe
una constante positiva k < 1 tal que |g 0 (x)| < k para x ∈ (a, b), entonces la sucesión generada por la iteración
de punto fijo converge al único punto fijo de g en [a, b]. Es más, |pn − p| ≤ k n |p0 − p|.

Demostración. Por el Teorema 4.1.3 sabemos que el punto fijo existe y es único. Por otro lado, el Teorema
1.1.9 garantiza que (para algunos ζn ∈ (a, b)):

|pn − p| = |g(pn−1 ) − g(p)| = |g 0 (ζn )| |pn−1 − p| ≤ k|pn−1 − p|,

con lo que por inducción |pn − p| ≤ k n |p0 − p|. Como k < 1, tenemos que |pn − p| converge a 0, lo que es
equivalente a que pn converge a p, como querı́amos.

El término de error contiene un factor desconocido (|p − p0 |), que obviamente podemos acotar por el máximo
de los valores b − p0 , p0 − a. También podemos acotar de forma concreta el error calculando una única iteración:

Corolario 4.1.5 En las hipótesis del Teorema 4.1.4, tenemos también

kn
|pn − p| ≤ |p1 − p0 |.
1−k

Demostración. Siguiendo el mismo procedimiento que en la demostración del Teorema 4.1.3, tenemos

|pn+1 − pn | = |g(pn ) − g(pn−1 )| ≤ k|pn − pn−1 | ≤ · · · ≤ k n |p1 − p0 |.

Por lo tanto, para m > n ≥ 1 tenemos:

|pm − pn | ≤|pm − pm−1 | + · · · + |pn+1 − pn |

≤k m−1 |p1 − p0 | + · · · + k n |p1 − p0 |
=k n (1 + · · · + k m−n−1 )|p1 − p0
1 − k m−n
=k n |p1 − p0 |,
1−k
y tomando el lı́mite cuando m → ∞ se obtiene el resultado deseado.


Ejemplo 4.1.6 Resolver la ecuación x2 − 2 = 0 equivale a encontrar un punto fijo de g(x) = x + 1 − x2 /2. Para
x ∈ [1,1, 1,9] se tiene que g(x) ∈ [1,1, 1,9] y además g 0 (x) = 1 − x tiene valor absoluto acotado por 1, con lo que
se cumplen las hipótesis del Teorema 4.1.3 con k = 0,9. Consecuentemente,
√ la iteración de punto fijo aplicada
a g(x) con punto inicial cualquiera en [1,1, 1,9] ha de converger a 2 y cumplir la cota del Corolario 4.1.5. Por
ejemplo, tomando x0 = 1,1 generamos la sucesión:

1,1 1,4950 1,3775

1,4288 1,4081 1,4167 1,4132 1,4146 1,4140 1,4143
√
que aparentemente converge en efecto a 2, como de hecho sabemos que sucede.
4.1. ECUACIONES DE UNA VARIABLE REAL 75

Dada una ecuación f (x) = 0, hay muchas maneras de escribir un problema de punto fijo g(p) = p tal que la
solución a la ecuación y la solución al problema de punto fijo sean equivalentes. De acuerdo con los resultados
vistos, lo óptimo serı́a elegir una función g tal que g 0 sea lo más pequeña (en valor absoluto) posible.
Una implementación en Matlab del método de punto fijo serı́a por ejemplo:

function p=puntofijo(g,p,epsilon,M)
% Dada la funcion g y el punto inicial p, realiza el metodo del punto fijo
% hasta alcanzar una precision epsilon. El programa realiza un maximo de M iteraciones
contador=1;
mierror=1;
while (mierror>epsilon) && (contador<=M)
contador=contador+1;
gp=g(p);
mierror=abs(gp-p);
p=gp;
contador=contador+1;
end
if contador==(M+1)
error('Se ha alcanzado el numero maximo de iteraciones');
end

4.1.3. Métodos de Newton y de la secante

El método de Newton (o Newton-Raphson) es otra de las grandes ideas que uno de los mayores genios de la
historia nos dejó. Para explicarlo, aparte de su conocida interpretación geométrica, supongamos que f ∈ C 2 [a, b]
y que x ∈ (a, b) con f (p) = 0. Por el Teorema 1.3.1, sabemos que para cada x0 ∈ [a, b] tenemos
f 00 (ζx )
f (x) = f (x0 ) + f 0 (x0 )(x − x0 ) + (x − x0 )2 ,
2
para algún ζx ∈ (a, b). Supongamos ahora que x0 está lo suficientemente cerca de p como para poder despreciar
(p − x0 )2 . Evaluando entonces la expresión anterior en x = p tenemos:
0 = f (p) ≈ f (x0 ) + f 0 (x0 )(p − x0 ),
de donde
f (x0 )
p ≈ x0 − .
f 0 (x0 )
Esto es, a partir de una aproximación x0 de la raı́z p, obtenemos una aproximación (siendo optimistas, una
mejor aproximación) de la raı́z p. El método de Newton consiste en repetir este procedimiento, de forma que se
genera la sucesión
f (xn )
xn+1 = xn − 0 .
f (xn )
Nótese que si f 0 (xn ) = 0 para algún n, no podemos definir la siguiente iteración. En ese caso concluimos que
hemos elegido mal x0 . Lo realmente sorprendente de este método es que en numerosas ocasiones converge a
una solución muy rápidamente, incluso si x0 está muy lejos de ella al principio. Lo que podemos afirmar es el
siguiente resultado.
Teorema 4.1.7 Sea f ∈ C 2 [a, b]. Si p ∈ [a, b] es tal que f (p) = 0 y f 0 (p) 6= 0, entonces existe δ > 0 tal que el
método de Newton genera una sucesión definida para todo n ≥ 0 y que converge a p, siempre que x0 se elija en
(p − δ, p + δ).
76 CAPÍTULO 4. LA RESOLUCIÓN DE ECUACIONES Y DE SISTEMAS DE ECUACIONES

Demostración. Consideremos g(x) = x − f (x)/f 0 (x). Como f 0 (p) 6= 0 y f 0 es continua, podemos suponer
(tras pasar a un intervalo [p − δ1 , p + δ1 ] lo bastante pequeño) que f 0 6= 0, con lo que g está bien definida en ese
intervalo. Además,
f (x)f 00 (x)
g 0 (x) = ,
f 0 (x)2

y como f ∈ C 2 [a, b] tenemos que g ∈ C 1 [p − δ1 , p + δ1 ]. Dado que f (p) = 0, tenemos g 0 (p) = 0. Sea k ∈ (0, 1)
un número cualquiera. Podemos encontrar δ < δ1 tal que |g 0 (x)| < k para x ∈ (p − δ, p + δ). Veamos ahora que
g([p − δ, p + δ]) ⊆ [p − δ, p + δ]. En efecto, si |x − p| < δ, por el Teorema 1.1.9 tenemos que

|g(x) − g(p)| = |g 0 (ζ)| |x − p| ≤ kδ < δ.

Todas las hipótesis del Teorema 4.1.3 se cumplen, con lo que existe un único punto fijo (una única solución
de f (x) = 0) en el intervalo (p − δ, p + δ), y la sucesión generada por el método de Newton converge a dicha
solución. El teorema queda ası́ demostrado.
Vemos una ilustración gráfica de este método en la Figura 4.2.

Figura 4.2: Ilustración de las dos primeras iteraciones del método de Newton.

Ejemplo 4.1.8 Utilizar las sucesivas iteraciones del método de Newton para resolver f (x) = x2 − 2 con x0 = 2.
Solución: La fórmula de la iteración es:
2
f (x) x2 − 2 x+ x
x 7→ x − 0
= x − = .
f (x) 2x 2

El cálculo de tres iteraciones produce:

2
1,5 1,4167 1,4142,
√
alcanzándose rápidamente la convergencia a 2. Nótese que
4.1. ECUACIONES DE UNA VARIABLE REAL 77

Ejemplo 4.1.9 Calcular las primeras tres iteraciones del método de Newton para resolver f (x) = 0 donde
f (x) = x6 − x − 1. Tomar como punto inicial de la iteración x0 = 2. Solución: la iteración de Newton tiene la
forma
f (xi ) x6 − xi − 1 5x6 + 1
xi+1 = xi − 0 = xi − i 5 = 5 .
f (xi ) 6xi − 1 6x − 1
Con ello tenemos:
x0 = 2, x1 = 1,6806, x2 = 1,4307, x3 = 1,2549,
que era lo que se pedı́a.

Una implementación en Matlab del método de Newton serı́a por ejemplo:

function x1=newton(f,fp,c,epsilon,M)
% Calcula una solucion aproximada de f(x)=0 con el metodo de Newton
% se comienza la iteracion en x 0=c y se pone un limite de iteraciones
% Tambien se pide al usuario la precision epsilon.
% se pide fp la derivada de f
x1=c;
mierror=1;
contador=1;
while (mierror>epsilon) && (contador<=M)
x0=x1;
x1=x1-f(x1)/fp(x1);
mierror=abs(x0-x1);
contador=contador+1;
end
if contador==(M+1)
error('Se ha alcanzado el numero maximo de iteraciones');
end

Pese a todas las virtudes del método de Newton, en ocasiones su uso resulta complicado. El motivo es que
necesitamos conocer f 0 (x) y no siempre tenemos una fórmula para esta función (basta considerar algunos de
los ejemplos mencionados en la introducción del capı́tulo). Una manera de evitar este problema es sustituir en
el cálculo el término f 0 (xn ) por su aproximación

f (xn ) − f (xn−1 )
f 0 (xn ) ≈ ,
xn − xn−1

lo que transforma la iteración de Newton dando la nueva fórmula

f (xn )(xn − xn−1 )

xn+1 = xn − .
f (xn ) − f (xn−1 )

Ası́ se define el llamado método de la secante, que requiere para ser inicializado dos valores aproximados de la
raı́z x0 , x1 (o bien x0 y f 0 (x0 ) para poder calcular x1 usando el método de Newton).
Vemos una ilustración gráfica de este método en la Figura 4.3.
Ejemplo 4.1.10 Las sucesivas iteraciones del método de la secante para resolver f (x) = x2 −2 con x0 = 2, x1 =
1 producen:
2 1 1,333 1,4286 1,4138 1,4142,
√
alcanzándose rápidamente (aunque algo más lentamente que con el método de Newton) la convergencia a 2.
78 CAPÍTULO 4. LA RESOLUCIÓN DE ECUACIONES Y DE SISTEMAS DE ECUACIONES

Figura 4.3: Ilustración de las dos primeras iteraciones del método de la secante.

Una posible implementación en Matlab de este método serı́a:

function x1=secante(f,c0,c1,epsilon,M)
% Calcula una solucion aproximada de f(x)=0 con el metodo de la secande
% se comienza la iteracion en x 0=c0, x 1=c1 y se pone un limite de iteracione
% Tambien se pide al usuario la precision epsilon.
x0=c0;
fx0=f(x0);
x1=c1;
fx1=f(x1);
mierror=1;
contador=1;
while (mierror>epsilon) && (contador<=M)
xaux=x1;
x1=x1-fx1*(x1-x0)/(fx1-fx0)
mierror=abs(x0-x1);
contador=contador+1;
x0=xaux;
fx0=fx1;
fx1=f(x1);
end
if contador==(M+1)
error('Se ha alcanzado el numero maximo de iteraciones');
end

Otros métodos inspirados en estos que hemos mencionado (por ejemplo el método de Regula Falsi) son de uso
menos frecuente y no los estudiaremos en este curso.
Finalmente señalamos que, aunque el método de Newton y el método de la secante han sido explicados para
funciones de una variable real, también son válidos para funciones de una variable compleja.
4.2. LA RESOLUCIÓN DE SISTEMAS DE ECUACIONES NO–LINEALES 79

4.2. La resolución de sistemas de ecuaciones no–lineales

En muchas ocasiones tenemos que resolver sistemas de ecuaciones no–lineales. Se trata ésta de una tarea
nada sencilla en según qué casos. Además, lo más probable es que el estudiante no se haya enfrentado a ella
anteriormente más que en algunos casos extremadamente sencillos. Pongamos que conocemos la trayectoria de
dos cuerpos en el sistema solar y que, tras simplificar enormemente el problema, concluimos que uno de ellos (la
Tierra) sigue una trayectoria aproximadamente esférica que en las unidades apropiadas y fijándonos únicamente
en el plano de la elı́ptica podemos escribir como

x2 + y 2 = 1, esto es, una circunferencia de radio 1UA.

Ahora consideremos un cometa en trayectoria parabólica cuya órbita obedece a la ecuación

y 2 + x2 + 2xy − y = 1, que es una parábola como se puede ver al dibujarla

La pregunta que nos hacemos es muy simple: ¿en qué punto el cometa atraviesa la órbita de la Tierra?
La pregunta que acabamos de plantear nos demanda simplemente resolver el sistema de ecuaciones
(
x2 + y 2 − 1 = 0
y 2 + x2 + 2xy − y − 1 = 0

esto es resolver f (v) = 0 donde f : R2 → R2 , un sistema de dos ecuaciones con dos variables.
Otro ejemplo, menos atractivo desde el punto de vista filosófico pero mucho más desde el punto de vista
industrial, viene dado por el movimiento de los brazos robóticos en las grandes fábricas. En el problema de los
brazos robóticos (o en el mismo problema de la intersección de las órbitas si lo consideramos en tres coordenadas
espaciales) fácilmente nos requiere resolver sistemas del orden de 20, 30 o 100 ecuaciones no–lineales. ¿Qué
métodos podemos usar?
La teorı́a en este caso es más compleja que en el caso de una variable. Obviamente son necesarios los conceptos
de continuidad y derivabilidad para funciones de varias variables reales, y además las demostraciones de los
resultados requieren un esfuerzo mucho mayor que en el caso anterior. Por ello indicamos simplemente algunos
métodos que, combinados, son suficientes para resolver una considerable cantidad de problemas, sin centrarnos
demasiado en su justificación teórica.

4.2.1. El método de Newton para varias variables

Sea f : Rn → Rn (o f : Cn → Cn ). El método de Newton consiste en tomar v0 ∈ Rn (una aproximación inicial
a la solución) y considerar la iteración:

vn+1 = vn − Jf (vn )−1 f (vn ),

siendo Jf (vn ) la matriz Jacobiana de f en vn . Al igual que en el caso de una variable, tenemos un teorema que
nos garantiza la convergencia bajo ciertas condiciones. Mencionamos aquı́ una versión sencilla de dicho teorema,
debido a Kantorovich.

Teorema 4.2.1 Si f : Rn → Rn es de clase C 2 y si v ∈ Rn con f (v) = 0, det(Jf (v)) 6= 0, entonces existe

un entorno de v tal que para todo v0 en dicho entorno la secuencia de Newton está bien definida y converge
cuadráticamente a v.
80 CAPÍTULO 4. LA RESOLUCIÓN DE ECUACIONES Y DE SISTEMAS DE ECUACIONES

Lamentablemente es difı́cil encontrar un punto v0 para el que este teorema (o sus hermanos mayores) nos
garanticen la convergencia. Afortunadamente, el método de Newton funciona en muchı́simos casos prácticos (más
allá de las garantı́as del Teorema, esto es, para “aproximaciones” v0 bastante groseras). Una implementación
posible para el método de Newton para sistemas de ecuaciones serı́a como sigue.

function x1=newtonvariasvariables(f,Jf,c,epsilon,M)
% Calcula una solucion aproximada de f(x)=0 con el metodo de Newton
% se comienza la iteracion en x 0=c y se pone un limite de iteraciones
% Tambien se pide al usuario la precision epsilon.
% se pide Jf que es la jacobiana de f.
x1=c;
mierror=1;
contador=1;
while (mierror>epsilon) && (contador<=M)
x0=x1;
x1=x1-Jf(x1)\f(x1);
mierror=norm(x0-x1);
contador=contador+1;
end
if contador==(M+1)
error('Se ha alcanzado el numero maximo de iteraciones');
end

Ejemplo 4.2.2 Calcular la fórmula que tiene la iteración de Newton para resolver el sistema de ecuaciones
(
x4 + y 3 = 0
xy = 1

Solución: sea    4
x + y3

x
f =
y xy − 1
La iteración de Newton tiene la fórmula
 −1      3 −1 
3y 2 x4 + y 3
    
x x x x x 4x
7→ − Df f = − .
y y y y y y x xy − 1

La inversa de la matriz es:

−1
4x3 3y 2 −3y 2
  
1 x
= 4 ,
y x 4x − 3y 3 −y 4x3

con lo que la iteración tiene la expresión

 4
−3y 2 x + y3
     
x x 1 x
7→ − 4 =
y y 4x − 3y 3 −y 4x3 xy − 1
 5
x + xy 3 − 3xy 3 + 3y 2
  
x 1
− 4 .
y 4x − 3y 3 −x4 y − y 4 + 4x4 y − 4x3
4.2. LA RESOLUCIÓN DE SISTEMAS DE ECUACIONES NO–LINEALES 81

4.2.2. Método del gradiente o del descenso más rápido

Si bien el método de Newton funciona con mucha frecuencia, en muchos problemas sucede que la convergencia
es solo posible cuando uno elige muy bien la aproximación inicial. Es por ello que a veces conviene usar un
método más lento pero más seguro, al menos que proporcione alguna información sobre la función. El método
del gradiente cumple a la perfección con este papel: no siempre nos devuelve un cero del sistema, pero sı́
que nos proporciona una secuencia de iteraciones x0 , x1 , . . . con la propidedad de que la norma de f (xn ) es
decreciente, y en el contexto teórico se detiene cuando esta norma alcanza un mı́nimo local. Obviamente, un
mı́nimo local de kf (x)k2 no tiene por qué ser una solución de f (x) = 0, pero el recı́proco es cierto: toda
solución de f (x) = 0 es un mı́nimo global de kf (x)k2 . Además, el método del gradiente puede ser utilizado para
minimizar (o maximizar) funciones cualesquiera, independientemende de que correspondan a los ceros de un
sistema o no. Este problema (la minimización de funciones complicadas) es de importancia esencial en Fı́sica,
habiendo llevado a la comprensión de muchos problemas, por ejemplo, el de la estructura cristalina de algunos
materiales.
Comenzamos discutiendo en general el método del gradiente. Dada g : Rn → R diferenciable, y dado v0 ∈ Rn ,
podemos preguntarnos cuál es la dirección de Rn tal que f decrece con mayor rapidez en dicha dirección, si
es que existe una (y si es que es única). Un poco de cálculo de varias variables demuestra rápidamente que la
dirección, de existir, es única, y viene dada precisamente por el opuesto del gradiente en el punto. Esto es, la
dirección de más rápido descenso es precisamente −∇g(v0 ). Un metodo apropiado para minimizar g(v) parece
ser por tanto el tomar v0 inicial (próximo al mı́nimo, si es posible) y mover v0 en la dirección de ∇g(v0 ) una
cantidad apropiada, repitiendo el proceso después. Por lo tanto la iteración ha de ser:

vn+1 = vn − α∇g(vn ), para un valor apropiado de α > 0.

La manera exacta de definir α depende del autor, pero una forma muy apropiada es primero definir h(α) =
g(vn − α∇g(vn )), encontrar por tanteo α∗ ∈ [0, 1] tal que h(α) < h(0) = g(vn ), interpolar la función h por un
polinomio cuadrático en α = 0, α∗ /2, α∗ , y calcular el α que minimiza dicho polinomio, bien en [0, α∗ ], bien
más allá si el valor de h disminuye más. En otras palabras, después de interpolar nos quedamos o bien con el
valor en que se anula la derivada del polinomio cuadrático o bien con el α∗ que hemos encontrado. Todos estos
procesos han de ser controlados con la tolerancia numérica deseada por el usuario.
Un vistazo al capı́tulo 2 nos da una fórmula concreta para hallar α una vez que hayamos encontrado α∗ : la
tabla de valores que ha de cumplirse es la mostrada en la Tabla 4.2.

α 0 α∗ /2 α∗
h(α) g(vn ) g(vn + (α∗ /2)∇g(vn )) g(vn + α∗ ∇g(vn ))

Cuadro 4.2: Tabla de valores que cumple el polinomio de grado 2 que interpola a h(α) en los 3 nodos descritos.

Por lo tanto tenemos el polinomio p(α) interpolante de h:

(α − α∗ /2)(α − α∗ ) α(α − α∗ ) α(α − α∗ )

p(α) = g(vn ) 2
+ g(vn + (α∗ /2)∇g(vn ) 2
+ g(vn + α∗ ∇g(vn ) .
α∗ /2 −α∗ /4 α∗2 /2

La derivada de este polinomio se anula cuando p0 (α) = 0, lo que viene a ser:

g(vn + α∗ ∇g(vn ) − 4g(vn + (α∗ /2)∇g(vn ) + 3g(vn )

α = α∗ .
4g(vn + α∗ ∇g(vn ) − 8g(vn + (α∗ /2)∇g(vn ) + 4g(vn )
82 CAPÍTULO 4. LA RESOLUCIÓN DE ECUACIONES Y DE SISTEMAS DE ECUACIONES

El α que tomamos como óptimo es, de entre este último valor y α∗ , en el que h sea menor. Con todo ello,
podemos afirmar lo siguiente.
Teorema 4.2.3 El método del gradiente tal y como ha sido descrito más arriba, tomando como entradas una
función g : Rn → R diferenciable y un vector v0 ∈ Rn , genera una sucesión v1 , v2 , . . . con la propiedad de que
g(vn+1 ) ≤ g(vn ) para n ≥ 1, deteniéndose (con la consideración de la precisión) en un mı́nimo local de g, o
divergiendo a −∞.
Al igual que en el caso de Newton, puede darse el caso de que calcular el gradiente sea muy complicado.
En ese caso, sustituiremos el valor del gradiente por el de su aproximación estudiada en el capı́tulo 3. Una
implementación posible del método del gradiente (con esta última variación incluida) se muestra a continuación.

function v=gradiente(g,v0,epsilon,M)
% Encuentra un minimo local de g utilizando el metodo del gradiente
% comenzando en v0 con maximo numero de pasos M
% y usando epsilon para la precision
n=length(v0);
v=v0;
mierror=1;
contador=1;
hcero=g(v);
while (mierror>epsilon) && (contador<=M)
gradiente=zeros(n,1);
for j=1:n
ej=zeros(n,1);
ej(j)=1;
gradiente(j)=(g(v+epsilon*ej)-g(v-epsilon*ej))/2/epsilon;
end
mierror=norm(gradiente);
gradiente=gradiente/mierror;
% Con eso queda definido el gradiente (unitario) en el punto actual v
% Ahora buscamos alpha * con la propiedad descrita arriba.
alphastar=1;
halphastar=g(v-gradiente);
otrocontador=1;
while (halphastar>hcero) && (otrocontador<=1/epsilon)
alphastar=alphastar/2;
halphastar=g(v-alphastar*gradiente);
otrocontador=otrocontador+1;
end
% Con eso queda definido alpha *. Ahora buscamos el valor optimo de alpha:
numerador=halphastar-4*g(v-alphastar/2*gradiente)+3*hcero;
denominador=4*halphastar-8*g(v-alphastar/2*gradiente)+4*hcero;
alpha1=halphastar*numerador/denominador;
halpha1=g(v+alpha1*gradiente);
if halpha1<halphastar
alphastar=alpha1;
hcero=halpha1;
else
hcero=halphastar;
end
% Con esto, alphastar es el valor optimo para el tamanno del paso
contador=contador+1;
v=v-alphastar*gradiente;
hcero
end
4.3. EXERCISES. SOLVING F (X) = 0 83

Ahora, dado f : Rn → Rn , podemos usar el método del gradiente para buscar mı́nimos locales de g(v) = kf (v)k2 ,
con la esperanza de que de hecho lo que encontremos sean mı́nimos globales (esto es, soluciones de f (v) = 0) y
la respuesta del gradiente puede ser mejorada grandemente con el método de Newton mediante 2 ó 3 iteraciones
finales.

4.3. Exercises. Solving f (x) = 0

Exercise 4.1
Use the bisection method to find solutions accurate to within 10−2 for x3 − 7x2 + 14x − 6 = 0 on each of the
intervals [0, 1], [1, 3,2] and [3,2, 4].

Exercise 4.2
Use the bisection method to find a solution to within 10−2 for x = tan x + 1.

Exercise 4.3
Let f (x) = (x + 2)(x + 1)2 x(x − 1)3 (x − 2). To which zero of f does the bisection method converge when applied
to the following intervals?
[−1,5, 2,5], [−0,5, 2,4], [−0,5, 3], [−3, −0,5].

If instead of the bisection method we consider the secant method, do we get to the same or to different zeros?

Exercise 4.4
If we want to solve f (x) = 0 and we know f (−1) = −4, f (1) = 10, how many iterations of the bisection method
are needed to get a precision of 10−6 ?

Exercise 4.5
√
Find an approximation to 3 to within 10−3 using the bisection method (hint: consider f (x) = x2 − 3). Do the
same but this time using Newton’s method.

Exercise 4.6
√
Write down a general formula for the Newton iteration for approximating c for c > 0.

Exercise 4.7
The criterion for stopping Newton’s method is sometimes taken as “the difference between two iterations is
sufficiently small”. Can it happen that for a sequence pn of real numbers, |pn − pn+1 | → 0 (that is, the difference
becomes arbitrarily small) but still pn does not converge?

Exercise 4.8
84 CAPÍTULO 4. LA RESOLUCIÓN DE ECUACIONES Y DE SISTEMAS DE ECUACIONES

Use bisection and Newton’s method to find approximate solutions of the following problems in the given intervals:
x3 − 2x2 − 5 = 0, [1, 4]
x = cos x, [0, π/2]
3 2
x + 3x = 1, [−3, −2]
x − 0,8 = 0,2 sin x, [0, π/2]
2 −x −2x
x − 2xe +e = 0, [0, 1]
6x 2 2x 4x 3
e + 3(log 2) e − (log 8)e − (log 2) = 0, [−1, 0]

Exercise 4.9
Find approximations to all the solutions of the given equations in the given intervals, if possible.
x2 sin(1/x) = 0, [0, π]
2x cos x − (x − 2)2 = 0, [0, 4]
(x − 2)2 = log x, [1, 5]
x 2
e = 3x , [0, 5]
sin x = ex , [0, 7]
x3 − 3x2 = −4, [−3, 3]
25x3 − 300x2 + 900x = 2, [−1, 7]
100x3 + 900x2 = −27, [−10, 2]
230x4 + 18x3 + 9x2 − 221x − 9 = 0, [−5, 5].
Compare your results to those of Matlab’s fsolve and fzero.

Exercise 4.10
Produce a Matlab function closest which on input (a, b) ∈ R2 produces the point(s) in the graph of y = x2
which is closest to (a, b).

Exercise 4.11
2 0,4x
The function f (x) = log(x + 1) − e cos(πx) has an infinite number of zeros. Approximate the only negative
zero, the four smallest positive zeros and the 25–th smallest positive zero of f .

Exercise 4.12
Write down a generic program extremalzeros which on input f and [a, b] attempts to find the greatest and
the smallest zero of f in [a, b].

Exercise 4.13
Write down a Matlab function Newtonseveralvariables for Newton’s method in several variables. Use your
program to compute as many solutions as possible for each of the following systems:
( (
−x(x + 1) + 2y = 18 5x2 − y 2 = 0
, ,
2 2
(x − 1) + (y − 6) = 25 y = sin x+cos
4
y
4.3. EXERCISES. SOLVING F (X) = 0 85


2 2
x + 2y − y − 2z = 0
 (
2 2
sin(4πxy) − 2y − x = 0
x − 8y + 10z = 0 , 4π−1

,

 2 4π e2x − e + 4ey 2 − 2ex = 0
x − 7yz = 0

Exercise 4.14
Write down a program which receives a system of equations and an initial point and performs the Gradient
Descent method to find a solution of the system. Apply it to the systems of Exercise (4.13). In matlab, a similar
function is fminunc. Learn how to use it and compare the results with the ones you got.

Exercise 4.15
Find out the angle and speed of a parabolic trajectory from the center of the bottom of a 5m deep and 2m
wide pit, if the impact with the floor must happen at 1m distance from the pit, if an impact angle of −2π/5 is
wanted. Answer: speed= 10,7944m/s, angle= 1,4476 rad.

Exercise 4.16
A bullet has been shot with unknown speed and shooting angle, but we have checked that the impact against
a wall situated 100 away from the shooting point took place at a height of 34,5 meters, and with an angle of
π/10. The weapon was found to be one whose cannon is 20cm long. Assuming no air resistance or wind, and
knowing from the kind of bullet that the initial speed was 800m/s, give an approximate value for the height
of the shooter. You can assume that the shooter was standing up and the arm whose hand pulled the trigger
was straight. You can also assume that the dimensions of the shooter’s body are those of the famous Da Vinci’s
Man of Vitrubio. Answer: 1,877m.

Exercise 4.17
A robot consisting of 3 arms of lengths 3, 2, 1m joined by rotors and based on the origin is asked to point with
the extreme of its smallest (and final) arm to the point (0, 4), while keeping the x component of its center of
mass as close as zero as possible. Which positions should the arms have? Can you get it to be actually at 0?
86 CAPÍTULO 4. LA RESOLUCIÓN DE ECUACIONES Y DE SISTEMAS DE ECUACIONES