Resolución de ecuaciones

Sea, por ejemplo, la ecuación no lineal:

$x = e^{-x}$

La representación gráfica permite saber que la ecuación tiene solución, pudiéndose establecer un intervalo dentro del cual se encuentra ésta (o incluso, directamente, estimarla).

Si se escribe:

$ \begin{array}{c} x - e^{-x} = 0 \\[1ex] f(x) = x - e^{-x} \end{array} $

Entonces un valor de $x$ que haga $f(x) = 0$ resuelve la ecuación planteada.

Así que el problema de dada una función $f$ hallar $x$ tal que $f(x) = 0$ es el que se presenta cuando se quieren resolver ecuaciones. En general, no es resoluble explícitamente, i.e. de forma analítica mediante métodos algebraicos, por lo que se hacen necesarios métodos de resolución aproximados.

Algunos métodos:

1) Método de la bisección:

Sea $f:[a,b] \to \mathbb{R}$ continua y $f(a)$, $f(b)$ con distinto signo, entonces $\mathrel{\exists} x$, $a < x <b$, tal que $f(x) = 0$.

Siendo, siguiendo con la misma función que se tomó como ejemplo, que:

$ \begin{array}{l} f(0) = -1 < 0 \\[1ex] f(1) = 1 - e^{-1} > 0 \end{array} $

Por tanto, la solución, $x$ que hace $f(x)$ cero, está dentro del intervalo (0,1).

A continuación se calcularía el valor de la función en el medio del intervalo. Esto es:

$f(0{,}5) = \text{? } \Biggl\{$ Si $f(0{,}5) > 0 \Rightarrow$ Solución entre $0$ y $0{,}5$. Si $f(0{,}5) < 0 \Rightarrow$ Solución entre $0{,}5$ y $1$.

Se genera, con el nuevo punto, otro intervalo, de la mitad de tamaño, de tal manera que se mantenga el signo opuesto para la función en los extremos.

Siendo:

$f(0{,}5) = -0{,}1065$

Por tanto la solución $x$ es tal que $0{,}5 < x < 1$.

Se repite el procedimiento sucesivamente:

$ \begin{array}{l} f(0{,}750) = 0{,}2776 \Rightarrow 0{,}5 < x < 0{,}75 \\[1ex] f(0{,}625) = 0{,}0897 \Rightarrow 0{,}5 < x < 0{,}625 \\[1ex] f(0{,}5625) = -0{,}0073 \Rightarrow 0{,}5625 < x < 0{,}625 \\[1ex] f(0{,}59375) = 0{,}0415 \Rightarrow 0{,}5625 < x < 0{,}59375 \\[1ex] f(0{,}578125) = 0{,}0172 \Rightarrow 0{,}5625 < x < 0{,}578125 \\[1ex] f(0{,}5703125) = 0{,}0050 \Rightarrow 0{,}5625 < x < 0{,}5703125 \\ \vdots \end{array} $

En cada nueva iteración el punto medio es una mejor aproximación a la solución, ya que ésta se encuentra dentro del intervalo que cada vez se hace más pequeño. Si se escoge el punto medio como solución aproximada, la distancia (el error) entre él y la solución como máximo puede ser la mitad del tamaño del intervalo. También puede darse el caso que el punto medio en una de las iteraciones haga la función cero, siendo la solución.

2) Método de Newton-Raphson:

Sean $f(x)$ una función doblemente derivable, $x_n$ una aproximación de una raíz $\alpha$ de dicha función, i.e. $f(\alpha) = 0$, y el desarrollo en serie de Taylor en torno a $x_n$ tal que:

$f(x) = f(x_n) + f'(x_n)(x-x_n) + \dfrac{f''(\xi(x))}{2}(x-x_n)^2 \quad \xi(x) \in (x,x_n)$

Si $\alpha = x_n + h$, entonces:

$0 = f(\alpha) = f(x_n + h) = f(x_n) + f'(x_n)h + \dfrac{f''(\xi(x_n+h))}{2} h^2$

Si se supone que $h$ es muy pequeño puede despreciarse el último término:

$ \begin{array}{c} 0 \approx f(x_n) + f'(x_n) h \\[1ex] h \approx h_n = -\dfrac{f(x_n)}{f'(x_n)} \end{array} $

