Discussion :

[Heimdall]

Bonjour,

J'aimerais comprendre prendre en compte les conditions limites (Surtout Neumann homogène et inhomogènes) avec une méthode des relaxation type Gauss-Seidel pour la résolution d'un problème elliptique ?

Soit une grille 2D G(i,j) avec i = 0, 1,.., nx et j=0, 1, ..., ny
Disons que j'ai une condition limite dirichlet f(i,j) sur mes bords (i=0,nx; j=0,ny). Je me donne une grille avec f(i,j) sur mes bords et 0 ou un autre initial guess, et je relaxe *dans* la grille jusqu'à avoir une convergence satisfaisante. Les bords "influencent" mon calcul et me font converger vers ma solution.

Est-ce la bonne façon de faire ? J'ai pas vraiment l'impression...

1/ Comment faire pour une condition de dirichlet homogène (Solution=0 aux bords) pour résoudre l'équation de laplace par exemple ? Avec quoi initier la relaxation ?

2/ Comment faire pour assurer une condition de Neumann ? homogène ? inhomogène ?


Merci, j'espère être clair.

[FR119492]

Salut!
Si je te comprends bien, tu cherches à résoudre une équation elliptique sur un domaine rectangulaire par la méthode desdifférences finies. Ton équation peut être homogène (Laplace) ou non homogène (Poisson), ça ne change pas grand chose.
Je suppose que ta grille est à pas constant h. Sur chaque noeud intérieur, tu as une équation de la forme:
4*f(i,j)-f(i-1,j)-f(i+1,j)-f(i,j-1)-f(i,j+1)=h^2*Q(i,j)
Tu remarques que j'ai changé tous les signes; c'est pour que la matrice soit définie positive. Le terme Q(i,j) est nul dans le cas de l'équation de Laplace.

Regardons maintenant les conditions aux limites:

Dans le cas d'une condition de Dirichlet (homogène ou non), la valeur de la fonction est connue; ce n'est donc pas une inconnue (La Palice dixit!). Si tu connais la valeur F(0,j), alors l'équation en (1,j) devient:
4*f(1,j)-f(2,j)-f(1,j-1)-f(1,j+1)=h^2*Q(1,j)+F(0,j)

Dans le cas d'une condition de Neumann homogène, la valeur de la fonction est inconnue; c'est donc une inconnue (La Palice dixit!). Mais tout se passe comme si ta frontière était un plan de symétrie; on aurait donc f(-1,j)=f(1,j). Ton équation sur la frontière devient:
4*f(0,j)-f(1,j)-f(1,j)-f(0,j-1)-f(0,j+1)=h^2*Q(i,j)
ou en simplifiant
4*f(0,j)-2*f(1,j)-f(0,j-1)-f(0,j+1)=h^2*Q(i,j)

Le cas d'une condition de Neumann non homogène est un peu plus compliqué, mais on ne le rencontre pratiquement jamais.

Remarque finale: Tu constateras que je n'ai pas parlé de la méthode de résolution, car il faut toujours séparer très distinctement la modélisation du problème et sa résolution numérique.

Jean-Marc Blanc

[Heimdall]

salut,

Ok merci de ta réponse, c'est à peu près ce que j'avais pensé pour diriclet homogène ou non, en fixant la valeur sur la frontière puis en ne travaillant qu'à l'intérieur du domaine ça roule.

A ceci près que mon équation n'est pas laplace ou poisson mais (toujours la même) :

en 1D : u (x) - d^2 d^2u/dx^2 = f(x)
en 2D = u(x,y) - d^2*laplacien(u) = f(x,y)

(avec d une constante)

j'ai d'abord essayé en 1D avec l'équation :

u(x) - u''(x) = cos(ax) avec a=cte

soit :

u(i) - 1/dx^2*(u(i+1) - 2u(i) + u(i-1)) = cos(a*i*dx)


Comment appliquer une méthode de relaxation gauss-seidel pour cette équation, en se donnant par exemple u(0)=u(X) = 0.0 ?

[FR119492]

Alors ton équation devient simplement quelque-chose comme
(1+4*D)*f(i,j)-D*f(i-1,j)-D*f(i+1,j)-D*f(i,j-1)-D*f(i,j+1)=h^2*Q(i,j)
A tout hasard, vérifie les signes.

Jean-Marc Blanc

PS: A part ça, j'éprouve toujours les mêmes doûtes concernant la validité de ton équation aux dérivées partielles.

[Heimdall]

Alors ton équation devient simplement quelque-chose comme
(1+4*D)*f(i,j)-D*f(i-1,j)-D*f(i+1,j)-D*f(i,j-1)-D*f(i,j+1)=h^2*Q(i,j)
A tout hasard, vérifie les signes.

