En este tema vamos a introducir la técnica de las diferencias divididas para hallar soluciones aproximadas de ecuaciones en derivadas parciales (EDPs).

Una ecuación en derivadas parciales en el plano \(\mathbb{R}^2\) es una ecuación del tipo \[ F\left(x,y,u(x,y),\frac{\partial u}{\partial x},\frac{\partial u}{\partial y},\ldots,\frac{\partial^{i+j} u}{\partial x^i\partial y^j}\right)=0, \] donde \(u(x,y)\) es una función definida en un dominio \(\mathbf{D}\subseteq\mathbb{R}^2\), \(F\) es una función de diversas variables que depende de las variables \(x,y,u\) y las parciales de la función \(u\) y \(i+j\) sería el llamado orden de la ecuación en derivadas parciales.

De la misma manera, una ecuación en derivadas parciales en el espacio \(\mathbb{R}^3\) es una ecuación del tipo \[ F\left(x,y,z,u(x,y,z),\frac{\partial u}{\partial x},\frac{\partial u}{\partial y},\frac{\partial u}{\partial z},\ldots,\frac{\partial^{i+j+k} u}{\partial x^i\partial y^j\partial z^k}\right)=0, \] donde \(u(x,y,z)\) es una función definida en un dominio \(\mathbf{D}\subseteq\mathbb{R}^3\), \(F\) es una función de diversas variables que depende de las variables \(x,y,z,u\) y las parciales de la función \(u\) y \(i+j+k\) sería el llamado orden de la ecuación en derivadas parciales.

Toda ecuación en derivadas parciales necesita unas restricciones para que tenga solución única.

Dichas restricciones suelen ser condiciones asociadas a la frontera del dominio \(\mathbf{D}\) en el sentido que se tiene algún conocimiento de la función \(u\) o de sus parciales en dicho dominio.

Hallar soluciones de ecuaciones en derivadas parciales es uno de los problemas más complejos en el campo del cálculo diferencial y tiene aplicaciones en una gran diversidad de campos tales como física atmosférica, mecánica, ingeniería incluso en machine learning o en big data.

Además, es uno de los campos más activos de investigación hoy en día.

En este tema, vamos a hallar soluciones aproximadas de las ecuaciones en derivadas parciales más conocidas en la literatura basadas en diferencias finitas ya introducidas en el capítulo de Problemas de frontera para ecuaciones diferenciales ordinarias.

Mas concretamente, nos restringiremos en ecuaciones lineales en derivadas parciales de segundo orden.

La forma más general de una EDP lineal de segundo orden es de la forma siguiente: \[ A\frac{\partial^2 u}{\partial x^2}+B\frac{\partial^2 u}{\partial x\partial y}+C\frac{\partial^2 u}{\partial y^2}+D\frac{\partial u}{\partial x}+E\frac{\partial u}{\partial y}+F u=f(x,y), \] donde \(A,B,C,D,E\) y \(F\) determinan el tipo de EDP lineal de segundo orden.

Dichas EDPs se clasifican de la forma siguiente:

• Si \(B^2-4AC>0\), la EDP se denomina elíptica.
• Si \(B^2-4AC=0\), la EDP se denomina parabólica.
• Si \(B^2-4AC<0\), la EDP se denomina hiperbólica.

La EDP elíptica de la que vamos a hallar una solución aproximada es la llamada ecuación de Poisson: \[ \frac{\partial^2 u}{\partial x^2}(x,y)+\frac{\partial^2 u}{\partial y^2}(x,y)=f(x,y). \] Dicha ecuación sirve para modelar varios problemas físicos dependiendo del tiempo tales como

• la distribución en el equilibrio del calor de una región plana \(\mathbf{R}\). En este caso, tenemos que suponer que \(f(x,y)\equiv 0\), dando lugar a la denominada ecuación de Laplace: \[ \frac{\partial^2 u}{\partial x^2}(x,y)+\frac{\partial^2 u}{\partial y^2}(x,y)=0. \] La función \(u(x,y)\) representa la temperatura en el punto \((x,y)\).

• la energía potencial en un punto del plano donde actúan fuerzas gravitacionales y
• problemas bidimensionales en el equilibrio de fluidos incomprimibles.

Para que la ecuación de Poisson tenga solución única, se tienen que imponer unas restricciones.

Las restricciones se denominan condiciones de frontera de Dirichlet y consisten en la suposición de una función \(g(x,y)\) en todo punto de la frontera \(\mathbf{\partial D}\) del dominio \(\mathbf{D}\). \[ u(x,y)=g(x,y),\ (x,y)\in\mathbf{\partial D} \]

La gráfica siguiente muestra un dominio \(\mathbf{D}\) rectangular.

Para ilustrar cómo se aproxima la ecuación de Poisson vamos a suponer que el dominio \(\mathbf{D}\) es un rectángulo en el plano \(\mathbb{R}^2\) de la forma \([a,b]\times [c,d]\) tal como muestra el gráfico anterior.

El primer paso es seleccionar un mallado en el rectángulo anterior, es decir, elegir dos enteros positivos \(n\) y \(m\) y a partir de éstos, considerar dos pasos \(h\) y \(k\) con respecto a los ejes \(X\) e \(Y\), respectivamente: \[ h=\frac{b-a}{n},\quad k=\frac{d-c}{m}. \]

Los puntos del mallado serán de la forma \((x_i,y_j)\), donde \[ x_i=a+ih,\ i=0,1,\ldots,n,\quad y_j=b+jk,\ j=0,1,\ldots,m. \]

Las rectas verticales \(x=x_i\) y las rectas horizontales \(y=y_j\) son las rectas del mallado y sus intersecciones son precisamente los puntos del mallado.

La aproximación de la función \(u(x,y)\) que vamos hallar será en los puntos del mallado \((x_i,y_j)\), \(i=0,1,\ldots,n\), \(j=0,1,\ldots,m\), \(u(x_i,y_j)\).

Dicha aproximación estará basada en aproximar las parciales de segundo orden con respecto \(x\) e \(y\) por sus respectivas diferencias divididas que vimos en el capítulo de Diferenciación e Integración Numérica del curso Métodos numéricos con Python. Cálculo Numérico: \[ \begin{align*} \frac{\partial^2 u}{\partial x^2}(x_i,y_j)= & \frac{u(x_{i+1},y_j)-2u(x_i,y_j)+u(x_{i-1},y_j)}{h^2}-\frac{h^2}{12}\frac{\partial^4 u}{\partial x^4}(\xi_i,y_j),\\ \frac{\partial^2 u}{\partial y^2}(x_i,y_j)= & \frac{u(x_{i},y_{j+1})-2u(x_i,y_j)+u(x_{i},y_{j-1})}{k^2}-\frac{k^2}{12}\frac{\partial^4 u}{\partial y^4}(x_i,\nu_j), \end{align*} \] donde \(\xi_i\in (x_{i-1},x_{i+1})\) y \(\nu_j\in (y_{j-1},y_{j+1})\).

Los índices \(i\) van desde \(i=1\) hasta \(i=n-1\) y los índices \(j\) van desde \(j=1\) hasta \(j=m-1\).

Para \(i=0\) e \(i=n\) no podemos aplicar las fórmulas anteriores ya que no conocemos \(u(x_{i-1},y_j)=u(x_{-1},y_j)\) en el caso \(i=0\) y \(u(x_{i+1},y_j)=u(x_{n+1},y_j)\) en el caso \(i=n\).

De la misma manera, no podemos aplicar las fórmulas anteriores para \(j=0\) y \(j=m\) ya que no conocemos \(u(x_{i},y_{j-1})=u(x_{i},y_{-1})\) en el caso \(j=0\) y \(u(x_{i},y_{j+1})=u(x_{i},y_{m+1})\) en el caso \(j=m\).

Sin embargo, las consideraciones anteriores no representan ningún problema ya que los valores \(u(x_0,y_j),\ u(x_n,y_j),\ u(x_i,y_0),\ u(x_i,y_m)\) vienen dados por las restricciones: recordemos que \(u(x,y)=g(x,y),\) si \((x,y)\in\partial\mathbf{D}\), (frontera del rectángulo).

Precisamente los puntos donde evaluamos la función \(u\), \((x_0,y_j)=(a,y_j),\ (x_n,y_j)=(b,y_j),\ (x_i,y_0)=(x_i,c),\ (x_i,y_m)=(x_i,d)\) son los valores de la frontera del rectángulo.

Es decir, \[ \begin{align*} u(x_0,y_j)= & u(a,y_j)=g(x_0,y_j)=g(a,y_j),\\ u(x_n,y_j)= & u(b,y_j)=g(x_n,y_j)=g(b,y_j),\\ u(x_i,y_0)= & u(x_i,c)=g(x_i,y_0)=g(x_i,c),\\ u(x_i,y_m)= & u(x_i,d)=g(x_i,y_m)=g(x_i,d), \end{align*} \] para \(i=0,1,\ldots,n\) y \(j=0,1,\ldots,m\).

Usando las aproximaciones anteriores de las parciales de segundo orden en la ecuación en derivadas parciales, obtenemos el siguiente sistema de ecuaciones lineal donde las incógnitas son las aproximaciones de la función \(u\) en el mallado, \(u(x_i,y_j)\), \(i=0,1,\ldots,n\), \(j=0,1,\ldots,m\): \[ \begin{align*} & \frac{u(x_{i+1},y_j)-2u(x_i,y_j)+u(x_{i-1},y_j)}{h^2}+\frac{u(x_{i},y_{j+1})-2u(x_i,y_j)+u(x_{i},y_{j-1})}{k^2}\\ & =f(x_i,y_j), \end{align*} \] para \(i=1,\ldots,n-1\) y \(j=1,\ldots,m-1\).

Para simplificar la notación, llamamos \(\hat{u}_{ij}=u(x_i,y_j)\) al valor de la aproximación de la función \(u\) en el punto del mallado \((x_i,y_j)\) y \(f(x_i,y_j)=f_{ij}\).

