Método Crank-Nicolson

Método de diferencias finitas para resolver numéricamente ecuaciones diferenciales parabólicas

En análisis numérico , el método de Crank-Nicolson es un método de diferencias finitas que se utiliza para resolver numéricamente la ecuación de calor y ecuaciones diferenciales parciales similares . ^[1] Es un método de segundo orden en el tiempo. Está implícito en el tiempo, puede escribirse como un método implícito de Runge-Kutta y es numéricamente estable . El método fue desarrollado por John Crank y Phyllis Nicolson a mediados del siglo XX. ^[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 (decrecientes) si la relación entre el paso de tiempo multiplicado por la difusividad térmica y el cuadrado del paso espacial, 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 grandes pasos de tiempo o una alta resolución espacial, se suele utilizar el método de Euler hacia atrás, menos preciso, que es a la vez estable e inmune a las oscilaciones. ^{[ cita necesaria ]} ${\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 del trapecio , que proporciona una convergencia de segundo orden en el tiempo. Para ecuaciones lineales, la regla trapezoidal es equivalente al método implícito del punto medio ^{[ cita necesaria ]}  , el ejemplo más simple de un método implícito de Runge-Kutta 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).}$$

Dejando y evaluado para y , la ecuación para el método de Crank-Nicolson es una combinación del método de Euler directo en y el 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 el método inverso La ecuación de Euler 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}$

Tenga en cuenta que este es un método implícito : para obtener el "siguiente" valor de en el tiempo, se debe resolver un sistema de ecuaciones algebraicas. Si la ecuación diferencial parcial es no lineal, la discretización también lo será, 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 de difusión lineal, el problema algebraico es tridiagonal y puede resolverse eficientemente con el algoritmo de matriz tridiagonal , que da una solución directa rápida, a diferencia de lo habitual para una matriz completa, en el que se 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 Crank-Nicolson se aplica a menudo a problemas de difusión . Como ejemplo, para difusión lineal,

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

aplicando una discretización espacial en 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 eficientemente 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 . Otras veces, puede ser posible estimar 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 múltiples conexiones de canales

Esta es una solución que generalmente se emplea para muchos propósitos cuando hay un problema de contaminación en arroyos o ríos en condiciones de flujo constante, pero la información se brinda en una sola dimensión. A menudo, el problema se puede simplificar a un problema unidimensional y aun así generar información útil.

Aquí modelamos la concentración de un soluto contaminante en el agua. Este problema se compone de tres partes: la ecuación de difusión conocida ( elegida como constante), una componente advectiva (que significa que el sistema está evolucionando en el espacio debido a un campo de velocidades), 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 del PDE 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 sustituya ( 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, ahora debemos ver que primero se deben dar condiciones de frontera al inicio 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,

${\displaystyle C_{N0}^{j}}$: condición inicial para el 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 sólo 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). Expresamos esto 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 arreglos formados por cuatro subarreglos 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 los siguientes arrays, 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 aquí se utiliza para mantener las condiciones de contorno. En este ejemplo es un vector de 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

Cuando se extiende a dos dimensiones en una cuadrícula cartesiana uniforme , la derivación es similar y los resultados pueden conducir a un sistema de ecuaciones diagonales de banda en lugar de tridiagonales . La ecuación del 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 reordenando los términos y usando el número CFL ${\displaystyle \Delta x=\Delta y}$

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

Para el esquema numérico de Crank-Nicolson, no se requiere un número CFL bajo para la estabilidad; sin embargo, sí se requiere 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 alterna para resolver la PDE numérica, mediante el cual una dimensión se trata implícitamente y otra dimensión explícitamente durante la mitad del paso de tiempo asignado y, a la inversa, durante 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 resolverlo. La diferencia entre la verdadera solución de Crank-Nicolson y la solución aproximada de ADI tiene un orden de precisión de y, por lo tanto, puede ignorarse en 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 exactamente. En cambio, se debe utilizar una técnica iterativa para converger hacia la solución. Una opción es utilizar el método de Newton para converger en 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 resultar inviable calcular este jacobiano.

Una alternativa libre de jacobiano 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 ampliada, 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],}$

¿Dónde está la estimación actual y 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.

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

Aplicación en matemáticas financieras

Debido a que se pueden modelar muchos otros fenómenos con la ecuación del calor (a menudo llamada ecuación de difusión en matemáticas financieras ), el método Crank-Nicolson también se ha aplicado a esas áreas. ^[5] En particular, la ecuación diferencial del modelo de valoración 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 valoración 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 pueden resolverse de forma cerrada, pero pueden resolverse utilizando este método. Sin embargo, tenga en cuenta que para condiciones finales no uniformes (lo que ocurre con 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 resulta en una oscilación en el valor gamma alrededor del precio de ejercicio . Por lo tanto, son necesarios pasos especiales de inicialización de la amortiguación (por ejemplo, método de diferencias finitas totalmente implícito).

Ver 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. Manivela, 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. Fil. Soc . 43 (1): 50–67. Código Bib : 1947PCPS...43...50C. doi :10.1017/S0305004100023197. S2CID  16676040.
3. Thomas, JW (1995). Ecuaciones diferenciales parciales numéricas: métodos de diferencias finitas . Textos en Matemática Aplicada. 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 Informática . 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 el estudiante . Universidad de Cambridge. Prensa. ISBN 0-521-49789-2. Las matemáticas de los derivados financieros Wilmott.

enlaces externos

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