Ce petit "notebook" a pour but de montrer quelques propriétés de résolution d'un systeme linéaire de type "Elliptique contraint". Le langage utilisé est Julia, l'interface utilisée -le notebook- est Ipython avec le package IJulia. Les graphiques sont réalisés à l'aide du package PyPlot.
Les différents tests (benchmarks) seront menés sur un système d'équation provenant de la discrétisation différence finie (schéma de MAC) du problème se Stockes suivant :
$$ \begin{cases} -\mu \Delta \vec{u} + \nabla p = \vec{f} &\text{ dans } \Omega\\ \text{div}(\vec{u}) = 0&\text{ dans } \Omega \\ \vec{u}=\vec{g}& \text{ sur } \partial \Omega \quad \text{ avec } \int_{\partial \Omega} \vec{g}.\vec{n} d\sigma =0 \end{cases} $$On se donne une grille uniforme sur $]0,1[^2$ avec un pas $h=\frac{1}{n+1}$. Dans un premier temps prenons ne nous soucions pas des conditions de Dirichlet ($\vec{g}$).
On note
On cherche
Le schéma de MAC donne le système suivant : $$ \begin{cases} - \mu \Delta_h u^1_{i,j} + \frac{1}{2h}\left(p_{i+\frac{1}{2},j+\frac{1}{2}}+p_{i+\frac{1}{2},j-\frac{1}{2}}-p_{i-\frac{1}{2},j+\frac{1}{2}}-p_{i-\frac{1}{2},j-\frac{1}{2}} \right) = f^1_{i,j} & i,j=0,..,n+1\\- \mu \Delta_h u^2_{i,j} + \frac{1}{2h}\left(p_{i+\frac{1}{2},j+\frac{1}{2}}+p_{i-\frac{1}{2},j+\frac{1}{2}}-p_{i+\frac{1}{2},j-\frac{1}{2}}-p_{i-\frac{1}{2},j-\frac{1}{2}} \right) = f^2_{i,j} & i,j=0,..,n+1 \\\frac{1}{2h}\left(u^1_{i,j}+u^1_{i,j+1}-u^1_{i+1,j}-u^1_{i+1,j+1} \right) + \frac{1}{2h}\left(u^2_{i,j}+u^2_{i+1,j}-u^2_{i,j+1}-u^2_{i+1,j+1} \right)=0 & i,j=0,..,n \end{cases}$$
avec $-\Delta_h u_{i,j}=\frac{1}{h^2}\left(-u_{i-1,j} - u_{i,j-1} +4u_{i,j} -u_{i+1,j} -u_{i,j+1}\right)$ et les valeurs $p_{-\frac{1}{2},..}=p_{..,-\frac{1}{2}}=p_{n+\frac{3}{2},..}=p_{..,n+\frac{3}{2}}=0$.
Ceci donne un système matriciel de la forme : $$ \left(\begin{array}{c c c}A_1 & 0 &B_1^T\\0 & A_2 &B_2^T \\B_1& B_2& 0 \end{array}\right)\left(\begin{array}{c}u^1\\u^2\\p\end{array}\right)=\left(\begin{array}{c}f_1\\f_2\\0\end{array}\right) $$
avec
In [1]:
function Laplace1d(n)
h=1./(n+1)
d1=[-1/h^2 for i=1:n+1];
d2=[2/h^2 for i=1:n+2];
d3=[-1/h^2 for i=1:n+1];
spdiagm((d1,d2,d3),(-1,0,1));
end
function Laplace2d(n,m)
kron(speye(m+2),Laplace1d(n))+kron(Laplace1d(m),speye(n+2))
end
# différentiation 1d (forward)
function Dx1d(n)
h=1./(n+1)
d1=[-1/h for i=1:n+1];
d2=[1/h for i=1:n+1];
spdiagm((d1,d2),(-1,0))';
end
# différentiation 2d
function Dx2d(n,m)
d1=ones(m+1);
d2=ones(m+1);
tmp=spdiagm((d1/2,d2/2),(-1,0))
kron(tmp,Dx1d(n))
end
function Dy2d(n,m)
tmp=spdiagm((ones(n+1)/2,ones(n+1)/2),(-1,0))
kron(Dx1d(m),tmp)
end
Out[1]:
Pour l'instant nous n'avons pas pris en compte les conditions aux limites autrement dit il n'est pas tenu compte des éléments faisant parti de la frontière.
Supposons que l'on ai numéroté les sommets de tel sorte que le système s'écrive $$\left(\begin{array}{rr}A_{11} & A_{12}\\A_{21} & A_{22}\end{array}\right) \left(\begin{array}{l}V\\V_D\end{array}\right)=\left(\begin{array}{l}b\\b_D\end{array}\right)$$ avec $V_D$ les valeurs de $V$ sur la frontière donc connus ! Ce ne sont pas des inconnues on peut écrire formellement $V_D=g$
Atrement dit on doit résoudre le système : $$ A_{11} V=b-A_{12}g$$
Une autre solution est de travailler par pénalisation, on garde le premiès système et on ajoute 2 termes : $$\left(\begin{array}{rr}A_{11} & A_{12}\\A_{21} & A_{22}+\frac{1}{\varepsilon}Id\end{array}\right) \left(\begin{array}{l}V\\V_D\end{array}\right)=\left(\begin{array}{l}b\\b_D+\frac{1}{\varepsilon}g\end{array}\right)$$
En prenant $\varepsilon<<1$ on fixe numériquement la valeur de $V_D$ car on a de manière équivallente $$\begin{cases} A_{11}V + A_{12}V_D = b\\ \left(Id+{\varepsilon}A_{22}\right)V_D=g+ {\varepsilon}\left(A_{21}V+b_D\right)\end{cases}$$
Et donc $V_D=g+o(\varepsilon)$, le choix de $\varepsilon$ va dépendre du rayon spectral de $A_{11}$.
En prenant $\varepsilon=\rho(A_{11})^{-1}\times10^{-20}$ donne $\rho\left(Id+{\varepsilon}A_{22}\right)\leq 1+10^{-20}$...
In [2]:
n=3; #
mu=1e-2;
A1=mu*Laplace2d(n,n) #Matrice de diffusion pour la première composante de la vitesse
A2=mu*Laplace2d(n,n) #Matrice de diffusion pour la seconde composante de la vitesse
B1=Dx2d(n,n)
B2=Dy2d(n,n)
N=(n+2)^2;
A=[A1,spzeros(N,N),B1;spzeros(N,N), A2 , B2 ; B1' , B2' , spzeros((n+1)^2,(n+1)^2)]
Dans cet exemple, on impose une vitesse nulle partout sur le bord du carré $]0,1[^2$ sauf sur le bord supérieur $y = 1$ où la vitesse est horizontale.
Plus précisément, on choisit
In [8]:
full(Dx1d(N))
Out[8]:
In [4]:
full(B1[2,1:12])
Out[4]:
In [10]:
size(B2)
Out[10]:
In [58]:
size(A1)
Out[58]:
In [ ]: