Oscilación extraña al resolver la ecuación de advección por diferencia finita con condiciones de contorno de Neumann completamente cerradas (reflexión en los límites)

Estoy tratando de resolver la ecuación de advección, pero aparece una extraña oscilación en la solución cuando la onda se refleja desde los límites. ¡Si alguien ha visto este artefacto antes, me interesaría saber la causa y cómo evitarlo!

Este es un gif animado, se abre en una ventana separada para ver la animación (¡solo se reproducirá una vez o no una vez que se haya almacenado en caché!)

Observe que la propagación parece altamente estable hasta que la onda comienza a reflejarse desde el primer límite. ¿Qué crees que podría estar pasando aquí? He pasado unos días revisando mi código y no puedo encontrar ningún error. Es extraño porque parece haber dos soluciones de propagación: una positiva y otra negativa; después de la reflexión desde el primer límite. Las soluciones parecen estar viajando a lo largo de puntos de malla adyacentes.

Los detalles de implementación siguen.

La ecuación de advección,

$\frac{\partial u}{\partial t} = \boldsymbol{v}\frac{\partial u}{\partial x}$

donde es la velocidad de propagación.$\boldsymbol{v}$

Crank-Nicolson es una discretización estable incondicional (enlace pdf) para la ecuación de advección, siempre que varíe lentamente en el espacio (solo contiene componentes de bajas frecuencias cuando Fourier se transforma).$u(x)$

La discretización que he aplicado es,

$\frac{\phi_{j}^{n+1} - \phi_{j}^{n}}{\Delta t} = \boldsymbol{v} \left[ \frac{1-\beta}{2\Delta x} \left( \phi_{j+1}^{n} - \phi_{j-1}^{n} \right) + \frac{\beta}{2\Delta x} \left( \phi_{j+1}^{n+1} - \phi_{j-1}^{n+1} \right) \right]$

Poner las incógnitas en el lado derecho permite que esto se escriba en forma lineal,

$\beta r\phi_{j-1}^{n+1} + \phi_{j}^{n+1} -\beta r\phi_{j+1}^{n+1} = -(1-\beta)r\phi_{j-1}^{n} + \phi_{j}^{n} + (1-\beta)r\phi_{j+1}^{n}$

donde (para tomar el tiempo promedio ponderado de manera uniforme entre el punto presente y futuro) y .r = v Δ $\beta=0.5$$r=\boldsymbol{v}\frac{\Delta t}{2\Delta x}$

Este conjunto de ecuaciones tiene la forma de matriz , donde,$A\cdot u^{n+1} = M\cdot u^n$

$\boldsymbol{A} = \left( \begin{matrix} 1 & -\beta r & & & 0 \\ \beta r & 1 & -\beta r & & \\ & \ddots & \ddots & \ddots & \\ & & \beta r & 1 & -\beta r \\ 0 & & & \beta r & 1 \\ \end{matrix} \right)$

$\boldsymbol{M} = \left( \begin{matrix} 1 & (1 - \beta)r & & & 0 \\ -(1 - \beta)r & 1 & (1 - \beta)r & & \\ & \ddots & \ddots & \ddots & \\ & & -(1 - \beta)r & 1 & (1 - \beta)r \\ 0 & & &-(1 - \beta)r & 1 \\ \end{matrix} \right)$

Los vectores y son los conocidos y desconocidos de la cantidad que queremos resolver.$u^n$$u^{n+1}$

Luego aplico condiciones de contorno de Neumann cerradas en los límites izquierdo y derecho. Por límites cerrados quiero decir en ambas interfaces. Para los límites cerrados resulta que (no mostraré mi trabajo aquí) solo necesitamos resolver la ecuación matricial anterior. Como señaló @DavidKetcheson, las ecuaciones matriciales anteriores en realidad describen las condiciones de contorno de Dirichlet . Para las condiciones de contorno de Neumann,$\frac{\partial u}{\partial x} = 0$

$\boldsymbol{A} = \left( \begin{matrix} 1 & 0 & & & 0 \\ \beta r & 1 & -\beta r & & \\ & \ddots & \ddots & \ddots & \\ & & \beta r & 1 & -\beta r \\ 0 & & & 0 & 1 \\ \end{matrix} \right)$

Actualizar

El comportamiento parece bastante independiente de la elección de las constantes que uso, pero estos son los valores para la gráfica que ves arriba:

• $\boldsymbol{v}$ = 2
• dx = 0.2
• dt = 0.005
• $\sigma$ = 2 (gaussiano hwhm)
• $\beta$ = 0.5

Actualización II

Una simulación con un coeficiente de difusión distinto de cero, (ver comentarios más abajo), la oscilación desaparece, pero la onda ya no se refleja. No entiendo porque$D=1$

La ecuación que está resolviendo no permite soluciones correctas, por lo que no existe tal cosa como un condición límite que refleje esta ecuación. Si considera las características, se dará cuenta de que solo puede imponer una condición de límite en el límite correcto. Está tratando de imponer una condición de límite de Dirichlet homogénea en el límite izquierdo, que es matemáticamente inválido.

$x-\nu t = C$$C$

A diferencia de la ecuación, su esquema numérico no admite soluciones que van derecho. Los modos de marcha a la derecha se denominan modos parásitos e implican frecuencias muy altas. Observe que la onda que va hacia la derecha es un paquete de ondas de diente de sierra, asociado con las frecuencias más altas que se pueden representar en su cuadrícula. Esa ola es puramente un artefacto numérico, creado por su discretización.

Para enfatizar: no ha escrito el problema completo del valor del límite inicial que está tratando de resolver. Si lo hace, quedará claro que no es un problema matemáticamente bien planteado.

Sin embargo, me alegro de que haya publicado esto aquí, ya que es una hermosa ilustración de lo que puede suceder cuando discrepa un problema que no está bien planteado, y del fenómeno de los modos parasitarios. Un gran +1 para tu pregunta de mí.

Aprendí mucho de las respuestas anteriores. Quiero incluir esta respuesta porque creo que ofrece diferentes puntos de vista sobre el problema.

$u(x,t)=f(x-ct)$

$u(x,t_0)=p(x)$$x \in (-\infty, \infty)$$p[x-c(t-t_0)]$ . No hay límite y esa es la solución.

$x \in [a,b]$$a$$b$

$t_0$

Esto produce un pulso que corre hacia la derecha hasta que desaparece en el borde derecho.

presione aquí para animación en Dirichlet en el límite izquierdo

Todavía recibo algo de ruido que no puedo entender (¿alguien podría ayudarme aquí, por favor?)

La otra opción es imponer condiciones de contorno periódicas. Es decir, en lugar de imponer la condición de límite de Dirichlet a la izquierda, podemos imponer el paquete de ondas que es consistente con el límite a la izquierda. Es decir:

$a-c(t-t_0) <a$$t>t_0$$c>0$$[a,b]$$b-a$$a-c(t-t_0)+b-a=a+b(t-t_0)$$u(a,t_0)=p[b-c(t-t_0]$

Este enlace muestra lo que yo llamaría condiciones de contorno periódicas.

Hice las animaciones en Python y el esquema es un esquema de Crank-Nicholson como se indica en la pregunta aquí.

Todavía estoy luchando con el patrón de ruido después de que la ola golpea el primer límite (derecho).