Así, el sistema anterior será: \[ \begin{align*} \frac{\hat{u}_{i+1,j}-2\hat{u}_{ij}+\hat{u}_{i-1,j}}{h^2}+\frac{\hat{u}_{i,j+1}-2\hat{u}_{ij}+\hat{u}_{i,j-1}}{k^2} = & f_{ij},\ \Rightarrow \\ \hat{u}_{i+1,j}-2\hat{u}_{ij}+\hat{u}_{i-1,j}+\left(\frac{h}{k}\right)^2(\hat{u}_{i,j+1}-2\hat{u}_{ij}+\hat{u}_{i,j-1}) = & h^2 f_{ij},\ \Rightarrow\\ -2\left(1+\left(\frac{h}{k}\right)^2\right)\hat{u}_{ij}+\hat{u}_{i+1,j}+\hat{u}_{i-1,j}+\left(\frac{h}{k}\right)^2(\hat{u}_{i,j+1}+\hat{u}_{i,j-1}) = & h^2 f_{ij}, \ \Rightarrow\\ 2\left(1+\left(\frac{h}{k}\right)^2\right)\hat{u}_{ij}-(\hat{u}_{i+1,j}+\hat{u}_{i-1,j})-\left(\frac{h}{k}\right)^2(\hat{u}_{i,j+1}+\hat{u}_{i,j-1}) = & -h^2 f_{ij}, \end{align*} \] para \(i=1,\ldots,n-1\) y \(j=1,\ldots,m-1\).

Nos queda un sistema lineal de la forma \(\mathbf{A}\mathbf{\hat{u}}=\mathbf{b}\) a resolver donde \(\mathbf{A}\) es una matriz de dimensiones \((n-1)\cdot (m-1)\times (n-1)\cdot (m-1)\), \(\hat{u}\), el vector de incógnitas, tiene \((n-1)\cdot (m-1)\) componentes y el término independiente \(\mathbf{b}\) también tiene \((n-1)\cdot (m-1)\) componentes.

En el sistema anterior son conocidos los valores siguientes: \[ \begin{align*} \hat{u}_{0j}=& g(x_0,y_j)=g(a,y_j),\quad \hat{u}_{nj}=g(x_n,y_j)=g(b,y_j),\\ \hat{u}_{i0}=& g(x_i,y_0)=g(x_i,c),\quad \hat{u}_{im}=g(x_i,y_m)=g(x_i,d), \end{align*} \] para \(i=0,1,\ldots,n\) y \(j=0,1,\ldots,m\).

La observación anterior significa que cada fila de la matriz del sistema anterior \(\mathbf{A}\) tiene sólo \(5\) componentes no nulas como máximo (fijémonos que para \(i=j=1\) sólo hay \(3\) o \(4\) componentes no nulas ya que siempre hay \(1\) o \(2\) puntos del mallado que son de la frontera y se suponen conocidas. Por tanto, no forman parte de la matriz \(\mathbf{A}\)).

El número de componentes no nulas de la matriz \(\mathbf{A}\) será como máximo \(5\cdot (n-1)\cdot (m-1)\) de entre \((n-1)^2\cdot (m-1)^2\) componentes en total, es decir, \[ \frac{5\cdot (n-1)\cdot (m-1)}{(n-1)^2\cdot (m-1)^2}=\frac{5}{(n-1)\cdot (m-1)}, \] o el porcentaje de ceros vale como mínimo: \[ 100\cdot\left(1-\frac{5}{(n-1)\cdot (m-1)}\right)\%. \]

Por ejemplo, en el ejemplo considerado (\(n=12\) y \(m=6\)), la matriz \(\mathbf{A}\) tiene como mínimo: \(100\cdot\left(1-\frac{5}{11\cdot 5}\right)\% =90.91\%\) de ceros.

Esto significa que la matriz \(\mathbf{A}\) para \(n\) y \(m\) relativamente grandes es una matriz denominada sparse, es decir, una matriz con muchos ceros.

De cara a obtener una matriz del sistema \(\mathbf{A}\) lo más óptima posible en el sentido que los ceros tengan una estructura lo más definida posible, se propone el etiquetado siguiente de los puntos del mallado donde se muestra el etiquetado antiguo de la forma \((i,j)\) y el propuesto de la forma \(P_k\):

Con este etiquetado y usando el ejemplo con \(n=12\) y \(m=6\), la fila \(13\) de la matriz \(\mathbf{A}\) tendría las componentes \(2,12,13,14\) y \(24\) distintas de cero.

La aplicación que nos pasa de los índices del mallado antiguo al mallado propuesto es la siguiente: \[ \begin{align*} \mathrm{Mallado\ antiguo} & \longrightarrow\mathrm{\ Mallado\ propuesto}\\ (i,j) & \longrightarrow k=i+(m-1-j)\cdot (n-1),\\ \end{align*} \] siendo su inversa: \[ \begin{align*} & \mathrm{Mallado\ propuesto} \longrightarrow\mathrm{\ Mallado\ antiguo}\\ & k \longrightarrow\left(i= k-\left(\left\lceil\frac{k}{n-1}\right\rceil-1\right)(n-1),j=m-\left\lceil\frac{k}{n-1}\right\rceil\right). \end{align*} \]

La fila correspondiente al punto del mallado \(P_k\) tienen los elementos \(k-1,k+1,k-(n-1)\) y \(k+n+1\) diferentes de cero siempre que \(k\pm 1,k\pm (n-1)\in \{1,\ldots,(n-1)\cdot (m-1)\}\).

Los valores de dichos elementos son los siguientes: \[ \begin{align*} a_{k,k-(n-1)}= & -\left(\frac{h}{k}\right)^2,\ a_{k,k-1}=-1,\ a_{k,k}=2\left(1+\left(\frac{h}{k}\right)^2\right),\\ a_{k,k+1}= & -1,\ a_{k,k+(n-1)}=-\left(\frac{h}{k}\right)^2,\ k=1,\ldots,(n-1)\cdot (m-1). \end{align*} \] Por tanto, la matriz del sistema lineal a resolver es simétrica.

Para resolver el sistema lineal \(\mathbf{A}\mathbf{\hat{u}}=\mathbf{b}\), debido a la estructura de la matriz \(\mathbf{A}\), podemos usar tanto un método directo o un método iterativo: (ver los capítulos correspondientes a la resolución de sistemas lineales por métodos directos o iterativos del curso Métodos numéricos con Python. Álgebra Lineal Numérica)

En el caso de usar un método iterativo, se puede usar los métodos clásicos de Jacobi, Gauss-Seidel o SOR.

Si se usa el método de Jacobi, se puede demostrar que el radio espectral de la matriz de iteración \(\mathbf{M_J}\) vale: \[ \rho(\mathbf{M_J})=\frac{1}{2}\left(\cos\left(\frac{\pi}{m}\right)+\cos\left(\frac{\pi}{n}\right)\right). \] Como \(\rho(\mathbf{M_J})<1\) ya que: \[ 0<\rho(\mathbf{M_J})\leq \frac{1}{2}(1+1)=1, \] el método de Jacobi siempre es convergente.

Si se usa el método SOR, sabemos (ver el curso anterior) que el valor del parámetro \(w\) óptimo vale: \[ w_{\mathrm{óptimo}}=\frac{2}{1+\sqrt{1-\rho(\mathbf{M_J})^2}}=\frac{4}{2+\sqrt{4-\left(\cos\left(\frac{\pi}{m}\right)+\cos\left(\frac{\pi}{n}\right)\right)^2}}. \]

• INPUT \(m\), \(n\), \(a\), \(b\), \(c\), \(d\), rutina que da \(f(x,y)\) donde \((x,y)\in\mathbf{D}\) es del dominio, rutina que da \(g(x,y)\) donde \((x,y)\in\partial\mathbf{D}\) es de la frontera del dominio.
• Set \(h=\frac{b-a}{n}\). (Hallamos los pasos para \(x\) y para \(y\))
• Set \(k=\frac{d-c}{m}\).
• Set \(\mathbf{x}=\mathbf{0}\) (En los vectores \(\mathbf{x}\) y \(\mathbf{y}\) guardaremos el mallado)
• Set \(\mathbf{y}=\mathbf{0}\)
• For i=1,...,n-1 (Hallamos el mallado)
  • Set \(x_i=a+i\cdot h\).
• For j=1,...,m-1
  • Set \(y_j=c+j\cdot k\).

• Set \(\lambda = \frac{h^2}{k^2}\). (Vamos a definir la matriz del sistema \(\mathbf{A}\) y el término independiente \(\mathbf{b}\))
• Set \(\mu = 2\cdot (1+\lambda)\).
• Set \(\mathbf{A}=\mathbf{0}\).
• Set \(\mathbf{b}=\mathbf{0}\).
• For l=1,...,(n-1)(m-1)
  • Set \(i=l-\left(\lceil\frac{l}{n-1}\rceil-1\right)\cdot (n-1)\).
  • Set \(j=m-\lceil\frac{l}{n-1}\rceil\).
  • Set \(b_l=-h^2\cdot f(x_i,y_j)\)

• • Set \(a_{ll}=\mu\).
  • If \((i-1)\geq 1\)
    • Set \(a_{l,l-1}=-1\).
  • If \((i+1)\leq n-1\)
    • Set \(a_{l,l+1}=-1\).
  • If \((j-1) >= 1\)
    • Set \(a_{l,l+(n-1)}=-\lambda\).
  • If \((j+1) <= m-1\)
    • Set \(a_{l,l-(n-1)}=-\lambda\).

• • If \(i=1\)
    • Set \(b_l=b_l+g(a,y_j)\).
  • If \(i=n-1\)
    • Set \(b_l=b_l+g(b,y_j)\).
  • If \(j=1\)
    • Set \(b_l=b_l+\lambda\cdot g(x_i,c)\).
  • If \(j=m-1\)
    • Set \(b_l=b_l+\lambda\cdot g(x_i,d)\).
• Solve \(\mathbf{A}{\mathbf{\hat{u}}}=\mathbf{b}\). (Resolvemos el sistema correspondiente)
• Print \(\mathbf{\hat{u}}\). (Damos la solución)

Para ilustrar cómo hallar una aproximación de una EDP parabólica, vamos a considerar la ecuación siguiente: \[ \frac{\partial u}{\partial t}(x,t)-\alpha^2\frac{\partial^2 u}{\partial x^2}(x,t)=0, \] que modela el flujo de calor a lo largo de una varilla de longitud \(l\), es decir \(u(x,t)\) representa la temperatura de la varilla en la posición \(x\) en el tiempo \(t\) y \(\alpha\) es un parámetro que depende de la varilla:

