Método de manivela-Nicolson

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 pueden contener oscilaciones espurias (en desintegración) si la relación entre el paso de tiempo Δ t multiplicado por la difusividad térmica y el cuadrado del paso del espacio, Δ x ² , es grande (normalmente, mayor que 1/2 por Análisis de estabilidad de Von Neumann ). Por esta razón, siempre que se necesitan grandes intervalos de tiempo o una alta resolución espacial, a menudo se utiliza el método de Euler hacia atrás menos preciso , que es estable e inmune a las oscilaciones. ^{[ cita requerida ]}

El método

La plantilla Crank-Nicolson para un problema 1D

El método de Crank-Nicolson se basa en la regla trapezoidal , lo que proporciona una convergencia de segundo orden en el tiempo. Para las 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 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, suponga 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 ${\displaystyle u(i\Delta x,n\Delta t)=u_{i}^{n}}$ y ${\displaystyle F_{i}^{n}=F}$ evaluado para ${\displaystyle i,n}$ y ${\displaystyle u_{i}^{n}}$, la ecuación para el método de Crank-Nicolson es una combinación del método de Euler directo en${\displaystyle n}$y el método de Euler hacia atrás en n  + 1 (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 hacia atrás tiene una dependencia implícita 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 {\text{(forward Euler)}},}$

${\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 {\text{(backward Euler)}},}$

${\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 {\text{(Crank–Nicolson)}}.}$

Tenga en cuenta que este es un método implícito : para obtener el "siguiente" valor de u 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 será no lineal, por lo 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 en la difusión lineal, el problema algebraico es tridiagonal y puede resolverse eficientemente con el algoritmo de matriz tridiagonal , que da una rápida${\displaystyle {\mathcal {O}}(N)}$ solución directa, a diferencia de la habitual ${\displaystyle {\mathcal {O}}(N^{3})}$ para una matriz completa, en la que ${\displaystyle N}$ indica el tamaño de la matriz.

Ejemplo: difusión 1D

El método Crank-Nicolson se aplica a menudo a problemas de difusión . Como 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 se