Ok merci, j'avais bien compris comment discrétiser en 2D avec dx=dy. Ma question portait en fait sur comment résoudre ça par relaxation.

J'ai lu dans numerical recipes l'exemple (chapitre 19 sur la relaxation) assez physique. Soit une équation Lu=f, avec u la solution, f le terme source et L l'opérateur elliptique. On construit une équation de diffusion :

du/dt = Lu-f

En discretisant le tout, on obtient u(i)_(k+1) en fonction de u(i)_k, L(i) et f(i), où 'k' est l'indice temporel, c'est à dire l'indice d'iteration de la relaxation. Ils font ça avec une équation de poisson, j'ai fait la même manip avec mon équation, en 1D pour commencer :

(u(i)_(k+1) - u(i)_k)/dt = u(i) - D/dx^2 * (u(i+1) - 2u(i) + u(i-1)) - f(i)

En prenant le plus grand pas de temps possible dt = dx^2/2, on obtient :

(u(i)_(k+1) - u(i)_k) = u(i)_k*dx^2/2 - D/2 (u(i+1)_k - 2u(i)_k + u(i-1)_k) - f(i)dx^2/2

soit :

u(i)_(k+1) = 2*D*u(i)_k + u(i)_k*dx^2/2 - D/2 (u(i+1)_k + u(i-1)_k) - f(i)dx^2/2

En factorisant par u(i)_k le membre de droite :

u(i)_(k+1) = u(i)_k*(1 + D+ dx^2/2) - D/2 (u(i+1)_k + u(i-1)_k) - f(i)dx^2/2


Est-ce correct ?

voici le code que j'ai écrit pour résoudre cette équation 1D :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
 
 
 
int main(void)
{
 
int i, i1, i2, nx;
double x, dx ,dx2,a, xm, D;
int nsteps, istep;
 
FILE *fp;
 
xm  =  1.0;
nx  =  64;
D   = 1.0;
dx  =  xm/(double)nx;
 
a = 1.0;
nsteps = 1000;
istep  = 0; 
 
double *u, *r, *f;
 
u = calloc((si.nx+1), sizeof(*u));
f = calloc((si.nx+1), sizeof(*f));
 
 
/*Dirichlet boundary conditions*/
x     = 0;
u[0]  = 0;
x     = dx*(si.nx+1);
u[si.nx] = 0;
 
 
/*source term*/
for (i = 0; i < si.nx+1; i++)
    {
    x = i*dx;
    f[i] = sin(a*2*M_PI*x/si.xm); 
    }
 
 
/*Relaxation*/
dx2 = dx*dx;
while (istep ++ < nsteps) 
      {
      for (i = 1; i < si.nx; i++)
          {
          i1  = i-1;
          i2  = i+1;
          u[i] =  u[i]*(1+D+dx2*0.5) - 0.5*D*(u[i2] + u[i1])  - f[i]*dx2*0.5;
          }
      };
 
 
 
printf("%d iterations effectuees\n",istep-1);
 
fp = fopen("solution.txt", "w");
for (i = 0; i < si.nx+1; i++)
    {
    x = i*dx;
    fprintf(fp,"%lg %lg\n",x,u[i]); 
    }
fclose(fp);
 
 
 
 
free(u);
free(f);
 
return 0;
}

[Heimdall]

ne pas tenir compte de mon post précédent, j'ai été induit en erreur par une page de numerical recipes... ce que je fais (et qui marche) :

u(i) = (u(i-1)+u(i+1) - f(i)*dx^2)/(dx^2 + 2)

[Heimdall]

Le cas d'une condition de Neumann non homogène est un peu plus compliqué, mais on ne le rencontre pratiquement jamais.

Pourrais-tu m'expliquer ? Il se pourrait que j'ai à m'en servir.

Merci

[FR119492]

Salut!

Sur la frontière:

• Dirichlet homogène: la fonction est nulle;
• Dirichlet non homogène: la fonction est donnée non nulle;
• Neumann homogène: la dérivée normale de la fonction est nulle;
• Neumann non homogène: la dérivée normale de la fonction est donnée non nulle

Jean-Marc Blanc

[Ehouarn]

Bonjour,

Pour une condition de Neumann, je recommanderai plutôt d'utiliser un stencil décentré (sur trois points, pour rester d'ordre 2) pour établir la relation à utiliser pour calculer la valeur sur la frontière.
Par exemple, si la frontière est en x_0, pour des points équidistants espacés de dx, la dérivée première x'_0 est telle que:
x'_0=(-3 x_0 + 4 x_1 - x_2)/(2dx)
ce qui permet d'écrire x_0 en fonction de x_1, x_2 et de la condition au bord x'_0, sans avoir à introduire de point fictif, ni avoir à distinguer entre les cas x'_0 nul ou non-nul.