Vamos a considerar las restricciones siguientes:

• Caso 1:
  • Distribución del flujo inicial para \(t=0\): \(u(x,0)=f(x)\) donde \(f\) se supone una función conocida.
  • Distribución del flujo en los extremos de la varilla: \(u(0,t)=u_1\) y \(u(l,t)=u_2\), donde \(u_1\) y \(u_2\) son valores conocidos. Además, suponemos que para \(t\to\infty\), la temperatura tiende a la función siguiente: \[ \lim_{t\to\infty} u(x,t)=u_1+\frac{(u_2-u_1)x}{l}, \] es decir la temperatura aumenta de forma lineal desde \(u=u_1\) en la posición \(x=0\) hasta \(u=u_2\) en la posición \(x=l\).

• Caso 2:
  • Distribución del flujo inicial para \(t=0\): \(u(x,0)=f(x)\) donde \(f\) se supone una función conocida.
  • Suponemos que no hay flujo de calor en los extremos. En este caso: \[ \frac{\partial u}{\partial x}(0,t)=0,\quad \frac{\partial u}{\partial x}(l,t)=0. \] En este caso, el calor nunca escapa de la varilla y en el caso límite para \(t\to\infty\), la temperatura de la varilla se mantendría constante.

Las EDPs parabólicas son fundamentales en el estudio de difusión de gases.

Vamos a ilustrar cómo hallar una aproximación de la EDP parabólica que estudia la distribución de la temperatura \(u(x,t)\) en una varilla en el caso 1 suponiendo que la temperatura es cero en los extremos.

Es decir, vamos a hallar una aproximación de la EDP siguiente: \[ \frac{\partial u}{\partial t}(x,t)=\alpha^2\frac{\partial^2 u}{\partial x^2}(x,t),\ 0<x<l,\ t>0, \] con restricciones: \[ u(0,t)=u(l,t)=0,\ t>0,\quad u(x,0)=f(x),\ 0\leq x\leq l. \]

Para seleccionar el paso espacial o para la variable \(x\), elegimos un entero \(m>0\) y elegimos \(h=\frac{l}{m}\).

A continuación, seleccionamos un paso temporal \(k\).

El mallado consiste en los puntos \((x_i,t_j)\), donde \(x_i=ih\), \(i=0,1,\ldots,m\) y \(t_j=jk\), \(j=0,1,\ldots\)

El objetivo es hallar una aproximación \(\hat{u}(x_i,t_j)\), \(i=0,1,\ldots,m\), \(j=0,1,\ldots\) en los puntos del mallado.

Observar que los valores \(\hat{u}(x_i,t_j)\) para \(i=0,l\) y para \(j=0\) son conocidos a partir de las restricciones: \[ \hat{u}(x_0,t_j)=\hat{u}(0,t_j)=0,\quad \hat{u}(x_m,t_j)=\hat{u}(l,t_j)=0,\quad \hat{u}(x_i,0)=f(x_i), \] para \(i=0,1,\ldots,m,\ j=0,1,\ldots\)

Para hallar la aproximación \(\hat{u}(x_i,t_j)\) en el mallado consideramos las aproximaciones siguientes vistas en el capítulo de Derivación e Integración Numérica del curso Métodos numéricos con Python. Cálculo Numérico: \[ \begin{align*} \frac{\partial u}{\partial t}(x_i,t_j) = & \frac{u(x_i,t_j+k)-u(x_i,t_j)}{k}-\frac{k}{2}\frac{\partial^2 u}{\partial t^2}(x_i,\mu_j),\\ \frac{\partial^2 u}{\partial x^2}(x_i,t_j) = & \frac{u(x_i+h,t_j)-2u(x_i,t_j)+u(x_i-h,t_j)}{h^2}-\frac{h^2}{12}\frac{\partial^4 u}{\partial t^4}(\xi_i,t_j), \end{align*} \] donde \(\mu_j\in (t_j,t_{j+1})\) y \(\xi_i\in (x_{i-1},x_{i+1})\).

Llamando \(\hat{u}_{ij}=\hat{u}(x_i,t_j)\), \(i=0,1,\ldots,m\) y \(j=0,1,\ldots\) y sustituyendo las aproximaciones anteriores en la EDP del flujo de calor obtenemos: \[ \begin{align*} & \frac{\hat{u}_{i,j+1}-\hat{u}_{ij}}{k}-\alpha^2\frac{\hat{u}_{i+1,j}-2\hat{u}_{ij}+\hat{u}_{i-1,j}}{h^2}= 0,\ \Rightarrow\\ & \hat{u}_{i,j+1}-\hat{u}_{ij}-\frac{k\alpha^2}{h^2}(\hat{u}_{i+1,j}-2\hat{u}_{ij}+\hat{u}_{i-1,j}) = 0,\ \Rightarrow\\ & \hat{u}_{i,j+1} = \left(1-\frac{2k\alpha^2}{h^2}\right)\hat{u}_{ij}+\frac{k\alpha^2}{h^2}(\hat{u}_{i+1,j}+\hat{u}_{i-1,j}),\\ & \hat{u}_{i,j+1} = \left(1-2\lambda\right)\hat{u}_{ij}+\lambda (\hat{u}_{i+1,j}+\hat{u}_{i-1,j}), \end{align*} \] donde \(\lambda = \frac{k\alpha^2}{h^2}\), para \(i=1,2,\ldots,m-1\) y \(j=1,2,\ldots\)

La ecuación anterior nos da el vector de aproximaciones en el tiempo \(t_{j+1}\), \((\hat{u}_{1,j+1},\ldots,\hat{u}_{m-1,j+1})\) en función del vector de aproximaciones en el tiempo \(t_j\), \((\hat{u}_{1,j},\ldots,\hat{u}_{m-1,j})\).

Como conocemos el vector de aproximaciones en el tiempo inicial \(t=0\): \[ (\hat{u}_{1,0},\ldots,\hat{u}_{m-1,0})=(f(x_1),\ldots,f(x_{m-1})), \] basta aplicar la recurrencia anterior para hallar el vector de aproximaciones en cualquier valor \(t_j\) del mallado.

Llamando \(\hat{u}^{(j)}\) al vector de aproximaciones en el tiempo \(t_j\), \(\hat{u}^{(j)}=(\hat{u}_{1,j},\ldots,\hat{u}_{m-1,j})^\top\), la recurrencia hallada se puede escribir matricialmente de la forma siguiente: \[ \hat{u}^{(j+1)}=\mathbf{A}\hat{u}^{(j)},\ \hat{u}^{(0)}=(f(x_1),\ldots,f(x_{m-1}))^\top, \] donde \(\mathbf{A}\) es la siguiente matriz tridiagonal:

\[ \mathbf{A}=\begin{bmatrix}(1-2\lambda)&\lambda&0&0&\ldots&0 \\\lambda&(1-2\lambda)&\lambda&0&\ldots&0 \\0&\lambda&(1-2\lambda)&\lambda&\ldots&0 \\\vdots&\vdots&\vdots&\ddots&\vdots&\vdots \\0&\ldots&0&\lambda&(1-2\lambda)&\lambda \\0&\ldots&\ldots&0&\lambda&(1-2\lambda) \\\end{bmatrix}, \]

Pensemos que las restricciones imponen que \(\hat{u}_{0j}=\hat{u}_{mj}=0\).

Así, por ejemplo: \[ \begin{align*} \hat{u}^{(1)}= & \begin{bmatrix}\hat{u}_{11}\\\hat{u}_{21}\\\vdots\\\hat{u}_{m-1,1}\end{bmatrix}=\mathbf{A}\hat{u}^{(0)}=\begin{bmatrix}(1-2\lambda)&\lambda&0&0&\ldots&0 \\\lambda&(1-2\lambda)&\lambda&0&\ldots&0 \\0&\lambda&(1-2\lambda)&\lambda&\ldots&0 \\\vdots&\vdots&\vdots&\ddots&\vdots&\vdots \\0&\ldots&0&\lambda&(1-2\lambda)&\lambda \\0&\ldots&\ldots&0&\lambda&(1-2\lambda) \\\end{bmatrix}\begin{bmatrix}\hat{u}_{10}\\\hat{u}_{20}\\\vdots\\\hat{u}_{m-1,0}\end{bmatrix},\\ = & \begin{bmatrix}(1-2\lambda)&\lambda&0&0&\ldots&0 \\\lambda&(1-2\lambda)&\lambda&0&\ldots&0 \\0&\lambda&(1-2\lambda)&\lambda&\ldots&0 \\\vdots&\vdots&\vdots&\ddots&\vdots&\vdots \\0&\ldots&0&\lambda&(1-2\lambda)&\lambda \\0&\ldots&\ldots&0&\lambda&(1-2\lambda) \\\end{bmatrix}\begin{bmatrix}f(x_1)\\ f(x_2)\\\vdots\\ f(x_{m-1})\end{bmatrix} \end{align*} \]

En componentes, \[ \begin{align*} \hat{u}_{11} = & (1-2\lambda) f(x_1)+\lambda f(x_2),\\ \hat{u}_{21} = & \lambda f(x_1)+(1-2\lambda) f(x_2)+\lambda f(x_3),\\ &\vdots \\ \hat{u}_{m-2,1} = & \lambda f(x_{m-3})+(1-2\lambda) f(x_{m-2})+\lambda f(x_{m-1}),\\ \hat{u}_{m-1,1} = & \lambda f(x_{m-2})+(1-2\lambda) f(x_{m-1}), \end{align*} \] y así sucesivamente.

El método descrito anteriormente se denomina método de diferencias divididas hacia adelante debido a la aproximación de \(\frac{\partial u}{\partial t}\) usada.

• INPUT \(m\), \(k\), \(l\), rutina que da \(f(x)\) donde \(x\in [0,l]\), \(\alpha\), \(N\) (\(N\) es el número de iteraciones temporales a realizar)
• Set \(h=\frac{l}{m}\). (Hallamos el paso espacial)
• Set \(\mathbf{x}=\mathbf{0}\) (En el vector \(\mathbf{x}\) guardaremos el mallado espacial)
• For i=1,...,m-1 (Hallamos el mallado espacial)
  • Set \(x_i=i\cdot h\).
