Hasta ahora, hemos estudiado dos propiedades fundamentales de un método numérico para resolver un problema de valores iniciales \(y'=f(t,y),\ a\leq t\leq b,\ y(a)=y_0\):
Dicho en otras palabras, un método es consistente cuando la ecuación en diferencias que define el método “tiende” a la ecuación diferencial del problema de valores iniciales cuando \(h\) tiende a cero.
En todo el estudio de las propiedades anteriores, no hemos tenido en cuenta el error de redondeo que se comete en las operaciones básicas (sumas, restas, multiplicaciones y divisiones) ni el posible error que tengamos en la condición inicial \(y_0\) o en la evaluación de la función \(f\).
Pensemos que, en la práctica, las operaciones no son exactas ni tampoco lo son los valores numéricos con los que trabajamos.
Por dicho motivo, hemos de estudiar la estabilidad del método numérico en cuestión, es decir, ¿cómo afectan los errores de redondeo en los valores iniciales y en las operaciones básicas a los valores aproximados \(\hat{y}(t_i)\) en los puntos del mallado?
Dicho de otra forma, si tenemos pequeñas perturbaciones en los valores iniciales, ¿qué resultado producen dichas perturbaciones en los valores aproximados, \(\hat{y}(t_i)\) en los puntos del mallado?
Entonces, un método numérico de un paso o multipaso es estable si el algoritmo para calcular los valores aproximados, \(\hat{y}(t_i)\) en los puntos del mallado es estable (ver el tema de errores en la parte I del curso de Métodos Numéricos).
Es decir, pequeñas perturbaciones en las condiciones iniciales provocan pequeños cambios en los resultados finales que son en nuestro caso, los valores aproximados, \(\hat{y}(t_i)\) en los puntos del mallado.
El siguiente resultado nos da una condición suficiente para que un método numérico de un paso sea estable dando además un resultado adicional sobre la consistencia y convergencia:
Consideramos un problema de valores iniciales \(y'=f(t,y),\ a\leq t\leq b,\ y(a)=y_0\). Consideramos el método explícito siguiente para hallar la solución aproximada \(\hat{y}(t)\) fijado un paso \(h\) en los puntos del mallado \(t_i=a+ih\), \(i=0,\ldots,n\), donde la relación entre \(h\) y \(n\) es \(h=\frac{b-a}{n}\): \[ \hat{y}_{i+1}=\hat{y}_i+h\varphi(t_i,\hat{y}_i,h), \] donde \(\hat{y}_i=\hat{y}(t_i)\) y \(\hat{y}_0=y_0\).
Supongamos que existe un valor \(h_0>0\) tal que la función del método numérico \(\varphi(t,\hat{y},h)\) es continua y verifica la condición de Lipschitz con respecto la variable \(\hat{y}\) con constante de Lipschitz \(L\) en el conjunto: \[ D=\{(t,\hat{y},h)\ |\ a\leq t\leq b,\ -\infty <\hat{y}<\infty,\ 0\leq h\leq h_0\}. \]
Entonces el método es estable.
Recordemos el método de Euler modificado: \[ \hat{y}_{i+1}=\hat{y}_i+h f\left(t_i+\frac{h}{2},\hat{y}_i+\frac{h}{2}f(t_i,\hat{y}_i)\right),\ \hat{y}_0 =y_0. \] Veamos que es estable suponiendo que la función \(f(t,y)\) del problema de valores iniciales verifica la condición de Lipschitz respecto la variable \(y\) con constante de Lipschitz \(L\), es decir, para todo \((t,y), (t,y'),\ a\leq t\leq b,\ -\infty <y,y'<\infty\): \(|f(t,y)-f(t,y')|\leq L|y-y'|\).
Para ello aplicaremos el Teorema anterior.
En este caso, la función \(\varphi\) vale \(\varphi(t,\hat{y},h)=f\left(t+\frac{h}{2},\hat{y}+\frac{h}{2}f(t,\hat{y})\right)\). Veamos que verifica la condición de Lipschitz respecto la variable \(\hat{y}\): \[ \begin{align*} |\varphi(t,\hat{y})-\varphi(t,\hat{y}')|= & \left|f\left(t+\frac{h}{2},\hat{y}+\frac{h}{2}f(t,\hat{y})\right)-f\left(t+\frac{h}{2},\hat{y}'+\frac{h}{2}f(t,\hat{y}')\right)\right|\\\leq & L\left|\hat{y}+\frac{h}{2}f(t,\hat{y})-\left(\hat{y}'+\frac{h}{2}f(t,\hat{y}')\right)\right|\leq L\left(|\hat{y}-\hat{y}'|+\frac{h}{2}\left|f(t,\hat{y})-f(t,\hat{y}')\right|\right)\\\leq & L\left(|\hat{y}-\hat{y}'|+\frac{h_0}{2}L|\hat{y}-\hat{y}'|\right)=L\left(1+\frac{1}{2}h_0 L\right)|\hat{y}-\hat{y}'|. \end{align*} \]
En conclusión, la función \(\varphi(t,\hat{y},h)=f\left(t+\frac{h}{2},\hat{y}+\frac{h}{2}f(t,\hat{y})\right)\) verifica la condición de Lipschitz con respecto la variable \(\hat{y}\) con constante de Lipschitz \(L\left(1+\frac{1}{2}h_0 L\right)\).
Entonces, usando el teorema anterior, podemos afirmar que el método de Euler modificado es estable.
Recordemos que un método multipaso para resolver el problema de valores iniciales \(y'=f(t,y),\ a\leq t\leq b,\ y(a)=y_0\) tenía como ecuación en diferencias:
\[ \begin{align*} \hat{y}_{i+1} & = a_{m-1}\hat{y}_i +a_{m-2}\hat{y}_{i-1}+\cdots +a_0 \hat{y}_{i+1-m}+\\ &\quad h(b_m f(t_{i+1},\hat{y}_{i+1})+b_{m-1}f(t_i,\hat{y}_i)+\cdots +b_0 f(t_{i+1-m},\hat{y}_{i+1-m})), \end{align*} \] para \(i=m-1,m,\ldots,n-1\), donde se tienen que especificar los valores iniciales siguientes: \[ \hat{y}_0=y_0,\ \hat{y}_1=\tilde{y}_1,\ldots,\hat{y}_{m-1}=\tilde{y}_{m-1}. \]
Llamemos parte lineal a \(a_{m-1}\hat{y}_i +a_{m-2}\hat{y}_{i-1}+\cdots +a_0 \hat{y}_{i+1-m}\) del método multipaso y \[ \begin{align*} & h G(t_i,h,\hat{y}_{i+1},\hat{y}_i,\ldots,\hat{y}_{i+1-m}) := \\ &\quad h(b_m f(t_{i+1},\hat{y}_{i+1})+b_{m-1}f(t_i,\hat{y}_i)+\cdots +b_0 f(t_{i+1-m},\hat{y}_{i+1-m})), \end{align*} \] a la parte no lineal.
Para estudiar la estabilidad de un método multipaso, vamos a considerar la aplicación del método multipaso al problema trivial de valores iniciales \(y'=0,\ a\leq t\leq b,\ y(a)=y_0\) con solución \(y(t)=y_0\) para todo \(t\in [a,b]\).
El método multipaso aplicado al problema anterior será: \[ \hat{y}_{i+1} = a_{m-1}\hat{y}_i +a_{m-2}\hat{y}_{i-1}+\cdots +a_0 \hat{y}_{i+1-m}, \] ya que \(G\equiv 0\) al valer \(f\equiv 0\).
El resultado siguiente nos da la solución de la ecuación en dferencias anterior:
Sea \[ \hat{y}_{i+1} = a_{m-1}\hat{y}_i +a_{m-2}\hat{y}_{i-1}+\cdots +a_0 \hat{y}_{i+1-m}, \] una ecuación en diferencias lineal de orden \(m\). Sea \[ P(\lambda)=\lambda^m-a_{m-1}\lambda^{m-1}-a_{m-2}\lambda^{m-2}-\cdots -a_1\lambda-a_0, \] el polinomio característico asociado a la ecuación anterior.
Sean \(\lambda_1,\ldots,\lambda_m\) los ceros del polinomio característico, es decir \(P(\lambda_i)=0\), \(i=1,\ldots,m\). Supongamos que los ceros anteriores son distintos dos a dos.
Entonces la solución general de la ecuación en diferencias es la siguiente: \[ \hat{y}_k=\sum_{j=1}^m \alpha_j \lambda_j^k, \] donde \(\alpha_1,\ldots,\alpha_m\) son constantes libres que se calculan a partir de las condiciones iniciales.
Como la ecuación en diferencias es lineal, la demostración de la proposición anterior se reduce a demostrar que las funciones \(\hat{y}_k =\lambda_j^k\) son soluciones de la ecuación en diferencias del método multipaso, siendo \(\lambda_j\) un cero del polinomio característico \(P(\lambda)\): \(P(\lambda_j)=0\).
Hemos de demostrar por tanto que: \[ \begin{align*} \hat{y}_{i+1}-(a_{m-1}\hat{y}_i +a_{m-2}\hat{y}_{i-1}+\cdots +a_0 \hat{y}_{i+1-m})= & 0\\ \lambda_j^{i+1}-(a_{m-1}\lambda_j^i +a_{m-2}\lambda_j^{i-1}+\cdots +a_0 \lambda_j^{i+1-m})= & 0,\\ \lambda_j^{i+1-m}\left(\lambda_j^m-(a_{m-1}\lambda_j^{m-1} +a_{m-2}\lambda_j^{m-2}+\cdots +a_0)\right)= & 0,\\ P(\lambda_j) = &0. \end{align*} \] La última igualdad es cierta ya que \(\lambda_j\) es un cero del polinomio característico \(P(\lambda)\).
Recordemos que la solución del problema de valores iniciales trivial \(y'=0,\ a\leq t\leq b,\ y(a)=y_0\) era \(y(t)=y_0\).
Entonces, la solución \(\hat{y}_n=y_0\) debe ser una solución de la ecuación en diferencias eligiendo como valores iniciales \(\hat{y}_0 =\hat{y}_1 =\cdots =\hat{y}_{m-1}=y_0\): \[ \begin{align*} \hat{y}_{i+1} = & a_{m-1}\hat{y}_i +a_{m-2}\hat{y}_{i-1}+\cdots +a_0 \hat{y}_{i+1-m},\ \Rightarrow\\ y_0 = & a_{m-1} y_0 +a_{m-2} y_0+\cdots +a_0 y_0, \ \Rightarrow\\ 1 = & a_{m-1} +a_{m-2} +\cdots +a_0,\ \Rightarrow \\ 1-(a_{m-1} +a_{m-2} +\cdots +a_0) = & 0,\ \Rightarrow P(1)=0, \end{align*} \] lo que significa que un cero del polinomio característico es siempre el valor \(\lambda =1\).
Las soluciones de la ecuación en diferencias se escriben de la forma: \[ \hat{y}_k = y_0+\sum_{j=2}^m \alpha_j \lambda_j^k, \] donde \(\alpha_2,\ldots,\alpha_k\) son constantes libres que deben ser \(0\) para la solución de nuestro problema trivial de valores iniciales.
Ahora bien, en la práctica, las constantes \(\alpha_j\), \(j=2,\ldots,m\) no serán cero debido a los errores de redondeo en \(y_0\) y en las operaciones básicas.
Fijémonos que, en la expresión anterior, si una constante \(\alpha_j\neq 0\) y \(|\lambda_j| >1\), el valor de \(\hat{y}_k\) crecerá exponencialmente hacia infinito haciendo que el método multipaso no sea estable.
Para que el método multipaso sea estable, necesitamos que todos ceros del polinomio característico \(P(\lambda)\), \(\lambda_j\) cumplan que \(|\lambda_j|\leq 1\), \(j=2,\ldots,m\) ya que, en caso contrario, como hemos comentado, se rompe la estabilidad del método.
Hemos tratado el caso en que todos los ceros \(\lambda_1,\lambda_2,\ldots,\lambda_m\) son distintos dos a dos.
En el caso en que haya un cero de multiplicidad \(p\), \(\lambda_l=\lambda_{l+1}=\cdots =\lambda_{l+p}\), la “parte” de la solución de la ecuación en diferencias correspondiente a este cero vale: \[ \alpha_l \lambda_l^k +\alpha_{l+1}k\lambda_l^{k-1} +\alpha_{l+2}k(k-1)\lambda_l^{k-2}+\cdots +\alpha_{l+p}k(k-1)\cdots (k-p+1)\lambda_l^{k-p}, \] es decir, suponiendo que los demás ceros son distintos, la solución general de la ecuación en diferencias vale: \[ \hat{y}_k =\sum_{j\neq l}\alpha_j\lambda_j^k +\alpha_l\lambda_l^k+\sum_{j=1}^{p}\alpha_{l+j} k(k-1)\cdots (k-j+1)\lambda_l^{k-j}. \]
Fijémonos que el número de parámetros libres vale \(m\): \(m-(p+1)+p+1=m\).
Veamos la demostración de la expresión anterior. Antes necesitamos el lema siguiente:
Sea un valor \(\lambda\) un cero de multiplicidad \(p\) del polinomio \[ P(\lambda)=\lambda^m-a_{m-1}\lambda^{m-1}-a_{m-2}\lambda^{m-2}-\cdots -a_1\lambda-a_0. \] Entonces, el valor \(\lambda\) es un cero de multiplicidad \(p\) del polinomio: \[ P_i(\lambda)=\lambda^{i+1-m} P(\lambda)=\lambda^{i+1}-a_{m-1}\lambda^{i}-a_{m-2}\lambda^{i-1}-\cdots -a_1\lambda^{i+2-m}-a_0\lambda^{i+1-m}. \]
Para demostrar el lema, aplicaremos la fórmula de Leibnitz para calcular la derivada \(j\)-ésima de un producto de funciones: \[ (u(\lambda)\cdot v(\lambda))^{(j)}=\sum_{q=0}^j \binom{j}{q} (u(\lambda))^{(j-q)}\cdot (v(\lambda))^{(q)}. \] Supongamos que \(\lambda\) es un cero de multiplicidad \(p\) de \(P(\lambda)\). Entonces \(P(\lambda)=P'(\lambda)=\cdots = P^{(p-1)}(\lambda)=0\).
Si hacemos la derivada \(j\)-ésima del polinomio \(P_i(\lambda)=\lambda^{i+1-m}P(\lambda)\) con \(j\leq p-1\), aplicando la fórmula de Leibnitz con \(u(\lambda)=\lambda^{i+1-m}\) y \(v(\lambda)=P(\lambda)\), obtenemos: \[ P_i^{(j)}(\lambda)=(\lambda^{i+1-m}P(\lambda))^{(j)}=\sum_{q=0}^j \binom{j}{q} (\lambda^{i+1-m})^{(j-q)}\cdot (P(\lambda))^{(q)} \] Como \((P(\lambda))^{(q)}=0\), para \(q=0,\ldots,j\) debido a que \(\lambda\) es un cero de multiplicidad \(p\) del polinoimo \(P(\lambda)\), tendremos que \(P_i^{(j)}(\lambda)=0\), \(j=0,\ldots,p-1\) y concluimos que \(\lambda\) es un cero de multiplicidad \(p\) del polinomio \(P_i(\lambda)\).
Lo único que tenemos que demostrar es que si \(\lambda_l\) es un cero de multiplicidad \(p\) del polinomio característico \(P(\lambda)\), es decir \(P(\lambda_l)=P'(\lambda_l)=\cdots = P^{(p-1)}(\lambda_l)=0\), entonces, las funciones siguientes: \[ \hat{y}_{k,1} =k\lambda_l^{k-1},\ \hat{y}_{k,2}=k(k-1)\lambda_l^{k-2},\ldots,\hat{y}_{k,p-1}=k(k-1)\cdots (k-p+1)\lambda_l^{k-p}, \] son soluciones de la ecuación en diferencias del método multipaso.
Decir que \(\hat{y}_{k,r} =k(k-1)\cdots (k-r+1)\lambda_l^{k-r}\), \(r=1,\ldots,p-1\) es solución de la ecuación en diferencias del método multipaso es equivalente a: \[ \begin{align*} & \hat{y}_{i+1,r}-(a_{m-1}\hat{y}_{i,r} +a_{m-2}\hat{y}_{i-1,r}+\cdots +a_0 \hat{y}_{i+1-m,r})= 0,\\ & (i+1)\cdots (i-r+2)\lambda_l^{i+1-r}-(a_{m-1}i\cdots (i-r+1)\lambda_l^{i-r}+a_{m-2}(i-1)\cdots (i-r)\lambda_l^{i-r-1}\\ & \quad +\cdots +a_0 (i+1-m)\cdots (i-m-r+2)\lambda_l^{i-m-r+1})= 0. \end{align*} \]
La última expresión es la derivada \(r\)-ésima del polinomio \[ P_i(\lambda)=\lambda^{i+1-m} P(\lambda)=\lambda^{i+1}-a_{m-1}\lambda^{i}-a_{m-2}\lambda^{i-1}-\cdots -a_1\lambda^{i+2-m}-a_0\lambda^{i+1-m}, \] evaluada en \(\lambda =\lambda_l\) ya que: \[ \begin{align*} & P_i^{(r)}(\lambda)=(i+1)\cdots (i-r+2)\lambda^{i+1-r}-(a_{m-1}i\cdots (i-r+1)\lambda^{i-r}+a_{m-2}(i-1)\cdots (i-r)\lambda^{i-r-1}\\ & \quad +\cdots +a_0 (i+1-m)\cdots (i-m-r+2)\lambda^{i-m-r+1}). \end{align*} \] Usando el lema anterior, podemos afirmar que \(\lambda_l\) es un cero de multiplicidad \(r\) del polinomio \(P_i(\lambda)\). Por tanto \(P_i^{(r)}(\lambda_l)=0\), \(r=1,\ldots,p-1\), tal como queríamos demostrar.
Entonces en el caso de ceros múltiples, si el cero \(\lambda_l\) verifica que \(|\lambda_l|\geq 1\), la solución de la ecuación en diferencias crecerá en forma polinomial y la solución tenderá a infinito ya que: \[ \lim_{k\to \infty}k (k-1)\cdots (k-p+1)|\lambda_l|^{k-p}=\infty. \] En resumen, para que un método multipaso sea estable, necesitamos que el polinomio característico verifique la denominada condición de raíz:
Sea \(P(\lambda)\) un polinomio de grado \(m\) y sean \(\lambda_1,\ldots,\lambda_m\) las (no necesariamente distintas) raíces de dicho polinomio: \(P(\lambda_j)=0,\ j=1,\ldots,m\).
Diremos que dicho polinomio verifica la condición de raíz si \(|\lambda_j|\leq 1\), \(j=1,\ldots,m\) y si existe alguna raíz \(|\lambda_l|\) de módulo unidad, \(|\lambda_l|=1\), dicha raíz debe ser simple.
Existen tres tipos de estabilidad para un método multipaso:
Estudiemos si el método de Adams-Bashford de orden \(4\) es estable: \[ \hat{y}_{i+1} = \hat{y}_i +\frac{h}{24}(55f(t_i,\hat{y}_i)-59 f(t_{i-1},\hat{y}_{i-1})+37 f(t_{i-2},\hat{y}_{i-2})-9f(t_{i-3},\hat{y}_{i-3})) \] La parte lineal del método es \(\hat{y}_i\). Al ser un método multipaso de orden \(4\), el polinomio característico vale: \(P(\lambda)=\lambda^4-\lambda^3\).
Las raíces del polinomio anterior son \(\lambda=0\) de multiplicidad \(3\) y \(\lambda=1\).
Dicho polinomio verifica la condición de raíz y además todas las raíces excepta la trivial \(\lambda=1\) tienen módulo menor que \(1\). Por tanto, es fuertemente estable.
Como la parte lineal del método de Adams-Moulton de orden \(3\) coincide con la parte lineal del método de Adams-Bashford de orden \(4\), podemos decir también que el método de Adams-Moulton es fuertemente estable.
En todos los métodos numéricos para resolver problemas de valores iniciales \(y'=f(t,y),\ a\leq t\leq b,\ y(a)=y_0\), en las expresiones de los errores cometidos tanto en el error local de truncamiento como en el error respecto la solución exacta de la ecuación diferencial aparecen derivadas de órdenes superiores de la solución exacta \(y(t)\).
Si dichas derivadas están controladas o están acotadas, podemos decir que el método numérico funciona ya que el término del error está acotado o controlado.
Sin embargo, existen casos en que la solución exacta tiende a cero o se hace muy pequeña pero en cambio sus derivadas crecen exponencialmente de tal forma que éstas son dominantes en los cálculos.
Cuando ocurre el fenómeno anterior, diremos que nos encontramos en una ecuación diferencial o un sistema diferencial tipo Stiff.
Este tipo de ecuaciones o sistemas aparecen en el estudio de vibraciones, reacciones químicas o circuitos eléctricos.
Los sistemas diferenciales tipo Stiff se caracterizan por el hecho de que en una de sus soluciones exactas aparecen términos del tipo \(\mathrm{e}^{\lambda t}\) con \(\lambda <0\).
Dicho término se denomina término transitorio y no tiene importancia en el cálculo de la solución exacta ya que decrece relativamente rápido a \(0\) a medida que \(t\) aumenta.
El término importante en la solución exacta es el llamado término estacionario.
Sin embargo, aunque el término transitorio \(\mathrm{e}^{\lambda t}\) tiene a cero, sus derivadas no lo hacen con la misma velocidad.
La derivada \(n\)-ésima de dicho término vale \(\lambda^n\mathrm{e}^{\lambda t}\) que no tiende a cero tan rápidamente como el término transitorio y son las responsables de los problemas al aplicar un método numérico.
Además, como las derivadas que aparecen en los términos del error del método numérico en cuestión están evaluadas no en \(t\) sino en un valor entre \(0\) y \(t\), pueden crecer muy rápidamente.
Para ver cómo se comporta un sistema diferencial tipo Stiff consideremos el problema de valores iniciales siguiente: \[ \left. \begin{align*} y_1' & =32y_1+66y_2+\frac{2}{3}t+\frac{2}{3},\\ y_2' & =-66y_1-133y_2-\frac{1}{3}t-\frac{1}{3}, \end{align*} \right\} \] con condiciones iniciales \(y_1(0)=y_2(0)=\frac{1}{3}\), \(0\leq t\leq 0.5\) y con solución exacta: \[ \begin{align*} y^{(1)}(t) &= \frac{2}{3}t+\frac{2}{3}\mathrm{e}^{-t}-\frac{1}{3}\mathrm{e}^{-100t},\\ y^{(2)}(t) &= -\frac{1}{3}t-\frac{1}{3}\mathrm{e}^{-t}+\frac{2}{3}\mathrm{e}^{-100t}. \end{align*} \]
El término transitorio en las funciones anteriores sería \(\frac{2}{3}\mathrm{e}^{-t}-\frac{1}{3}\mathrm{e}^{-100t}, -\frac{1}{3}\mathrm{e}^{-t}+\frac{2}{3}\mathrm{e}^{-100t}\) y el término estacionario, \(\frac{2}{3}t, -\frac{1}{3}t\).
Vamos a aplicarle el método de Runge-Kutta 4 con \(h=0.1\).
La siguiente tabla muestra los valores de las soluciones aproximada y exacta en los puntos del mallado:
| \(t_i\) | \(\hat{y}_i^{(1)}\) | \(\hat{y}_i^{(2)}\) | \(y_{i}^{(1)}\) | \(y_{i}^{(2)}\) | \(\|y_i-\hat{y}_i\|_2\) |
|---|---|---|---|---|---|
| 0 | 0.3333333 | 0.3333333 | 0.3333333 | 0.3333333 | 0 |
| 0.1 | -96.3301083 | 193.6650542 | 0.6698765 | -0.3349155 | 216.8985599 |
| 0.2 | -28226.3208461 | 56453.660423 | 0.6791538 | -0.3395769 | 63117.4908008 |
| 0.3 | -8214056.3061211 | 16428113.6530605 | 0.6938788 | -0.3469394 | 18367189.8230579 |
| 0.4 | -2390290586.28645 | 4780581173.64323 | 0.7135467 | -0.3567733 | 5344852238.50988 |
| 0.5 | -695574560816.264 | 1391149121633.63 | 0.7376871 | -0.3688436 | 1555352001406.38 |
Observad los errores cometidos.
Si en lugar de considerar \(h=0.1\), consideramos \(h=0.05\) el efecto es prácticamente el mismo:
| \(t_i\) | \(\hat{y}_i^{(1)}\) | \(\hat{y}_i^{(2)}\) | \(y_{i}^{(1)}\) | \(y_{i}^{(2)}\) | \(\|y_i-\hat{y}_i\|_2\) |
|---|---|---|---|---|---|
| 0 | 0.3333333 | 0.3333333 | 0.3333333 | 0.3333333 | 0 |
| 0.05 | -3.9019582 | 8.8051457 | 0.6652403 | -0.3292512 | 10.2125662 |
| 0.1 | -61.969576 | 124.9439894 | 0.6698765 | -0.3349155 | 140.0660738 |
| 0.15 | -858.0088963 | 1717.0285005 | 0.6738052 | -0.3369025 | 1920.0728916 |
| 0.2 | -11770.429547 | 23541.8778248 | 0.6791538 | -0.3395769 | 26320.9992256 |
| 0.25 | -161361.595907 | 322724.2206147 | 0.6858672 | -0.3429336 | 360817.0310515 |
| \(t_i\) | \(\hat{y}_i^{(1)}\) | \(\hat{y}_i^{(2)}\) | \(y_{i}^{(1)}\) | \(y_{i}^{(2)}\) | \(\|y_i-\hat{y}_i\|_2\) |
|---|---|---|---|---|---|
| 0.3 | -2212007.2521086 | 4424015.5450354 | 0.6938788 | -0.3469394 | 4946200.1339975 |
| 0.35 | -30322941.556452 | 60645884.167592 | 0.7031254 | -0.3515627 | 67804160.1702161 |
| 0.4 | -415676999.42816 | 831353999.926639 | 0.7135467 | -0.3567733 | 929482029.000045 |
| 0.45 | -5698238876.21747 | 11396477753.5226 | 0.7250854 | -0.3625427 | 12741649480.8756 |
| 0.5 | -78113357937.3498 | 156226715875.806 | 0.7376871 | -0.3688436 | 174666778300.337 |
Si consideramos \(h=0.025\), es cuando los errores quedan controlados
| \(t_i\) | \(\hat{y}_i^{(1)}\) | \(\hat{y}_i^{(2)}\) | \(y_{i}^{(1)}\) | \(y_{i}^{(2)}\) | \(\|y_i-\hat{y}_i\|_2\) |
|---|---|---|---|---|---|
| 0 | 0.3333333 | 0.3333333 | 0.3333333 | 0.3333333 | 0 |
| 0.025 | 0.4507274 | 0.098855 | 0.6395116 | -0.2787133 | 0.4221342 |
| 0.05 | 0.5273292 | -0.053429 | 0.6652403 | -0.3292512 | 0.3083786 |
| 0.075 | 0.5776126 | -0.1524816 | 0.6683113 | -0.3338791 | 0.2028085 |
| 0.1 | 0.6109596 | -0.2170818 | 0.6698765 | -0.3349155 | 0.1317421 |
| 0.125 | 0.6334509 | -0.2594049 | 0.6716634 | -0.3358298 | 0.0854457 |
| 0.15 | 0.6490261 | -0.2873442 | 0.6738052 | -0.3369025 | 0.0554078 |
| \(t_i\) | \(\hat{y}_i^{(1)}\) | \(\hat{y}_i^{(2)}\) | \(y_{i}^{(1)}\) | \(y_{i}^{(2)}\) | \(\|y_i-\hat{y}_i\|_2\) |
|---|---|---|---|---|---|
| 0.175 | 0.6602369 | -0.3060168 | 0.6763047 | -0.3381523 | 0.0359286 |
| 0.2 | 0.6687349 | -0.318739 | 0.6791538 | -0.3395769 | 0.0232975 |
| 0.225 | 0.6755881 | -0.32766 | 0.6823441 | -0.3411721 | 0.015107 |
| 0.25 | 0.6814863 | -0.3341719 | 0.6858672 | -0.3429336 | 0.0097959 |
| 0.275 | 0.686874 | -0.3391759 | 0.6897147 | -0.3448574 | 0.006352 |
| 0.3 | 0.6920368 | -0.3432554 | 0.6938788 | -0.3469394 | 0.0041189 |
| 0.325 | 0.6971571 | -0.3467869 | 0.6983516 | -0.3491758 | 0.0026708 |
| \(t_i\) | \(\hat{y}_i^{(1)}\) | \(\hat{y}_i^{(2)}\) | \(y_{i}^{(1)}\) | \(y_{i}^{(2)}\) | \(\|y_i-\hat{y}_i\|_2\) |
|---|---|---|---|---|---|
| 0.35 | 0.7023509 | -0.3500137 | 0.7031254 | -0.3515627 | 0.0017319 |
| 0.375 | 0.7076906 | -0.353092 | 0.7081929 | -0.3540964 | 0.001123 |
| 0.4 | 0.713221 | -0.356122 | 0.7135467 | -0.3567733 | 0.0007282 |
| 0.425 | 0.7189687 | -0.3591676 | 0.7191799 | -0.3595899 | 0.0004722 |
| 0.45 | 0.7249485 | -0.3622689 | 0.7250854 | -0.3625427 | 0.0003062 |
| 0.475 | 0.7311679 | -0.3654508 | 0.7312567 | -0.3656284 | 0.0001985 |
| 0.5 | 0.7376295 | -0.3687284 | 0.7376871 | -0.3688436 | 0.0001287 |
En el ejemplo anterior, vimos que cuando nos encontramos con un sistema diferencial tipo Stiff, el paso \(h\) tiene que estar controlado en el sentido que existe un umbral \(h_0\) de forma que \(h\) debe cumplir \(h\leq h_0\) para que los errores cometidos por el método numérico en cuestión estén controlados.
El estudio de la aplicación de un método numérico a un sistema diferencial tipo Stiff se puede predecir estudiando el error cometido cuando el método es aplicado a la llamado ecuación test: \[ y'=\lambda y,\ y(0)=y_0, \] cuando \(\lambda <0\) y solución exacta \(y(t)=y_0\mathrm{e}^{\lambda t}\).
Notemos que \(\displaystyle\lim_{t\to\infty}y(t)=\lim_{t\to\infty}y_0\mathrm{e}^{\lambda t}=0.\)
Ilustremos dicho estudio usando el conocido método de Euler: \[ \hat{y}_{i+1}=\hat{y}_i +hf(t_i,\hat{y}_i),\ \hat{y}_0=y_0. \]
Si aplicamos el método de Euler a la ecuación test, obtenemos la ecuación en diferencias siguiente: \[ \hat{y}_{i+1}=\hat{y}_i +h\lambda\hat{y}_i =(1+h\lambda)\hat{y}_i,\ \hat{y}_0=y_0. \] La solución de la ecuación en diferencias anterior es: \[ \hat{y}_{i+1}=(1+h\lambda)^{i+1}y_0,\ i=0,\ldots,n-1. \] Los errores cometidos serán, pues: \[ |y(t_{i})-\hat{y}_i|=\left|y_0\mathrm{e}^{\lambda t_i}-(1+h\lambda)^{i}y_0\right|=\left|\mathrm{e}^{\lambda h i}-(1+h\lambda)^{i}\right|\cdot |y_0|. \]
Como \(\displaystyle\lim_{i\to\infty}\mathrm{e}^{\lambda h i}=0\), los errores cometidos estarán controlados si \(\displaystyle\lim_{i\to\infty} \left|(1+h\lambda)^{i}\right|=0\), es decir, \(|1+h\lambda|<1\).
La condición anterior es equivalente a: \[ |1+h\lambda|<1,\ \Leftrightarrow -1<1+h\lambda <1,\ \Leftrightarrow -2<h\lambda <0,\ \Leftrightarrow h<\frac{2}{|\lambda|}. \] Entonces el valor \(h_0\) para el método de Euler sería \(\frac{2}{|\lambda|}\).
Otra manera de ver el estudio realizado es comprobar que el método de Euler es estable en la ecuación test si \(h<\frac{2}{|\lambda|}\). Veámoslo.
Supongamos que cometemos un error de redondeo \(\delta_0\) en el valor inicial \(y_0\), es decir \(\hat{y}_0 =y_0+\delta_0\).
Entonces el valor “real” del valor de la solución aproximada \(\hat{y}_i\) valdrá: \[ \hat{y}_i =(1+h\lambda)^i (y_0+\delta_0)=(1+h\lambda)^i y_0 +(1+h\lambda)^i \delta_0. \] Por tanto, el error de redondeo cometido al evaluar \(\hat{y}_i\) será \((1+h\lambda)^i \delta_0\) que tenderá a cero si \(\frac{2}{|\lambda|}\), o, dicho en otras palabras, el método de Euler será estable.
En general, si aplicamos la ecuación en diferencias de un método de un paso a la ecuación test, obtendremos una expresión del tipo: \[ \hat{y}_{i+1}=Q(h\lambda)\hat{y}_i, \] donde \(Q(h\lambda)\) es una función que depende del método tratado.
La solución de la ecuación en diferencias anterior será: \(\hat{y}_i = Q(h\lambda)^i y_0\).
Entonces, usando el mismo razonamiento que en el método de Euler, para que los errores cometidos estén controlados, el paso \(h\) debe cumplir \(|Q(h\lambda)|<1\).
Por ejemplo, si consideramos un método de Taylor de orden \(k\), \[ \hat{y}_{i+1}= \hat{y}_i+hf(t_i,\hat{y}_i)+\frac{h^2}{2}f'(t_i,\hat{y}_i)+\cdots +\frac{h^k}{k!}f^{(k-1)}(t_i,\hat{y}_i), \] podemos probar, usando un proceso inductivo, que \[ f(t,\lambda)=\lambda y,\ f'(t,\lambda)=\lambda y'=\lambda^2 y,\ldots, f^{(k-1)}(t,\lambda)=\lambda^k y. \] Entonces, podemos escribir la ecuación en diferencias de la manera siguiente: \[ \begin{align*} \hat{y}_{i+1}= & \hat{y}_i+h\lambda \hat{y}_i+\frac{h^2}{2}\lambda^2\hat{y}_i+\cdots +\frac{h^k}{k!}\lambda^k\hat{y}_i,\\ \hat{y}_{i+1} = & \left(1+h\lambda+\frac{h^2\lambda^2}{2}+\cdots +\frac{h^k\lambda^k}{k!}\right)\hat{y}_i =Q(h\lambda)\hat{y}_i. \end{align*} \]
La función \(Q(h\lambda)\) será en este caso: \[ Q(h\lambda)=\left(1+h\lambda+\frac{h^2\lambda^2}{2}+\cdots +\frac{h^k\lambda^k}{k!}\right). \] Los errores cometidos al aplicar un método de Taylor de orden \(k\) a un sistema diferencial Stiff estarán controlados si \(h\) verifica: \[ \left|1+h\lambda+\frac{h^2\lambda^2}{2}+\cdots +\frac{h^k\lambda^k}{k!}\right|<1. \]
Recordemos que un método multipaso para resolver el problema de valores iniciales \(y'=f(t,y),\ a\leq t\leq b,\ y(a)=y_0\) tenía como ecuación en diferencias: \[ \begin{align*} \hat{y}_{i+1} & = a_{m-1}\hat{y}_i +a_{m-2}\hat{y}_{i-1}+\cdots +a_0 \hat{y}_{i+1-m}+\\ &\quad h(b_m f(t_{i+1},\hat{y}_{i+1})+b_{m-1}f(t_i,\hat{y}_i)+\cdots +b_0 f(t_{i+1-m},\hat{y}_{i+1-m})), \end{align*} \] para \(i=m-1,m,\ldots,n-1\), donde se tienen que especificar los valores iniciales siguientes: \[ \hat{y}_0=y_0,\ \hat{y}_1=\tilde{y}_1,\ldots,\hat{y}_{m-1}=\tilde{y}_{m-1}. \]
Si aplicamos dicho método a la ecuación test \(y'=\lambda y,\ t\geq 0,\ y(a)=y_0\), obtenemos: \[ \begin{align*} & \hat{y}_{i+1} = a_{m-1}\hat{y}_i +a_{m-2}\hat{y}_{i-1}+\cdots +a_0 \hat{y}_{i+1-m}+\\ &\quad h(b_m \lambda\hat{y}_{i+1}+b_{m-1}\lambda\hat{y}_i+\cdots +b_0 \lambda\hat{y}_{i+1-m}),\\ & (1-h\lambda b_m)\hat{y}_{i+1}-(a_{m-1}+h\lambda b_{m-1})\hat{y}_i -(a_{m-2}+h\lambda b_{m-2})\hat{y}_{i-1}\\ & \quad -\cdots - (a_0+h\lambda b_0)\hat{y}_{i+1-m}=0. \end{align*} \] Obtenemos una ecuación en diferencias lineal.
Recordemos que las soluciones de una ecuación en diferencias lineal vienen dadas en función del llamado polinomio característico que, en este caso, será: \[ \begin{align*} Q(z,h\lambda)= & (1-h\lambda b_m)z^m-(a_{m-1}+h\lambda b_{m-1})z^{m-1} -(a_{m-2}+h\lambda b_{m-2})z^{m-2} \\ & \quad -\cdots - (a_0+h\lambda b_0). \end{align*} \]
Concretamente, sean \(\beta_1(h\lambda),\ldots,\beta_m(h\lambda)\) las raíces del polinomio característico. Si las \(\beta_i(h\lambda)\) son distintas dos a dos, la solución de la ecuación en diferencias anterior será: \[ \hat{y}_i = \sum_{j=1}^m\alpha_j \beta_j(h\lambda)^i,\ i=0,\ldots,n, \] donde las \(\alpha_j\), \(j=1,\ldots,m\) se obtienen a partir de las condiciones iniciales del método multipaso \[ \hat{y}_0=y_0,\ \hat{y}_1=\tilde{y}_1,\ldots,\hat{y}_{m-1}=\tilde{y}_{m-1}. \]
Entonces, los errores cometidos del método multipaso en cuestión estarán controlados si las raíces \(\beta_j(h\lambda)\) están acotadas por \(1\): \(|\beta_j(h\lambda)|<1\), \(j=1,\ldots,m\).
Por tanto, el paso \(h\) debe elegirse de tal forma que se cumplan las \(m\) condiciones siguientes: \(|\beta_j(h\lambda)|<1\), \(j=1,\ldots,m\).
En el caso en que haya una raiz de multiplicidad \(p\), \(\beta_l(h\lambda)\), los términos correspondientes a dicha raíz serán: \[ \begin{align*} & \alpha_l\beta_l(h\lambda)^i+\alpha_{l+1}i\beta_l(h\lambda)^{i-1}+\alpha_{l+2}i(i-1)\beta_l(h\lambda)^{i-2}\\ & \quad +\cdots +\alpha_{l+p}i(i-1)\cdots (i-p+1)\beta_l(h\lambda)^{i-p}. \end{align*} \] En este caso, fijémonos que si \(|\beta_l(h\lambda)|<1\), \(j=1\ldots,m\), los términos anteriores también tenderán a cero cuando \(i\to\infty\).
Por tanto, la condición \(|\beta_j(h\lambda)|<1\) es suficiente también para ceros múltiples.
La región de estabilidad absoluta es una región en el plano complejo de tal forma que si \(h\lambda\) está en dicha región, para \(h>0\), el método se puede aplicar a un sistema Stiff:
La región de estabilidad absoluta para un método numérico de un paso vale: \[ R=\{h\lambda\in\mathbb{C}\ |\ |Q(h\lambda)|<1\}, \] y para un método numérico multipaso vale: \[ R=\{h\lambda\in\mathbb{C}\ |\ |\beta_j(h\lambda)|<1\}, \] para todos los ceros \(\beta_k(h\lambda)\) del polinomio característico del método \(Q(z,h\lambda)\).
La región anterior dice que \(h\lambda\) tiene que estar en \(R\) para que la solución aproximada del problema de valores iniciales tienda a cero cuando se aplica a un sistema Stiff.
Los métodos más vulnerables a los sistemas Stiff son los métodos con control de paso ya que aunque los errores de truncamiento nos permitan aumentar el paso \(h\), la región de estabilidad absoluta dice que el paso \(h\) no puede ser mayor que un cierto umbral \(h_0\).
La idea es hallar métodos numéricos tanto de un paso como multipaso con una región de estabilidad absoluta lo más grande posible.
En este sentido, un método numérico se denomina A-estable si la región de estabilidad absoluta contiene el plano complejo izquierdo: \[ \{z\in\mathbb{C}\ |\ \mathrm{Re\ }z\leq 0\}\subseteq R. \] En este caso, podemos aplicar cualquier valor de \(h>0\) a un sistema diferencial Stiff para tener los errores cometidos controlados ya que como \(h\lambda\in \{z\in\mathbb{C}\ |\ \mathrm{Re\ }z\leq 0\}\subseteq R\), se verificará \(h\lambda\in R\), para cualquier \(h>0\).
Hallemos la región de estabilidad absoluta para el método del trapecio: \[ \hat{y}_{i+1} =\hat{y}_i+ \frac{h}{2}(f(t_i,\hat{y}_i)+f(t_{i+1},\hat{y}_{i+1})). \] Recordemos que se trata de un método implícito.
Si aplicamos dicho método a la ecuación test obtenemos: \[ \hat{y}_{i+1} =\hat{y}_i+ \frac{h}{2}(\lambda\hat{y}_i+\lambda \hat{y}_{i+1}),\ \Rightarrow \left(1-\frac{h\lambda}{2}\right)\hat{y}_{i+1}-\left(1+\frac{h\lambda}{2}\right)\hat{y}_i,\ \Rightarrow \hat{y}_{i+1}=\frac{1+\frac{h\lambda}{2}}{1-\frac{h\lambda}{2}}\hat{y}_i =\frac{2+h\lambda}{2-h\lambda}\hat{y}_i. \] La función \(Q(\lambda h)\) vale: \[ Q(\lambda h)=\frac{2+h\lambda}{2-h\lambda}. \] Para que los errores estén controlados, el valor \(h\) debe cumplir que \[ |Q(h\lambda)|=\left|\frac{2+h\lambda}{2-h\lambda}\right|<1,\ \Rightarrow |2+h\lambda| < |2-h\lambda|. \]
Sea \(h\lambda = a+b\mathrm{i}\) un número complejo que pertenece a \(R\). Entonces: \[ |2+a+b\mathrm{i}|< |2-(a+b\mathrm{i})|,\ \Rightarrow (a+2)^2+b^2 < (a-2)^2+b^2,\ \Rightarrow 4a< 0,\ \Rightarrow a<0. \] La región \(R\) será: \[ R=\{h\lambda\in\mathbb{C}\ |\ \mathrm{Re}(h\lambda)<0\}. \] El método es A-estable porque la región \(R\) de estabilidad absoluta es la parte izquierda del plano complejo.´