Ehouarn

PS: Le cas d'une condition de Neumann non nulle n'a rien de rare. On rencontre fréquemment des situations physiques où un flux est imposé à la frontière du domaine considéré.

[Heimdall]

Salut,


Pourrais-tu m'expliquer pourquoi tu me conseilles plutot un schéma décentré ? Il se trouve que j'utilise actuellement un schéma centré d'ordre deux (donc pour u'(0) je fais u(-1) = u(1)) et j'ai beaucoup d'imprécision sur ma condition de Neumann. Quelle en est la cause ? Cela rejoint-il ta remarque ?

exemple sur le bord droit : http://nico.aunai.free.fr/by.png
SOlution analytique en bleu, solution numérique en point rouge. (ne pas tenir compte de la droite bleue, mon plot est une projection dans le plan (y,By(x,y)) cette droite relie juste le dernier point de cette ligne x=cte de la ligne x=cte+1.

On voit clairement a l'oeil nu que mon dernier point n'est pas exactement sur ma solution analytique.

Voici un graphe montrant l'erreur "analytique-numérique"

http://nico.aunai.free.fr/ErrBy.png

on voit clairement un probleme sur le bord droit avec une erreur de 2.5e-5 environ au max. Ne pas croire qu'il n'y a pas de probleme au bord gauche, c'est juste que mon terme source est naturellement à dérivée nulle sur ce bord, la condition est donc je pense facilement respectée. J'ai d'ailleurs testé en décalant le terme source afin de lui virer sa tangente horizontale a gauche, et j'ai le meme type d'erreur. Il semblerait donc que je n'arrive pas a satisfaire très précisément la condition de neumann avec les différences finies.

Autre probleme, j'ai également des problèmes aux bords en ce qui concerne mes conditions dirichlet. Ok les points sur la frontière sont parfaitement sur la courbe vu que je les y met "à la main", mais l'erreur 'analytique-numérique' augmente légèrement à l'intérieur.

voici le graphe solution numérique (rouge) et solution analytique (vert) :
http://nico.aunai.free.fr/bx.png

et voici l'erreur (lun moins l'autre) :

http://nico.aunai.free.fr/errBx.png


On remarque directement que l'erreur grossi aux endroits ou la dérivée est grande (pente abrupte) pour la solution. Il y a forcément un lien mais comment gérer ce probleme sans augmenter la résolution en augmentant le nombre de points à mort ? Est-ce lié a la discrétisation des différences finies ?


PS : Notez que je ne comprends pas comment ceci est possible car sur le graphe on voit une erreur de 2e-5 au maximum alors que mon solver me dit :

cycle No 0, err = 34.556
cycle No 1, err = 1.24569
cycle No 2, err = 0.0360623
cycle No 3, err = 0.00455891
cycle No 4, err = 0.0031296
cycle No 5, err = 0.000371798
cycle No 6, err = 0.000202667
cycle No 7, err = 1.66485e-07
cycle No 8, err = 1.51832e-05
cycle No 9, err = 1.4456e-07
cycle No 10, err = 1.14143e-06
cycle No 11, err = 3.28487e-10
cycle No 12, err = 8.61014e-08
cycle No 13, err = 2.53432e-11
cycle No 14, err = 6.50848e-09
cycle No 15, err = 5.10001e-13
converge au cycle No 15, err = 5.10001e-13

