Método de Crank-Nicolson

Método de diferencias finitas para la solución numérica de ecuaciones diferenciales parabólicas

En el análisis numérico , el método de Crank-Nicolson es un método de diferencias finitas utilizado para resolver numéricamente la ecuación del calor y ecuaciones diferenciales parciales similares . ^[1] Es un método de segundo orden en el tiempo. Es implícito en el tiempo, se puede escribir como un método Runge-Kutta implícito y es numéricamente estable . El método fue desarrollado por John Crank y Phyllis Nicolson en la década de 1940. ^[2]

Para las ecuaciones de difusión (y muchas otras ecuaciones), se puede demostrar que el método de Crank-Nicolson es incondicionalmente estable . ^[3] Sin embargo, las soluciones aproximadas aún pueden contener oscilaciones espurias (en descomposición) si la relación entre el paso de tiempo por la difusividad térmica y el cuadrado del paso de espacio, , es grande (normalmente, mayor que 1/2 según el análisis de estabilidad de Von Neumann ). Por esta razón, siempre que se necesitan pasos de tiempo grandes o una alta resolución espacial, a menudo se utiliza el método de Euler inverso menos preciso , que es estable e inmune a las oscilaciones. ^{[ cita requerida ]} ${\displaystyle \Delta t}$ ${\displaystyle \Delta x^{2}}$

Principio

La plantilla de Crank-Nicolson para un problema 1D

El método de Crank-Nicolson se basa en la regla trapezoidal , que proporciona una convergencia de segundo orden en el tiempo. Para ecuaciones lineales, la regla trapezoidal es equivalente al método del punto medio implícito ^{[ cita requerida ]} —el ejemplo más simple de un método de Runge-Kutta implícito de Gauss-Legendre— que también tiene la propiedad de ser un integrador geométrico . Por ejemplo, en una dimensión, supongamos que la ecuación diferencial parcial 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).}$$

Si se evalúa para y , la ecuación para el método de Crank-Nicolson es una combinación del método de Euler directo en y del método de Euler inverso en (tenga en cuenta, sin embargo, que el método en sí no es simplemente el promedio de esos dos métodos, ya que la ecuación de Euler inverso tiene una dependencia implícita de la solución): ${\displaystyle u(i\Delta x,n\Delta t)=u_{i}^{n}}$ ${\displaystyle F_{i}^{n}=F}$ ${\displaystyle i,n}$ ${\displaystyle u_{i}^{n}}$ ${\displaystyle n}$ ${\displaystyle n+1}$

Obsérvese que este es un método implícito : para obtener el "próximo" valor de en el tiempo, se debe resolver un sistema de ecuaciones algebraicas. Si la ecuación diferencial parcial no es lineal, la discretización también será no lineal, de modo que avanzar en el tiempo implicará la solución de un sistema de ecuaciones algebraicas no lineales, aunque las linealizaciones son posibles. En muchos problemas, especialmente la difusión lineal, el problema algebraico es tridiagonal y puede resolverse de manera eficiente con el algoritmo de matriz tridiagonal , que proporciona una solución directa rápida , a diferencia de la habitual para una matriz completa, en la que indica el tamaño de la matriz. ${\displaystyle u}$ ${\displaystyle {\mathcal {O}}(N)}$ ${\displaystyle {\mathcal {O}}(N^{3})}$ ${\displaystyle N}$

Ejemplo: difusión 1D

El método de Crank-Nicolson se aplica a menudo a problemas de difusión . Por ejemplo, para la difusión lineal,

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

Aplicando una discretización espacial de diferencias finitas para el lado derecho, la discretización de Crank-Nicolson es entonces

