Conditions aux limites pour l'équation d'advection discrétisée par une méthode aux différences finies

14

J'essaie de trouver des ressources pour aider à expliquer comment choisir les conditions aux limites lors de l'utilisation de méthodes aux différences finies pour résoudre les PDE.

Les livres et notes auxquels j'ai actuellement accès disent tous des choses similaires:

Les règles générales régissant la stabilité en présence de frontières sont beaucoup trop compliquées pour un texte introductif; ils nécessitent des machines mathématiques sophistiquées

(A. Iserles Un premier cours dans l'analyse numérique des équations différentielles)

Par exemple, lorsque vous essayez d'implémenter la méthode saute-mouton en 2 étapes pour l'équation d'advection:

uin+1=uin1+μ(ui+1nui1n)

en utilisant MATLAB

M = 100; N = 100;

mu = 0.5;

c = [mu 0 -mu];
f = @(x)(exp(-100*(x-0.5).^2));

u  = zeros (M, N);
x = 1/(M+1) * (1:M);

u(:,1) = f(x);
u(:,2) = f(x + mu/(M+1));

for i = 3:N
    hold off;
    u(:,i) = conv(u(:,i-1),c,'same') + u(:,i-2);
    plot(x, u(:,i));
    axis( [ 0 1 0 2] )
    drawnow;
end

La solution se comporte bien jusqu'à ce qu'elle atteigne la limite, lorsqu'elle commence très soudainement à se comporter mal.

Où puis-je apprendre à gérer des conditions aux limites comme celle-ci?

Simon Morris
la source

Réponses:

12

La réponse de Sloede est très complète et correcte. Je voulais juste ajouter quelques points pour le rendre plus facile à saisir.

Fondamentalement, toute équation d'onde a une vitesse et une direction inhérentes. Pour une équation d'onde unidimensionnelle: la vitesse de l'onde est la constante qui détermine non seulement la vitesse à laquelle l'information se propage dans le domaine mais également sa direction. Si , les informations vont de gauche à droite et si c'est l'inverse.a a > 0 a < 0

ut+aux=0
aa>0a<0

Pour la méthode de saut de grenouille, lorsque vous discrétisez les équations, vous obtenez: ou: où . Dans votre cas,

uinuin22Δt+aui+1n1ui1n12Δx=0
uin=uin2+μ(ui+1n1ui1n1)
μ=aΔt/Δxμ>0ce qui se traduit par une vague allant vers la gauche. Maintenant, si vous y réfléchissez, une onde qui se déplace vers la gauche n'aura besoin que d'une condition aux limites à la droite car toutes les valeurs à gauche sont mises à jour via leurs voisins de droite. En fait, la spécification d'une valeur à la limite gauche est incompatible avec la nature du problème. Dans certaines méthodes, comme au vent simple, cela est pris en charge automatiquement car le schéma n'implique également que les bons voisins dans son gabarit. Dans d'autres méthodes, comme la grenouille bondissante, vous devez spécifier une valeur "correcte".

Cela se fait généralement par extrapolation à partir du domaine interne pour trouver la valeur manquante. Pour les problèmes multidimensionnels et non canoniques, cela implique de trouver tous les vecteurs propres du flux jacobien pour déterminer quelles parties de la frontière ont réellement besoin de conditions aux limites et quelles parties nécessitent une extrapolation.

GradGuy
la source
Physiquement, que signifierait l'utilisation de cette équation avec une condition aux limites à gauche et à droite?
Frank
5

Réponse générale
Votre problème est que vous ne définissez pas (ou même ne spécifiez pas) les conditions aux limites - votre problème numérique est mal défini.

En règle générale, il existe deux manières possibles de spécifier les conditions aux limites:

  1. Définissez les conditions aux limites en spécifiant et externe, par exemple via la solution exacte.u0u101
  2. Modifiez le gabarit numérique afin qu'il n'utilise que des informations intérieures à la limite.

La voie à suivre dépend fortement de la physique de votre problème. Pour les problèmes de type équation d'onde, on détermine généralement les valeurs propres du flux jacobien afin de décider si des conditions aux limites externes sont nécessaires ou si la solution intérieure doit être utilisée (cette méthode est communément appelée `` remontée '').



ui1nui+1nin+1i=1u0nu100n+1u101n

u1nu100n

Vous pouvez trouver une version modifiée de votre code source ci-dessous:

M = 100; N = 100;

mu = 0.5;

c = [mu 0 -mu];
f = @(x)(exp(-100*(x-0.5).^2));

u  = zeros (M, N);
x = 1/(M+1) * (1:M);

u(:,1) = f(x);
u(:,2) = f(x + mu/(M+1));

for i = 3:N
    hold off;
    %u(:,i) = conv(u(:,i-1),c,'same') + u(:,i-2);

    % Apply the numerical stencil to all interior points
    for j = 2:M-1
        u(j,i) = u(j,i-2) + mu*(u(j+1,i-1) - u(j-1,i-1));
    end

    % Set the boundary values by interpolating linearly from the interior
    u(1,i) = 2*u(2,i) - u(3,i);
    u(M,i) = 2*u(M-1,i) - u(M-2,i);

    plot(x, u(:,i));
    axis( [ 0 1 0 2] )
    drawnow;
end
Michael Schlottke-Lakemper
la source
Belle réponse, et bienvenue à Scicomp, Sloede. Une question, normalement je vois "remontée" définie comme l'utilisation d'un pochoir unilatéral où l'information est tirée d'une seule limite du domaine. Vouliez-vous dire cela dans votre réponse?
Aron Ahmadia
1
Oui en effet. Désolé, si ma réponse n'était pas assez claire. En règle générale, cependant, «remontée» signifie simplement que vous tenez compte du sens du flux d'informations. Cela ne signifie pas que vous jetez complètement un côté de la solution, cela signifie simplement que vous donnez une préférence à la partie de la solution qui se trouve dans le sens "au près".
Michael Schlottke-Lakemper
Si vous créez N = 1000et exécutez le code un peu plus longtemps, vous constatez qu'il ne se comporte pas tout à fait comme prévu.
Simon Morris
La raison en est que ma solution "quick fix" n'est pas physiquement saine, et en plus de cela plutôt sensible aux oscillations parasites dans la solution. Ne l'utilisez pas pour des calculs scientifiques réels!
Michael Schlottke-Lakemper
2

J'ai donc examiné cela plus en détail et il semble que cela (au moins dans les cas de base que je gère) dépend de la vitesse de groupe de la méthode.

La méthode leapfrog (par exemple) est:

ujen+1=ujen-1+μ(uje+1n-uje-1n)

ukn=eje(ζkΔX+ω(ζ)nΔt)

e2jeωΔt=1+μejeωΔt(ejeζΔX-e-jeζΔX)

péché(ωΔt)=μpéché(ζΔX)

dωdζ=cos(ζΔx)1μ2sin2(ζΔx)[1,1]

Maintenant, nous devons trouver la vitesse de groupe des conditions aux limites:

u1n+2=u1n+μu2n+1

Nous pouvons calculer la vitesse du groupe limite comme suit:

2isin(ωΔt)=μeiζΔx

donc pour trouver des vitesses de groupe que les limites permettent, nous devons trouver:

ω=cζ

cos(ζΔx)=0,μsin(ζΔx)=2sin(ζcΔt)

ζ=π2Δxμ=2sin(cμπ2)c[1,1]μ


u0n+1=u1n[1,1]

J'ai encore pas mal de choses à lire à ce sujet avant de le comprendre complètement. Je pense que les mots clés que je recherche sont la théorie GKS.

Source pour toutes ces notes A Iserles Part III


Un calcul plus clair de ce que j'ai fait peut être trouvé ici: http://people.maths.ox.ac.uk/trefethen/publication/PDF/1983_7.pdf

Simon Morris
la source
-2

Les gars, je suis très nouveau sur ce site. Peut-être que ce n'est pas l'endroit pour demander, mais s'il vous plaît pardonnez-moi car je suis très nouveau ici :) J'ai un problème extrêmement similaire, la seule différence étant la fonction de démarrage qui, dans mon cas, est une onde cosinusoïdale. Mon code est le suivant: effacer tout; clc; ferme tout;

M = 1000; N = 2100;

mu = 0,5;

c = [mu 0 -mu]; f = @ (x) 1- cos (20 * pi * x-0,025). ^ 2; u = zéros (M, N); x = 0: (1 / M): 0,05; u (1: longueur (x), 1) = f (x); u (1: longueur (x), 2) = f (x - mu / (M)); x = espace lins (0,1, M);

pour i = 3: N retient;

% Apply the numerical stencil to all interior points
for j = 2:M-1
    u(j,i) = u(j,i-2) - mu*(u(j+1,i-1) - u(j-1,i-1));
end

% Set the boundary values by interpolating linearly from the interior
u(M,i) =  2*u(M-1,i-1) - u(M-2,i-1);

graphique (x, u (:, i)); axe ([0 1,5 -0,5 2]) tracé; % pause fin

Il y a déjà ce code ici, mais pour une raison quelconque, probablement liée à la vague cosinus, mon code échoue: / toute aide serait appréciée :) merci!

John
la source
2
Bienvenue sur SciComp.SE! Vous devriez en faire une nouvelle question. (Les réponses ne sont destinées qu'aux réponses réelles.) Si vous utilisez le lien "posez votre propre question" en bas (il est jaune foncé sur jaune clair, certes un peu difficile à voir si vous ne le savez pas) , il liera automatiquement la question à celle-ci. (Vous pouvez également inclure un lien vers cette question dans la vôtre.)
Christian Clason