Método de Crank-Nicolson

En el campo del análisis numérico, el método de Crank-Nicolson es un método de diferencias finitas usado para la resolución numérica de ecuaciones en derivadas parciales, tales como la ecuación del calor.^[1]​ Se trata de un método de segundo orden en tiempo, implícito y numéricamente estable. El método fue desarrollado por John Crank y Phyllis Nicolson a mediados del siglo XX.^[2]​

Para ecuaciones difusivas (y para muchos otros tipos de ecuaciones), puede demostrarse que el método de Crank–Nicolson es incondicionalmente estable.^[3]​ Sin embargo, las soluciones aproximadas pueden aún contener algunas oscilaciones espurias (decrecientes) si el ratio entre el paso de tiempo y el cuadrado de la malla en espacio es grande (típicamente, mayor que 1/2). Por este motivo, siempre que sean necesarios pasos de tiempo grande o pequeñas mallas espaciales, puede considerarse el uso del método de Euler implícito, que es a la vez estable e inmune a oscilaciones (aunque es de menor orden).

Descripción del método de Crank-Nicolson

Esquema de Crank–Nicolson para un problema 1D.

El método de Crank–Nicolson se basa en diferencias centrales en espacio y en la Regla del trapecio en tiempo, resultando así en un método con convergencia de segundo orden en tiempo. Por ejemplo en una dimensión, si la ecuación en derivadas parciales es:

${\displaystyle {\frac {\partial u}{\partial t}}=F\left(u,x,t,{\frac {\partial u}{\partial x}},{\frac {\partial ^{2}u}{\partial x^{2}}}\right)}$

entonces, tomando ${\displaystyle u(i\Delta x,n\Delta t)=u_{i}^{n}\,}$, la ecuación para el método de Crank–Nicolson es una combinación del método de Euler implícito y el método de Euler explícito en la etapa de tiempo n + 1 (obsérvese que, sin embargo, el método en sí no es simplemente la media de estos dos métodos, puesto que la ecuación depende implícitamente de la solución):

${\displaystyle {\frac {u_{i}^{n+1}-u_{i}^{n}}{\Delta t}}=F_{i}^{n}\left(u,x,t,{\frac {\partial u}{\partial x}},{\frac {\partial ^{2}u}{\partial x^{2}}}\right)\qquad {\mbox{(Euler explícito)}}}$

${\displaystyle {\frac {u_{i}^{n+1}-u_{i}^{n}}{\Delta t}}=F_{i}^{n+1}\left(u,x,t,{\frac {\partial u}{\partial x}},{\frac {\partial ^{2}u}{\partial x^{2}}}\right)\qquad {\mbox{(Euler implícito)}}}$

${\displaystyle {\frac {u_{i}^{n+1}-u_{i}^{n}}{\Delta t}}={\frac {1}{2}}\left[F_{i}^{n+1}\left(u,x,t,{\frac {\partial u}{\partial x}},{\frac {\partial ^{2}u}{\partial x^{2}}}\right)+F_{i}^{n}\left(u,x,t,{\frac {\partial u}{\partial x}},{\frac {\partial ^{2}u}{\partial x^{2}}}\right)\right]\qquad {\mbox{(Crank-Nicolson)}}}$

La función F debe ser discretizada espacialmente mediante diferencias centrales.

Como se puede ver, se trata de un método implícito: para obtener el valor de u en el "siguiente instante de tiempo", debe resolverse un sistema de ecuaciones algebraicas. Si la ecuación en derivadas parciales no es lineal, la discretización tampoco lo será, de forma que para avanzar en tiempo resultará necesario resolver un sistema de ecuaciones algebraicas no lineales, siendo necesario emplear para ello algún tipo de método numérico. En muchos casos, especialmente en el de difusión lineal, el sistema de ecuaciones algebraicas tiene asociada una matriz tridiagonal y puede ser resuelto eficientemente mediante algoritmos adaptados a este tipo de matrices, que son de orden ${\displaystyle {\mathcal {O}}(n)}$, mientras que los métodos habituales (para resolución de sistemas lineales con matrices llenas) son de orden ${\displaystyle {\mathcal {O}}(n^{3})}$.

Ejemplo: difusión 1D.

El método de Crank–Nicolson se usa a menudo en problemas difusivos. Por ejemplo, en el caso de difusión lineal,

${\displaystyle {\partial u \over \partial t}=a{\frac {\partial ^{2}u}{\partial x^{2}}}}$

cuya discretización mediante Crank–Nicolson es:

${\displaystyle {\frac {u_{i}^{n+1}-u_{i}^{n}}{\Delta t}}={\frac {a}{2(\Delta x)^{2}}}\left((u_{i+1}^{n+1}-2u_{i}^{n+1}+u_{i-1}^{n+1})+(u_{i+1}^{n}-2u_{i}^{n}+u_{i-1}^{n})\right)}$