• Set \(\lambda = \frac{k\alpha^2}{h^2}\).
• Set \(\mathbf{\hat{u}}^{(0)}=\mathbf{0}\) (Vamos a definir el vector \(\hat{u}^{(0)}\))
• For i=1,...,m-1
  • Set \(\hat{u}_i^{(0)} =f(x_i)\)

• Set \(\mathbf{\hat{u}}=\mathbf{0}\)
• For j=1,...,N (Vamos calculando los vectores de aproximación \(\hat{u}^{(j)}\))
  • For i=2,...,m-2
    • Set \(\hat{u}_i =\lambda {\hat{u}}^{(0)}_{i-1}+(1-2\lambda){\hat{u}}^{(0)}_i+\lambda\hat{u}^{(0)}_{i+1}\).
  • Set \(\hat{u}_1 = (1-2\lambda)\hat{u}_1^{(0)}+\lambda \hat{u}_2^{(0)}\).
  • Set \(\hat{u}_{m-1}=\lambda\hat{u}_{m-2}^{(0)}+(1-2\lambda)\hat{u}_{m-1}^{(0)}\).
  • Set \(\mathbf{\hat{u}}^{(0)}=\mathbf{\hat{u}}\).
• Print \(\mathbf{\hat{u}}\).

El método anterior de resolución de la EDP parabólica consiste en aplicar la recurrencia \(\mathbf{\hat{u}}^{(j+1)}=\mathbf{A}\mathbf{\hat{u}}^{(j)}\) a partir de un valor inicial \(\mathbf{\hat{u}}^{(0)}\).

En la práctica, dicho valor inicial \(\mathbf{\hat{u}}^{(0)}\) tendrá errores de redondeo.

Nos preguntamos cómo se van propagando dichos errores de redondeo a los diferentes valores \(\mathbf{\hat{u}}^{(j)}\), es decir, queremos estudiar si el método de aproximación anterior es estable. Ver el capítulo correspondiente a Errores del curso Métodos numéricos con Python. Cálculo Numérico para más detalles sobre estabilidad de métodos numéricos.

Más concretamente, el valor inicial que tenemos en la práctica es \(\mathbf{\hat{u}}^{(0)}+\mathbf{e}^{(0)}\), donde \(\mathbf{e}^{(0)}\) son los errores de redondeo de dicho valor inicial. Entonces el valor \(\mathbf{\hat{u}}^{(1)}\) que obtendremos será: \[ \mathbf{\hat{u}}^{(1)}=\mathbf{A}(\mathbf{\hat{u}}^{(0)}+\mathbf{e}^{(0)})=\mathbf{A}\mathbf{\hat{u}}^{(0)}+\mathbf{A}\mathbf{e}^{(0)}. \]

Entonces el error de redondeo correspondiente al vector de aproximación en el tiempo \(t_1\), \(\mathbf{\hat{u}}^{(1)}\) será: \(\mathbf{e}^{(1)}=\mathbf{A}\mathbf{e}^{(0)}.\)

En general, el error de redondeo correspondiente al vector de aproximación en el tiempo \(t_j\), \(\mathbf{\hat{u}}^{(j)}\) será: \(\mathbf{e}^{(j)}=\mathbf{A}^j\mathbf{e}^{(0)}.\)

Para ver si dicho error es grande, consideremos una norma vectorial (ver el capítulo correspondiente a Preliminares del curso Métodos numéricos con Python. Álgebra Lineal Numérica) y acotemos su valor de la forma siguiente: \[ \|\mathbf{e}^{(j)}\|=\|\mathbf{A}^j\mathbf{e}^{(0)}\|\leq \|\mathbf{A}\|^j\|\cdot \|\mathbf{e}^{(0)}||, \] donde \(\|\mathbf{A}\|\) es la norma matricial asociada a la norma vectorial usada.

Entonces la estabilidad del método viene dada por \(\|\mathbf{A}\|\) en el sentido de que tenemos que tener que \(\|\mathbf{A}\|\leq 1\).

La condición anterior es equivalente a afirmar que \(\rho(\mathbf{A}^j)=\rho(\mathbf{A})^j\leq 1\), donde \(\rho(\mathbf{A})\) representa el radio espectral de la matriz \(\mathbf{A}\). En el capítulo de Preliminares del curso anterior hay todos los detalles de dicho resultado.

En resumen, el método de diferencias divididas hacia adelante será estable si \(\rho(\mathbf{A})\leq 1\).

Como \(\rho(\mathbf{A})=\max_{\mu_i\mathrm{\ valor\ propio\ de\ }\mathbf{A}}|\mu_i|\), necesitamos calcular los valores propios de la matriz \(\mathbf{A}\).

Proposición. Valores y vectores propios de la matriz \(\mathbf{A}\).

Los valores propios de la matriz de iteración \(\mathbf{A}\) de dimensiones \((m-1)\times (m-1)\) correspondiente al método de diferencias divididas hacia adelante son: \[ \mu_i = 1-4\lambda\left(\sin\left(\frac{\mathrm{i}\pi}{2m}\right)\right)^2,\ i=1,2,\ldots,m-1. \] donde el correspondiente vector propio \(\mathbf{v}^{(i)}\) vale: \[ \mathbf{v}^{(i)}=\left(\sin\left(\frac{\mathrm{i}\pi}{m}\right),\sin\left(\frac{2\mathrm{i}\pi}{m}\right),\ldots,\sin\left(\frac{(m-1)\mathrm{i}\pi}{m}\right)\right), \] es decir, la componente \(k\)-ésima del vector propio \(\mathbf{v}^{(i)}\) vale \(v_k^{(i)}=\sin\left(\frac{k\mathrm{i}\pi}{m}\right)\).

Usando la proposición anterior, podemos decir que el método de las diferencias finitas hacia adelante es estable si para todo \(i\in\{1,\ldots,m-1\}\), \[ \begin{align*} & |\mu_i|\leq 1,\ \Leftrightarrow \left| 1-4\lambda\left(\sin\left(\frac{\mathrm{i}\pi}{2m}\right)\right)^2\right|\leq 1,\ \Leftrightarrow -1\leq 1-4\lambda\left(\sin\left(\frac{\mathrm{i}\pi}{2m}\right)\right)^2\leq 1,\\ & -2\leq -4\lambda\left(\sin\left(\frac{\mathrm{i}\pi}{2m}\right)\right)^2\leq 0,\ \Leftrightarrow 0\leq 4\lambda\left(\sin\left(\frac{\mathrm{i}\pi}{2m}\right)\right)^2\leq 2,\\ & 0\leq \lambda\left(\sin\left(\frac{\mathrm{i}\pi}{2m}\right)\right)^2\leq \frac{1}{2}. \end{align*} \] La condición \(\lambda\left(\sin\left(\frac{\mathrm{i}\pi}{2m}\right)\right)^2\geq 0\) siempre es cierta por definición de \(\lambda\), recordemos que \(\lambda = \frac{k\alpha^2}{h^2}>0\).

Respecto a la otra condición, podemos observar que para todo \(i\in\{1,\ldots,m-1\}\), \[ \left(\sin\left(\frac{\mathrm{i}\pi}{2m}\right)\right)^2 \leq \left(\sin\left(\frac{\mathrm{i(m-1)}\pi}{2m}\right)\right)^2\to 1, \mbox{ si }m\to\infty. \] En resumen, si \(0\leq \lambda\leq \frac{1}{2}\), se cumplirá para todo \(i\in\{1,\ldots,m-1\}\) la condición anterior.

Entonces, el método de las diferencias divididas hacia adelante es estable si: \[ 0\leq \lambda\leq \frac{1}{2},\ \Leftrightarrow 0\leq \frac{k\alpha^2}{h^2}\leq \frac{1}{2},\ \Leftrightarrow 0\leq k\leq \frac{h^2}{2\alpha^2}, \] condición que nos dice que el paso temporal debe ser mucho más pequeño que el paso espacial ya que el error que teníamos en la aproximación de \(\frac{\partial u}{\partial t}\) era lineal mientras que el error en la aproximación de \(\frac{\partial^2 u}{\partial x^2}\) era cuadrático.

En este caso se dice que el método de las diferencias divididas hacia adelante es condicionalmente estable.

Vamos a obtener un método incondicionalmente estable. Para ello, en lugar de considerar diferencias divididas hacia adelante en la aproximación de \(\frac{\partial u}{\partial t}\), consideramos diferencias divididas hacia atrás: \[ \frac{\partial u}{\partial t}(x_i,t_j)=\frac{u(x_i,t_j)-u(x_i,t_{j-1})}{k}+\frac{k}{2}\frac{\partial^2 u}{\partial t^2}(x_i,\mu_j), \] donde \(\mu_j\in (t_{j-1},t_j)\).

Como hacíamos en el método anterior, llamamos \(\hat{u}_{ij}=\hat{u}(u_i,t_j)\), \(i=0,1,\ldots,m\) y \(j=0,1,\ldots\) y sustituyendo la aproximación anterior junto con la aproximación de \(\frac{\partial^u}{\partial x^2}(x_i,t_j)\) vista en el método anterior en la EDP del flujo de calor, llamando \(\lambda = \frac{k\alpha^2}{h^2}\), obtenemos: \[ \begin{align*} & \frac{\hat{u}_{i,j}-\hat{u}_{i,j-1}}{k}-\alpha^2\frac{\hat{u}_{i+1,j}-2\hat{u}_{ij}+\hat{u}_{i-1,j}}{h^2}= 0,\ \Rightarrow\\ & \hat{u}_{i,j}-\hat{u}_{i,j-1}-\frac{k\alpha^2}{h^2}(\hat{u}_{i+1,j}-2\hat{u}_{ij}+\hat{u}_{i-1,j}) = 0,\ \Rightarrow\\ & (1+2\lambda)\hat{u}_{i,j} -\lambda \hat{u}_{i+1,j}-\lambda \hat{u}_{i-1,j} =\hat{u}_{i,j-1}, \end{align*} \] para \(i=1,2,\ldots,m-1\) y \(j=1,2,\ldots\)

Las ecuaciones anteriores dan un sistema de ecuaciones lineal para el vector de aproximaciones \(\mathbf{\hat{u}}^{(j)}\), supuesto conocido el vector de aproximaciones en el tiempo anterior del mallado \(\mathbf{\hat{u}}^{(j-1)}\).