Esto da lugar a otra aproximación mejor de la solución:

$ \begin{array}{c} \alpha \approx x_n + h_n = x_n - \dfrac{f(x_n)}{f'(x_n)} \\[1ex] x_{n+1} = x_n - \dfrac{f(x_n)}{f'(x_n)} \end{array} $

Así pues, puede construirse la sucesión:

$ \begin{array}{l} x_0 \\[1ex] x_1 = x_0 - \dfrac{f(x_0)}{f'(x_0)} \\[1ex] x_2 = x_1 - \dfrac{f(x_1)}{f'(x_1)} \\ \vdots \\ x_{n+1} = x_n - \dfrac{f(x_n)}{f'(x_n)} \end{array} $

Siendo que:

$\lim\limits_{n \to \infty} x_n = \alpha$

Usando, por ejemplo, de nuevo, la función anterior:

$ \begin{array}{c} f(x) = x - e^{-x} \\[1ex] f'(x) = 1 + e^{-x} \end{array} $

Entonces:

$ \begin{array}{l} x_0 = 0{,}55 \\[1ex] x_1 = 0{,}55 - \dfrac{f(0{,}55)}{f'(0{,}55)} = 0{,}567090 \\[1ex] x_2 = 0{,}567090 - \dfrac{f(0{,}567090)}{f'(0{,}567090)} = 0{,}567143 \end{array} $

Este método converge más rápidamente.

Es importante que el valor inicial sea muy cercano a la raíz, para hacer más propicio que el método converja, podría no hacerlo, hacia ella.

Para soluciones múltiples, como por ejemplo:

$f(x) = x^2 - 2x + 1 = (x-1)^2$

Donde:

$x = 1$, cero doble.

Si $\,f$ tiene un cero doble en $x = \alpha\,,$ resulta:

$f(\alpha) = 0\,, \, f'(\alpha) = 0\,, \,f''(\alpha) \neq 0$

En general, si $f$ posee un cero múltiple con multiplicidad $m$ (el número de ceros) en $x = \alpha\,,$ entonces:

$f(\alpha) = 0\,, \, f'(\alpha) = 0\,, \, \dotsc \,, \, f^{(m-1)}(\alpha) = 0\,, \, f^{(m)}(\alpha) \neq 0$

La convergencia en una solución o cero múltiple ($m \geq 2$) es mucho más lenta que en una solución o cero simple ($m = 1$).

El desarrollo en serie de Taylor de $f(x)$ en el entorno de la raíz $\alpha$, de multiplicidad $m > 1$ conocida, es:

$ \begin{array}{l} f(x) = (x - \alpha)^m \underbrace{\dfrac{f^{(m)}(\xi(x))}{m!}}_{h(x)} \,, \quad \xi(x) \in (x,\alpha) \\[1ex] f(x) = (x - \alpha)^m h(x) \end{array} $

Sea:

$g(x) = f(x)^{1{∕}m} = (x-\alpha) h(x)^{1{∕}m}$

Siendo pues $\alpha$ raíz simple de $g(x)$, si se aplica el método de Newton-Raphson sobre $g(x)$ la convergencia es más rápida:

$ \begin{array}{c} x_{n+1} = x_n - \dfrac{g(x_n)}{g'(x_n)} \\[1ex] x_{n+1} = x_n - \dfrac{f(x_n)^{1{∕}m}}{\dfrac{1}{m}\,f(x_n)^{1{∕}m-1}\,f'(x_n)} \\[1ex] \boxed{x_{n+1} = x_n - m \dfrac{f(x_n)}{f'(x_n)}} \end{array} $ Fórmula de Newton para soluciones múltiples.


Ejemplo:

$f(x) = x^2 e^{2x} - 2x e^x + 1 = 0$

La solución debe ser la misma que para $x = e^{-x}$, (i.e. 0,567143), aunque con multiplicidad doble, ya que:

$ \begin{array}{c} x^2 e^{2x} - 2xe^x + 1 = 0 \\[1ex] e^{2x} (x^2 - 2xe^{-x} + e^{-2x}) = 0 \\[1ex] e^{2x} (x - e^{-x})^2 = 0 \\[1ex] \end{array} $

La derivada:

$ \begin{align} f'(x) &= 2x^2 e^{2x}+ 2xe^{2x} - 2xe^x - 2e^x = \\[1ex] &= 2e^{2x} (x^2 + x - xe^{-x} - e^{-x}) = \\[1ex] &= 2e^{2x} (x - e^{-x})(x+1) \end{align} $

Siendo la fórmula de Newton-Raphson sin modificar:

$x_{n+1} = x_n - \dfrac{f(x_n)}{f'(x_n)}$

Usándola:

$n$ $x_n$
0 0,55
1 0,558693
2 0,562948
3 0,565053
4 0,566100
5 0,566622
6 0,566883
7 0,567013
$n$ $x_n$
8 0,567078
9 0,567111
10 0,567127
11 0,567135
12 0,567139
13 0,567141
14 0,567142
15 0,567143

Al estar elevado $(x - e^{-x})$ al cuadrado, pasar de solución simple a doble, la convergencia es más lenta.

El método de Newton-Raphson para solución doble:

$x_{n+1} = x_{n} - 2 \dfrac{f(x_n)}{f'(x_n)}$

Empleándola:

$n$ $x_n$
0 0,55
1 0,567387
2 0,567143

La convergencia es mucho más rápida.


En una función que está elevada al cuadrado, por tanto raíz doble, nunca se tendrán valores negativos, siempre es positiva o cero, y, por tanto, el método de la bisección no puede utilizarse. Esto es, por ejemplo:

$y = xe^x - 1$ $y = (xe^x - 1)^2$

3) Método de la secante:

Cuando en el método de Newton-Raphson la derivada es difícil de evaluar puede aproximarse. Esto es:

$f'(x_n) \approx \dfrac{f(x_n) - f(x_{n-1})}{x_n - x_{n-1}}$

Por tanto:

$x_{n+1} = x_n - \dfrac{f(x_n)}{f'(x_n)} \approx x_n - f(x_n) \dfrac{x_n - x_{n-1}}{f(x_n) - f(x_{n-1})}$

Se necesitan pues dos estimaciones iniciales distintas.

Por ejemplo, repitiendo con la misma función, $f(x) = x - e^{-x}$, que antes:

$ \begin{array}{l} x_0 = 0{,}55 \\[1ex] x_1 = 0{,}575 \\[1ex] x_2 = 0{,}575 - f(0{,}575) \dfrac{0{,}575 - 0{,}55}{f(0{,}575) - f(0{,}55)} = 0{,}567168 \\[1ex] x_3 = 0{,}567168 - f(0{,}567168) \dfrac{0{,}567168 - 0{,}575}{f(0{,}567168) - f(0{,}575)} = 0{,}567143 \end{array} $

4) Método iterativo simple:

Cada termino de la sucesión, que busca la solución, proviene del anterior. Esto es:

$ \begin{array}{l} x_0 \\[1ex] x_1 = \varphi(x_0) \\ \vdots \\ x_n = \varphi(x_{n-1}) \\ \vdots \\ \alpha \end{array} $

Donde, siendo $\alpha$ la solución:

$\varphi(\alpha) = \alpha$   Punto fijo.

Convergencia cuando existe:

$\alpha = \lim\limits_{n \to \infty} x_n$

Los errores:

$\varepsilon_n = x_n - \alpha\,, \quad n = 0,1,\dotsc$

Cuanto más rápido desaparezca el error, más rápida es la convergencia.

El punto fijo puede que no exista para una función, y para otra aun existiendo tener que hacer un gran número de sucesiones para llegar a él, o incluso no poder seguir la sucesión.

Por ejemplo, siguiendo con la misma función ya empleada:

$f(x) = x - e^{-x} = 0 \Rightarrow \Biggl\{$ $x = e^{-x} \Rightarrow$ ¿Punto fijo de $g(x) = e^{-x}$? $x = -\ln x \Rightarrow$ ¿Punto fijo de $h(x) = -\ln x\,$?

Entonces:

$x_{n+1} = g(x_n)$
$n$ $x_n$
0 0,55
1 0,576950
2 0,561609
3 0,570291
4 0,565361
$\vdots$$\vphantom{\text{5}}$ $\vdots$$\vphantom{\text{0,895393}}$
18 0,567142
19 0,567144
20 0,567143
$x_{n+1} = h(x_n)$
$n$ $x_n$
0 0,55
1 0,597837
2 0,514437
3 0,664682
4 0,408447
$\vphantom{\vdots}$$\text{5}$ $\vphantom{\vdots}$$\text{0,895393}$
6 0,110493
7 2,202803
8 -0,789731

A pesar de ser equivalentes, en el método iterativo simple una converge y la otra no.


Ejemplo:

$f(x) = e^x - x - 2$

Siendo pues:

$ f'(x) = e^x - 1 \, \left\{ \begin{array}{l} > 0 \text{ si } x > 0 \enspace \text{(Creciente).} \\[1ex] < 0 \text{ si } x < 0 \enspace \text{(Decreciente).} \\[1ex] = 0 \text{ si } x = 0 \enspace \text{(Mínimo).} \end{array} \right. $

Representando:

También:

$ \begin{array}{l} \lim\limits_{x \to +\infty} f(x) = +\infty \\[1ex] \lim\limits_{x \to -\infty} f(x) = +\infty \end{array} $

Centrándose en la solución más pequeña de las dos, por bisección:

$ \begin{array}{ll} f(-1) = e^{-1} - 1 < 0 & \\[1ex] f(-2) = e^{-2} > 0 & -2 < \alpha_1 < -1 \\[1ex] f(-1{,}5) = -0{,}277 & -2 < \alpha_1 < -1{,}5 \\[1ex] f(-1{,}75) = -0{,}076 & -2 < \alpha_1 < -1{,}75 \end{array} $

Una vez localizada más o menos la solución por el método de la bisección, hay que proceder con el método de Newton o el método de la secante. Esto es:

- Newton:

$ \begin{array}{l} x_0 = \dfrac{-2 - 1{,}75}{2} = -1{,}875 \\[1ex] x_1 = -1{,}875 - \dfrac{f(-1{,}875)}{f'(-1{,}875)} = -1{,}84151 \\[1ex] x_2 = -1{,}84151 - \dfrac{f(-1{,}84151)}{f'(-1{,}84151)} = -1{,}84141 \end{array} $

- Secante:

$ \begin{array}{l} x_0 = -1{,}875 \\[1ex] x_1 = -1{,}75 \\[1ex] x_2 = -1{,}75 - f(-1{,}75) \dfrac{-1{,}75 - (-1{,}875)}{f(-1{,}75) - f(-1{,}875)} = -1{,}84111 \\[1ex] x_3 = -1{,}84111 - f(-1{,}84111) \dfrac{-1{,}84111 - (-1{,}75)}{f(-1{,}84111) - f(-1{,}75)} = -1{,}84141 \end{array} $

Tambíen aquí, el punto fijo:

$x = e^x - 2 \Rightarrow g(x) = e^x - 2$

En el método iterativo simple:

$x_n = g(x_{n-1})$

¿Estimación inicial $x_0$?

Según sea el valor inicial, puede suceder que los valores salgan divergentes. Esto es, sucesión:

$x_0, x_1, x_2, \dotsc, x_n \to +\infty$

Dibujando:

Siendo:

$e^x - 2 = x \, \Rightarrow \,$ Intersección $ \left\{ \begin{array}{l} y = x \\[1ex] y = e^x - 2 \end{array} \right. $

También:

$ \begin{array}{l} \lim\limits_{x \to +\infty} e^x - 2 = +\infty \\[1ex] \lim\limits_{x \to -\infty} e^x - 2 = -2 \end{array} $

Observando la gráfica, según donde esté situado el valor inicial, la tendencia de los valores al iterar:

Esto quiere decir que la solución más pequeña se hallará siempre que $x_0 < \alpha_2$, pero $\alpha_2$ no se hallará nunca.

Así pues:

$ \begin{array}{l} x_0 = -1{,}875 \\[1ex] x_1 = g(-1{,}875) = -1{,}84665 \\[1ex] x_2 = g(-1{,}84665) = -1{,}84224 \\[1ex] x_3 = g(-1{,}84224) = -1{,}84154 \\[1ex] x_4 = g(-1{,}84154) = -1{,}84143 \\[1ex] x_5 = g(-1{,}84143) = -1{,}84141 \end{array} $


Si $g(x)$ es una función continua en el intervalo $[a,b]$, tal que $g(x) \in [a,b] \mathrel{\forall}$ $x \in [a,b]$, y derivable en $(a,b)$, tal que $|g'(x)| \leq k < 1 \mathrel{\forall}$ $x \in (a,b)$, entonces existe un único punto fijo de $g$ en $[a,b]$, hacia el cual converge la sucesión:

$ \left\{ \begin{array}{l} x_0 \in [a,b] \\[1ex] x_{n+1} = g(x_n)\,, \enspace n \geq 0 \end{array} \right. $


Ejemplo:

$g(x) = e^x - 2$

Entonces:

$g'(x) = e^x < 1$, si $x < 0$

También:

$g(x) \leq 0$, si $x \leq 0$

Por lo que $g$, según el teorema anterior, con valor igual o inferior a cero sólo tiene un único punto fijo, convergiendo en él, empleando cualquier valor inicial que también sea igual o menor de cero, el metodo iterativo.


Ejemplo:

Sea una mezcla equimolecular de mónoxido de carbono y oxígeno que reaccionan alcanzando el equilibrio, a una temperatura de $\pu{300 K}$ y $\pu{5 atm}$ de presión:

$ \begin{array}{cccccc} \ce{ &CO &\mspace{-18mu}+ &\mspace{-18mu}1/2 O2 &\mspace{-18mu}<=> &\mspace{-12mu}CO2 } \\[1ex] n_{i,0} &a &\mspace{-18mu} &\mspace{-18mu}a &\mspace{-18mu} &\mspace{-12mu}\text{-} \\[1ex] n_{i,\rm eq} &a-\xi &\mspace{-18mu} &\mspace{-18mu}a-\frac{1}{2}\xi &\mspace{-18mu} &\mspace{-12mu}\xi \end{array} $

Si se define:

$\dfrac{n_{\ce{CO},\rm eq}}{n_{\ce{CO},0}} = \dfrac{a-\xi}{a} = x\,,\,$ fracción de $\ce{CO}$ que queda.

Entonces:

$ \begin{array}{l} n_{\ce{CO},\rm eq} = a - \xi = ax \\[1ex] n_{\ce{O2},\rm eq} = a - \dfrac{1}{2} (a - ax) = \dfrac{a}{2}(1 + x) \\[1ex] n_{\ce{CO2},\rm eq} = a - ax = a(1-x) \end{array} $

Por lo que la reacción también puede escribirse como:

$\ce{CO + O2 -> x CO + 1/2 $(1+x)$ O2 + $(1\mathrel{-}x)$ CO2}\,, \quad 0 < x < 1$

La fracción molar de cada especie:

$y_i = \dfrac{n_i}{\displaystyle \sum_i n_i}$

Aquí, en el equilibrio:

$ \begin{align} \sum_i n_i &= n_{\ce{CO},\rm eq} + n_{\ce{O2},\rm eq} + n_{\ce{CO2},\rm eq} = \\[1ex] &= ax + \dfrac{a}{2}(1+x) + a(1-x) = \\[1ex] &= \dfrac{a}{2} (3+x) \end{align} $

Por tanto:

$ \begin{array}{l} y_{\ce{CO}} = \dfrac{\cancel{a}\!x}{\dfrac{\cancel{a}}{2}(3+x)} = \dfrac{2x}{3+x} \\[1ex] y_{\ce{O2}} = \dfrac{\cancel{\dfrac{a}{2}}\!(1+x)}{\cancel{\dfrac{a}{2}}\!(3+x)} = \dfrac{1+x}{3+x} \\[1ex] y_{\ce{CO2}} = \dfrac{\cancel{a}\!(1-x)}{\dfrac{\cancel{a}}{2}(3+x)} = \dfrac{2(1-x)}{3+x} \end{array} $

Entonces, la constante de equilibrio:

$ \begin{align} K_P &= \dfrac{P_{\ce{CO2}}}{P_{\ce{CO}} \cdot P_{\ce{O2}}^{1/2}} \mspace{-12mu} \overset{ \begin{subarray}{c} \text{Dalton} \\ \downarrow \end{subarray} }{=} \mspace{-12 mu} \dfrac{ y_{\ce{CO2}} \mspace{-1mu} \cancel{P} }{ y_{\ce{CO}} \mspace{-2mu} \cancel{P} \! \cdot y_{\ce{O_2}}^{1/2} P^{1/2} } = \\[1ex] &= \dfrac{ \cancel{2}\!(1-x) {∕} \!\cancel{(3+x)} }{ \cancel{2}\!x {∕} \!\cancel{(3+x)}\!\! \cdot ((1+x) {∕} (3+x))^{1/2} P^{1/2} } = \\[1ex] &= \dfrac{(1-x)(3+x)^{1/2}}{x(1+x)^{1/2} P^{1/2}} \end{align} $

Lo que permite obtener una ecuación con la que calcular $x$, si se conoce $K_P$ para unas condiciones de presión y temperatura determinadas.

Siendo $K_P = 3{,}06$ a $\pu{300 K}$ y $\pu{5 atm}$, se tiene pues:

$ \begin{array}{c} 3{,}06 = \dfrac{(1-x)(3+x)^{1/2}}{x(1+x)^{1/2} 5^{1/2}} \\[1ex] 6{,}84237 = \dfrac{(1-x)(3+x)^{1/2}}{x(1+x)^{1/2}} \end{array} $

Que equivale:

$f(x) = 6{,}84237 x(1+x)^{1/2} - (1-x)(3+x)^{1/2} = 0$

Su representación gráfica:

Haciendo:

a) Newton-Raphson:

$x_{n+1} = x_n - \dfrac{f(x_n)}{f'(x_n)}$

Se necesita la primera derivada:

$f'(x) = 6{,}84237 ((1+x)^{1/2} + \dfrac{x}{2}(1+x)^{-1/2}) + (3+x)^{1/2} - \dfrac{1}{2} (1-x)(3+x)^{-1/2}$

Escogiendo $x_0 = 0\,$:

$n$ $x_n$
0 0
1 0,209040
2 0,193049
3 0,192962

b) Secante:

$x_{n+1} = x_n - f(x_n) \dfrac{x_n - x_{n-1}}{f(x_n) - f(x_{n-1})}$

Eligiendo $x_0 = 0$ y $x_1 = 1\,$:

$n$ $x_n$
0 0
1 1
2 0,151819
3 0,184755
4 0,193081
5 0,192961
6 0,192962

c) Siendo:

$ \begin{array}{c} (a-b) = 0 \\[1ex] (a-b)(a+b) = 0 \\[1ex] a^2 - b^2 = 0 \end{array} $

Por tanto, otra ecuación equivalente para hallar la solución:

$g(x) = 46{,}818 x^2 (1+x) - (1-x)^2 (3+x) = 0$

Si se aplica en el método de Newton-Raphson:

$x_{n+1} = x_n - \dfrac{g(x_n)}{g'(x_n)}$

La derivada:

$g'(x) = 46{,}818(2x(1+x) + x^2) + 2(1-x)(3+x) - (1-x)^2$

Entonces:

$n$ $x_n$
0 0
1 0,6
2 0,358908
3 0,236327
4 0,197142
5 0,193007
6 0,192962

Converge más lentamente que antes.

d) Punto fijo:

$h(x) = 0{,}146148 \dfrac{(1-x)(3+x)^{1/2}}{(1+x)^{1/2}} = x$

Por tanto:

$x_{n+1} = h(x_n)$

$n$ $x_n$
0 0
1 0,253136
2 0,175868
3 0,197944
4 0,191520
5 0,193379
$n$ $x_n$
6 0,192841
7 0,192996
8 0,192952
9 0,192964
10 0,192961
11 0,192962

- Orden de convergencia:

Sea la sucesión:

$x_0, x_1, x_2, \dotsc, x_n, \dotsc$, con límite $\alpha$.

Siendo el error:

$\varepsilon_n = x_n - \alpha$

De existir, con $p$ y $C$ constantes positivas, el siguiente límite:

$\lim\limits_{n \to \infty} \dfrac{|\varepsilon_{n+1}|}{|\varepsilon_n|^p} = C \neq 0,\infty$

Entonces, si $n$ es grande $|\varepsilon_{n+1}|$ es aproximadamente proporcional a $|\varepsilon_n|^p$, donde $p$, conocido como el orden de convergencia, es una medida de la rapidez, velocidad con que el valor de $x$ se acerca, converge a $\alpha$, i.e. de cuantos pasos (iteraciones) se han de hacer para que $\varepsilon_n$ sea menor que una cierta cantidad.

Casos: