\documentclass{TDPhysique}
\partie{Informatique pour tous}
\chapitre{Analyse numérique}
\annee{2020-2021}


\begin{document}

\begin{TD}
\begin{center}
\begin{LARGE}
\textcolor{DodgerBlue!50!black}{Résolution numérique de l'équation de Laplace}
\end{LARGE}
\end{center}
\end{TD}

Les équations de Maxwell permettent de montrer que, en régime stationnaire et dans le vide, le potentiel électrostatique $V$ 
vérifie l'équation de Laplace, c'est-à-dire
\[\Delta V = 0\]

où $\Delta$ est l'opérateur laplacien scalaire. 

Cette équation se retrouve dans de nombreux domaines : la température de
tout système vérifie $\Delta T = 0$ en régime permanent et sans source de production, mais on la rencontre également en gravitation avec le potentiel associé au champ gravitationnel, en mécanique des fluides, pour le potentiel des vitesses, en mécanique quantique, etc. 

Des solutions analytiques existent dans des cas simples, mais dans le cas général, résoudre cette équation demande d'avoir recours à des méthodes
numériques : c'est le but du travail suit. 

Nous nous limiterons à deux dimensions, ce qui revient physiquement à supposer une invariance par translation dans la troisième dimension.

Vous disposez du script Python \lstinline$Eq_laplace_script-depart.py$, qui contient quelques fonctions qui nous serviront à l'initialisation des calculs et à l'affichage des résultats. 

Ouvrir ce fichier et l'\textit{enregistrer sous} un nom différent et dans \textit{vos} documents.

\section{Principe de la résolution par méthode de relaxation}

\subsection{Schéma général}
Il est évidemment impossible de résoudre numériquement l'équation de Laplace en tout point de l'espace : on se restreint à un domaine $\mathcal{D}$ fini, discrétisé sous forme d'une grille de dimensions $N_x \times N_y$.

On peut démontrer mathématiquement que l'équation de Laplace admet une unique solution dans un domaine fermé $\mathcal{D}$ si la valeur de $V$ est fixée sur les bords $B$ du domaine (on parle en mathématiques de « problème de Dirichlet »). 

Dans notre cas, les bords correspondent au contour de la grille, mais peuvent aussi inclure d'autres points à l'intérieur de celle-ci sur lesquels le potentiel serait imposé, par exemple les armatures d'un condensateur.

Les conditions aux limites étant fixées, on utilise ensuite une méthode de résolution par différences finies. 

L'idée est exactement la même que pour résoudre une équation différentielle par la méthode d'Euler : les dérivées spatiales 
du potentiel sont approximées par des différences entre les valeurs du potentiel entre deux points voisins de la grille.

La nouveauté par rapport à la méthode d'Euler est la résolution itérative : partant d'une « condition initiale » arbitraire (valeur de $V$ fixée arbitrairement et « aléatoirement » en tout point de la grille), on améliore par récurrence la précision du résultat pour qu'il soit conforme à l'équation de Laplace.
 
\subsection{Discrétisation de l'équation de Laplace}

On rappelle qu'on se place sur une grille à deux dimensions de taille $N_x \times N_y$, schématisée figure 1. Le pas spatial
de la grille est noté $\delta$, si bien que le point $(i, j)$ de la grille a pour coordonnées $(x_i, y_j) = (i\delta, j\delta)$. Les indices sont donc écrits dans tout ce travail dans une convention de type « abscisse, ordonnée » plutôt que « ligne, colonne », comme
ils le seraient pour un élément de matrice.

\begin{center}
\includegraphics[width=8 cm]{Fig1_Grille.png}
\caption{Grille de discrétisation}
\end{center}

À deux dimensions cartésiennes, le laplacien s'écrit

\[\Delta V = \frac{\partial^2 V}{\partial x^2}+\frac{\partial^2 V}{\partial y^2}\]

Il faut donc discrétiser la dérivée seconde sous forme de différence finie, ce qui se fait comme toujours grâce à des
développements limités. 

En transposant la formule de Taylor,

\[f(a + h) = f(a) + h f'(a) + \frac{h^2}{2} f''(a) + O(h^3) \] 

on peut écrire : 

\[V(x_i+\delta,y_j)=V(x_i,y_j)+\delta\frac{\partial V}{\partial x}(x_i,y_j)+\frac{\delta^2}{2}\frac{\partial^2 V}{\partial x^2}(x_i,y_j)+\mathrm{O}(\delta^3)\]

et de même

\[V(x_i-\delta,y_j)=V(x_i,y_j)-\delta\frac{\partial V}{\partial x}(x_i,y_j)+\frac{\delta^2}{2}\frac{\partial^2 V}{\partial x^2}(x_i,y_j)+\mathrm{O}(\delta^3)\]


Le potentiel $V$ sera stocké dans une variable globale sous forme d'un tableau \lstinline$numpy$ \lstinline$V$ tel que \lstinline$V[i,j]$ donne la valeur du potentiel $V(x_i, y_j)$. 

Par simplicité de notation, on notera donc dans la suite $V[i,j]=V(x_i,y_j)$

\quest 1 - Montrer que

\[\frac{\partial^2V}{\partial x^2}[i,j]\simeq\frac{V[i+1,j]+V[i-1,j]-2V[i,j]}{\delta^2} 
\]

puis que

\[\Delta V [i, j] \simeq \frac{V [i + 1, j] + V [i - 1, j] + V [i, j + 1] + V [i, j - 1] - 4V [i, j]}{\delta^2}.\]

\vspace{0.5cm}

\begin{correction}
On fait les sommes : 
$$V(x_i+ \delta,y_j)+V(x_i-\delta,y_j)= 2V(x_i,y_j)+\delta^2\left(\frac{\partial^2 V}{\partial x^2}\right)(x_i,y_j) + O(\delta^3)$$

D'où le résultat avec les notations proposées. 

\[\frac{\partial^2V}{\partial x^2}[i,j]\simeq\frac{V[i+1,j]+V[i-1,j]-2V[i,j]}{\delta^2} 
\]

\end{correction}

\begin{correction}
De même pour la dérivée seconde par rapport à $y$.

On fait les sommes : 
$$V(x_i,y_j+ \delta)+V(x_i,y_j-\delta)= 2V(x_i,y_j)+\delta^2\left(\frac{\partial^2 V}{\partial y^2}\right)(x_i,y_j) + O(\delta^3)$$

D'où le résultat avec les notations proposées. 

\[\frac{\partial^2V}{\partial y^2}[i,j]\simeq\frac{V[i,j+1]+V[i,j-1]-2V[i,j]}{\delta^2} 
\]

En sommant ces deux derniers résultats, soit les quatre termes calculés avec la formule de Taylor, on obtient, comme  $\Delta V = \frac{\partial^2 V}{\partial x^2}+\frac{\partial^2 V}{\partial y^2}$: 

\[\frac{\partial^2V}{\partial x^2}[i,j]+\frac{\partial^2V}{\partial y^2}[i,j]\simeq\frac{V[i+1,j]+V[i-1,j] + V[i,j+1]+V[i,j-1]-4V[i,j]}{\delta^2} 
\]
\end{correction}


Ainsi, l'équation de Laplace discrétisée s'écrit, pour tout point $(i, j)$ du domaine de résolution,

\[\Delta V [i, j] = 0\]
Soit 

\[\frac{V [i + 1, j] + V [i - 1, j] + V [i, j + 1] + V [i, j - 1] - 4V [i, j]}{\delta^2} = 0\]

ce que l'on peut encore mettre sous la forme

\begin{equation}
V [i, j] = \frac{V [i + 1, j] + V [i - 1, j] + V [i, j + 1] + V [i, j - 1]}{4} 
\label{EqLaplace}
\end{equation}


Lorsqu'il vérifie l'équation de Laplace, le potentiel en un point est donc une moyenne des potentiels aux points voisins : c'est ainsi qu'il sera possible de construire par récurrence les valeurs de $V$ en tout point de la grille en
partant des bords où les conditions aux limites sont connues.


\subsection{Définition des conditions aux limites}
Rappelons que l'on appelle ici « bord » du domaine de résolution l'ensemble des points où la valeur du potentiel est imposée. Pour que la méthode utilisée s'applique, le bord doit nécessairement inclure le contour géométrique de
la grille mais peut également inclure d'autres points intérieurs.

Les points appartenant au bord du domaine sont stockés sous forme d'un tableau de booléens \lstinline$B$, de taille $N_x \times N_y$.

Ce tableau est construit de la façon suivante :
\begin{itemize}
	\item si le point $(i, j)$ appartient au bord, alors \lstinline$B[i,j]$ vaut \lstinline$True$ : le potentiel en ce point est imposé, sa valeur ne doit
jamais être modifiée. 

	\item si le point $(i, j)$ n'appartient pas au bord, alors \lstinline$B[i,j]$ vaut \lstinline$False$ : le potentiel en ce point est inconnu \textit{a priori}
mais doit vérifier l'équation (\ref{EqLaplace}).

\end{itemize}

Les tableaux décrivant le bord et le potentiel sont respectivement initialisés à des \lstinline$False$ et à des zéros.

\quest Écrire une fonction \lstinline$initialisation_contour(V0)$ qui prend comme argument un flottant \lstinline$V0$. Cette fonction initialise les points du contour de la grille : elle leur affecte la valeur \lstinline$True$ dans \lstinline$B$ et la valeur \lstinline$V0$ dans \lstinline$V$ sans toucher
aux autres points. Les variables \lstinline$B$ et \lstinline$V$ étant globales, la fonction les modifie mais ne retourne rien.

\begin{correction}
\begin{lstlisting}
def initialisation_contour(V0):
    """ Initialise les bords de la grille a  la valeur V0 """
    V[:,0] = V0
    V[:,-1] = V0
    V[0,:] = V0
    V[-1,:] = V0
    B[:,0] = True
    B[:,-1] = True
    B[0,:] = True
    B[-1,:] = True
\end{lstlisting}
\end{correction}


Exécuter cette fonction et, pour vérifier son bon fonctionnement, afficher les bords du domaine grâce à la fonction \lstinline$graphe_bords()$ que l'on trouve dans la rubrique des outils graphiques, au début du script fourni, et qui affiche en noir les points appartenant au bord du domaine et en blanc les autres (avec un dégradé grisé entre les deux). Choisir par exemple $N_x = N_y = 60$ pour obtenir un résultat visuel.

\begin{correction}
Dans la corps de programme principal

\begin{lstlisting}
# Affichage du bord 
initialisation_contour(0)
graphe_bords()
\end{lstlisting}
\end{correction}

\section{Algorithme de Jacobi}

\subsection{Principe}

L'algorithme de Jacobi permet une résolution itérative de l'équation (\ref{EqLaplace}) : la valeur du potentiel $V_n[i,j]$ au
point $[i,j]$ à l'itération $n$ se déduit de la valeur des potentiels voisins à l'itération $n-1$ par

\begin{equation}
V_n[i,j]=\frac{V_{n-1}[i+1,j]+V_{n-1}[i-1,j]+V_{n-1}[i,j+1] + V_{n-1}[i,j-1]}{4} 
\label{EqIterative}
\end{equation}

Ainsi, la valeur du potentiel en tout point est recalculée à chaque itération en fonction de ce qu'elle aurait dû être si l'itération précédente vérifiait l'équation de Laplace. La convergence de la méthode est démontrable mathématiquement : au bout d'un nombre finies d'itérations, on est assuré d'approcher de la solution exacte et les nouvelles itérations ne font que très peu  évoluer le potentiel.

Le nombre d'itérations dépend de la précision souhaitée. Dans ce travail, on utilise un critère basé sur l'écart quadratique moyen

\begin{equation}
e=\sqrt{\frac{1}{N_xN_y}\sum_{i,j}(V_{n+1}(i,j))-V_n(i,j))^2} 
\label{EqEcartMoyen}
\end{equation}

qui représente l'écart entre la nouvelle et l'ancienne valeur du potentiel moyenné sur toute la grille. On choisit de
stopper la simulation lorsque $e$ devient inférieur à une valeur $\varepsilon$ définie au début de la simulation, en fonction d'un
compromis entre précision des résultats et temps de calcul.


	\subsection{Mise en pratique}
	
L'algorithme de Jacobi est le suivant :
\begin{description}
	\item[Initialisation] :
	
	$\rightarrow$ créer le tableau \lstinline$B$ et le remplir avec des \lstinline$True$ et des \lstinline$False$ pour décrire les bords du domaine
	
	$\rightarrow$ créer le tableau \lstinline$V$ et l'initialiser avec les valeurs voulues sur les bords du domaine.
	
	\item[Itérations] :	actualiser le tableau \lstinline$V$ : pour tout point $(i,j)$ n'appartenant pas au bord la nouvelle valeur est calculée selon la relation de récurrence (\ref{EqIterative}).
	
	\item[Terminaison] : calculer après chaque itération l'écart $e$ à la précédente avec l'équation (\ref{EqEcartMoyen}) et cesser la procédure lorsque $e$ devient inférieur à $\varepsilon$.
\end{description}

\quest Écrire une fonction \lstinline$ecart(V1,V2)$ prenant en argument deux tableaux \lstinline$numpy$ et renvoyant l'écart entre ces deux tableaux au sens de l'équation (\ref{EqIterative}).

\begin{correction}
\begin{lstlisting}
def ecart(V1,V2):
    """ les tableaux ne sont pas forcément carrés
    """
    Nx = V1.shape[0]
    Ny = V1.shape[1]
    diff_carre = 0
    for i in range(Nx):
        for j in range(Ny):
            diff_carre += (V1[i,j]-V2[i,j])**2
    return (diff_carre/(Nx*Ny))**(1/2)

\end{lstlisting}
\end{correction}

\quest Écrire une fonction \lstinline$iteration_jacobi()$ qui effectue une itération de l'algorithme ci-dessus. Cette fonction ne
prend aucun argument et doit renvoyer l'écart $e$ entre les potentiels avant et après itération.

\begin{correction}
\begin{lstlisting}
def iteration_jacobi():
    '''renvoie le nouveau potentiel actualise si le point fait parti
       du domaine D a l exclusion du bord '''
    V_copie = V.copy()
    for i in range(Nx):
        for j in range(Ny):
            if B[i,j]: # si [i,j] est un point du bord ... ne rien faire
                V[i,j] = V_copie[i,j]
            else:
                V[i,j] = (V_copie[i+1,j]+V_copie[i-1,j]+V_copie[i,j+1]+V_copie[i,j-1])/4
    return ecart(V_copie,V)

\end{lstlisting}
\end{correction}


\textbf{Rappels utiles :}

Pour effectuer une copie d'un tableau \lstinline$numpy$ on utilise \lstinline$V_copie = V.copy()$ mais pas \lstinline$V_copie = V$,
sans quoi les deux sont liés et toute modification de \lstinline$V$ entraîne une modification de \lstinline$V_copie$.

Avoir choisi un tableau de booléens pour \lstinline$B$ permet d'écrire directement les tests sous la forme \lstinline$if B[i,j]:$, qui est exactement équivalent à écrire \lstinline$if B[i,j] == True:$, et même \lstinline$if B[i,j] == 1:$ 
car en Python les booléens sont équivalents à des 0 et des 1.

\quest Écrire une fonction \lstinline$jacobi(eps)$ qui prend en argument un flottant \lstinline$eps$ et qui itère l'algorithme tant que l'écart
est supérieur à \lstinline$eps$. À des fins de comparaison avec la méthode suivante, votre fonction doit renvoyer le nombre
d'itérations. Pour suivre la convergence en temps réel, on pourra ajouter un \lstinline$print(e)$ à un endroit bien choisi.

\begin{correction}
\begin{lstlisting}
def jacobi(eps):
    ''' actualise le tableau V donc inutile de le renvoyer
    '''
    compteur = 0
    while iteration_jacobi() > eps:
        iteration_jacobi()
        compteur += 1
        print(iteration_jacobi() - eps)
    return compteur

\end{lstlisting}
\end{correction}


	\subsection{Application au condensateur plan de taille fini}

\quest Représenter les surfaces équipotentielles d'un condensateur plan de taille finie, en procédant de la façon suivante :

\begin{enumerate}
	\item Choisir une grille de taille $60 \times 60$ et créer les tableaux \lstinline$B$ et \lstinline$V$ dont toutes les valeurs sont nulles.
	\item Initialiser le contour de la grille à un potentiel nul
	\item Ajouter un condensateur dont les armatures sont aux  potentiels $\pm 10$ V avec la fonction \lstinline$initialisation_condensateur(V1,V2)$ qui prend en argument \lstinline$(V1,V2)$ respectivement les tensions de la plaque inférieure et supérieure.
	\item Résoudre l'équation de Laplace par l'algorithme de Jacobi en choisissant comme critère de terminaison $\varepsilon=1\cdot 10^{-3}$.
	\item 	Afficher les surfaces équipotentielles avec la fonction \lstinline$graphe_equipot()$.
	
\end{enumerate}

\begin{correction}
Dans le corps principal du programme. 
N'oubliez pas de mettre en commentaire (CTRL 3 Commenter/Décommenter) la partie du programme principal qui n'est pas utile.
 
\begin{lstlisting}
initialisation_contour(0)
initialisation_condensateur(-10,10)
print(jacobi(1e-3))
graphe_equipot()
\end{lstlisting}
\end{correction}


\quest Identifier les effets de bord et la distance sur laquelle ils se manifestent. Comme vous le savez les équipotentielles d'un condensateur plan infini sont équidistantes, parallèles entre elles, et parallèles aux armatures du condensateur.
\begin{correction}
\begin{center}
\includegraphics[width=10 cm]{Fig_condensateur.pdf}
\caption{Équipotentielles d'un condensateur}
\end{center}

Les effets de bord se font sentir, près des bords où les équipotentielles s'écartent. 
\end{correction}



\section{Algorithme de Gauss-Seidel adaptatif}

On peut montrer que la complexité de l'algorithme de Jacobi pour une grille de taille $N \times N$ est en $O(N^3)$, ce qui est néfaste pour les temps de calcul et devient rapidement gênant, même pour des géométries simples. L'algorithme de Gauss-Seidel adaptatif est une amélioration possible, mais moins immédiatement compréhensible.

	\subsection{Principe}
	
La méthode consiste essentiellement à remplacer la relation de récurrence  (\ref{EqIterative}) par

\begin{equation}
V_n[i,j]=(1-\omega)V_{n-1}[i,j]+\omega\times \frac{V_{n-1}[i+1,j]+V_n[i-1,j]  V_{n-1}[i,j+1]+V_n[i,j-1]}{4}
\label{EqGaussSeidel}
\end{equation}


Il y a donc deux différences majeures par rapport à l'algorithme précédent.

D'une part, on utilise la valeur $V_{n-1}[i, j]$ de $V$ au point $(i,j)$ à l'itération précédente, pondérée avec un poids $\omega$.

On peut montrer que pour une grille carrée de taille $N \times N$ il existe une valeur optimale $$ \omega_{r opt}=\frac{2}{1+\frac{\pi}{N}}$$

D'autre part, le calcul fait intervenir les valeurs du potentiel aux points voisins à l'itération précédente $n-1$ comme c'était déjà le cas, mais aussi certaines valeurs à l'itération $n$ elle-même. Ceci est rendu possible par le fait que la grille est parcourue à $i$ et $j$ croissants, ces valeurs ont donc déjà été actualisées lors de l'itération $n$ : on
comprend ainsi que cela améliore la vitesse de convergence. En pratique, il ne faut plus utiliser un tableau copie mais
directement itérer $V$ avec l'équation (\ref{EqGaussSeidel}).

\quest Écrire les fonctions \lstinline$iteration_gauss_seidel()$ et \lstinline$gauss_seidel(eps)$ qui jouent un rôle analogue aux fonctions du même nom écrites pour l'algorithme de Jacobi.

\begin{correction}
\begin{lstlisting}
def iteration_gauss_seidel():
    '''renvoie le nouveau potentiel actualise si le point fait parti
       du domaine D a l exclusion du bord '''
    V_copie = V.copy()
    omega = 2/(1+np.pi/Nx)
    for i in range(Nx):
        for j in range(Ny):
            if B[i,j]: # si [i,j] est un point du bord ... ne rien faire
                V[i,j] = V_copie[i,j]
            else:
                V[i,j] = (1-omega)*V_copie[i,j]+omega*(V_copie[i+1,j]+V[i-1,j]+V_copie[i,j+1]+V[i,j-1])/4
    return ecart(V_copie,V)
\end{lstlisting}
\end{correction}
\begin{correction}
\begin{lstlisting}
def gauss_seidel(eps):
    ''' actualise le tableau V donc inutile de le renvoyer
    '''
    compteur = 0
    while iteration_gauss_seidel() > eps:
        iteration_gauss_seidel()
        compteur += 1
        print(iteration_gauss_seidel() - eps)
    return compteur
\end{lstlisting}
\end{correction}


\quest Tester ces fonctions sur l'exemple précédent du condensateur. Vérifier que le nombre d'itérations est inférieur.

Pour que la comparaison ait un sens, penser à réinitialiser la matrice $V$ entre les deux simulations.

\begin{correction}
Pour l'algorithme de Jacobi, on trouve 166 itérations contre 27 pour Gauss-Seidel.

\begin{lstlisting}
initialisation_contour(0)
#graphe_bords()
initialisation_condensateur(-10,10)
print(gauss_seidel(1e-3))
graphe_equipot()
\end{lstlisting}

\end{correction}


	\subsection{Taille du domaine de calcul}
	
La taille du domaine de calcul est un point crucial dans les simulations numériques de ce type. Il est nécessaire, pour la simulation, de fixer la valeur du potentiel sur tout le contour de la grille. Cette nécessité de la modélisation ne doit pas ou peu influer sur le résultat de la simulation dans la zone intéressante.

On étudie cet effet dans le cas de domaines carrés de taille $N \times N$, les différentes valeurs étudiées étant stockées dans une liste \lstinline$lst_N$. En guise d'illustration, on étudie l'influence de la valeur de $N$ sur le potentiel en un point de référence choisi arbitrairement, légèrement à l'extérieur du condensateur. Les valeurs du potentiel en ce point pour les différentes valeurs de $N$ sont stockées dans une liste \lstinline$Vref$. On ne pourra prendre que des valeurs $N \geq  30$ compte tenu du point choisi.

\quest Compléter la fonction \lstinline$influence_taille_domaine(lst_N)$, qui renvoie la liste \lstinline$Vref$. On notera que quelques précautions s'imposent avec les variables choisies comme étant globales dans les simulations précédentes et dont cette fonction modifie les valeurs.

\begin{correction}
Comme nous devons redéfinir les grilles à chaque fois, les variables ont été précisées en entête de la fonction. La priorité sera donc donné aux valeurs définies dans la fonction. 
\begin{lstlisting}
def influence_taille_domaine(lst_N):
    global B, V, Nx, Ny
    Vref = []
    for N in lst_N:
        print(N)
        Nx, Ny = N, N
        V = np.zeros((Nx,Ny))
        B = np.zeros((Nx,Ny), dtype=bool)
        initialisation_contour(0)
        initialisation_condensateur(-10,10)
        gauss_seidel(1e-3)
        Vref.append(V[Nx//2+2,Ny//2+14])
    return Vref
\end{lstlisting}
\end{correction}


\quest Représenter le potentiel au point de référence en fonction de la taille du domaine de calcul.

\begin{correction}
\begin{lstlisting}
lst_N=[30,40,50,60,70,80,90,100,110,120,150]
Pot = influence_taille_domaine(lst_N)

plt.figure()
plt.plot(lst_N,Pot,marker = '+')
plt.xlabel('N : taille du domaine')
plt.ylabel('Potentiel au même point')
plt.show()
\end{lstlisting}
\end{correction}

\begin{correction}
\begin{center}
\includegraphics[width=12 cm]{Fig_Influence_N.pdf}
\caption{Influence de la taille du domaine sur la convergence du potentiel en un point}
\end{center}
\end{correction}

	\subsection{Application à l'effet de pointe}

		\textbf{Introduction}\\
		
En l'absence de perturbations, le potentiel électrostatique augmente quasi-linéaire dans l'atmosphère entre le sol -- sol qui est à un potentiel \lstinline$Vsol$ que l'on prendra égal à 0 V) -- et l'ionosphère. Dans une atmosphère non perturbée, par temps calme, le gradient de potentiel est de 100 V $\cdot$ m$^{-1}$ environ.

L'ajout d'une tige verticale conductrice, un paratonnerre par exemple, vient perturber cette situation : ses propriétés
conductrices font qu'elle est au même potentiel en tout point, égal à celui du sol auquel elle est reliée. 

Ceci va
déformer les surfaces équipotentielles, et donc modifier le champ électrique. On s'attend à ce que ce dernier augmente
vers la pointe : c'est ce que l'on nomme l'effet de pointe, et qui explique que la foudre tombe préférentiellement sur
les objets pointus.


On simule ici une tige verticale de hauteur $H$ et de largeur $D$, dont la pointe est une demi-sphère de diamètre $D$,
voir figure (\ref{FigParatonerre}). 

Plus précisément, notre programme étant bidimensionnel, nous simulons en fait un plan de hauteur $H$, épaisseur $D$, de longueur très grande dans la direction perpendiculaire au plan de la simulation.

\begin{center}
\schema{14cm}{Fig_Paratonerre}
\caption{Allure qualitative des équipotentielles en présence d’une pointe.}
\label{FigParatonerre}
\end{center}

L'objectif est d'obtenir numériquement une relation entre la valeur $E_{\rm max}$ du champ électrique et le diamètre $D$
de la pointe de la tige.\\

		\textbf{Calcul du potentiel}\\
		
La fonction \lstinline$initialisation_effet_pointe(Vmin,Vmax,h,d)$ fournie dans le script de départ permet d'initialiser un domaine avec : 

\begin{enumerate}
	\item un contour dont le bord bas est au potentiel $V_{\rm min}$, le bord haut au potentiel $V_{\rm max}$ et où le potentiel augmente
linéairement de $V_{\rm min}$ à $V_{\rm max}$ sur les bords droit et gauche. 
	\item une tige de hauteur $h$ et épaisseur $d$ (en nombre de points), au potentiel fixé $V_{\rm min}$. Attention, le nombre $d$ doit être impair !
\end{enumerate}


On raisonnera sur un pas de grille $\delta = 3$ cm, et on prendra : 

\begin{enumerate}
	\item  $h = 50$ et $d = 13$, ce qui correspond à une tige de hauteur $H = 1,5$ m et de largeur $D = d \times \delta = 30$ cm 
	\item un domaine de taille 120 $\times$ 120 points, ce qui correspond à 3,6 cm $\times$ 3,6m
	\item un potentiel $V_{\rm min} = 0$ V au niveau du sol et $V_{\rm max} =  N_x \times \delta \times 100 $ V en haut du cadre, ce qui donne bien un gradient au repos de 100 Vm$^{-1}$.
\end{enumerate}


\quest Initialiser le problème à l'aide de cette fonction. Vérifier que tout est correct à l'aide de \lstinline$graphe_bords()$.

\begin{correction}
\begin{lstlisting}
delta = 3
Nx = 120
Ny = 120
h = 50
d = 13

Vmin =0
Vmax = Nx * delta * 100

V = np.zeros((Nx,Ny))
B = np.zeros((Nx,Ny), dtype=bool)

initialisation_effet_pointe(Vmin,Vmax,h,d)
graphe_bords()
\end{lstlisting}
\end{correction}

\quest Utiliser ensuite l'algorithme de Gauss-Seidel pour obtenir le potentiel. On prendra encore $\varepsilon = 10^{-3}$. Visualiser
la solution à l'aide de \lstinline$graphe_equipot()$.\\

\begin{correction}
\begin{lstlisting}
initialisation_effet_pointe(Vmin,Vmax,h,d)
#graphe_bords()
gauss_seidel(1e-3)
graphe_equipot()
\end{lstlisting}
\end{correction}

\begin{correction}
\begin{center}
\includegraphics[width=10 cm]{Fig_EffetPointe.pdf}
\end{center}
\end{correction}

		\textbf{Calcul du champ électrique}\\

\quest Sur les équipotentielles précédentes, identifier la zone où le champ est le plus intense.\\


Pour calculer le champ électrique, on évaluera les dérivées premières par des relations du type 

\[\frac{\partial V}{\partial x}(x_i,y_j)=\frac{V(x_i+\delta,y_j)-V(x_i-\delta,y_j)}{2\delta} \]


qui donne une précision sur la dérivée en $O(\delta^2)$, d'un ordre meilleure que celle du schéma d'Euler habituel.

\quest - Écrire une fonction \lstinline$calcul_champ(delta)$ qui prend comme argument le pas de la grille et renvoie un tableau
\lstinline$norme_E$ de taille $N_x \times N_y$ qui contient la norme du champ électrique. On ne cherchera pas à calculer le champ sur
les bords du domaine.

\begin{correction}
\begin{lstlisting}
def calcul_E(delta):
    Ex = np.zeros((Nx,Ny))
    Ey = np.zeros((Nx,Ny))
    norme_E = np.zeros((Nx,Ny))
    for i in range(1,Nx-1): # on évite les bords
        for j in range(1,Ny-1): # on évite les bords 
            Ex[i,j] = -(V[i+1,j]-V[i-1,j])/(2.*delta)
            Ey[i,j] = -(V[i,j+1]-V[i,j-1])/(2.*delta)
            norme_E = (Ex**2+Ey**2)**0.5
    return Ex, Ey, norme_E
\end{lstlisting}
\end{correction}

\quest Appliquer cette fonction dans le cas précédent : par combien la présence de l'obstacle multiplie-t-elle le champ ?

\begin{correction}
\begin{lstlisting}
Ex, Ey, normeE = calcul_E(delta)
Emax = np.max(normeE)

# on trouve 530 .. donc fois 5.3
\end{lstlisting}

\end{correction}


\quest Reproduire la procédure pour quelques valeurs de $d$, par exemple avec $d \in [7, 13, 17, 23]$ et représenter le maximum de la norme
de $\overrightarrow{E}$, $E_{\rm max}$ en fonction de $d$. Ceci va-t-il dans le sens de l'effet de pointe mentionné précédemment ?

\begin{correction}
\begin{center}
\includegraphics[width=12 cm]{Fig_Max_E.pdf}
\caption{$E_{\rm max}$ en fonction de $d$, diamètre de la pointe}
\end{center}

\begin{center}
\includegraphics[width=12cm]{Fig_E_V_EffetPointe.pdf}
\caption{Lignes de champ et équipotentielles autour d'une pointe}
\end{center}
\end{correction}



\end{document}