o, tomando ${\displaystyle r={\frac {a\Delta t}{(\Delta x)^{2}}}}$:

${\displaystyle -ru_{i+1}^{n+1}+(2+2r)u_{i}^{n+1}-ru_{i-1}^{n+1}=ru_{i+1}^{n}+(2-2r)u_{i}^{n}+ru_{i-1}^{n}\,}$,

que es un sistema con matriz tridiagonal, por lo que, como se ha comentado, ${\displaystyle u_{i}^{n+1}\,}$ puede ser resuelto mediante métodos mucho más eficientes.

Sin embargo, si introducimos una ecuación no sea lineal como la siguiente (una ecuación cuasi-lineal):

${\displaystyle {\frac {\partial u}{\partial t}}=a(u){\frac {\partial ^{2}u}{\partial x^{2}}}}$

aun tratándose de un ejemplo minimalista y no general, éste conduce a un sistema de ecuaciones algebraicas no líneas que no puede ser resuelto tan fácilmente como el anterior. Aun así, en algunos casos es posible "linealizar" el problema, usando el valor de ${\displaystyle a}$ en el instante de tiempo anterior, es decir, ${\displaystyle a_{i}^{n}(u)\,}$ en lugar de ${\displaystyle a_{i}^{n+1}(u)\,}$. En otras ocasiones, puede ser imposible el estimar ${\displaystyle a_{i}^{n+1}(u)\,}$ mediante un método explícito y, a la vez, mantener la estabilidad.

Ejemplo: difusión 2D.

Cuando se extiende a dos dimensiones en una malla Cartesiana, el método se deduce de forma similar y el resultado conduce de nuevo a un sistema de ecuaciones algebraicas con matriz banda. Por ejemplo, la ecuación del calor bidimensional:

${\displaystyle {\frac {\partial u}{\partial t}}=a\left({\frac {\partial ^{2}u}{\partial x^{2}}}+{\frac {\partial ^{2}u}{\partial y^{2}}}\right)}$

puede resolverse mediante una discretización de Crank–Nicolson:

${\displaystyle {\begin{aligned}u_{i,j}^{n+1}&=u_{i,j}^{n}+{\frac {1}{2}}{\frac {a\Delta t}{(\Delta x)^{2}}}{\big [}(u_{i+1,j}^{n+1}+u_{i-1,j}^{n+1}+u_{i,j+1}^{n+1}+u_{i,j-1}^{n+1}-4u_{i,j}^{n+1})\\&\qquad {}+(u_{i+1,j}^{n}+u_{i-1,j}^{n}+u_{i,j+1}^{n}+u_{i,j-1}^{n}-4u_{i,j}^{n}){\big ]}\end{aligned}}}$

suponiendo que se usa una malla cuadrada en la que ${\displaystyle \Delta x=\Delta y}$. Estas ecuaciones pueden ser simplificadas algo, ordenando sus términos y usando el número de Courant:

${\displaystyle \mu ={\frac {a\Delta t}{(\Delta x)^{2}}}}$

En el esquema numérico de Crank–Nicolson, el tener un número de Courant pequeño no es un requisito para la estabilidad, aunque es necesario para la precisión numérica. Ahora podemos escribir el esquema como:

${\displaystyle {\begin{aligned}&(1+2\mu )u_{i,j}^{n+1}-{\frac {\mu }{2}}\left(u_{i+1,j}^{n+1}+u_{i-1,j}^{n+1}+u_{i,j+1}^{n+1}+u_{i,j-1}^{n+1}\right)\\&\quad =(1-2\mu )u_{i,j}^{n}+{\frac {\mu }{2}}\left(u_{i+1,j}^{n}+u_{i-1,j}^{n}+u_{i,j+1}^{n}+u_{i,j-1}^{n}\right)\end{aligned}}}$

Véase también

• Ecuaciones en derivadas parciales

Referencias

1. Tuncer Cebeci (2002). Convective Heat Transfer. Springer. ISBN 0966846141. 
2. Crank, J.; Nicolson, P. (1947). «A practical method for numerical evaluation of solutions of partial differential equations of the heat conduction type». Proc. Camb. Phil. Soc. 43 (1): 50-67. doi:10.1007/BF02127704. 
3. Thomas, J. W. (1995). Numerical Partial Differential Equations: Finite Difference Methods. Texts in Applied Mathematics 22. Berlin, New York: Springer-Verlag. ISBN 978-0-387-97999-1. . En el ejemplo 3.3.2 se demuestra que Crank–Nicolson es incondicionalmente estable cuando se aplica a la ecuación ${\displaystyle u_{t}=au_{xx}}$

Enlaces externos

• Module for Parabolic P.D.E.'s