6  Método de Crank Nicolson

El algoritmo introducido por J. Crank y P. Nicolson (Crank-Nicolson) en 1947 representa un esquema numérico ampliamente utilizado para resolver ecuaciones diferenciales parciales de tipo parabólico, como la ecuación de calor. Para comprender su fundamento, resulta ilustrativo considerar el problema de determinar la evolución temporal de la temperatura en una varilla metálica. Dado que la temperatura en cada punto varía de manera continua, es necesario discretizar el problema para su tratamiento computacional, lo que implica representar la varilla mediante un conjunto discreto de puntos y el tiempo mediante una secuencia de pasos finitos (Zill y Cullen 2008).

En este marco, las estrategias de solución numérica se clasifican principalmente en dos categorías. Por un lado, los métodos explícitos calculan el estado futuro del sistema a partir exclusivamente de información del estado presente, lo que los hace conceptualmente sencillos y computacionalmente eficientes por paso de tiempo. Sin embargo, presentan una limitación significativa: su estabilidad depende críticamente del tamaño del paso temporal. Si este paso excede un umbral crítico, la solución numérica puede volverse inestable, manifestando oscilaciones no físicas que divergen hacia infinito (Burden y Faires 2010).

Por otro lado, los métodos implícitos superan esta restricción de estabilidad al establecer una dependencia entre el estado futuro de un punto y el de sus vecinos en el mismo instante futuro. Esta característica garantiza estabilidad incondicional para una gama más amplia de parámetros, pero conlleva una mayor complejidad computacional, ya que requiere resolver un sistema de ecuaciones acoplado en cada paso de tiempo (Burden y Faires 2010).

El método de Crank-Nicolson surge como un esquema híbrido que sintetiza las ventajas de ambos enfoques. Se fundamenta en promediar la discretización espacial de la ecuación diferencial entre el instante de tiempo actual (\(n\)) y el futuro (\(n+1\)). Esta estrategia de promediado le confiere dos propiedades clave (Zill y Cullen 2008):

  1. Estabilidad Incondicional: A diferencia de los métodos explícitos, el esquema de Crank-Nicolson permanece estable para cualquier tamaño de paso temporal, lo que permite simulaciones más rápidas sin riesgo de divergencia.

  2. Precisión de Segundo Orden: Al ser un método de segundo orden en tiempo, el error de truncamiento local se reduce más ráticamente al disminuir el paso temporal, lo que se traduce en una mayor precisión global de la solución numérica comparado con métodos de primer orden.

En esencia, el algoritmo consiste en sustituir la segunda derivada parcial en \(c\frac{\partial^2 u}{\partial x^2} = \frac{\partial u}{\partial t}\) por el promedio de dos cocientes de diferencias centrales, uno evaluado en \(t\) y el otro en \(t+k\): \[ \begin{split} \dfrac{c}{2}& \left[\dfrac{u(x+h,t)-2u(x,t)+u(x-h,t)}{h^2}\right] + \\ + &\dfrac{c}{2} \left[\dfrac{u(x+h,t+k)-2u(x,t+k)+u(x-h,t+k)}{h^2} \right] \\ &= \dfrac{1}{k}[u(x,t+k)-u(x,t)] \end{split} \tag{6.1}\]

si se define \(\lambda = \frac{ck}{h^2}\) y \[ \begin{split} u(x+h,t) &=u_{i+1,j}, \quad u(x,t)=u_{ij}, \quad u(x-h,t)=u_{i-1,j}, \\ u(x+h,t+k) &=u_{i+1,j+1}, \quad u(x,t+k)=u_{i,j+1}, \quad u(x-h,t+k)=u_{i-1,j+1}, \end{split} \]

es posible reescribir a la Ec. 6.1 como: \[ -u_{i-1,j+1} + \alpha u_{i,j+1} - u_{i+1,j+1} = u_{i+1,j} - \beta u_{ij} + u_{i-1,j}, \tag{6.2}\]

donde \(\alpha=2(1+\frac{1}{\lambda})\), \(\ \beta=2(1-\frac{1}{\lambda})\), \(\ j=0,1,\dots, m-1, \ \ \ i=0,1,\dots, n-1\).

Para cada elección de \(j\), la ecuación diferencial Ec. 6.2 para \(i=0,1,\dots, n-1\) da \(n-1\) ecuaciones en \(n-1\) incógnitas \(u_{i,j+1}\). Debido a las condiciones de contorno preestablecidas, los valores de \(u_{i, j+1}\) se conocen para \(i=0\) y para \(i=n\). Por ejemplo, en el caso \(n=4\), el sistema de ecuaciones para determinar los valores aproximados de \(u\) en la línea de tiempo \((j+1)\) es:

\[ \begin{split} -u_{0,j+1} + \alpha u_{1,j+1} - u_{2,j+1} =& u_{2,j} - \beta u_{1,j} + u_{0,j} \\ -u_{1,j+1} + \alpha u_{2,j+1} - u_{3,j+1} =& u_{3,j} - \beta u_{2,j} + u_{1,j} \\ -u_{2,j+1} + \alpha u_{3,j+1} - u_{4,j+1} =& u_{4,j} - \beta u_{3,j} + u_{2,j} \end{split} \]

reordenando se llega a \[ \begin{matrix} \alpha u_{1,j+1} &-& u_{2,j+1}&& &= b_1 \\ -u_{1,j+1} &+& \alpha u_{2,j+1} &-& u_{3,j+1} &= b_2 \\ &-&u_{2,j+1} &+& \alpha u_{3,j+1} &= b_3 \end{matrix} \tag{6.3}\]

donde \[\begin{align*} b_1 &= u_{2,j} - \beta u_{1,j} + u_{0,j} + u_{0,j+1}, \\ b_2 &= u_{3,j} - \beta u_{2,j} + u_{1,j}, \\ b_3 &= u_{4,j} - \beta u_{3,j} + u_{2,j} + u_{4,j+1}. \end{align*}\]

En general, si utilizamos la ecuación diferencial Ec. 6.2 para determinar valores de \(u\) en la línea de tiempo \((j -1)\), es necesario resolver un sistema lineal \(\mathbf{AX=B}\), donde la matriz de coeficientes \(\mathbf{A}\) es una matriz tridiagonal, \[ A = \begin{pmatrix} \alpha & -1 & 0 & 0 & 0 & \cdots & 0 & 0 \\ -1 & \alpha & -1 & 0 & 0 & \cdots & 0 & 0 \\ 0 & -1 & \alpha & -1 & 0 & \cdots & 0 & 0 \\ 0 & 0 & -1 & \alpha & -1 & \cdots & 0 & 0 \\ \vdots & \vdots & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & 0 & 0 & 0 & \cdots & \alpha & -1 \\ 0 & 0 & 0 & 0 & 0 & \cdots & -1 & \alpha \end{pmatrix}, \]

y las componentes de la matriz columna \(\mathbf{B}\) son

\[\begin{align*} b_1 &= u_{2,j} - \beta u_{1,j} + u_{0,j} + u_{0,j+1}, \\ b_2 &= u_{3,j} - \beta u_{2,j} + u_{1,j}, \\ b_3 &= u_{4,j} - \beta u_{3,j} + u_{2,j}, \\ &\vdots \\ b_{n-1} &= u_{n,j} - \beta u_{n-1,j} + u_{n-2,j} + u_{n,j+1}. \end{align*}\]