El sistema lineal obtenido es de la forma \(\mathbf{A}\mathbf{\hat{u}}^{(j)}=\mathbf{\hat{u}}^{(j-1)}\), donde \(\mathbf{A}\) es la matriz del sistema tridiagonal simétrica siguiente:

\[ \mathbf{A}=\begin{bmatrix}(1+2\lambda)&-\lambda&0&0&\ldots&0 \\-\lambda&(1+2\lambda)&-\lambda&0&\ldots&0 \\0&-\lambda&(1+2\lambda)&-\lambda&\ldots&0 \\\vdots&\vdots&\vdots&\ddots&\vdots&\vdots \\0&\ldots&0&-\lambda&(1+2\lambda)&-\lambda \\0&\ldots&\ldots&0&-\lambda&(1+2\lambda) \\\end{bmatrix}. \] Para resolver el sistema \(\mathbf{A}\mathbf{\hat{u}}^{(j)}=\mathbf{\hat{u}}^{(j-1)}\), podemos usar el método directo de Crout o un método iterativo vistos en los capítulos de métodos directos e iterativos del curso Métodos numéricos con Python. Álgebra Lineal Numérica.

Para ir obteniendo las aproximaciones \(\mathbf{\hat{u}}^{(j)}\), seguimos los pasos siguientes:

• Paso 1: Obtenemos \(\mathbf{\hat{u}}^{(0)}=(f(x_1),f(x_2),\ldots,f(x_{m-1}))^\top\).
• Paso 2: Obtenemos \(\mathbf{\hat{u}}^{(1)}\) resolviendo el sistema lineal \(\mathbf{A}\mathbf{\hat{u}}^{(1)}=\mathbf{\hat{u}}^{(0)}\).
• Paso 3: Obtenemos \(\mathbf{\hat{u}}^{(2)}\) resolviendo el sistema lineal \(\mathbf{A}\mathbf{\hat{u}}^{(2)}=\mathbf{\hat{u}}^{(1)}\).
• Y así sucesivamente hasta llegar a \(\mathbf{\hat{u}}^{(j)}\) o la aproximación \(\hat{u}(x_i,t_j)\), \(i=1,\ldots,m-1\).

El método anterior se denomina método de diferencias divididas hacia atrás.

Observación.

Si las condiciones frontera en los extremos de la varilla no fuesen \(u(0,t)=u(l,t)=0\), \(t>0\) sino otras, hay que cambiar el término independiente del sistema anterior de la forma siguiente: \(\mathbf{A}\mathbf{\hat{u}}^{(j)}=\mathbf{b}\), siendo \(\mathbf{b}\): \[ \mathbf{b}=\mathbf{\hat{u}}^{(j-1)}+(\lambda u_{0j},0,\ldots,0,\lambda u_{mj})^\top, \] ya que la primera ecuación del sistema \(\mathbf{A}\mathbf{\hat{u}}^{(j)}=\mathbf{b}\) se arregla de la forma siguiente: \[ (1+2\lambda) \hat{u}_{1j}-\lambda \hat{u}_{2j}-\lambda \hat{u}_{0j}=\hat{u}_{1,j-1},\ \Rightarrow (1+2\lambda) \hat{u}_{1j}-\lambda \hat{u}_{2j}=\hat{u}_{1,j-1}+\lambda \hat{u}_{0j}, \] y la última, de la forma siguiente: \[ \begin{align*} (1+2\lambda) \hat{u}_{m-1,j}-\lambda \hat{u}_{m,j}-\lambda \hat{u}_{m-2,j}= & \hat{u}_{m-1,j-1},\ \Rightarrow \\ (1+2\lambda) \hat{u}_{m-1,j}-\lambda \hat{u}_{m-2,j}= & \hat{u}_{m-1,j-1}+\lambda \hat{u}_{m,j}. \end{align*} \]

• INPUT \(m\), \(k\), \(l\), rutina que da \(f(x)\) donde \(x\in [0,l]\), \(\alpha\), \(N\) (\(N\) es el número de iteraciones temporales a realizar)
• Set \(h=\frac{l}{m}\). (Hallamos el paso espacial)
• Set \(\mathbf{x}=\mathbf{0}\) (En el vector \(\mathbf{x}\) guardaremos el mallado espacial)
• For i=1,...,m-1 (Hallamos el mallado espacial)
  • Set \(x_i=i\cdot h\).
• Set \(\lambda = \frac{k\alpha^2}{h^2}\).
• Set \(\mathbf{\hat{u}}^{(0)}=\mathbf{0}\) (Vamos a definir el vector \(\hat{u}^{(0)}\))
• For i=1,...,m-1
  • Set \(\hat{u}_i^{(0)} =f(x_i)\)

• Set \(\mathbf{\hat{u}}=\mathbf{0}\)
• Set \(\mathbf{A}=\mathbf{0}\) (Definimos la matriz del sistema \(\mathbf{A}\))
• For l=1,...,m-1
  • Set \(a_{ll}=1+2\lambda\).
  • If \((l-1)\geq 1\)
    • Set \(a_{l,l-1}=-\lambda\).
  • If \((l+1)\leq m-1\)
    • Set \(a_{l,l+1}=-\lambda\).

• For j=1,...,N (Vamos calculando los vectores de aproximación \(\hat{u}^{(j)}\))
  • Solve \(\mathbf{A}\mathbf{\hat{u}}=\mathbf{\hat{u}}^{(0)}\).
  • Set \(\mathbf{\hat{u}}^{(0)}=\mathbf{\hat{u}}\).
• Print \(\mathbf{\hat{u}}\).

El método de las diferencias divididas hacia atrás puede considerarse desde el punto de vista teórico como un método recurrente tal como ocurría con el método de las diferencias divididas hacia adelante que calcula el vector de aproximación en el tiempo \(t_j\), \(\mathbf{\hat{u}}^{(j)}\) en función del vector de aproximación en el tiempo \(t_{j-1}\), \(\mathbf{\hat{u}}^{(j-1)}\): \[ \mathbf{\hat{u}}^{(j)}=\mathbf{A}\mathbf{\hat{u}}^{(j-1)},\ \Rightarrow \mathbf{\hat{u}}^{(j)}=\mathbf{A}^{-1}\mathbf{\hat{u}}^{(j-1)}, \] donde la matriz de recurrencia o de iteración vale en este caso \(\mathbf{A}^{-1}\).

En general, podemos escribir, \[ \mathbf{\hat{u}}^{(j)}=\left(\mathbf{A}^{-1}\right)^j\mathbf{\hat{u}}^{(0)}. \] La estabilidad del método de las diferencias divididas hacia atrás vendrá dado por el radio espectral de la matriz \(\mathbf{A}^{-1}\) en el sentido de que si \(\rho(\mathbf{A}^{-1})\leq 1\), método será estable y si \(\rho(\mathbf{A}^{-1})> 1\), el método será inestable.

El radio espectral de la matriz \(\mathbf{A}^{-1}\) valdrá: \[ \rho(\mathbf{A}^{-1})=\max_{\mu_i \mathrm{\ valor\ propio\ de\ }\mathbf{A}^{-1}}|\mu_i|. \] Ahora bien, \(\mu_i\) es un valor propio de \(\mathbf{A}^{-1}\) si, y sólo si, \(\frac{1}{\mu_i}\) es un valor propio de \(\mathbf{A}\) (ver el capítulo de Preliminares del curso Métodos numéricos con Python. Álgebra Lineal Numérica). Entonces, \[ \rho(\mathbf{A}^{-1})=\max_{\mu_i \mathrm{\ valor\ propio\ de\ }\mathbf{A}^{-1}}|\mu_i|=\max_{\beta_i \mathrm{\ valor\ propio\ de\ }\mathbf{A}}\left|\frac{1}{\beta_i}\right|=\frac{1}{\min_{\beta_i \mathrm{\ valor\ propio\ de\ }\mathbf{A}}|\beta_i|}. \]

Los valores propios de la matriz \(\mathbf{A}\) se pueden obtener a partir de la proposición que nos daba los valores propios de la matriz de recurrencia o de iteración correspondiente al método de las diferencias divididas hacia adelante ya que la matriz del sistema \(\mathbf{A}\) correspondiente al método de las diferencias divididas hacia atrás se obtiene a partir de la matriz de recurrencia o de iteración correspondiente al método de las diferencias divididas hacia adelante cambiando \(\lambda\) por \(-\lambda\).

Entonces, los valores propios correspondientes a la matriz del sistema del método de las diferencias divididas hacia atrás serán: \[ \beta_i = 1+4\lambda\left(\sin\left(\frac{\mathrm{i}\pi}{2m}\right)\right)^2,\ i=1,2,\ldots,m-1. \] Como \(\lambda>0\), tendremos que \(\beta_i > 1\) para \(i=1,2,\ldots,m-1\) y por tanto, \(\min_{i}\beta_i >1\).

En conclusión, \[ 0<\rho(\mathbf{A}^{-1}) =\frac{1}{\min_i \beta_i}<1, \] y el método de las diferencias divididas hacia atrás es incondicionalmente estable lo que significa que somos libres de elegir el paso temporal \(k\) y el paso espacial \(h\) que el método siempre será estable.

Aunque el método sea incondicionalmente estable, no debemos olvidar que el error en la aproximación de \(\frac{\partial u}{\partial t}\) es lineal mientras que el error en la aproximación de \(\frac{\partial^2 u}{\partial x^2}\) es cuadrático, lo que significa que el paso temporal debe ser siempre mucho menor que el paso espacial si queremos tener aproximaciones fiables.

Vamos a ver un método incondicionalmente estable que además tiene error cuadrático en las aproximaciones de \(\frac{\partial u}{\partial t}\) y \(\frac{\partial^2 u}{\partial x^2}\).