${\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 dejar ${\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 , de modo que${\displaystyle u_{i}^{n+1}}$puede resolverse de manera eficiente utilizando el algoritmo de matriz tridiagonal a favor de una inversión de matriz mucho más costosa .

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 antes; Sin embargo, en algunos casos es posible linealizar el problema utilizando el valor anterior para${\displaystyle a}$, es decir, ${\displaystyle a_{i}^{n}(u)}$ en vez de ${\displaystyle a_{i}^{n+1}(u)}$. Otras veces, es posible estimar${\displaystyle a_{i}^{n+1}(u)}$ utilizando un método explícito y mantener la estabilidad.

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

Esta es una solución que generalmente se emplea para muchos propósitos cuando existe un problema de contaminación en arroyos o ríos en condiciones de flujo constante, pero la información se proporciona en una sola dimensión. A menudo, el problema se puede simplificar en un problema unidimensional y aún así producir 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 (${\displaystyle D_{x}}$elegida como constante), una componente advectiva (lo que significa que el sistema está evolucionando en el espacio debido a un campo de velocidad), que elegimos como una constante Ux , y una interacción lateral entre canales longitudinales ( k ):

${\displaystyle {\frac {\partial C}{\partial t}}=D_{x}{\frac {\partial ^{2}C}{\partial x^{2}}}-U_{x}{\frac {\partial C}{\partial x}}-k(C-C_{N})-k(C-C_{M}),}$

( 1 )

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

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

${\displaystyle {\frac {\partial C}{\partial t}}\Rightarrow {\frac {C_{i}^{j+1}-C_{i}^{j}}{\Delta t}},}$

( 2 )

${\displaystyle {\frac {\partial ^{2}C}{\partial x^{2}}}\Rightarrow {\frac {1}{2(\Delta x)^{2}}}\left((C_{i+1}^{j+1}-2C_{i}^{j+1}+C_{i-1}^{j+1})+(C_{i+1}^{j}-2C_{i}^{j}+C_{i-1}^{j})\right),}$

( 3 )

${\displaystyle {\frac {\partial C}{\partial x}}\Rightarrow {\frac {1}{2}}\left({\frac {(C_{i+1}^{j+1}-C_{i-1}^{j+1})}{2(\Delta x)}}+{\frac {(C_{i+1}^{j}-C_{i-1}^{j})}{2(\Delta x)}}\right),}$

( 4 )

${\displaystyle C\Rightarrow {\frac {1}{2}}(C_{i}^{j+1}+C_{i}^{j}),}$

( 5 )

${\displaystyle C_{N}\Rightarrow {\frac {1}{2}}(C_{Ni}^{j+1}+C_{Ni}^{j}),}$

( 6 )

${\displaystyle C_{M}\Rightarrow {\frac {1}{2}}(C_{Mi}^{j+1}+C_{Mi}^{j}).}$

( 7 )

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 ( j  + 1) y los términos de tiempo presente a la derecha ( j ) para obtener

${\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 de que solo puede estar en contacto con el siguiente canal ( M ), por lo que la expresión se simplifica a

${\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 forma, para modelar el último canal, nos damos cuenta de que solo puede estar en contacto con el canal anterior ( N ), por lo que la expresión se simplifica a

${\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 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 intervalo de tiempo actual,

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

Para la última celda de los canales ( z ), la condición más conveniente se convierte en adiabática, por lo que

${\displaystyle {\frac {\partial C}{\partial x}}_{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 de matriz) para el caso de 3 canales y 5 nodos (incluida la condición de frontera 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 hechos de 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 las siguientes matrices, y 4 × 4 adicionales llenos 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 usa 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 en 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 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)}$

puede resolverse 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}}}$

asumiendo que se usa una cuadrícula cuadrada, de modo que ${\displaystyle \Delta x=\Delta y}$. Esta ecuación se puede simplificar un poco reordenando los términos y usando el número CFL

${\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, 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 no es práctico debido a la complejidad de tiempo extremadamente alta de resolver un sistema lineal por medio de la eliminación gaussiana o incluso el algoritmo de Strassen . Por tanto , se puede implementar un método implícito de dirección alterna para resolver el PDE numérico, mediante el cual una dimensión se trata implícitamente y otra dimensión explícitamente para la mitad del paso de tiempo asignado y, a la inversa, 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 verdadera solución de Crank-Nicolson y la solución aproximada de ADI tiene un orden de precisión de${\displaystyle O(\Delta t^{4})}$y por lo tanto puede ignorarse con un paso de tiempo suficientemente pequeño. ^[4]

Crank-Nicolson para problemas no lineales

Debido a que el método de Crank-Nicolson es implícito , generalmente es imposible resolver el futuro predicho cuando la ecuación diferencial no es lineal. En cambio, se debe utilizar una técnica iterativa para converger a la predicció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 ser inviable calcular este jacobiano.

Una alternativa es definir un mapa, para el cual la predicción de Crank-Nicolson es un punto fijo. Las iteraciones inteligentes de este mapa podrían converger en una predicción válida. Tenga en cuenta que para los problemas no lineales, no se garantiza que la predicción exista o que sea única. Para tamaños de paso lo suficientemente pequeños, estas técnicas iterativas deberían tener éxito. Dejar${\displaystyle {\dot {x}}=f(x),}$ dónde ${\displaystyle x\in \mathbb {R} ^{n}}$ y ${\displaystyle f:\mathbb {R} ^{n}\rightarrow \mathbb {R} ^{n}}$. La predicción de Crank-Nicolson será un punto fijo del mapa${\displaystyle \Phi (x)=x_{0}+{\frac {h}{2}}\left[f(x_{0})+f(x)\right].}$ Si la iteración del mapa ${\displaystyle x^{(i+1)}=\Phi (x^{(i)})}$ no converge, el mapa parametrizado ${\displaystyle \Theta (x,\alpha )=\alpha x+(1-\alpha )\Phi (x)}$, con ${\displaystyle \alpha \in (0,1)}$, puede comportarse mejor.

"> Reproducir medios

Una simulación dimensional de 128x128 de vórtices en una caja 2D, donde la integración numérica se realiza mediante Crank-Nicolson iterativo. Se necesitaron ~ 30 iteraciones en cada paso de tiempo para converger.