TP 5 Algorithmes mathématiques: applications.¶
Commençons par importer tous les packages dont nous aurons besoin dans la suite.
import numpy as np
import numpy.random as npr
from matplotlib import pyplot as plt
import scipy.stats as st
import scipy.linalg as alg
from scipy.integrate import quad
Exercice 1¶
1) Ecrire une fonction python utilisant la méthode de rejet (avec $g$ densité de la loi $\mathcal{N}(0,1)$) et retournant un nombre aléatoire de densité¶
$$ f(x) = \frac{1}{K}\frac{e^{-x^2/2}}{1+x^2},\qquad x\in \mathbb{R}\qquad\textrm{avec}\qquad K = \int \frac{e^{-x^2/2}}{1+x^2} dx.$$
Pour calculer $K$ on pourra utiliser la fonction quad du sous module scipy.integrate.¶
Dans cette question nous implémentons le pseudo-code du cours dans une fonction python, avec comme agrument la fonction
$$ q(x) = \frac{f(x)}{cg(x)} \qquad x\in\mathbb{R}, $$
où les densités de probabilités $f$ et $g$ vérifient l'inégalité $f \leq c g$, avec $c=\sqrt{2\pi}/K$
def reject():
U = npr.uniform(0, 1)
Y = npr.normal(0, 1)
while U > q(Y):
U = npr.uniform(0, 1)
Y = npr.normal(0, 1)
return Y
On pourra remarquer que dans notre cas ce ratio se simplifie sous la forme
$$ q(x) = \frac{f(x)}{cg(x)} = \frac{K}{\sqrt{2\pi}}\frac{ e^{-x^2/2}}{K(1+x^2)}\frac{\sqrt{2\pi}}{ e^{-x^2/2}} = \frac{1}{1+x^2}.$$
def q(x):
return 1/(1+x**2)
2) A l'aide de la fonction précédente générer n=1000 nombres aléatoires de densité $f$.¶
Cette définition simplifier va être utilisée dans la fonction ci-dessous générant une échantillon de taille $n$ dont la densité est donnée par $f$.
n = 1000
X = np.array([reject() for _ in range(n)])
# On prend un numpy array pour avoir les methodes .min() et .max()
Une autre possibilité est la fonction suivante
def sample(n):
count = 0
q = lambda x : 1/(1+x**2)
X = np.array([])
while count < n:
X = np.append(X, reject())
count += 1
return X
Une autre possibilité est la suivante. Cependant cette solution aloue directement la mémoire pour $X$, et si $n$ est très grand cela pourrait ralentir le code dès le début.
def sample_bis(n):
count = 0
q = lambda x : 1/(1+x**2)
X = np.zeros(n)
while count < n:
X[count] = reject(q)
count += 1
return X
3) Réaliser un histogramme normalisé avec les nombres aléatoires obtenus à la question précédente.¶
Traçons l'histogramme normalisé avec $50$ classes
nbins = 50
plt.hist(X, nbins, density = True)
plt.xlabel(r'$X$')
plt.title('histogramme')
plt.show()
4) Ajouter sur cette histogramme la densité $f$. Que remarquez vous?¶
Rajoutons la courbe $f$ à l'histogramme précédent
nbins = 50
plt.hist(X, nbins, density = True, label = 'histogramme')
plt.xlabel(r'$X$')
plt.ylabel('densité de probabilité')
plt.title('histogramme vs densité')
f = lambda x: np.exp(-x**2/2)/(1 + x**2)
K, e = quad(f, -20, 20)
x = np.linspace(X.min(), X.max(), 100)
plt.plot(x, f(x)/K, label = 'densité')
plt.legend()
plt.show()
Si on augemente la taille $n$ de l'échantillon, on obtient:
n = 10000
X = sample(n)
nbins = 50
plt.hist(X, nbins, density = True, label = 'histogramme')
plt.xlabel(r'$X$')
plt.ylabel('densité de probabilité')
plt.title('histogramme vs densité')
x = np.linspace(X.min(), X.max(), 100)
plt.plot(x, f(x)/K, label = 'densité')
plt.legend()
plt.show()
et l'histogramme approche bien la densité de probabilité $f$ pour laquelle on voulait simuler des nombres aléatoires.
Exercice 2¶
1)¶
On considère une matrice $M$ symétrique de taille $n\times n$ dont les entrées $m_{jl}$ sont des variables aléatoires de loi $\mathcal{N}(0,1/(4n))$ pour $j<l$ et $\mathcal{N}(0,1/(2n))$ sur la diagonale. Toutes ces variables aléatoires sont indépendantes. La matrice étant symétrique elle est diagonalisable, et on ordonne les valeurs propres $(\lambda^{(n)}_j)_{j\in\{1,\dots,n\}}$ dans un ordre décroissant.
On s'intéresse à la répartition de ces valeurs propres. Pour cela, on fera un histogramme normalisé avec les valeurs propres de $M$ pour n=1000, qu'on comparera avec la densité de la loi semi-circulaire de rayon 1. Essayer avec différents tirages de $M$.
On utilisera les commandes eigvals du sous module scipy.linalg pour déterminer les valeurs propres et semicircular.pdf du sous module scipy.stats pour la densité de la loi semi-circulaire.
Quelle hypothèse peut-on émettre de ces observations? (c.f. théorème de Wigner).¶
Commençons par génerer la matrice $M$. On peut le faire sous forme de fonction.
def M_rnd(n):
M = npr.normal(0,1/np.sqrt(2 * n), (n, n))
return 0.5*(M + M.T)
Faisons les deux remarques suivantes:
La variance en argument de la fonction npr.normal est $\sigma$, et non $\sigma^2$. Vous pouvez aller voir la documentation en ligne pour plus de rensignement.
La fonction ci-dessus est suffisante. En effet, c'est tout d'abord direct de voir que la matrice produite est symmétrique. Ici la méthode $M.T$ renvoie la transposé de $M$. Ensuite, il faut se rappeler que si $X$ et $Y$ sont deux variables aléatoires Gaussiennes $N(0,\sigma^2)$ indépendantes, alors $X+Y$ est aussi une variable aléatoire Gaussienne d'espérance $$ \mathbb{E}[X+Y] = \mathbb{E}[X]+ \mathbb{E}[Y] = 0 $$ et de variance $$ Var[X+Y] = Var[X] + Var[Y] = 2\sigma^2 $$ et donc $$ Var\Big[\frac{X+Y}{2}\Big] = \frac{1}{4} Var[X+Y] = \frac{\sigma^2}{2}.$$ Donc les termes diagonaux de la matrice produite dans la fonction ci-dessus sont bien de variance $\sigma^2 = 1/(2n)$ et hors de la diagonale $\sigma^2/2 = 1/(4n)$.
Pour obtenir les valeurs propres il suffit d'éxécuter la commande suivante
n = 1000
vp_M = alg.eigvals(M_rnd(n))
Mais là on se rend compte que les résultats sont de type complexe, même si la partie imaginaire est nulle.
print(vp_M[0:4])
[-0.99770535+0.j -0.98007166+0.j -0.96999853+0.j -0.96350273+0.j]
Pour retirer la partie imagniaire on rajoute la méthode real à vp_M
print(vp_M.real[0:4])
[-0.99770535 -0.98007166 -0.96999853 -0.96350273]
Une remarque est qu'elles ne sont pas nécessairement ordonnées. On peut ainsi faire
vp_M = vp_M.real
vp_M.sort()
vp_M = vp_M[::-1]
print(vp_M[0:4])
[0.98704472 0.97973292 0.97065481 0.96335827]
On peut maintenant afficher l'histogramme de vp_M illustrant la répartition des valeurs propres de $M$, avec 50 classes.
n_bins = 50
plt.hist(vp_M.real, nbins, density =True, label = 'histogramme des vp de M')
plt.xlabel(r'$\lambda^{(n)}_j$')
plt.title('histogramme: répartition des valeurs propres de M')
plt.show()
Superposons maintenant sur cette histogramme la densité de la loi semi-circulaire
n_bins = 50
plt.hist(vp_M.real, nbins, density =True, label = 'histogramme des vp de M')
x = np.linspace(-1, 1, 100)
plt.plot(x, st.semicircular.pdf(x), label = 'densité loi semi-circulaire de rayon 1')
plt.xlabel(r'$\lambda^{(n)}_j$')
plt.ylabel('densité de probabilité')
plt.title('histogramme: répartition des valeurs propres de M')
plt.legend(loc = 'lower center')
plt.show()
Si on recommence avec une nouvelle matrice $M$, on retrouve approximativement la même répartition des valeurs propres.
n = 1000
vp_M = alg.eigvals(M_rnd(n))
n_bins = 50
plt.hist(vp_M.real, nbins, density =True, label = 'histogramme des vp de M')
x = np.linspace(-1, 1, 100)
plt.plot(x, st.semicircular.pdf(x), label = 'densité loi semi-circulaire de rayon 1')
plt.xlabel(r'$\lambda^{(n)}_j$')
plt.ylabel('densité de probabilité')
plt.title('histogramme: répartition des valeurs propres de M')
plt.legend(loc = 'lower center')
plt.show()
C'est ce que traduit le théorème de Wigner. La répartition des valeurs propres de grande matrice aléatoire peut être décrit de manière déterministe (universelle) à l'aide de la densité de la loi semi-circulaire.
2)¶
On considère une matrice $V$ de taille $n\,\times\,2n$ dont toutes les entrées sont des variables aléatoires $\mathcal{N}(0,1)$ indépendantes, et on pose $M = \frac{1}{2} V V^t$. Cette matrice de taille $n\times n$ est symétrique, et donc diagonalisable. On note $O$ la matrice orthogonale des vecteurs propres de $M$ dans la base canonique, classés suivant l'ordre croissant des valeurs propres. On considère $x\in\mathbb{R}^n$ un vecteur normalisé, $y = O^t x$, et la fonction $$ X_n(t) = \sqrt{\frac{n}{2}} \sum_{j = 1}^{[n t]}\Big(y_j^2 - \frac{1}{n}\Big),\qquad t\in[0,1].$$
Illustrer $10$ trajectoires de cette fonction pour différentes matrices $V$. On remarquera que $X_n(0) =X_n(1)= 0$. Est-ce que ce résultat change si on change de $x$ (toujours normalisé)?¶
On utilisera la commandes eig du sous module scipy.linalg pour déterminer les valeurs propres et vecteurs propres de $M$.
Pour générer $V$ et $M$ il suffit de faire
n = 500
V = npr.normal(0, 1, (n, 2*n))
M = 0.5 * np.dot(V , V.T)
Pour obtenir les valeurs propres et vecteurs propres de $M$ il suffit d'appliquer la commande
vp_M, O = alg.eig(M)
les valeurs propres sont dans vp_M et les vecteurs propres dans $O$. Mais comme précédement, les valeurs propres de $M$ ne sont pas forcément ordonnées. Pour obtenir les vecteurs propres suivant l'odre croissant des valeurs propres il faut connaitre l'ordre de permutation des valeurs propres pour les ordonner. Pour cela on utilise la fonction argsort de numpy.
permut = np.argsort(vp_M.real)
vp_M[permut[0:4]].real # les 5 premières valeurs propres
array([44.13246286, 47.80576464, 48.6542788 , 51.30007539])
Ainsi, pour obtenir les vecteurs propres dans l'ordre croissant des valeurs propres on applique permut sur les colonnes de $O$.
O = O[:, permut]
Prenons un vecteur $x\in \mathbb{R}^n$ normalisé quelconque. Choisissons le de manière aléatoire. Pour le normaliser on peut utiliser la fonction norm de scipy.linalg.
u = npr.uniform(-1,1,n)
x = u/alg.norm(u) # par défaut il s'agit de la norme Euclidienne
Définissons maintenant l'intervalle $[0,1]$ et la fonction $X_n$.
T = np.linspace(0, 1, 300)
y = np.dot(O.T, x)
Xn = np.array([np.sqrt(n/2)*np.sum(y[:int(n*t)]**2 - 1/n) for t in T])
On obtient le graphe d'une trajectoire.
plt.plot(T, Xn, 'b', linewidth = 0.4, label = r'$X_n(t)$')
plt.plot(T, [0]*len(T) , 'k', label = 'axe abscisse')
plt.xlabel(r'$t$')
plt.ylabel(r'$X_n(t)$')
plt.title(r'graphe de $X_n$')
plt.legend()
plt.show()
Maintenant, pour obtenir 10 trajectoires on peut procéder à l'aide d'une boucle de la façon suivante.
n = 500
u = npr.uniform(-1,1,n)
x = u/alg.norm(u)
N = 10
for j in range(N):
V = npr.normal(0, 1, (n, 2*n))
M = 0.5 * np.dot(V , V.T)
vp_M, O = alg.eig(M)
permut = np.argsort(vp_M.real)
O = O[:, permut]
y = np.dot(O.T, x)
Xn = np.array([np.sqrt(n/2)*np.sum(y[:int(n*t)]**2 - 1/n) for t in T])
if j=0:
plt.plot(T, Xn, 'b', linewidth = 0.4, label = r'$X_n(t)$')
else:
plt.plot(T, Xn, 'b', linewidth = 0.4, label = r'$X_n(t)$')
plt.plot(T, [0]*len(T) , 'k', label = 'axe abscisse')
plt.xlabel(r'$t$')
plt.ylabel(r'$X_n(t)$')
plt.title(r'graphe de $X_n$')
plt.legend()
plt.show()
Pour voir si $x$ influence le fait que toutes courbes commencent et finissent en $0$, il suffit de refaire tourner la cellule précédente. On voit que ceci n'y change rien.
En effet, pour $n$ grand ces courbes sont proches de celles de ce qu'on appelle un pont Brownien.
Exercice 3¶
On cherche ici à approximer numériquement la solution de l'équation de Poisson $$ (P):\,\left\{\begin{array}{ccc} \displaystyle -u''(x) = f(x)& \textrm{pour}& x\in]0,1[ \\ && \\ f(0)=f(1)=0. & & \end{array}\right. $$ Cette équation modélise une corde attachée aux points $0$ et $1$ et sur laquelle on exerce une force $f$. On va approcher $u$ au point $$ x_j^N = \frac{j}{N+1} \in[0,1],\qquad j\in\{0,\dots,N+1\},$$ et on va utiliser le fait que $$ u''(x) = \frac{1}{h^2}(u(x-h)-2u(x)+u(x+h)) + o(1)\qquad \textrm{pour}\quad h=\frac{1}{N+1}\ll 1. $$ La méthode qu'on considère est la suivante. On approche l'équation $(P)$ par le système $$ \frac{1}{h^2}(-u(x_{j-1}^N)+2u(x_j^N)-u(x_{j+1}^N)) + o(1)= f(x_j^N)\qquad j\in\{1,\dots,N\} $$ qui nous amène à résoudre le système $$(P_N):\quad A^N U^N = F^N$$ avec $$ A^N = \frac{1}{h^2}\begin{pmatrix} 2 & -1 & 0 & 0 &\dots& 0 \\ -1& 2 & -1 & 0&\dots &0 \\ 0& -1 &2 &-1&\dots &0 \\ \vdots & & \ddots &\ddots&\ddots&\vdots \\ 0& \dots & 0&-1&2& -1 \\ 0&\dots&\dots&0&-1&2 \end{pmatrix}\qquad F^N = \begin{pmatrix} f(x_1^N)\\ \vdots \\ f(x_N^N) \end{pmatrix} $$ où $U^N \in \mathbb{R}^N$ est l'inconnue approchant $u$ au point $x^N_j$: $$ U^N =\begin{pmatrix} U^N_1\\ \vdots \\ U^N_N) \end{pmatrix} \simeq \begin{pmatrix} u(x_1^N)\\ \vdots \\ u(x_N^N) \end{pmatrix}. $$ Ici, on a $U^N_0 = U^N_{N+1}=0$. Cette approximation peut se justifier rigoureusement!
1.a) Crée la matrice $A^N$ pour $N=100$. On pourra utiliser la fonction diag de numpy. On utilisera l'argument 1 et -1 avec cette fonction pour crée les diagonales supérieur et inférieur de $A^N$.¶
Dans cette question on génére la matrice $A^N$ à l'aide de la fonction diag de numpy. Pour la diagonale on va procéder ainsi
np.diag([2]*3)
array([[2, 0, 0],
[0, 2, 0],
[0, 0, 2]])
pour la sous diagonale
np.diag([-1]*2,-1)
array([[ 0, 0, 0],
[-1, 0, 0],
[ 0, -1, 0]])
et pour la sur diagonale
np.diag([-1]*2,1)
array([[ 0, -1, 0],
[ 0, 0, -1],
[ 0, 0, 0]])
Ainsi, pour $N=3$, en sommant les 3 matrices précédentes on a
np.diag([2]*3) + np.diag([-1]*2,-1) + np.diag([-1]*2,1)
array([[ 2, -1, 0],
[-1, 2, -1],
[ 0, -1, 2]])
On obtient alors $A^N$ pour $N=100$, de la façon suivanteet on a
N = 100
h = 1/(N+1)
A_N = ( np.diag([2]*N) + np.diag([-1]*(N-1),-1) + np.diag([-1]*(N-1),1) )/h**2
1.b) Résoudre le système $(P_N)$ pour $f\equiv -1$ avec $N=100$. On utilisera la fonction \texttt{solve} du sous module \texttt{scipy.linalg}. Dans ce cas la solution de $(P)$ peut être calculer explicitement et est $u(x) = x(x-1)/2$. Déterminer l'erreur uniforme commise:¶
$$ \max_{j\in\{1,\dots,N\}}|u(x_j^N) -U^N_j|. $$
Dans cette question on a
F = - np.ones(N)
On peut alors résoudre le système
U = alg.solve(A_N, F)
ajouter le fait que $U_0 = U_{N+1} = 0$
U = np.append(0, U)
U = np.append(U, 0)
et observer graphiquement la solution calculer numériquement via notre système $(P_N)$
x = np.linspace(0, 1, N+2) # on génère les points j/(N+1) pour j = 0, ..., N+1
plt.plot(x, U, label = 'sol simu')
plt.xlabel(r'$x$')
plt.ylabel(r'$U^N$')
plt.title('solution numérique de $-u\'\' = f$, avec $f\equiv -1$')
plt.show()
On peut alors superposer la solution exacte de cette equation : $u(x) = x(x-1)/2$ pour observer l'erreur d'approximation.
x = np.linspace(0, 1, N+2)
plt.plot(x, U, label = 'sol simu')
plt.plot(x, x*(x-1)/2, label = 'sol exacte')
plt.xlabel(r'$x$')
plt.ylabel(r'$U^N$')
plt.title('solution numérique et exacte de $-u\'\' = f$, avec $f\equiv -1$')
plt.legend()
plt.show()
On ne voit aucune distinction, l'approximation est très bonne comme le confirme le calcul de l'erreur uniforme $$\max_{j\in \{0,\dots,N+1 \}} |u(x^N_j)- U^N_j|$$
alg.norm( x*(x-1)/2 - U ,np.inf) # np.inf comme argument pour avoir la norme infinie
8.7430063189231078e-16
1.c) Essayer ensuite différentes forces $f$ et observer le résultat.¶
On va anticiper sur la question suivante en introduisant la fonction $f$ définie en question 2.d adaptée au cas de $[0,1]$.
def f(X,X0,r):
F = np.array([])
for x in X:
F = np.append(alg.norm(x-X0) <= r,F)
return - F/r
r = 0.05
X0 = 0.5
F = f(x[1:N+1],X0,r)
U = alg.solve(A_N, F)
U = np.append(0, U)
U = np.append(U, 0)
plt.plot(x, U, label = 'sol simu')
plt.xlabel(r'$x$')
plt.ylabel(r'$U^N$')
plt.title('solution numérique de $-u\'\' = f$')
plt.show()
Avec deux points de poussé maintenant,
r = 0.05
X0 = 0.25
X1 = 0.75
F = f(x[1:N+1], X0, r) + f(x[1:N+1], X1, r)
U = alg.solve(A_N, F)
U = np.append(0, U)
U = np.append(U, 0)
plt.plot(x, U, label = 'sol simu')
plt.xlabel(r'$x$')
plt.ylabel(r'$U^N$')
plt.title('solution numérique de $-u\'\' = f$')
plt.show()
Q.2¶
On cherche ici à approximer numériquement la solution de l'équation de Poisson $$ (P):\,\left\{\begin{array}{ccl} \displaystyle - \Delta u(x,y) = f(x,y)& \textrm{pour}& x\in]0,1[^2 \\ && \\ f(x,y)=0 &\textrm{pour} & \partial [0,1]^2 \qquad \textrm{(sur le bord du carré)} \end{array}\right.$$ où $\Delta = \frac{\partial ^2}{\partial x^2}+\frac{\partial ^2}{\partial y^2}$ est le Laplacien 2D. Cette équation modélise une membrane carrée sur laquelle on exerce une force $f$. On va approcher $u$ aux points $$ x_j^N = y_j^N = \frac{j}{N+1} \in[0,1],\qquad j\in\{0,\dots,N+1\},$$ et on va utiliser le fait que $$\Delta u(x,y) = \frac{1}{h^2}(u(x-h,y)+u(x+h,y)-4u(x,y)+u(x,y-h)+u(x,y+h)) + o(1). $$ On va donc considérer le système $$ -\frac{1}{h^2}(U^N_{j-1,l}+U^N_{j+1,l}-4U^N_{j,l}+U^N_{j,l-1}+U^N_{j,l+1}) = f(x^N_j,y^N_l),\qquad (j,l)\in\{1,\dots,n\}^2, $$ d'inconnues $U^N=(U^N_{j,l})_{(j,l)\in\{1,\dots,n\}^2}$ approchant $u$ au point $(x^N_j,y^N_l)$: $U^N_{j,l}\simeq u(x^N_j,y^N_l)$. On peut réécrire le système précédent sous forme matricielle $$(P_N):\quad \mathbf{A}^N \mathbf{U}^N = \mathbf{F}^N$$ avec $$\mathbf{U}^N = {}^t (U^N_{1,1},U^N_{2,1},\dots,U^N_{N,1},U^N_{1,2},\dots,U^N_{1,N},\dots, U^N_{N,N})\in \mathbb{R}^{N^2},$$ $\mathbf{F}^N$ définie de la même manière avec les valeurs $f(x^N_j,y^N_l)$, et $$\mathbf{A}^N = \frac{1}{h^2} \begin{pmatrix} S_N & -I_N & O &\dots & O \\ - I_N & S_N & -I_N &\ddots& \vdots \\ O & \ddots &\ddots &\ddots &O \\ \vdots &\ddots& -I_N& S_N& I_N \\ O &\dots& O & -I_N & S_N \end{pmatrix} $$ avec $$S_N = \begin{pmatrix} 4 & -1 & 0 & 0 &\dots& 0 \\ -1& 4 & -1 & 0&\dots &0 \\ 0& -1 &4 &-1&\dots &0 \\ \vdots & & \ddots &\ddots&\ddots&\vdots \\ 0& \dots & 0&-1&4& -1 \\ 0&\dots&\dots&0&-1&4 \end{pmatrix}, $$ et $I_N$ la matrice identité de taille $N\times N$. Ici aussi, on a considéré $U^N_{0,l}=U^N_{N+1,l}= U^N_{j,N+1} = U^N_{j,N+1}=0$.
a) Crée la matrice $S_N$ comme dans la question précédente, puis $\mathbf{A}^N$.¶
Pour créer la matrice $S_N$ on procède comme dans la question précédente. Mais pour comprendre la construction de $\mathbf{A}^N$ commençons par considérer $N=3$
N = 3
S_N = ( np.diag([4]*N) + np.diag([-1]*(N-1),-1) + np.diag([-1]*(N-1),1) )
Pour créer la matrice $\mathbf{A}^N$ il faut réfléchire par bloc. Rappelons que cette matrice est de taille $N^2 \times N^2$. Donc ici de taille $9\times 9$.
bA_N = np.zeros((N**2, N**2))
for j in range(N):
bA_N[j*N:(j+1)*N, j*N:(j+1)*N] = S_N
if j <= N-2:
bA_N[j*N:(j+1)*N, (j+1)*N:(j+2)*N] = - np.eye(N,N)
if j>=1:
bA_N[j*N:(j+1)*N, (j-1)*N:j*N] = - np.eye(N,N)
print(bA_N)
[[ 4. -1. 0. -1. -0. -0. 0. 0. 0.] [-1. 4. -1. -0. -1. -0. 0. 0. 0.] [ 0. -1. 4. -0. -0. -1. 0. 0. 0.] [-1. -0. -0. 4. -1. 0. -1. -0. -0.] [-0. -1. -0. -1. 4. -1. -0. -1. -0.] [-0. -0. -1. 0. -1. 4. -0. -0. -1.] [ 0. 0. 0. -1. -0. -0. 4. -1. 0.] [ 0. 0. 0. -0. -1. -0. -1. 4. -1.] [ 0. 0. 0. -0. -0. -1. 0. -1. 4.]]
On observe bien les blocs $S_N$ sur la diagonale, et les bloc $-I_N$ pour les sur et sous diagonales de $\mathbf{A}^N$. Passons au cas $N=100$ pour construire $\mathbf{A}^N$.
N = 50
S_N = ( np.diag([4]*N) + np.diag([-1]*(N-1),-1) + np.diag([-1]*(N-1),1) )
bA_N = np.zeros((N**2, N**2))
for j in range(N):
bA_N[j*N:(j+1)*N, j*N:(j+1)*N] = S_N
if j <= N-2:
bA_N[j*N:(j+1)*N, (j+1)*N:(j+2)*N] = - np.eye(N,N)
if j>=1:
bA_N[j*N:(j+1)*N, (j-1)*N:j*N] = - np.eye(N,N)
h = 1/(N+1)
bA_N = bA_N/h
b) Résoudre le système $(P_N)$ pour $N=100$ (commencer avec moins pour que votre machine ne chauffe pas trop) pour $f\equiv-1$.¶
Construisons $F^N$
F = - np.ones(N**2)
et calculons $U^N$
U = alg.solve(bA_N, F)
c) Pour faire les représentations graphiques, remettre tout d'abord $\mathbf{U}^N$ sous la forme d'une matrice $N\times N$ (on pourra utiliser la méthode reshape(N,N)), et pour les graphiques eux-mêmes les fonctions matshow ou imshow du sous module matplotlib.pyplot.¶
On redimensionne $U^N$ pour correspondre au carré $[0,1]^2$
U = U.reshape((N, N))
et rajoutons les conditions de nullité aux bords
U = np.concatenate( (np.zeros((1,N)), U, np.zeros((1,N)) ), axis = 0)
U = np.concatenate( (np.zeros((N+2,1)), U, np.zeros((N+2,1)) ), axis = 1)
plt.imshow(U)
plt.show()
d) Résoudre le système pour¶
$$f(x,y)=\left\{\begin{array}{ccc} - 1/r^2 & \textrm{si}& \|(x,y)-(x_0,y_0)\|_2 < r \\ 0 &\textrm{sinon,} \end{array}\right.$$
avec $r=0.05$ et $(x_0,y_0)$ que vous choisirez. Cette fonction modélise l'action d'un doigt sur la membrane. On pourra placer différentes pressions à différents endroits simultanément.¶
On va maintenant réituliser notre fonction $f$
def f(X,Y,X0,r):
F = np.array([])
for x in X:
for y in Y:
F = np.append(alg.norm(np.array([x,y])-X0) <= r,F)
return - F/r
et on recalcule $U^N$ comme précédemment
r = 0.05
X0 = np.array([0.5, 0.5])
x = y = np.linspace(0,1,N)
F = f(x,y,X0,r)
U = alg.solve(bA_N, F)
U = U.reshape((N, N))
U = np.concatenate( (np.zeros((1,N)), U, np.zeros((1,N)) ), axis = 0)
U = np.concatenate( (np.zeros((N+2,1)), U, np.zeros((N+2,1)) ), axis = 1)
Ce qui nous donne
plt.imshow(U)
plt.show()
Maintenant, on va considérer un force à différents points pour modéliser l'action des 5 doigts de la main gauche sur la membrane.
F = 0
for X0 in [[0.3, 0.1], [0.7, 0.2], [0.8, 0.4], [0.8, 0.6], [0.7, 0.8]]:
F += f(x,y,np.array(X0),r)
On résoud notre système comme précédemment
U = alg.solve(bA_N, F)
U = U.reshape((N, N))
U = np.concatenate( (np.zeros((1,N)), U, np.zeros((1,N)) ), axis = 0)
U = np.concatenate( (np.zeros((N+2,1)), U, np.zeros((N+2,1)) ), axis = 1)
Ce qui nous donne
plt.imshow(U)
plt.show()
On peut aussi visualiser la pression de la main sur la membrane d'une autre manière.
from mpl_toolkits.mplot3d import axes3d
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
x = y = np.linspace(0,1,N+2)
X, Y = np.meshgrid(x, y)
ax.plot_wireframe(X, Y, U, rstride=1, cstride=1)
<mpl_toolkits.mplot3d.art3d.Line3DCollection at 0x7f3172774c18>