La idea del método es considerar los dos métodos vistos hasta ahora (diferencias divididas hacia adelante y hacia atrás) en el instante del tiempo \(t_{j+1}\) (¡ojo! que el hacia atrás está desarrollado en \(t=t_j\)) \[ \begin{align*} \frac{\hat{u}_{i,j+1}-\hat{u}_{i,j}}{k}-\alpha^2\frac{\hat{u}_{i+1,j}-2\hat{u}_{i,j}+\hat{u}_{i-1,j}}{h^2}=& 0,\quad (\mbox{hacia adelante en } t_{j+1}),\\ \frac{\hat{u}_{i,j+1}-\hat{u}_{i,j}}{k}-\alpha^2\frac{\hat{u}_{i+1,j+1}-2\hat{u}_{i,j+1}+\hat{u}_{i-1,j+1}}{h^2}=& 0,\quad (\mbox{hacia atrás en } t_{j+1}),\\ \end{align*} \] con errores de truncamiento \(\frac{k}{2}\frac{\partial^u}{\partial t^2}(x_i,\mu_j)+O(h^2)\) (hacia adelante) y \(-\frac{k}{2}\frac{\partial^u}{\partial t^2}(x_i,\tilde{\mu}_j)+O(h^2)\) (hacia atrás) donde \(\mu_j\in (t_j,t_j+k)\) y \(\tilde{\mu}_j\in (t_j-k,t_j)\) y considerar la suma de las expresiones partido por dos:

\[ \begin{align*} & \frac{1}{2}\left(\frac{\hat{u}_{i,j+1}-\hat{u}_{i,j}}{k}+\frac{\hat{u}_{i,j+1}-\hat{u}_{i,j}}{k}\right)\\ & -\frac{\alpha^2}{2}\left(\frac{\hat{u}_{i+1,j}-2\hat{u}_{i,j}+\hat{u}_{i-1,j}}{h^2}+\frac{\hat{u}_{i+1,j+1}-2\hat{u}_{i,j+1}+\hat{u}_{i-1,j+1}}{h^2}\right)=0,\\ & \frac{\hat{u}_{i,j+1}-\hat{u}_{i,j}}{k}-\frac{\alpha^2}{2}\left(\frac{\hat{u}_{i+1,j+1}-2\hat{u}_{i,j+1}+\hat{u}_{i-1,j+1}+\hat{u}_{i+1,j}-2\hat{u}_{i,j}+\hat{u}_{i-1,j}}{h^2}\right)=0,\\ &\hat{u}_{i,j+1}-\hat{u}_{i,j}-\frac{\lambda}{2}\left(\hat{u}_{i+1,j+1}-2\hat{u}_{i,j+1}+\hat{u}_{i-1,j+1}+\hat{u}_{i+1,j}-2\hat{u}_{i,j}+\hat{u}_{i-1,j}\right)=0,\\ & -\frac{\lambda}{2}\hat{u}_{i+1,j+1}+(1+\lambda)\hat{u}_{i,j+1}-\frac{\lambda}{2}\hat{u}_{i-1,j+1}=\frac{\lambda}{2}\hat{u}_{i+1,j}+(1-\lambda)\hat{u}_{i,j}+\frac{\lambda}{2}\hat{u}_{i-1,j}, \end{align*} \] para \(i=1,\ldots,m-1\) y \(j=0,1,\ldots.\)

Si suponemos que \(\frac{\partial^2 u}{\partial t^2}(x_i,\mu_j)\approx \frac{\partial^2 u}{\partial t^2}(x_i,\tilde{\mu}_j)\), el error cometido en la aproximación anterior será cuadrático, \(O(k^2)\) ya que dicho error vale: \[ \frac{1}{2}\left(\frac{k}{2}\frac{\partial^u}{\partial t^2}(x_i,\mu_j)-\frac{k}{2}\frac{\partial^u}{\partial t^2}(x_i,\tilde{\mu}_j)\right)+O(h^2)\approx O(k^2)+O(h^2). \]

Matricialmente, els sistema de ecuaciones anterior se puede escribir de la forma siguiente: \[ \mathbf{A}\mathbf{\hat{u}}^{(j+1)}=\mathbf{B}\mathbf{\hat{u}}^{(j)}, \] donde las matrices \(\mathbf{A}\) y \(\mathbf{B}\) son las siguientes:

\[ \begin{align*} \mathbf{A}= & \begin{bmatrix}(1+\lambda)&-\frac{\lambda}{2}&0&0&\ldots&0 \\-\frac{\lambda}{2}&(1+\lambda)&-\frac{\lambda}{2}&0&\ldots&0 \\0&-\frac{\lambda}{2}&(1+\lambda)&-\frac{\lambda}{2}&\ldots&0 \\\vdots&\vdots&\vdots&\ddots&\vdots&\vdots \\0&\ldots&0&-\frac{\lambda}{2}&(1+\lambda)&-\frac{\lambda}{2} \\0&\ldots&\ldots&0&-\frac{\lambda}{2}&(1+\lambda) \\\end{bmatrix},\\ \mathbf{B}= & \begin{bmatrix}(1-\lambda)&\frac{\lambda}{2}&0&0&\ldots&0 \\\frac{\lambda}{2}&(1-\lambda)&\frac{\lambda}{2}&0&\ldots&0 \\0&\frac{\lambda}{2}&(1-\lambda)&\frac{\lambda}{2}&\ldots&0 \\\vdots&\vdots&\vdots&\ddots&\vdots&\vdots \\0&\ldots&0&\frac{\lambda}{2}&(1-\lambda)&\frac{\lambda}{2} \\0&\ldots&\ldots&0&\frac{\lambda}{2}&(1-\lambda) \\\end{bmatrix}. \end{align*} \]

Para ir obteniendo las aproximaciones \(\mathbf{\hat{u}}^{(j)}\), seguimos los pasos siguientes:

• Paso 1: Obtenemos \(\mathbf{\hat{u}}^{(0)}=(f(x_1),f(x_2),\ldots,f(x_{m-1}))^\top\).
• Paso 2: Obtenemos \(\mathbf{\hat{u}}^{(1)}\) resolviendo el sistema lineal \(\mathbf{A}\mathbf{\hat{u}}^{(1)}=\mathbf{B}\mathbf{\hat{u}}^{(0)}\).
• Paso 3: Obtenemos \(\mathbf{\hat{u}}^{(2)}\) resolviendo el sistema lineal \(\mathbf{A}\mathbf{\hat{u}}^{(2)}=\mathbf{B}\mathbf{\hat{u}}^{(1)}\).
• Y así sucesivamente hasta llegar a \(\mathbf{\hat{u}}^{(j)}\) o la aproximación \(\hat{u}(x_i,t_j)\), \(i=1,\ldots,m-1\).

El método anterior se denomina método de Crank-Nicolson.

• INPUT \(m\), \(k\), \(l\), rutina que da \(f(x)\) donde \(x\in [0,l]\), \(\alpha\), \(N\) (\(N\) es el número de iteraciones temporales a realizar)
• Set \(h=\frac{l}{m}\). (Hallamos el paso espacial)
• Set \(\mathbf{x}=\mathbf{0}\) (En el vector \(\mathbf{x}\) guardaremos el mallado espacial)
• For i=1,...,m-1 (Hallamos el mallado espacial)
  • Set \(x_i=i\cdot h\).
• Set \(\lambda = \frac{k\alpha^2}{h^2}\).
• Set \(\mathbf{\hat{u}}^{(0)}=\mathbf{0}\) (Vamos a definir el vector \(\hat{u}^{(0)}\))
• For i=1,...,m-1
  • Set \(\hat{u}_i^{(0)} =f(x_i)\)

• Set \(\mathbf{\hat{u}}=\mathbf{0}\)
• Set \(\mathbf{A}=\mathbf{0}\) (Definimos la matriz del sistema \(\mathbf{A}\))
• For l=1,...,m-1
  • Set \(a_{ll}=1+\lambda\).
  • If \((l-1)\geq 1\)
    • Set \(a_{l,l-1}=-\frac{\lambda}{2}\).
  • If \((l+1)\leq m-1\)
    • Set \(a_{l,l+1}=-\frac{\lambda}{2}\).

• Set \(\mathbf{b}=\mathbf{0}\) (En el vector \(\mathbf{b}\) calcularemos \(\mathbf{B}\mathbf{\hat{u}}^{(j-1)}\))
• For j=1,...,N (Vamos calculando los vectores de aproximación \(\hat{u}^{(j)}\))
  • For i=2,...,m-2
    • Set \(b_i =\frac{\lambda}{2} {\hat{u}}^{(0)}_{i-1}+(1-\lambda){\hat{u}}^{(0)}_i+\frac{\lambda}{2}\hat{u}^{(0)}_{i+1}\).
  • Set \(b_1 = (1-\lambda)\hat{u}_1^{(0)}+\frac{\lambda}{2} \hat{u}_2^{(0)}\).
  • Set \(b_{m-1}=\frac{\lambda}{2}\hat{u}_{m-2}^{(0)}+(1-\lambda)\hat{u}_{m-1}^{(0)}\).
  • Solve \(\mathbf{A}\mathbf{\hat{u}}=\mathbf{b}\).
  • Set \(\mathbf{\hat{u}}^{(0)}=\mathbf{\hat{u}}\).
• Print \(\mathbf{\hat{u}}\).

Usando técnicas avanzadas de álgebra lineal numérica, se puede demostrar que el método de Crank-Nicolson es incondicionalmente estable por lo que es el más usado al tener error de truncamiento cuadrático tanto a nivel temporal como espacial y ser, tal como hemos comentado, incondicionalmente estable.

El valor \(u(x,t)\) seria la “altura” de la cuerda en el lugar \(x\) y en el instante \(t>0\).

Las restricciones son las siguientes:

• Se supone conocidos la posición y la velocidad de la cuerda en el instante inicial \(t=0\): \[ u(x,0)=f(x),\ \frac{\partial u}{\partial t}(x,0)=g(x),\ 0\leq x\leq l. \]
• Se supone que la cuerda está fija en los extremos: \[ u(0,t)=u(l,t)=0. \]

El mallado se selecciona de la misma manera que hicimos en la ecuación parabólica del flujo del calor.

Es decir, recordemos que para seleccionar el paso espacial o para la variable \(x\), elegíamos un entero \(m>0\) y elegíamos \(h=\frac{l}{m}\).

A continuación, seleccionábamos un paso temporal \(k\).

El mallado consiste en los puntos \((x_i,t_j)\), donde \(x_i=ih\), \(i=0,1,\ldots,m\) y \(t_j=jk\), \(j=0,1,\ldots\)

