Étrange oscillation lors de la résolution de l'équation d'advection par différence finie avec des conditions aux limites de Neumann complètement fermées (réflexion aux limites)

J'essaie de résoudre l'équation d'advection mais j'ai une étrange oscillation apparaissant dans la solution lorsque l'onde se réfléchit depuis les limites. Si quelqu'un a déjà vu cet artefact, je serais intéressé de connaître la cause et comment l'éviter!

Ceci est un gif animé, ouvert dans une fenêtre séparée pour voir l'animation (il ne jouera qu'une fois ou pas à la fois il a été mis en cache!)

Notez que la propagation semble très stable jusqu'à ce que l'onde commence à se refléter depuis la première frontière. Que pensez-vous qu'il pourrait se passer ici? J'ai passé quelques jours à vérifier mon code et je ne trouve aucune erreur. C'est étrange car il semble y avoir deux solutions de propagation: une positive et une négative; après la réflexion de la première frontière. Les solutions semblent se déplacer le long des points de maillage adjacents.

Les détails de l'implémentation suivent.

L'équation d'advection,

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

où est la vitesse de propagation.$\boldsymbol{v}$

Le Crank-Nicolson est une discrétisation stable inconditionnellement (lien pdf) pour l'équation d'advection à condition que varie lentement dans l'espace (ne contient que des composantes basses fréquences lors de la transformation de Fourier).$u(x)$

La discrétisation que j'ai appliquée est,

$\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]$

Mettre les inconnues sur le côté droit permet de l'écrire sous la forme linéaire,

$\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}$

où (pour prendre la moyenne temporelle uniformément pondérée entre le point présent et futur) et .r = v Δ $\beta=0.5$$r=\boldsymbol{v}\frac{\Delta t}{2\Delta x}$

Ces équations ont la forme matricielle , où,$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)$

Les vecteurs et sont connus et inconnus de la quantité que nous voulons résoudre.u n + $u^n$$u^{n+1}$

J'applique ensuite les conditions aux limites fermées de Neumann sur les limites gauche et droite. Par frontières fermées, j'entends sur les deux interfaces. Pour les frontières fermées, il s'avère que (je ne montrerai pas mon travail ici) nous avons juste besoin de résoudre l'équation de matrice ci-dessus. Comme l'a souligné @DavidKetcheson, les équations matricielles ci-dessus décrivent en fait les conditions aux limites de Dirichlet . Pour les conditions aux limites 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)$

Mise à jour

Le comportement semble assez indépendant du choix des constantes que j'utilise, mais ce sont les valeurs du graphique que vous voyez ci-dessus:

• $\boldsymbol{v}$ = 2
• dx = 0,2
• dt = 0,005
• $\sigma$ = 2 (hwhm gaussien)
• $\beta$ = 0,5

Mise à jour II

Une simulation avec un coefficient de diffusion non nul, (voir commentaires ci-dessous), l'oscillation disparaît, mais l'onde ne réfléchit plus!? Je ne comprends pas pourquoi?$D=1$

L'équation que vous résolvez ne permet pas de solutions à droite, il n'y a donc pas de condition aux limites réfléchissantes pour cette équation. Si vous considérez les caractéristiques, vous vous rendrez compte que vous ne pouvez imposer une condition aux limites qu'à la droite. Vous essayez d'imposer une condition aux limites de Dirichlet homogène à la limite de gauche, qui est mathématiquement invalide.

Pour réitérer: la méthode des caractéristiques dit que la solution doit être constante le long d' une ligne de la forme pour toute constante . Ainsi, la solution le long de la frontière gauche est déterminée par la solution à des moments antérieurs à l'intérieur de votre domaine de problème; vous ne pouvez pas y imposer de solution.$x-\nu t = C$$C$

Contrairement à l'équation, votre schéma numérique n'admet des solutions à droite en cours. Les modes à droite sont appelés modes parasites et impliquent des fréquences très élevées. Notez que l'onde de droite est un paquet d'ondes en dents de scie, associé aux fréquences les plus élevées qui peuvent être représentées sur votre grille. Cette onde est purement un artefact numérique, créé par votre discrétisation.

Pour mettre l'accent: vous n'avez pas noté le problème complet de la valeur de limite initiale que vous essayez de résoudre. Si vous le faites, il sera clair que ce n'est pas un problème mathématique bien posé.

Je suis heureux que vous ayez posté ceci ici, car c'est une belle illustration de ce qui peut se produire lorsque vous discrétisez un problème qui n'est pas bien posé et du phénomène des modes parasites. Un gros +1 pour votre question de ma part.

J'ai beaucoup appris des réponses ci-dessus. Je veux inclure cette réponse parce que je pense qu'elle offre des perspectives différentes sur le problème.

Considérons l'équation avec une vitesse d'onde constante .

Sans conditions initiales et aux limites, cette équation a une solution de la forme . (une impulsion se déplaçant vers la droite)$u(x,t)=f(x-ct)$

Si nous imposons une condition initiale , alors la solution de l'équation dans l'intervalle est . Il n'y a pas de frontière et c'est la solution.$u(x,t_0)=p(x)$$x \in (-\infty, \infty)$$p[x-c(t-t_0)]$

Supposons maintenant que nous définissions un domaine borné , après que tous les ordinateurs aient une mémoire limitée. Ensuite, nous devons spécifier les valeurs en et , sinon nous sommes bloqués sur le plan du calcul.$x \in [a,b]$$a$$b$

Une façon de définir les conditions aux limites consiste à utiliser Dirichlet sur la limite de gauche et une condition cohérente avec la solution en cours de propagation. Autrement dit, nous pouvons définir (supposer l'instant initial ) $t_0$

Cela produit une impulsion allant vers la droite jusqu'à ce qu'elle disparaisse sur le bord droit.

Poussez ici pour l'animation sur Dirichlet sur la limite gauche

Je reçois toujours du bruit que je ne peux pas comprendre (n'importe qui pourrait aider ici s'il vous plaît?)

L'autre option consiste à imposer des conditions aux limites périodiques. C'est au lieu d'imposer la condition aux limites de Dirichlet à gauche, nous pouvons imposer le paquet d'onde qui est cohérent avec la limite à gauche. C'est:

Cependant pour et , et comme nous devons mettre des données à l'intérieur de l'intervalle nous ajoutons une "période" de longueur , puis nous trouvons , de sorte que la condition à gauche serait (pareil!) et nous aurions une impulsion sortant sur à droite et en entrant à gauche.$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]$

Ce lien montre ce que j'appellerais des conditions aux limites périodiques.

J'ai fait les animations en python et le schéma est un schéma de Crank-Nicholson comme indiqué dans la question ici.

Je lutte toujours avec le modèle de bruit après que l'onde ait atteint la première limite (à droite).