Normalement ce sont les même grandeurs, je m'explique. Sur le graphe on a la soustraction de ma solution numérique par la solution analytique, l'erreur quoi. Ce que mon solveur affiche, c'est la correction que l'ont fait à la fin de chaque cycle V multigrille : Solution = Solution + erreur, ou 'erreur' est trouvée par la relaxation sur les multiples grilles. Normalement ce sont théoriquement les deux mêmes quantités. Pourtant on solveur me dit que le max de la correction (le max de l'erreur que je fais donc) est 5e-13, or quand je compare effectivement j'ai du 2e-5. Comment l'expliquer ?

[Ehouarn]

Bonjour,

Je conseille le schéma décentré parce qu'il permet de gérer n'importe quelle condition de Neumann alors que la méthode du point fictif ne permet de résoudre que le cas d'une dérivée nulle.
On pourrait effectivement s'amuser à regarder laquelle des deux méthodes donne un meilleur résultat, mais cela ne serait valable que sur ton cas test et ne permettrai pas de généraliser. De même en ce qui concerne les écarts entre solutions analytique et numérique que tu observes, les lieux ou tu obtiens les plus grands écarts dépendent évidement très fortement du problème particulier qui est considéré.
On peut par contre s'attendre à ce que les plus fortes erreurs se trouvent là ou la solution analytique varie le plus dans la mesure où les différences finies d'ordre 2 sont en gros une représentation (en chaque point) d'un développement de Taylor à l'ordre 2 de la solution.

Si tu souhaites augmenter la précision sans augmenter le nombre de points (ni changer de méthode), je ne vois pas d'autre solution que de passer à des différences finies d'ordre plus élevé. En passant simplement à l'ordre 4 tu devrais obtenir une sérieuse amélioration.

Pour ce qui est du fait d'obtenir une erreur (entre solution analytique et numérique) plus grande que les "erreurs" de relaxation d'une grille à la suivante, je pense que tu confonds deux choses: il y a d'une part les relaxations d'une grille à l'autre concernent la convergence vers la solution numérique d'une grille plus grossière vers la grille finale (la plus fine), et d'autre part l'erreur (c.-à-d. l'écart) entre cette solution numérique et la solution analytique.

Ehouarn

[Heimdall]

Salut, merci de ta réponse !

Bonjour,

Je conseille le schéma décentré parce qu'il permet de gérer n'importe quelle condition de Neumann alors que la méthode du point fictif ne permet de résoudre que le cas d'une dérivée nulle.

u(-1) - u(1) = 2dx f'(0)

Je peux pa remplacer dans mon équation discrétisée en 0, u(-1) par u(1)+2dxf'(0) ?

On pourrait effectivement s'amuser à regarder laquelle des deux méthodes donne un meilleur résultat, mais cela ne serait valable que sur ton cas test et ne permettrai pas de généraliser. De même en ce qui concerne les écarts entre solutions analytique et numérique que tu observes, les lieux ou tu obtiens les plus grands écarts dépendent évidement très fortement du problème particulier qui est considéré.
On peut par contre s'attendre à ce que les plus fortes erreurs se trouvent là ou la solution analytique varie le plus dans la mesure où les différences finies d'ordre 2 sont en gros une représentation (en chaque point) d'un développement de Taylor à l'ordre 2 de la solution.

Oui c'est bien ce qu'on voit dans le cas bx que j'ai montré. On voit clairement que l'erreur se produit près des deux bords là ou la dérivée de la solution est immense pour satisfaire dirichlet homogène. Ca me semble bien logique aussi.

Si tu souhaites augmenter la précision sans augmenter le nombre de points (ni changer de méthode), je ne vois pas d'autre solution que de passer à des différences finies d'ordre plus élevé. En passant simplement à l'ordre 4 tu devrais obtenir une sérieuse amélioration.

Ok why not

Pour ce qui est du fait d'obtenir une erreur (entre solution analytique et numérique) plus grande que les "erreurs" de relaxation d'une grille à la suivante, je pense que tu confonds deux choses: il y a d'une part les relaxations d'une grille à l'autre concernent la convergence vers la solution numérique d'une grille plus grossière vers la grille finale (la plus fine), et d'autre part l'erreur (c.-à-d. l'écart) entre cette solution numérique et la solution analytique.

Hum, je suis toujours pas sûr de comprendre alors.

Soit L un opérateur linéaire, x la solution de Lx=b.

Lorsque j'utilise une méthode itérative comme les multigrilles, j'ai à l'itération k, le vecteur x_k, approchant ma solution x.

J'ai donc L.x_k = b + r_k

Avec r_k le résidu de mon équation à l'itération k.

Par linéarité je peux écrire L(x_k - x) = r_k, autrement dit, l'erreur 'e_k' que je fais à l'itération k satisfait L.e_k = r_k.

Cette erreur il s'agit bien de x_k-x, autrement dit, elle quantifie la "proximité" de ma solution analytique à la vraie solution x.

Bien, maintenant, en multigrille, on smooth quelque fois sur la grille fine l'équation L. x_k = b. Ayant relaxé que quelques itérations, on a un gros résidu sur cette équation, r_k. Ce résidu est Restreint sur une grille plus grossière et devient le second membre de l'équation qu'on va relaxer là bas. Cette équation est donc L.e_k = r_k. Restons sur la méthode "two grid", j'ai disons une solution exacte pour e_k sur cette grille plus grossière, je la Prolonge sur ma grille fine, et je corrige x_k :

x_k+1 = x_k + e_k

On peut donc voir e_k comme étant la "distance" séparant le cycle V (k+1) du cycle V(k), mais par définition de e, c'est aussi et surtout l'erreur faite sur e_k par rapport à la solution exacte non ?

Si maintenant à chaque fin de cycle V, je regarde le Max(e_k), j'ai donc par définition de e_k, |X-X_k| < 1e-13, autrement dit, ma solution numérique est "proche" de la solution exacte à au plus 1e-13. Si par chance j'ai la solution analytique, donc exacte de mon équation. Pourquoi je ne pourrais pas comparer directement X-X_k "à la main" en faisant une simple soustraction et obtenir une valeur maximale de cette différence de 1e-13 ?

L'erreur e_k dont je parle n'est pas comme tu dis "l'erreur d'un grille grossière vers la grille plus fine", c'est l'erreur SUR la grille fine, que l'on a trouvé par relaxation/corrections successives dans le cycle V.

Où est-ce que j'ai faux ?

[Heimdall]

Je n'arrive pas à trouver la discrétisation du laplacian à l'ordre 4 pour dx différent de dy. J'ai bien trouvé l'opérateur 9 points du laplacien, mais pour dx=dy=h...

Quelqu'un aurait-il ça en main ?

(Si en plus vous savez comment dériver l'expression, j'avoue ne même pas trop savoir comment faire pour dériver l'opérateur 9 points pour dx=dy=h)


Merci.

[Ehouarn]

Salut,

u(-1) - u(1) = 2dx f'(0)
Je peux pa remplacer dans mon équation discrétisée en 0, u(-1) par u(1)+2dxf'(0) ?

Oui, tu as raison, on peux effectivement faire cela. J'aurais tendance à penser que ce sera un peux moins précis puisqu'on utiliserai alors moins d'informations (seulement u(1) et u'(0)) qu'avec un stencil décentré (où on utilise u'(0), u(1) et u(2)). Mais il faudrait une analyse plus poussée pour le savoir.

Pour ce qui est des "erreurs":

Par linéarité je peux écrire L(x_k - x) = r_k, autrement dit, l'erreur 'e_k' que je fais à l'itération k satisfait L.e_k = r_k.
Cette erreur il s'agit bien de x_k-x, autrement dit, elle quantifie la "proximité" de ma solution analytique à la vraie solution x.

Euh, non, x est ta solution numérique (obtenue sur un maillage donné avec une discrétisation donnée) et ce 'x' tends vers la solution analytique lorsque par exemple tu raffine ton maillage et/ou ta discrétisation. Lorsque le résidu e_k est (presque) nul, c'est simplement que ton processus itératif a convergé et que donc x_k=x, à e_k (très petit) près. Cet e_k n'a rien à voir avec la solution analytique.

L'expression du laplacien sur 9 points s'obtient juste en combinant les expressions suivant chaque dimension (sur 5 points pour être à l'ordre 4) que l'on retrouve soit en utilisant des développements de Taylor, soit en considérant le polynôme d'interpolation de Lagrange passant par ces mêmes points. J'ai quelques infos et résultats là dessus ici dans les parties "dérivation numérique" et "équations différentielles".

[Heimdall]

Euh, non, x est ta solution numérique (obtenue sur un maillage donné avec une discrétisation donnée) et ce 'x' tends vers la solution analytique lorsque par exemple tu raffine ton maillage et/ou ta discrétisation. Lorsque le résidu e_k est (presque) nul, c'est simplement que ton processus itératif a convergé et que donc x_k=x, à e_k (très petit) près. Cet e_k n'a rien à voir avec la solution analytique.

ok, très bien merci. Donc typiquement quand je dis que je vois ma condition de Neumann pas bien respectée sur le bord en comparant avec la solution analytique je dis :

Mon algorithme a convergé vers la solution numérique la plus précise posible AVEC cette discrétisation et ce maillage.

L'expression du laplacien sur 9 points s'obtient juste en combinant les expressions suivant chaque dimension (sur 5 points pour être à l'ordre 4) que l'on retrouve soit en utilisant des développements de Taylor, soit en considérant le polynôme d'interpolation de Lagrange passant par ces mêmes points. J'ai quelques infos et résultats là dessus ici dans les parties "dérivation numérique" et "équations différentielles".

ok

en fait en réfléchissant j'ai l'impression que je ne vais pas résoudre mon probleme en passant a l'ordre 4. regarde : http://nico.aunai.free.fr/logdivB.png

c'est la divergence (en log décimal). Au milieu j'ai bien 1e-14 comme div(S), par contre sur les bords Y=0 et Y=Ly ça montre quasi à 1e0 ! c'est pas en passant à l'ordre 4 qua je vais gagner 14 ordres de grandeur sur ma divergence.

Je suis a cours d'idée. Que pourriez-vous me proposer ?

Pensez vous déjà que ma condition limite est adaptée (de manière théorique) ?


Merci bcp

[Ehouarn]

Salut,

Au milieu j'ai bien 1e-14 comme div(S),

J'en déduis que tu parles là du système avec grilles décalées mentionné dans ce fil et effectivement une si mauvaise divergence n'est très certainement pas simplement une affaire de précision du solver. Mais avant de se pencher plus avant sur ce cas précis, je te proposerais plutôt de faire un test sur un cas simple 2D où tu résouts:
u(x,y) - d^2*laplacien(u) = f(x,y)
avec un f(x,y) choisi tel que tu ais une solution u(x,y) que tu connaisse analytiquement et avec les conditions aux limites de ton choix.
Comme cela tu pourras déterminer si la divergence non nulle que tu observes est due au solver où à la formulation du problème complet que tu souhaites résoudre.

Bon courage.

[Heimdall]

Salut,


J'en déduis que tu parles là du système avec grilles décalées mentionné dans ce fil et effectivement une si mauvaise divergence n'est très certainement pas simplement une affaire de précision du solver.

Mais avant de se pencher plus avant sur ce cas précis, je te proposerais plutôt de faire un test sur un cas simple 2D où tu résouts:
u(x,y) - d^2*laplacien(u) = f(x,y)
avec un f(x,y) choisi tel que tu ais une solution u(x,y) que tu connaisse analytiquement et avec les conditions aux limites de ton choix.
Comme cela tu pourras déterminer si la divergence non nulle que tu observes est due au solver où à la formulation du problème complet que tu souhaites résoudre.

Salut, j'ai déjà tenté quelque chose dans ce style, en rentrant un terme source f(x,y) à divergence nulle (numériquement 1e-7), la solution u(x,y) que j'ai trouvé avait également un probleme de divergence au bord y=cte.


Bon plutot qu'un long discours, je vais poster le code de ma fonction smooth() en expliquant la discrétisation, ça se trouve c'est ma discrétisation qui n'est pas bonne (je veux dire que j'introduit peut-etre sans le savoir une erreur sur ma condition de Neumann):

J'ai donc l'équation (E) :


u(i,j) - mu*/dx^2*(u(i-1,j)-2u(i,j) + u(i+1,j)) - mu/dy^2*(u(i,j-1)-2u(i,j) + u(i,j-1)) = f(i,j)

tout est périodique en i=0 et i=nx
ux = 1 en j = ny et ux=-1 en j=0
duy/dy = 0 en j=0 et j=ny
uz = 0 en j=0 et j=ny

J'écris l'équation (E) sous la forme :

u(i,j) = tout le reste


et je la résoud pour ux,uy et uz.
Techniquement, pour ux, je fais ça dans une boucle allant de i=0 à nx et j=1 à ny-1, car la ligne du bas et celle du haut sont les conditions de dirichlet qu'on ne touche pas. Idem pour uz.

Pour uy, je fais la boucle de i=0 à i=nx et j=0 à j=ny.

en j=0 et j=ny, les points en j-1 et j+1 respectivements sont hors du domaine, je me sert donc de la formule d'ordre deux suvante :

(u(i,ny+1) - u(i,ny-1))/2dy^2 = u'(i,ny) = 0
et
(u(i,1) - u(i,-1))/2dy^2 = u'(i,0) = 0


Dans les faits, je numérote les 4 points de mon stencil dans l'ordre suivant :

ij1 : (i-1), j
ij2 : i,j+1
ij3 : i+1, j
ij4 : i, j-1

pour j=0 et pour uy je fais donc

ij4 = ij2 de sorte a avoir uy[ij4] == uy[ij2]

pour j = ny je fais donc :

ij2 = ij4, de sorte à avoir uy[ij2] == uy[ij4]


ce qui me permet de n'avoir qu'une seule fois la relation

u(i,j) = ...


voici le code :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
 
typedef struct _mgs
{
double *bx, *by, *bz;
double *bxc, *byc, *bzc;
double *sx, *sy, *sz;
double  *rx, *ry, *rz;
int nx, ny;
double dx, dy;
double mu;
int k;
}MGS;
 
 
/**
EQUATION : 
 
S = B - MU*LAPLACIAN(B)
 
**/
 
 
void mgs_smooth(MGS *mgs, int m)
{
int i, j, ij, ij1, ij2, ij3, ij4;
double dx2,dy2;
int istep, nsteps;
 
nsteps = m;
istep = 0;
 
dx2 = mgs->dx*mgs->dx;
dy2 = mgs->dy*mgs->dy;
 
while (istep ++ < nsteps)
      {
      for (i = 0; i < mgs->nx+1; i++) /*PERIODIC (X) - NEUMANN (Y)*/
	      {
		  for (j = 0; j < mgs->ny+1; j++)
		      {
		      ij  = i    +  j*(mgs->nx+1);	
		      ij1 = i-1  + (j  )*(mgs->nx+1);
		      ij2 = i    + (j+1)*(mgs->nx+1);
		      ij3 = i+1  + (j  )*(mgs->nx+1);
		      ij4 = i    + (j-1)*(mgs->nx+1);
 
                      if (i == 0)
                         ij1 = mgs->nx-1 + j*(mgs->nx+1);
                      else if (i == mgs->nx)
                         ij3 = 1 + j*(mgs->nx+1);
 
                      if (j == 0)
                         ij4 = ij2;
                      else if (j == mgs->ny)
                         ij2 = ij4;
 
		      mgs->by[ij] = ((dx2*dy2)*mgs->sy[ij] + mgs->mu*dy2*(mgs->by[ij1] + mgs->by[ij3]) 
		                   + mgs->mu*dx2*(mgs->by[ij2] + mgs->by[ij4]))/(dx2*dy2 + 2*mgs->mu*(dx2+dy2));                   
		      }	
		  }   
	  /* __ Periodic(X) - Dirichlet(Y) __ */	
      for (i = 0; i < mgs->nx+1; i++)
		  {
		  for (j = 1; j < mgs->ny; j++)
		      {
		      ij  = i    +  j*(mgs->nx+1);	
		      ij1 = i-1  + (j  )*(mgs->nx+1);
		      ij2 = i    + (j+1)*(mgs->nx+1);
		      ij3 = i+1  + (j  )*(mgs->nx+1);
		      ij4 = i    + (j-1)*(mgs->nx+1);
 
		      if (i == 0)
		         ij1 = mgs->nx-1 + j*(mgs->nx+1);
		      else if (i == mgs->nx)
		         ij3 = 1 + j*(mgs->nx+1);
 
                      mgs->bx[ij] = ((dx2*dy2)*mgs->sx[ij] + mgs->mu*dy2*(mgs->bx[ij1] + mgs->bx[ij3]) 
		                   + mgs->mu*dx2*(mgs->bx[ij2] + mgs->bx[ij4]))/(dx2*dy2 + 2*mgs->mu*(dx2+dy2));
 
		      mgs->bz[ij] = ((dx2*dy2)*mgs->sz[ij] + mgs->mu*dy2*(mgs->bz[ij1] + mgs->bz[ij3])
		                  + mgs->mu*dx2*(mgs->bz[ij2] + mgs->bz[ij4]))/(dx2*dy2 + 2*mgs->mu*(dx2+dy2));		      
		      }	
		  }   			  
       };  
}

[Ehouarn]

Salut,

Salut, j'ai déjà tenté quelque chose dans ce style, en rentrant un terme source f(x,y) à divergence nulle (numériquement 1e-7), la solution u(x,y) que j'ai trouvé avait également un probleme de divergence au bord y=cte.

Ce qui est effectivement très inquiétant. J'ai quand même un doute sur ton cas test, est-ce qu'il suffit d'avoir f à divergence nulle pour que du coup u (tel que u-mu.laplacien(u)=f) le soit aussi?
Au cas où, tu pourrais tenter un cas test ou tu choisis u, à divergence nulle, puis tu en déduit le terme source f qu'il faudrait pour récupérer ton équation différentielle. Comme cela tu seras sûr de tenir un cas test où il faudrait obtenir une solution u à divergence nulle.

en j=0 et j=ny, les points en j-1 et j+1 respectivements sont hors du domaine, je me sert donc de la formule d'ordre deux suvante :
(u(i,ny+1) - u(i,ny-1))/2dy^2 = u'(i,ny) = 0
et
(u(i,1) - u(i,-1))/2dy^2 = u'(i,0) = 0

C'est plutôt
(u(i,ny+1) - u(i,ny-1))/2dy = u'(i,ny) = 0 et (u(i,1) - u(i,-1))/2dy = u'(i,0) = 0
mais sans conséquences puisqu'au final tu te ramène à u(i,-1)=u(i,1) et u(i,ny+1)=u(i,ny-1).

Je n'ai pas regardé ton code en détail, mais ce que tu présente là est juste une itération à la Gauss-Seidel de la résolution, n'est-ce pas? Et dans ce cas la méthode que tu présentes a l'air tout à fait correcte.

[Heimdall]

Salut,


Ce qui est effectivement très inquiétant. J'ai quand même un doute sur ton cas test, est-ce qu'il suffit d'avoir f à divergence nulle pour que du coup u (tel que u-mu.laplacien(u)=f) le soit aussi?
Au cas où, tu pourrais tenter un cas test ou tu choisis u, à divergence nulle, puis tu en déduit le terme source f qu'il faudrait pour récupérer ton équation différentielle. Comme cela tu seras sûr de tenir un cas test où il faudrait obtenir une solution u à divergence nulle.

Ok, supposons que le champ magnétique initial de mon code (B) soit à divergence nulle (!!) Que dire pour div(S) ?

DivS = divB - mu*div(laplacien(B))

par hypothèse le premier term est nul donc

divS = -mu*div(laplacien(B))
= -mu*div(rot(rot(B)) - grad(div(B))

Le premier terme est toujours nul (div(rot(f)) = 0
Le second est nul par hypothèse.

On a donc divB=0 => divS=0.


Au cours de ma simulation, j'update S dans le code par dS/dt = rot(f) avec f un autre champ vectoriel. donc numériquement si S est a divergence nulle au début de la simulation, il l'est pour toujours, du moment que mon schéma numérique n'introduit pas de la divergence.

Ok supposons maintenant divS = 0, je dois inverser mon équation elliptique, qu'advient-il de la divergence de B ?


divS = div(B) -mu*div(laplacien(B))
= div(B) -mu*div(rot(rot(B)) - grad(div(B))

comme tout a l'heure, div(rot)=0 donc il reste :

divS = 0 = divB + mu div(grad(div(B)))

Appelons f le scalaire tel que f = div(B), alors :

f + mu*grad(div(f)) = 0

Equation différentielle du second ordre, linéaire en f, c'est à dire en div(B).

Autrement dit, je pense qu'en imposant dans mon équation de départ, div(B)=0 sur mon contour, alors div(B) sera nul partout car la seule variation linéaire spatiale qu'il peut avoir entre 0 et 0, c'est rester à 0 !

C'est plutôt
(u(i,ny+1) - u(i,ny-1))/2dy = u'(i,ny) = 0 et (u(i,1) - u(i,-1))/2dy = u'(i,0) = 0
mais sans conséquences puisqu'au final tu te ramène à u(i,-1)=u(i,1) et u(i,ny+1)=u(i,ny-1).

En effet, faute de frappe.

Je n'ai pas regardé ton code en détail, mais ce que tu présente là est juste une itération à la Gauss-Seidel de la résolution, n'est-ce pas? Et dans ce cas la méthode que tu présentes a l'air tout à fait correcte.

Oui c'est juste 'm' itérations de gauss-seidel, la routine est "normalement" destinée à être utilisée en pré/post smoothing pour les multigrilles. Mais là je l'utilise seule en regardant le résidu de l'équation, qui diminue jusqu'à 1e-15.


Alors je suis pas fort en différences finies, je débute, mais hier je me suis renseigné et j'ai tenté de calculer l'erreur de troncation, qui est en dx^2 et dy^2 en fait c'est quasi la même que celle du laplacien bien sûr. Idem pour ma condition limite c'est du second ordre, tout ça est consistant.

Les valeurs propres de mon opérateur (I - mu*Laplacien) sont 1 moins celles du laplacien en gros 1- mu/h^2(2cos(A) - 2), la valeur propre de l'inverse de l'opérateur symétrique étant l'inverse de la valeur propre de l'opérateur, je regarde la plus petite, qui sera celle qui va limiter la convergence. Elle est bien une constante indépendante de h, donc le schéma est stable. Stable et consistant on a donc convergent à l'ordre 2.

Ca c'est pour les maths... intuitivement je me dis que je n'ai pas beaucoup de points pour ma dérivée nulle aux bords, je dis juste que le point fantome = le point à l'intérieur... je trouve ça peut. Est-ce que ça pourrait être utile d'introduire une dérivée décentrée du style :

a*u(-1) + b*u(0) + c*u(1) + d*u(2) ?

Intuitivement je dirait que c'est plus acceptable car on peut exprimer le point fantome dans mon laplacien en fonction de plus de points à l'intérieur, et surtout aussi en fonction du point u(0) au bord, qu'il suffit donc de passer dans le membre de gauche dans gauss seidel.

Mathématiquement je vois pas pourquoi je devrais faire quelque chose à l'ordre 3 et pourquoi l'ordre 2 n'assure pas....

[Heimdall]

J'ajoute une question :

que vaut le résidu L(u)-S = r sur les bords ?
Sur les bords ou j'impose une solution de dirichlet je mettais 0 (car je me disais que 'u' vérifiant l'équation on a un résidu nul. Cependant l'opérateur L fait intervenir la dérivée seconde de 'u' en ce point du bord, et cette dérivée fait intervenir des points dans la grille, et donc varie au cours des itérations, donc je dirais maintenant que 'r' le résidu, doit être calculé sur ce bord aussi.

je regarde numerical recipes, et je vois page 906 du pdf (page 881 du document) qu'ils mettent le résidu à 0 sur tous les bords pour leur probleme de poisson avec dirichlet homogene.

Que penser ?

Si je donne une fausse valeur au résidu sur ce bord, le résidu servant de terme source a l'équation de l'erreur sur la grille plus grossière, je fausse mon calcul et j'engendre un probleme au bord.

Vous connaissez surement la réponse ?

merci