El objetivo es hallar una aproximación \(\hat{u}(x_i,t_j)\), \(i=0,1,\ldots,m\), \(j=0,1,\ldots\) en los puntos del mallado.

Observar que los valores \(\hat{u}(x_i,t_j)\) para \(i=0,l\) son conocidos a partir de las restricciones: \[ \hat{u}(x_0,t_j)=\hat{u}(0,t_j)=0,\quad \hat{u}(x_m,t_j)=\hat{u}(l,t_j)=0, \] para \(j=0,1,\ldots\)

De la misma manera, los valores \(\hat{u}(x_i,0)\) y \(\frac{\partial u}{\partial t}(x_i,0)\) son conocidos a partir de las restricciones: \[ \hat{u}(x_i,0)=f(x_i),\ \frac{\partial u}{\partial t}(x_i,0)=g(x_i), \] para \(i=0,1,\ldots,m\)

Para hallar la aproximación \(\hat{u}(x_i,t_j)\) en el mallado consideramos las aproximaciones siguientes vistas en el capítulo de Derivación e Integración Numérica del curso Métodos numéricos con Python. Cálculo Numérico: \[ \begin{align*} \frac{\partial u}{\partial t^2}(x_i,t_j) = & \frac{u(x_i,t_j+k)-2u(x_i,t_j)+u(x_i,t_{j}-k)}{k^2}-\frac{k^2}{12}\frac{\partial^4 u}{\partial t^4}(x_i,\mu_j),\\ \frac{\partial^2 u}{\partial x^2}(x_i,t_j) = & \frac{u(x_i+h,t_j)-2u(x_i,t_j)+u(x_i-h,t_j)}{h^2}-\frac{h^2}{12}\frac{\partial^4 u}{\partial t^4}(\xi_i,t_j), \end{align*} \] donde \(\mu_j\in (t_{j-1},t_{j+1})\) y \(\xi_i\in (x_{i-1},x_{i+1})\).

Llamando \(\hat{u}_{ij}=\hat{u}(x_i,t_j)\), \(i=0,1,\ldots,m\) y \(j=0,1,\ldots\) y sustituyendo las aproximaciones anteriores en la EDP de ecuación de ondas obtenemos: \[ \begin{align*} \frac{\hat{u}_{i,j+1}-2\hat{u}_{ij}+\hat{u}_{i,j-1}}{k^2}-\alpha^2\frac{\hat{u}_{i+1,j}-2\hat{u}_{ij}+\hat{u}_{i-1,j}}{h^2}= & 0,\\ \hat{u}_{i,j+1}-2\hat{u}_{ij}+\hat{u}_{i,j-1}-\frac{k^2\alpha^2}{h^2}(\hat{u}_{i+1,j}-2\hat{u}_{ij}+\hat{u}_{i-1,j})= &0. \end{align*} \] Llamando \(\lambda=\frac{\alpha k}{h}\), el sistema de ecuaciones que verifica el vector de aproximaciones \(\mathbf{\hat{u}}^{(j+1)}\) es el siguiente: \[ \begin{align*} & \hat{u}_{i,j+1}-2\hat{u}_{ij}+\hat{u}_{i,j-1}-\lambda^2 (\hat{u}_{i+1,j}-2\hat{u}_{ij}+\hat{u}_{i-1,j})= 0,\\ & \hat{u}_{i,j+1} = 2(1-\lambda^2)\hat{u}_{ij}+\lambda^2\hat{u}_{i+1,j}+\lambda^2\hat{u}_{i-1,j}-\hat{u}_{i,j-1}, \end{align*} \] para \(i=1,\ldots,m-1\) y \(j=0,1,\ldots\) suponiendo conocidos \(\hat{u}_{i0}=f(x_i)\), \(\hat{u}_{0j}=\hat{u}_{mj}=0\).

Matricialmente: \[ \mathbf{\hat{u}}^{(j+1)}=\begin{bmatrix}\hat{u}_{1,j+1}\\\hat{u}_{2,j+1}\\\vdots\\\hat{u}_{m-1,j+1}\end{bmatrix}=\begin{bmatrix}2(1-\lambda^2)&\lambda^2&0&0&\ldots&0 \\\lambda^2&2(1-\lambda^2)&\lambda^2&0&\ldots&0 \\0&\lambda^2&2(1-\lambda^2)&\lambda^2&\ldots&0 \\\vdots&\vdots&\vdots&\ddots&\vdots&\vdots \\0&\ldots&0&\lambda^2&2(1-\lambda^2)&\lambda^2 \\0&\ldots&\ldots&0&\lambda^2&2(1-\lambda^2) \\\end{bmatrix}\begin{bmatrix}\hat{u}_{1,j}\\\hat{u}_{2,j}\\\vdots\\\hat{u}_{m-1,j}\end{bmatrix}-\begin{bmatrix}\hat{u}_{1,j-1}\\\hat{u}_{2,j-1}\\\vdots\\\hat{u}_{m-1,j-1}\end{bmatrix} \] El sistema anterior tiene un problema: cuando intentamos calcular el vector de aproximaciones \(\mathbf{\hat{u}}^{(1)}\) para \(t=t_1\) que resulta de aplicar la recurrencia anterior para \(j=0\), necesitamos conocer \(\mathbf{\hat{u}}^{(-1)}\) cuyos valores no conocemos!

Observar que aún no hemos aplicado la restricción \(\frac{\partial g}{\partial t}(x_i,0)=g(x_i)\). Dicha restricción nos permitirá conocer el vector de aproximaciones \(\mathbf{\hat{u}}^{(1)}\) de la forma siguiente:

Aproximamos \(\frac{\partial g}{\partial t}(x_i,0)\) usando la diferencia dividida siguiente: \[ \begin{align*} \frac{\partial g}{\partial t}(x_i,0)= & \frac{u(x_i,t_1)-u(x_i,0)}{k}-\frac{k}{2}\frac{\partial^2 u}{\partial t^2}(x_i,\tilde{\mu}_i)\\ = & \frac{u(x_i,t_1)-f(x_i)}{k}-\frac{k}{2}\frac{\partial^2 u}{\partial t^2}(x_i,\tilde{\mu}_i)=g(x_i), \end{align*} \]

donde \(\tilde{\mu}_i\in (0,t_{1})\). Despreciando el término del error, podemos calcular \(\hat{u}_{i1}=\hat{u}(x_i,t_1)\): \[ \hat{u}_{i1}=f(x_i)+kg(x_i), \] para \(i=1,2,\ldots,m-1\).

El problema es que hemos usado una aproximación con error de orden lineal temporal \(O(k)\) mientras que las aproximaciones usadas para deducir el sistema de ecuaciones que verifica \(\mathbf{\hat{u}}^{(j+1)}\) tiene error de orden cuadrático tanto temporal como espacial, \(O(k^2)\), \(O(h^2)\).

Esto hará que nuestras aproximaciones futuras \(\mathbf{\hat{u}}^{(j+1)}\), \(j=1,2,\ldots\) se ven afectadas por dicho error.

Para mejorar la aproximación \(\hat{u}_{i1}=\hat{u}(x_i,t_1)\) desarrollamos por Taylor la función \(u(x_i,t_1)\) alrededor del punto del mallado \((x_i,0)\), para \(i=1,\ldots,m-1\): \[ \begin{align*} u(x_i,t_1)=& u(x_i,0)+k\frac{\partial u}{\partial t}(x_i,0)+\frac{k^2}{2}\frac{\partial^2 u}{\partial t^2}(x_i,0)+\frac{k^3}{6}\frac{\partial^3 u}{\partial t^3}(x_i,\hat{\mu}_i)\\ = & f(x_i)+k g(x_i)+\frac{k^2}{2}\frac{\partial^2 u}{\partial t^2}(x_i,0)+\frac{k^3}{6}\frac{\partial^3 u}{\partial t^3}(x_i,\hat{\mu}_i), \end{align*} \] donde \(\hat{\mu}_i\in (0,t_1)\).

Ahora bien, como se tiene que verificar la EDP de la ecuación de ondas, \[ \frac{\partial^2 u}{\partial t^2}(x,t)-\alpha^2\frac{\partial^2 u}{\partial x^2}(x,t)=0,\ 0<x<l,\ t>0, \] en particular, \[ \frac{\partial^2 u}{\partial t^2}(x_i,0)-\alpha^2\frac{\partial^2 u}{\partial x^2}(x_i,0)=0,\ i=1,\ldots,m-1, \] es decir: \[ \frac{\partial^2 u}{\partial t^2}(x_i,0)=\alpha^2 \frac{\partial^2 u}{\partial x^2}(x_i,0),\ i=1,\ldots,m-1. \] Como \(u(x,0)=f(x)\), \(0\leq x\leq l\) si suponemos que \(f\) tiene derivada segunda, \(\frac{\partial^2 u}{\partial x^2}(x,0)=f''(x)\) y, en particular \(\frac{\partial^2 u}{\partial x^2}(x_i,0)=f''(x_i)\).

En resumen, \[ \frac{\partial^2 u}{\partial t^2}(x_i,0)=\alpha^2 f''(x_i),\ i=1,\ldots,m-1, \] y la aproximación que usaremos para calcular \(\hat{u}^{(1)}\) es la siguiente: \[ \hat{u}_{i1}=\hat{u}(x_i,t_1)=f(x_i)+kg(x_i)+\frac{\alpha^2 k^2}{2}f''(x_i), \] para \(i=1,\ldots,m-1\) cuyo error de trucamiento es de orden \(O(k^3)\).

En el caso de que no exista la derivada segunda de la función \(f\), podemos aproximar \(f''(x_i)\) usando la siguiente diferencia dividida con error de truncamiento \(O(h^2)\): \[ f''(x_i)=\frac{f(x_{i+1})-2f(x_i)+f(x_{i-1})}{h^2}-\frac{h^2}{12}f^{(4)}(\tilde{\xi}), \] donde \(\tilde{\xi}\in (x_{i-1},x_{i+1})\). En este caso, la aproximación usada para calcular \(\hat{u}^{(1)}\) es la siguiente:

