Mathematics">
PDF Exam Dec 2008 - Copie
PDF Exam Dec 2008 - Copie
PDF Exam Dec 2008 - Copie
—
PHYSIQUE NUMÉRIQUE
Examen du 11 décembre 2008
I— On veut calculer le potentiel V (x, y) (le problème est ici à 2 dimensions) dû à une distribution
de charges ρ(x, y). Pour cela, on peut utiliser l’équation de Poisson :
∂2V ∂2V ρ
2
+ 2
=−
∂x ∂y ε0
où ε0 est la permittivité du vide, ρ(x, y) est connu et V (x, y) ce que l’on cherche à calculer. On choisira
d’utiliser un système d’unités dans lequel ε0 = 1.
1– On discrétise le problème selon les directions x et y. Montrer que l’on peut reécrire l’équation
de Poisson de la façon suivante :
Vi+1,j − 2Vi,j + Vi−1,j Vi,j+1 − 2Vi,j + Vi,j−1
2
+ = −ρi,j
δx δy2
3– On pose
V1 ρ1
.. ..
.
.
V= Vℓ et R= ρℓ
.. ..
. .
Vn x n y ρn x n y
On écrit alors l’équation de Poisson discrétisée sous forme matricielle :
D·V =R
1
contenant les éléments de D. On s’attachera à écrire le sous-programme complet, y compris les
déclarations. On prendra garde à ce que les indices des éléments de tableaux utilisés ne débordent
pas l’intervalle [1,nx*ny].
4– On dispose d’un sous-programme subroutine linsyst(a,b,n) où a est un tableau
real, dimension(n,n) :: a
b un tableau
real, dimension(n) :: b
et n un scalaire entier.
Ce sous-programme calcule X = A−1 · B où la matrice A est représentée par le tableau a et le
vecteur B par le tableau b. Il faut fournir en entrée à linsyst les tableaux a et b ainsi que n. En
sortie, linsyst fournit le vecteur X solution dans le tableau b.
A– On veut évidemment utiliser ce sous-programme pour résoudre l’équation de Poisson. Quelle
valeur doit-on donner à n ?
B– Il faudra créer un tableau real, dimension(nx,ny) :: rho dans lequel on mettra les va-
leurs de ρ(x, y). Pourra-t-on envoyer ce tableau comme deuxième argument dans linsyst ? Pourquoi ?
C – Ecrire un programme principal qui doit :
1. Créer et remplir le tableau rho. On prendra comme exemple un dipôle : rho est nul partout sauf
en (nx/2-10,ny/2) et (nx/2+10,ny/2) où il prend les valeurs −1 et +1 respectivement.
2. Copier ce tableau dans un autre de dimensions convenablement choisies pour pouvoir l’utiliser
comme deuxième argument de linsyst. On utilisera la relation (1).
3. Appeler linsyst et écrire le résultat dans un fichier de résultats.
II— Le programme ainsi conçu, pour nx = ny = 50, produit la figure ci-dessous, à gauche, ce qui
paraı̂t assez satisfaisant :
0.0003 0.00035
0.0002 0.0003
0.00025
1e-04
0.0002
0
0.00015
-0.0001
0.0001
1 1
-0.0002 0.9 5e-05 0.9
0.8 0.8
0.7 0.7
-0.0003 0 0.6 00 0.6
0.1 0.2 0.5 0.1 0.2 0.5
0.3 0.4 0.4 Y 0.3 0.4 0.4 Y
0.5 0.6 0.3 0.5 0.6 0.3
0.7 0.8 0.2 0.7 0.8 0.2
X 0.1 X 0.1
0.9 1 0 0.9 1 0
Potentiel obtenu pour la distribution de charges Potentiel obtenu pour la distribution de charges
d’un dipole. décrite à la question II-1
1– Toutefois, la figure de droite est obtenue en définissant la distribution de charge de la façon
suivante (toujours avec nx = ny = 50) :
do i = 1, nx
do j = 1, ny
r = sqrt( ((i-nx/2)*delta_x)**2 + ((j-ny/2)*delta_y)**2 )
rho(i,j) = exp(-r**2/(2*0.002**2))
enddo
enddo
A– Quelle est la symétrie de cette distribution ? Quelle est celle de la figure ? Les deux axes x
et y sont-ils équivalents ? Est-ce compatible ?
B– Si, par exemple, on choisit nx = ny = 3, la matrice D devient une matrice (9 × 9) qui a
l’allure ci-dessous :
2
⋆ ⋆ 0 ⋆ 0 0 0 0 0
⋆ ⋆ ⋆ 0 ⋆ 0 0 0 0
0 ⋆ ⋆ ⋆ 0 ⋆ 0 0 0
⋆ 0 ⋆ ⋆ ⋆ 0 ⋆ 0 0
D=
0 ⋆ 0 ⋆ ⋆ ⋆ 0 ⋆ 0
0 0 ⋆ 0 ⋆ ⋆ ⋆ 0 ⋆
0 0 0 ⋆ 0 ⋆ ⋆ ⋆ 0
0 0 0 0 0
⋆ ⋆ ⋆ ⋆
0 0 0 0 0 ⋆ 0 ⋆ ⋆
où le symbole ⋆ représente les éléments non nuls. Un élément du tableau est repéré par les indices ℓ
et ℓ′ .
a) Suivant l’équation (1), à quels indices (i, j, i′ , j ′ ) correspondent les trois premières lignes de
cet exemple.
Combien d’éléments de la matrice ont-ils été supprimés du fait du débordement de l’intervalle
[1, nx ny ] pour ℓ ou ℓ′ :
b) à cause des termes en i + 1, i − 1
c) à cause des termes en j + 1, j − 1 ?
C – Afin de rétablir la symétrie, on modifie le calcul de D de la façon suivante :
d = 0. ! initialisation
do i = 1, nx ! termes non-diagonaux
do j = 1, ny
l = i + (j-1)*nx
if ( i+1 <= nx ) d(l,l+1) = -1./delta_x**2
if ( i-1 >= 1 ) d(l,l-1) = -1./delta_x**2
if ( j+1 <= ny ) d(l,l+nx) = -1./delta_y**2
if ( j-1 >= 1 ) d(l,l-nx) = -1./delta_y**2
enddo
enddo
Les éléments de tableau d(l,m) pour lesquels m n’est pas dans l’intervalle [1, nx ny ] sont simplement
ignorés.
a) Qu’est-ce qui a changé par rapport 0.00035
0.00025
3-B où l’on vérifiait simplement que les
0.0002
indices composites étaient dans l’intervalle 0.00015
[1, nx ny ] ? 0.0001
1
b) Le résultat est montré sur la figure 5e-05
0.7
0.8
0.9
00 0.6
ci-contre : quel rapport cela peut-il avoir 0.1 0.2
0.3 0.4
0.3
0.4
0.5
Y
0.5 0.6 0.2
avec la modification apportée ? X
0.7 0.8
0.9 1 0
0.1
3
Corrigé
I—
1– Discrétisation en x : xi = iδx , où δx est le pas de discrétisation. De même en y : yj = jδy . On
pose alors :
V (xi , yj ) = Vi,j et ρ(xi , yj ) = ρi,j
2
La dérivée seconde ∂∂xV2 , prise en (x, y), est la dérivée première de la dérivée première que l’on peut
approximer comme suit :
∂V δx ∂V δx V (x+δx ,y)−V (x,y) V (x,y)−V (x−δx ,y) Vi+1,j −Vi,j Vi,j −Vi−1,j
∂2V ∂x (x + 2 , y) − ∂x (x − 2 , y) δx − δx δx − δx
(x, y) ∼ ∼ =
∂x2 δx δx δx
La même opération avec y donne l’expression cherchée.
2– Le couple (i + 1, j) donne ℓ + 1, alors que le couple (i, j + 1) donne ℓ + nx . Même chose avec le
signe −.
3–
2 2
A– Les éléments diagonaux de D sont Dℓ,ℓ = δx2 + δy2 . Les éléments Dℓ,ℓ±1 s’écrivent − δ12 et les
x
éléments Dℓ,ℓ±nx , − δ12 . Tous les autres éléments de la matrice sont nuls.
y
B–
d = 0. ! initialisation de l’ensemble.
do l = 1, nx*ny
! termes diagonaux
d(l,l) = 2./delta_x**2 + 2./delta_y**2
! termes non-diagonaux
k = l+1
if ( k <= nx*ny ) d(l,k) = -1./delta_x**2
k = l-1
if ( k >= 1 ) d(l,k) = -1./delta_x**2
k = l+nx
if ( k <= nx*ny ) d(l,k) = -1./delta_y**2
k = l-nx
if ( k >= 1 ) d(l,k) = -1./delta_y**2
enddo
end
Autre possibilité
d = 0. ! initialisation de l’ensemble.
do l = 1, nx*ny
! termes diagonaux
4
d(l,l) = 2./delta_x**2 + 2./delta_y**2
enddo
! termes non-diagonaux
do l = 2, nx*ny
d(i,i-1) = -1./delta_x**2
d(i-1,i) = -1./delta_x**2
enddo
do l = nx+1, nx*ny
d(i,i-nx) = -1./delta_y**2
d(i-nx,i) = -1./delta_y**2
enddo
4–
A– n=nx*ny
B– Non ! Le tableau b est de dimension (nx*ny) (un indice ∈ [1, nx ny ]) alors que rho est de
dimension (nx,ny) (deux indices).
C – Une possibilité (ici, tout est en double précision, par précaution) :
program poisson_simple
implicit none
integer :: nx, ny, i, j, l
real(8), dimension(:,:), allocatable :: d, rho
real(8), dimension(:), allocatable :: v
real(8) :: delta_x, delta_y
do i = 1, nx ; do j = 1, ny
l = i + (j-1)*nx ; v(l) = rho(i,j)
enddo ; enddo
! les 3 lignes ci-dessus peuvent etre remplacées par:
! v = pack(rho,.TRUE.)
call linsyst(d, v, nx*ny)
open(1,file=’poisson_simple.res’)
do i = 1, nx ; do j = 1, ny
write(1,*) i*delta_x, j*delta_y, v(i + (j-1)*nx)
enddo ; write(1,*) ; enddo
deallocate(d,rho,v)
end
II—
1–
5
A– La distribution de charges est à symétrie de révolution autour du centre, alors que le potentiel
obtenu n’a que des symétries d’ordre 2 par rapport aux axes, en particulier l’axe x et l’axe y ne sont
pas équivalents.
Ce n’est évidemment pas compatible avec la distribution de charges initiale.
B–
a) j ′ = 1. Les trois premières colonnes j = 1 et i, i′ = 1, 2, 3, les trois suivants j = 2, etc.
b) 2
c) 6
C–
a) Au lieu de vérifier globalement que les indices de éléments de matrice tombent bien dans
l’intervalle [1, nx ny ], on vérifie séparement pour chaque direction que les indices i − 1, i + 1, d’une part,
sont dans l’intervalle [1, nx ] et que les indices j − 1 et j + 1, de l’autre, sont dans l’intervalle [1, ny ].
La matrice D (dans l’exemple nx = 3, ny = 3) prend l’allure suivante :
⋆ ⋆ 0 ⋆ 0 0 0 0 0
⋆ ⋆ ⋆ 0 ⋆ 0 0 0 0
0 ⋆ ⋆ 0 0 ⋆ 0 0 0
⋆ 0 0 ⋆ ⋆ 0 ⋆ 0 0
D=
0 ⋆ 0 ⋆ ⋆ ⋆ 0 ⋆ 0
0 0 ⋆ 0 ⋆ ⋆ 0 0 ⋆
0 0 0 ⋆ 0 0 ⋆ ⋆ 0
0 0 0 0 0
⋆ ⋆ ⋆ ⋆
0 0 0 0 0 ⋆ 0 ⋆ ⋆
et qu’on se borne à tronquer la matrice correspondante aux bords de l’intervalle le potentiel est
implicitement considéré comme nul à l’extérieur du domaine : or il est considéré comme continu.
—
Pour info, voici un linsyst pas toujours très performant (pour (50 × 50), ça fait quelques dizaines
de secondes sur un PC) mais très simple et raisonnablement fiable :
subroutine linsyst(a, b, n)
implicit none
integer, intent(in) :: n
real(8), dimension(n,n), intent(inout) :: a
real(8), dimension(n) , intent(inout) :: b
integer :: i, j
! LU decomposition
do j = 1, n
do i = 2, j
6
a(i,j) = a(i,j) - sum(a(i,1:i-1)*a(1:i-1,j))
enddo
do i = j+1, n
a(i,j) = (a(i,j) - sum(a(i,1:j-1)*a(1:j-1,j)))/a(j,j)
enddo
enddo
! L.y = b (L(i,i) = 1)
do i = 2, n
b(i) = b(i) - sum(a(i,1:i-1)*b(1:i-1))
enddo
! U.x = y
b(n) = b(n)/a(n,n)
do i = n-1, 1, -1
b(i) = (b(i) - sum(a(i,i+1:n)*b(i+1:n)) )/a(i,i)
enddo
end