${\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, dejando , ${\displaystyle r={\frac {a\Delta t}{2(\Delta x)^{2}}}}$

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

Dado que se conocen los términos del lado derecho de la ecuación, este es un problema tridiagonal , por lo que puede resolverse de manera eficiente utilizando el algoritmo de matriz tridiagonal en lugar de la inversión de matriz mucho más costosa . ${\displaystyle u_{i}^{n+1}}$

Una ecuación cuasilineal, como (este es un ejemplo minimalista y no general)

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

conduciría a un sistema no lineal de ecuaciones algebraicas, que no podría resolverse fácilmente como se indicó anteriormente; sin embargo, en algunos casos es posible linealizar el problema utilizando el valor antiguo de , es decir, en lugar de . En otras ocasiones, puede ser posible realizar una estimación utilizando un método explícito y mantener la estabilidad. ${\displaystyle a}$ ${\displaystyle a_{i}^{n}(u)}$ ${\displaystyle a_{i}^{n+1}(u)}$ ${\displaystyle a_{i}^{n+1}(u)}$

Ejemplo: difusión 1D con advección para flujo constante, con conexiones de múltiples canales

Esta es una solución que se suele emplear para muchos propósitos cuando existe un problema de contaminación en arroyos o ríos en condiciones de caudal constante, pero la información se proporciona solo en una dimensión. A menudo, el problema se puede simplificar a un problema unidimensional y aun así proporcionar información útil.

Aquí modelamos la concentración de un contaminante soluto en el agua. Este problema se compone de tres partes: la ecuación de difusión conocida ( elegida como constante), un componente advectivo (que significa que el sistema está evolucionando en el espacio debido a un campo de velocidad), que elegimos como constante , y una interacción lateral entre canales longitudinales ( ): ${\displaystyle D_{x}}$ ${\displaystyle U_{x}}$ ${\displaystyle k}$

donde es la concentración del contaminante, y los subíndices y corresponden al canal anterior y siguiente . ${\displaystyle C}$ ${\displaystyle N}$ ${\displaystyle M}$

El método Crank-Nicolson (donde representa la posición y el tiempo) transforma cada componente de la EDP en lo siguiente: ${\displaystyle i}$ ${\displaystyle j}$

Ahora creamos las siguientes constantes para simplificar el álgebra:

${\displaystyle \lambda ={\frac {D_{x}\,\Delta t}{2\,\Delta x^{2}}},}$

${\displaystyle \alpha ={\frac {U_{x}\,\Delta t}{4\,\Delta x}},}$

${\displaystyle \beta ={\frac {k\,\Delta t}{2}},}$

y sustituimos ( 2 ), ( 3 ), ( 4 ), ( 5 ), ( 6 ), ( 7 ), y en ( 1 ). Luego colocamos los nuevos términos de tiempo a la izquierda ( ) y los términos de tiempo presente a la derecha ( ) para obtener ${\displaystyle \alpha }$ ${\displaystyle \beta }$ ${\displaystyle \lambda }$ ${\displaystyle j+1}$ ${\displaystyle j}$

${\displaystyle -\beta C_{Ni}^{j+1}-(\lambda +\alpha )C_{i-1}^{j+1}+(1+2\lambda +2\beta )C_{i}^{j+1}-(\lambda -\alpha )C_{i+1}^{j+1}-\beta C_{Mi}^{j+1}={}}$

${\displaystyle \qquad \beta C_{Ni}^{j}+(\lambda +\alpha )C_{i-1}^{j}+(1-2\lambda -2\beta )C_{i}^{j}+(\lambda -\alpha )C_{i+1}^{j}+\beta C_{Mi}^{j}.}$

Para modelar el primer canal, nos damos cuenta que solo puede estar en contacto con el siguiente canal ( ), por lo que la expresión se simplifica a ${\displaystyle M}$

${\displaystyle -(\lambda +\alpha )C_{i-1}^{j+1}+(1+2\lambda +\beta )C_{i}^{j+1}-(\lambda -\alpha )C_{i+1}^{j+1}-\beta C_{Mi}^{j+1}={}}$

${\displaystyle \qquad {}+(\lambda +\alpha )C_{i-1}^{j}+(1-2\lambda -\beta )C_{i}^{j}+(\lambda -\alpha )C_{i+1}^{j}+\beta C_{Mi}^{j}.}$

De la misma manera, para modelar el último canal, nos damos cuenta que solo puede estar en contacto con el canal anterior ( ), por lo que la expresión se simplifica a ${\displaystyle N}$

${\displaystyle -\beta C_{Ni}^{j+1}-(\lambda +\alpha )C_{i-1}^{j+1}+(1+2\lambda +\beta )C_{i}^{j+1}-(\lambda -\alpha )C_{i+1}^{j+1}={}}$

${\displaystyle \qquad \beta C_{Ni}^{j}+(\lambda +\alpha )C_{i-1}^{j}+(1-2\lambda -\beta )C_{i}^{j}+(\lambda -\alpha )C_{i+1}^{j}.}$

Para resolver este sistema lineal de ecuaciones, debemos ver ahora que las condiciones de contorno deben darse primero al comienzo de los canales:

${\displaystyle C_{0}^{j}}$: condición inicial para el canal en el paso de tiempo actual,

${\displaystyle C_{0}^{j+1}}$: condición inicial para el canal en el siguiente paso de tiempo,

${\displaystyle C_{N0}^{j}}$: condición inicial del canal anterior al analizado en el paso de tiempo actual,

${\displaystyle C_{M0}^{j}}$: condición inicial para el canal siguiente al analizado en el paso de tiempo actual.

Para la última celda de los canales ( ), la condición más conveniente pasa a ser adiabática, por lo que ${\displaystyle z}$

${\displaystyle \left.{\frac {\partial C}{\partial x}}\right|_{x=z}={\frac {C_{i+1}-C_{i-1}}{2\,\Delta x}}=0.}$

Esta condición se cumple si y solo si (independientemente de un valor nulo)

${\displaystyle C_{i+1}^{j+1}=C_{i-1}^{j+1}.}$

Resolvamos este problema (en forma matricial) para el caso de 3 canales y 5 nodos (incluida la condición de contorno inicial). Lo expresamos como un problema de sistema lineal:

${\displaystyle \mathbf {AA} \,\mathbf {C^{j+1}} =\mathbf {BB} \,\mathbf {C^{j}} +\mathbf {d} ,}$

dónde

${\displaystyle \mathbf {C^{j+1}} ={\begin{bmatrix}C_{11}^{j+1}\\C_{12}^{j+1}\\C_{13}^{j+1}\\C_{14}^{j+1}\\C_{21}^{j+1}\\C_{22}^{j+1}\\C_{23}^{j+1}\\C_{24}^{j+1}\\C_{31}^{j+1}\\C_{32}^{j+1}\\C_{33}^{j+1}\\C_{34}^{j+1}\end{bmatrix}},\quad \mathbf {C^{j}} ={\begin{bmatrix}C_{11}^{j}\\C_{12}^{j}\\C_{13}^{j}\\C_{14}^{j}\\C_{21}^{j}\\C_{22}^{j}\\C_{23}^{j}\\C_{24}^{j}\\C_{31}^{j}\\C_{32}^{j}\\C_{33}^{j}\\C_{34}^{j}\end{bmatrix}}.}$

Ahora debemos darnos cuenta de que AA y BB deben ser matrices compuestas de cuatro submatrices diferentes (recuerde que solo se consideran tres canales para este ejemplo, pero cubre la parte principal discutida anteriormente):

${\displaystyle \mathbf {AA} ={\begin{bmatrix}AA1&AA3&0\\AA3&AA2&AA3\\0&AA3&AA1\end{bmatrix}},\quad \mathbf {BB} ={\begin{bmatrix}BB1&-AA3&0\\-AA3&BB2&-AA3\\0&-AA3&BB1\end{bmatrix}},}$

donde los elementos mencionados anteriormente corresponden a las siguientes matrices y un 4×4 adicional lleno de ceros. Tenga en cuenta que los tamaños de AA y BB son 12×12:

${\displaystyle \mathbf {AA1} ={\begin{bmatrix}(1+2\lambda +\beta )&-(\lambda -\alpha )&0&0\\-(\lambda +\alpha )&(1+2\lambda +\beta )&-(\lambda -\alpha )&0\\0&-(\lambda +\alpha )&(1+2\lambda +\beta )&-(\lambda -\alpha )\\0&0&-2\lambda &(1+2\lambda +\beta )\end{bmatrix}},}$

${\displaystyle \mathbf {AA2} ={\begin{bmatrix}(1+2\lambda +2\beta )&-(\lambda -\alpha )&0&0\\-(\lambda +\alpha )&(1+2\lambda +2\beta )&-(\lambda -\alpha )&0\\0&-(\lambda +\alpha )&(1+2\lambda +2\beta )&-(\lambda -\alpha )\\0&0&-2\lambda &(1+2\lambda +2\beta )\end{bmatrix}},}$

${\displaystyle \mathbf {AA3} ={\begin{bmatrix}-\beta &0&0&0\\0&-\beta &0&0\\0&0&-\beta &0\\0&0&0&-\beta \end{bmatrix}},}$

${\displaystyle \mathbf {BB1} ={\begin{bmatrix}(1-2\lambda -\beta )&(\lambda -\alpha )&0&0\\(\lambda +\alpha )&(1-2\lambda -\beta )&(\lambda -\alpha )&0\\0&(\lambda +\alpha )&(1-2\lambda -\beta )&(\lambda -\alpha )\\0&0&2\lambda &(1-2\lambda -\beta )\end{bmatrix}},}$

${\displaystyle \mathbf {BB2} ={\begin{bmatrix}(1-2\lambda -2\beta )&(\lambda -\alpha )&0&0\\(\lambda +\alpha )&(1-2\lambda -2\beta )&(\lambda -\alpha )&0\\0&(\lambda +\alpha )&(1-2\lambda -2\beta )&(\lambda -\alpha )\\0&0&2\lambda &(1-2\lambda -2\beta )\end{bmatrix}}.}$

El vector d se utiliza aquí para contener las condiciones de contorno. En este ejemplo, es un vector 12×1:

${\displaystyle \mathbf {d} ={\begin{bmatrix}(\lambda +\alpha )(C_{10}^{j+1}+C_{10}^{j})\\0\\0\\0\\(\lambda +\alpha )(C_{20}^{j+1}+C_{20}^{j})\\0\\0\\0\\(\lambda +\alpha )(C_{30}^{j+1}+C_{30}^{j})\\0\\0\\0\end{bmatrix}}.}$

Para encontrar la concentración en cualquier momento, se debe iterar la siguiente ecuación:

${\displaystyle \mathbf {C^{j+1}} =\mathbf {AA} ^{-1}(\mathbf {BB} \,\mathbf {C^{j}} +\mathbf {d} ).}$

Ejemplo: difusión 2D

Al extenderse a dos dimensiones en una cuadrícula cartesiana uniforme , la derivación es similar y los resultados pueden conducir a un sistema de ecuaciones de banda diagonal en lugar de tridiagonales . La ecuación de calor bidimensional

${\displaystyle {\frac {\partial u}{\partial t}}=a\,\nabla ^{2}u,}$

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

se puede resolver con la discretización de Crank-Nicolson de

${\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})\\&+(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 utiliza una cuadrícula, de modo que . Esta ecuación se puede simplificar un poco reorganizando los términos y utilizando el número CFL ${\displaystyle \Delta x=\Delta y}$

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

En el esquema numérico de Crank-Nicolson, no se requiere un número CFL bajo para la estabilidad, pero sí para la precisión numérica. Ahora podemos escribir el esquema como

${\displaystyle (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)}$

${\displaystyle \qquad =(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).}$

Resolver un sistema lineal de este tipo es costoso. Por lo tanto, se puede implementar un método implícito de dirección alternada para resolver la EDP numérica, mediante el cual una dimensión se trata implícitamente y la otra dimensión explícitamente para la mitad del paso de tiempo asignado y viceversa para la mitad restante del paso de tiempo. El beneficio de esta estrategia es que el solucionador implícito solo requiere un algoritmo de matriz tridiagonal para ser resuelto. La diferencia entre la solución verdadera de Crank-Nicolson y la solución aproximada de ADI tiene un orden de precisión de y, por lo tanto, se puede ignorar con un paso de tiempo suficientemente pequeño. ^[4] ${\displaystyle O(\Delta t^{4})}$

Crank-Nicolson para problemas no lineales

Debido a que el método de Crank-Nicolson es implícito , generalmente es imposible resolverlo con exactitud. En cambio, se debe utilizar una técnica iterativa para converger a la solución. Una opción es utilizar el método de Newton para converger a la predicción, pero esto requiere el cálculo del jacobiano . Para un sistema de alta dimensión como los de la dinámica de fluidos computacional o la relatividad numérica , puede que no sea factible calcular este jacobiano.

Una alternativa libre de jacobianos es la iteración de punto fijo . Si es la velocidad del sistema, entonces la predicción de Crank-Nicolson será un punto fijo del mapa. Si la iteración del mapa no converge, el mapa parametrizado , con , puede comportarse mejor. En forma expandida, la fórmula de actualización es ${\displaystyle f}$ ${\displaystyle \Phi (x)=x_{0}+{\frac {h}{2}}\left[f(x_{0})+f(x)\right].}$ ${\displaystyle x^{(i+1)}=\Phi (x^{(i)})}$ ${\displaystyle \Theta (x,\alpha )=\alpha x+(1-\alpha )\Phi (x)}$ ${\displaystyle \alpha \in (0,1)}$

${\displaystyle x^{i+1}=\alpha x^{i}+(1-\alpha )\left[x_{0}+{\frac {h}{2}}\left(f(x_{0})+f(x^{i})\right)\right],}$

donde es la estimación actual y es el paso de tiempo anterior. ${\displaystyle x^{i}}$ ${\displaystyle x_{i-1}}$

Incluso para sistemas de alta dimensión, la iteración de este mapa puede converger sorprendentemente rápido.

Solución numérica de las ecuaciones de Navier-Stokes en forma de vorticidad. En este caso, era necesaria la iteración de punto fijo de Crank-Nicolson para que convergiera. ${\displaystyle \alpha =7/8}$

Aplicación en matemáticas financieras

Dado que se pueden modelar varios otros fenómenos con la ecuación de calor (a menudo llamada ecuación de difusión en matemáticas financieras ), el método de Crank-Nicolson también se ha aplicado a esas áreas. ^[5] En particular, la ecuación diferencial del modelo de fijación de precios de opciones de Black-Scholes se puede transformar en la ecuación de calor y, por lo tanto, se pueden obtener soluciones numéricas para la fijación de precios de opciones con el método de Crank-Nicolson.

La importancia de esto para las finanzas es que los problemas de fijación de precios de opciones, cuando se extienden más allá de los supuestos estándar (por ejemplo, incorporando dividendos cambiantes), no se pueden resolver en forma cerrada, pero se pueden resolver utilizando este método. Sin embargo, tenga en cuenta que para las condiciones finales no suaves (que ocurren para la mayoría de los instrumentos financieros), el método de Crank-Nicolson no es satisfactorio ya que las oscilaciones numéricas no se amortiguan. Para las opciones vainilla , esto da como resultado una oscilación en el valor gamma alrededor del precio de ejercicio . Por lo tanto, son necesarios pasos especiales de inicialización de amortiguación (por ejemplo, el método de diferencia finita completamente implícito).

Véase también

• Matemáticas financieras
• Regla trapezoidal

Referencias

1. Tuncer Cebeci (2002). Transferencia de calor por convección. Saltador. ISBN 0-9668461-4-1.
2. Crank, J.; Nicolson, P. (1947). "Un método práctico para la evaluación numérica de soluciones de ecuaciones diferenciales parciales del tipo de conducción de calor". Proc. Camb. Phil. Soc . 43 (1): 50–67. Bibcode :1947PCPS...43...50C. doi :10.1017/S0305004100023197. S2CID  16676040.
3. Thomas, JW (1995). Ecuaciones diferenciales parciales numéricas: métodos de diferencias finitas . Textos de matemáticas aplicadas. Vol. 22. Berlín, Nueva York: Springer-Verlag . ISBN. 978-0-387-97999-1.El ejemplo 3.3.2 muestra que Crank-Nicolson es incondicionalmente estable cuando se aplica a . ${\displaystyle u_{t}=au_{xx}}$
4. "Problemas parabólicos multidimensionales" (PDF) . Departamento de Ciencias de la Computación . RPI . Consultado el 29 de mayo de 2016 .
5. Wilmott, P.; Howison, S.; Dewynne, J. (1995). Las matemáticas de los derivados financieros: una introducción para estudiantes . Cambridge Univ. Press. ISBN 0-521-49789-2. Las matemáticas de los derivados financieros Wilmott.

Enlaces externos

• Técnicas de ecuaciones diferenciales numéricas para científicos e ingenieros, acceso abierto. Conferencias y códigos para ecuaciones diferenciales numéricas
• Un ejemplo de cómo aplicar e implementar el método Crank-Nicolson para la ecuación de advección