\[ \begin{align*} \hat{u}_{i1}= & \hat{u}(x_i,t_1)=f(x_i)+kg(x_i)+\frac{\alpha^2 k^2}{2h^2}(f(x_{i+1})-2f(x_i)+f(x_{i-1})),\\ = & f(x_i)+kg(x_i)+\frac{\lambda^2}{2}(f(x_{i+1})-2f(x_i)+f(x_{i-1}))\\ = & \left(1-\lambda^2\right)f(x_i)+\frac{\lambda^2}{2}(f(x_{i+1})+f(x_{i-1}))+kg(x_i), \end{align*} \] para \(i=1,\ldots,m-1\) con error de truncamento temporal \(O(k^3)\) y error de truncamiento espacial \(O(h^2)\).

• INPUT \(m\), \(k\), \(l\), rutina que da \(f(x)\), rutina que da \(g(x)\) donde \(x\in [0,l]\), \(\alpha\), \(N\) (\(N\) es el número de iteraciones temporales a realizar)
• Set \(h=\frac{l}{m}\). (Hallamos el paso espacial)
• Set \(\mathbf{x}=\mathbf{0}\) (En el vector \(\mathbf{x}\) guardaremos el mallado espacial)
• For i=0,...,m (Hallamos el mallado espacial)
  • Set \(x_i=i\cdot h\).
• Set \(\lambda=\frac{k\alpha}{h}\).
• Set \(\mathbf{\hat{u}}^{(0)}=\mathbf{0}\) (Vamos a definir el vector \(\hat{u}^{(0)}\))
• For i=1,...,m-1
  • Set \(\hat{u}_i^{(0)} =f(x_i)\)

• Set \(\mathbf{\hat{u}}^{(1)}=\mathbf{0}\) (Vamos a definir el vector \(\hat{u}^{(1)}\))
• For i=1,...,m-1
  • Set \(\hat{u}_i^{(1)} =(1-\lambda^2) f(x_i)+\frac{\lambda^2}{2}(f(x_{i+1})+f(x_{i-1})+kg(x_i)\).
• Set \(\mathbf{\hat{u}}=\mathbf{0}\)

• For j=2,...,N (Vamos calculando los vectores de aproximación \(\hat{u}^{(j)}\))
  • For i=2,...,m-2
    • Set \(\hat{u}_i =\lambda^2 {\hat{u}}^{(1)}_{i-1}+2(1-\lambda^2){\hat{u}}^{(1)}_i+\lambda^2\hat{u}^{(1)}_{i+1}-\hat{u}_{i}^{(0)}\).
  • Set \(\hat{u}_1 = 2(1-\lambda^2)\hat{u}_1^{(1)}+\lambda^2 \hat{u}_2^{(1)}-\hat{u}_1^{(0)}\).
  • Set \(\hat{u}_{m-1}=\lambda^2 \hat{u}_{m-2}^{(1)}+2(1-\lambda^2)\hat{u}_{m-1}^{(1)}-\hat{u}_{m-1}^{(0)}\).
  • Set \(\mathbf{\hat{u}}^{(0)}=\mathbf{\hat{u}}^{(1)}\).
  • Set \(\mathbf{\hat{u}}^{(1)}=\mathbf{\hat{u}}\).
• Print \(\mathbf{\hat{u}}\).

El estudio de la estabilidad en el caso de la resolución de la EDP de ecuación de ondas es muy complicado ya que el esquema recursivo es de la forma: \[ \mathbf{\hat{u}}^{(j+1)}=\mathbf{A}\mathbf{\hat{u}}^{(j)}-\mathbf{\hat{u}}^{(j-1)}, \] donde \(\mathbf{A}\) recordemos que es la matriz de iteración siguiente: \[ \mathbf{A}=\begin{bmatrix}2(1-\lambda^2)&\lambda^2&0&0&\ldots&0 \\\lambda^2&2(1-\lambda^2)&\lambda^2&0&\ldots&0 \\0&\lambda^2&2(1-\lambda^2)&\lambda^2&\ldots&0 \\\vdots&\vdots&\vdots&\ddots&\vdots&\vdots \\0&\ldots&0&\lambda^2&2(1-\lambda^2)&\lambda^2 \\0&\ldots&\ldots&0&\lambda^2&2(1-\lambda^2) \\\end{bmatrix} \]

Para realizar el estudio de la estabilidad deberíamos hallar las expresiones \(P_j(\mathbf{A})\) y \(Q_j(\mathbf{A})\) tal que \[ \mathbf{\hat{u}}^{(j)}=P_j(\mathbf{A})\mathbf{\hat{u}}^{(1)}+Q_j(\mathbf{A})\mathbf{\hat{u}}^{(0)}. \] Por ejemplo, los primeros valores \(P_j\) y \(Q_j\) serían:

• \(P_2(\mathbf{A})=\mathbf{A}\) y \(Q_2(\mathbf{A})=-\mathbf{Id}\) ya que \(\mathbf{\hat{u}}^{(2)}=\mathbf{A}\mathbf{\hat{u}}^{(1)}-\mathbf{\hat{u}}^{(0)}\).
• \(P_3(\mathbf{A})=(\mathbf{A}^2-\mathbf{Id})\) y \(Q_3(\mathbf{A})=-\mathbf{A}\) ya que: \[ \mathbf{\hat{u}}^{(3)}=\mathbf{A}\mathbf{\hat{u}}^{(2)}-\mathbf{\hat{u}}^{(1)}=\mathbf{A}(\mathbf{A}\mathbf{\hat{u}}^{(1)}-\mathbf{\hat{u}}^{(0)})-\mathbf{\hat{u}}^{(1)}=(\mathbf{A}^2-\mathbf{Id})\mathbf{\hat{u}}^{(1)}-\mathbf{A}\mathbf{\hat{u}}^{(0)}. \]

Como se puede observar las funciones \(P_j\) y \(Q_j\) son polinomios de grados \(j\) y \(j-2\) respectivamente con respecto la matriz \(\mathbf{A}\).

El siguiente paso sería hallar el radio espectral de \(P_j(\mathbf{A})\) y \(Q_j(\mathbf{A})\) y ver cuando dichos radios espectrales son menores o iguales que \(1\).

Usando la propiedad de que los valores propios de las matrices \(P_j(\mathbf{A})\) y \(Q_j(\mathbf{A})\) son \(P_j(\mu)\) y \(Q_j(\mu)\), donde \(\mu\) es un valor propio de la matriz de iteración \(\mathbf{A}\), hallar los radios espectrales anteriores no es tan complicado.

Los valores propios de la matriz \(\mathbf{A}\) son relativamente sencillos de calcular ya que dicha matriz se puede escribir de la forma siguiente: \(\mathbf{A}=2\mathbf{B}\), siendo \(\mathbf{B}\) la matriz siguiente:

\[ \mathbf{B}=\begin{bmatrix}1-\lambda^2&\frac{\lambda^2}{2}&0&0&\ldots&0 \\\frac{\lambda^2}{2}&1-\lambda^2&\frac{\lambda^2}{2}&0&\ldots&0 \\0&\frac{\lambda^2}{2}&1-\lambda^2&\frac{\lambda^2}{2}&\ldots&0 \\\vdots&\vdots&\vdots&\ddots&\vdots&\vdots \\0&\ldots&0&\frac{\lambda^2}{2}&1-\lambda^2&\frac{\lambda^2}{2} \\0&\ldots&\ldots&0&\frac{\lambda^2}{2}&1-\lambda^2 \\\end{bmatrix}. \]

Los valores propios de \(\mathbf{A}\) serán de la forma \(2\mu\), donde \(\mu\) es un valor propio de la matriz \(\mathbf{B}\).

Ahora bien, para calcular los valores propios de la matriz \(\mathbf{B}\) podemos usar la proposición que calculaba los valores propios de la matriz de iteración en el caso del método de diferencias divididas hacia adelante para resolver la EDP parabólica del flujo de calor.

Recordemos que la matriz de iteración era: \[ \begin{bmatrix}(1-2\lambda)&\lambda&0&0&\ldots&0 \\\lambda&(1-2\lambda)&\lambda&0&\ldots&0 \\0&\lambda&(1-2\lambda)&\lambda&\ldots&0 \\\vdots&\vdots&\vdots&\ddots&\vdots&\vdots \\0&\ldots&0&\lambda&(1-2\lambda)&\lambda \\0&\ldots&\ldots&0&\lambda&(1-2\lambda) \\\end{bmatrix}. \]

Si cambiamos los papeles de \(\lambda\) por \(\frac{\lambda^2}{2}\), tenemos la matriz \(\mathbf{B}\).

Como los valores propios de la matriz de iteración para resolver la EDP parabólica del flujo de calor eran: \[ 1-4\lambda\left(\sin\left(\frac{\mathrm{i}\pi}{2m}\right)\right)^2,\ i=1,2,\ldots,m-1, \] los valores propios de la matriz \(\mathbf{B}\) serán: \[ \mu_i = 1-2\lambda^2\left(\sin\left(\frac{\mathrm{i}\pi}{2m}\right)\right)^2,\ i=1,2,\ldots,m-1, \]

y los valores propios de la matriz de iteración \(\mathbf{A}\) en el caso de la resolución de la EDP hiperbólica de la ecuación de ondas serán: \[ 2\mu_i = 2\cdot\left(1-2\lambda^2\left(\sin\left(\frac{\mathrm{i}\pi}{2m}\right)\right)^2\right),\ i=1,2,\ldots,m-1. \]

Entonces el método de resolución de la EDP hiperbólica de la ecuación de ondas será estable si los máximos siguientes: \[ \begin{align*} & \max_{i=1,\ldots,m-1}\left|P_j\left(2\cdot\left(1-2\lambda^2\left(\sin\left(\frac{\mathrm{i}\pi}{2m}\right)\right)^2\right)\right)\right|,\\ &\max_{i=1,\ldots,m-1}\left|Q_j\left(2\cdot\left(1-2\lambda^2\left(\sin\left(\frac{\mathrm{i}\pi}{2m}\right)\right)^2\right)\right)\right| \end{align*} \] son menores o iguales que \(1\).

Fijarse que dichos máximos serían los radios espectrales de las matrices \(P_j(\mathbf{A})\) y \(Q_j(\mathbf{A})\).

El único resultado es que si \(\lambda >1\), el método de resolución de la EDP hiperbólica de la ecuación de ondas no es estable.

Por tanto, es necesario que \(\lambda\leq 1\) para que dicho método sea estable.