diff --git a/labo/F10/TemplateTPII.pdf b/labo/F10/TemplateTPII.pdf index 0f54d46..a4d0f1c 100644 Binary files a/labo/F10/TemplateTPII.pdf and b/labo/F10/TemplateTPII.pdf differ diff --git a/labo/F10/TemplateTPII.synctex.gz b/labo/F10/TemplateTPII.synctex.gz new file mode 100644 index 0000000..84e2fbb Binary files /dev/null and b/labo/F10/TemplateTPII.synctex.gz differ diff --git a/labo/F10/TemplateTPII.tex b/labo/F10/TemplateTPII.tex index f11b9b7..eb70edd 100644 --- a/labo/F10/TemplateTPII.tex +++ b/labo/F10/TemplateTPII.tex @@ -1,539 +1,539 @@ \documentclass[a4paper, 12pt,oneside]{article} %On peut changer "oneside" en "twoside" si on sait que le résultat sera recto-verso. %Cela influence les marges (pas ici car elles sont identiques à droite et à gauche) % pour l'inclusion de figures en eps,pdf,jpg,.... \usepackage{graphicx, caption} \usepackage{gensymb} \usepackage{subcaption} %Marges. Désactiver pour utiliser les valeurs LaTeX par défaut %\usepackage[top=2.5cm, bottom=1.5cm, left=2cm, right=2cm, showframe]{geometry} \usepackage[top=2.5cm, bottom=1.5cm, left=2cm, right=2cm]{geometry} % quelques symboles mathematiques en plus \usepackage{amsmath} % chemie \usepackage{mhchem} % quelques symboles mathematiques en plus \usepackage{amsmath} \usepackage{multirow} \usepackage{mathtools} \usepackage{gensymb} \usepackage{url} \usepackage{makecell} % absolute value \DeclarePairedDelimiter\abs{\lvert}{\rvert} % le tout en langue francaise %\usepackage[francais]{babel} % on peut ecrire directement les charactères avec l'accent \usepackage[T1]{fontenc} \usepackage[french]{babel} % a utiliser sur Linux/Windows %\usepackage[latin1]{inputenc} % a utiliser avec UTF8 \usepackage[utf8]{inputenc} %Très utiles pour les groupes mixtes mac/PC. Un fichier texte enregistré sous codage UTF-8 est lisible dans les deux environnement. %Plus de problème de caractères accentués et spéciaux qui ne s'affichent pas % a utiliser sur le Mac %\usepackage[applemac]{inputenc} % pour l'inclusion de liens dans le document (pdflatex) \usepackage[colorlinks,bookmarks=false,linkcolor=black,urlcolor=blue, citecolor=black]{hyperref} %Pour l'utilisation plus simple des unités et fractions \usepackage{siunitx} %Pour utiliser du time new roman... Comenter pour utiliser du ComputerModern %\usepackage{mathptmx} %graphs \usepackage{tikz} \usepackage{pgfplots} \usepackage{pgfplotstable} \usetikzlibrary{patterns} \usepgfplotslibrary{external} \pgfplotsset{compat=newest} \usepackage{float} %Pour du code non interprété \usepackage{verbatim} \usepackage{verbdef}% http://ctan.org/pkg/verbdef %Pour changer la taille des titres de section et subsection. Ajoutez manuellement les autres styles si besoin. \makeatletter \renewcommand{\section}{\@startsection {section}{1}{\z@}% {-3.5ex \@plus -1ex \@minus -.2ex}% {2.3ex \@plus.2ex}% {\normalfont\normalsize\bfseries}} \makeatother \makeatletter \renewcommand{\subsection}{\@startsection {subsection}{1}{\z@}% {-3.5ex \@plus -1ex \@minus -.2ex}% {2.3ex \@plus.2ex}% {\normalfont\normalsize\bfseries}} \makeatother %pgfplot setup \makeatletter \pgfplotsset{ /pgfplots/flexible xticklabels from table/.code n args={3}{% \pgfplotstableread[#3]{#1}\coordinate@table \pgfplotstablegetcolumn{#2}\of{\coordinate@table}\to\pgfplots@xticklabels \let\pgfplots@xticklabel=\pgfplots@user@ticklabel@list@x }, % layer definition layers/my layer set/.define layer set={ background, main, foreground }{ % you could state styles here which should be moved to % corresponding layers, but that is not necessary here. % That is why wo don't state anything here }, % activate the newly created layer set set layers=my layer set } \makeatother %Début du document \begin{document} \title{\normalsize{Lab Work Report - Group N$^\circ$\\ XX - Experiment}} \date{\normalsize{\today}} \author{\normalsize{Name} 1\and \normalsize{Name 2}} %Crée la page de titre %\maketitle %Ajoute la table des matières %\tableofcontents %Début du rapport à la page suivante %\newpage %De manière à ce que template latex ressemble au mieux au template word, on empêche latex de créer la page de titre et la créons à la main %En taille de police 12, la commande \large donne une taille de police 14 %On utilise la commande \sffamily pour créer des caractères sans-serif \begin{center} \large\textbf{\sffamily Experiment N$^\circ$ F10: Fourier optics}\\% \large\sffamily Group N$^\circ$16: Ancarola Raffaele, Cincotti Armando\\% \large\sffamily \today\qquad Tarrago Velez Santiago\\%TODO cambiare nome assistente \end{center} % Introduction \section{Introduction} Based on wave physics, Fourier optics is a much useful tool than geometrical optics to study light propagation phenomenons, as it takes in account all wave properties of light. It permits in fact to study perfectly complex phenomenon as diffraction, and it is at the basis for developping useful technique for image and signal treatments. In the following report, as an example of direct application of this field in physics, it is shown that using convergent lenses it is possible to obtain the Fourier transform of a complex light signal, and that the latter can be easily treated by filtering it directly on its Fourier transform. \section{A brief Theroretical basis} \subsection{Principles of geometrical optic}\label{prince} Optics of thin convergent lenses can be studied in a first approximation by usig two geometrical principles \cite{ref1}: \begin{enumerate} \item Each beam parallel to the main axis of the lens, will refract so that it passes really or virtually through the focal point. \item Each beam passing through the center of the lens, won't get deflected. \end{enumerate} \subsection{Obtaining a Fourier Transform with a lens} The optical systems that will be studied in the following report can be schematised as : \begin{minipage}{\textwidth} \vspace{0.3cm} \hspace{-0.8cm} \begin{minipage}{0.45\textwidth} An object with a certain transmittance function $f(x,y)$, emitting a light signal $U(x,y,z)$ which creates an image $g(x,y) = U(x,y,d)$ on a plane screen positionned at a given distance $z = d$. One can show that a convergent lens with a given focal plane at distance $f$ permit to obtain the Fourier Transform $F(\nu_x, \nu_y)$ of the signal emitted in $z=0$, $U(x,y,0)=f(x,y)$, where $\nu_x$ and $\nu_y$ are the spatial frequencies of the signal in directions $\vec{e}_x$ and $\vec{e}_y$. \end{minipage} %\hspace{0.005\textwidth} \begin{minipage}{0.65\textwidth} \includegraphics[width=\textwidth]{simple_transform.png} \captionof{figure}{Plane wave composed by one spatial frequency projected in one point of the Fourier plane.} \label{simple_trans} \end{minipage} \vspace{0.3cm} \end{minipage} In fact, the luminous intensity of the signal at a distance $f$ from the lens he passed trough, is given by : \begin{equation}\label{I} I(x,y) = \frac{1}{(\lambda f)^2} \abs*{F(\frac{x}{\lambda f}, \frac{y}{\lambda f})} \\ \end{equation} where $\lambda$ is the wave lenght of the signal. Equation (\ref{I}) proove in fact that the Fourier Transform can be ``projected'' on the focal plane of the lens (see exemple in Figure \ref{simple_trans}). One can also show that, considering the inverse Fourier Transform $f(x,y) = F^{-1}(\nu_x,\nu_y)$, if the screen is now in front of a second lens, on its focal plane, then $g(x,y)$ follows this relation : \begin{equation}\label{g} g(-x,-y) \propto F^{-1}(\frac{x}{\lambda f}, \frac{y}{\lambda f}) = f(x,y) \end{equation} This means that, appling a second convergent lens to a signal already transformed by one lens, creates a real inververse image (because of the $-$ sign) of the emitted signal. {\bf N.B.} the whole mathematical procedure for the computation of equations (\ref{I}) and (\ref{g}) can be found in the {\it notice} for this report \cite{ref2}. \paragraph{About the Fourier Plane.} A Fourier Transform qualitatively consist on the projection of a signal (here $f(x,y)$) on each of the vector from its Fourier basis. In fact, the Fourier Transform decomposes the signal in the frequencies that make it up. As just explained in this section, a lens consinst in this case in a sort of projection from the real space to the frequencies space. So at the focal plane, also called {\bf Fourier plane} because it's where the Fourier Transform {\it lays}, can be observed several points consisting on the set of frequencies that composes the signal $f(x,y)$. Each point on this plane is the projection of a certain frequency which can be measured using the following relation \cite{ref2} : \begin{equation}\label{nu} \nu_x = \frac{x}{\lambda f} = \frac{x^{'}}{\lambda f} \end{equation} where $x$ the distance in $\vec{e}_x$ direction between the point and the origin of the Fourier plane where $\nu_x = \nu_y = 0$, $x^{'}$ is the same distance but mesured on a zoomed image of the Fourier Plane, and where $\gamma$ is the extension coefficient which will be defined in the next section. The equation is similar for $\nu_y$. \subsection{Projection of a real image.} When a real image is projected on a plane in front of the system of lenses and luminous source, the projection can be smaller or bigger than the image from the source. Let $d$ be the distance between one lens and a given object, and $d^{'}$ the distance between the same lens and the projected image of this object, then an extension coefficient $\gamma$ is geometrically given as : \begin{equation}\label{gamma} \gamma = \frac{d^{'}}{d} \end{equation} Having $\gamma$, then the following relation is maintained \cite{ref1}: \begin{equation}\label{d} \frac{1}{d} + \frac{1}{d^{'}} = \frac{1}{f} \quad \Rightarrow \quad d^{'} = \frac{df}{d - f} \end{equation} where $f$ is the distance between the lens and its focal plane. \section{Experimental setup} To study an optical image and obtain its Fourier transform, one can start by set up the {\bf 4f configuration} illustrated in Figure \ref{4f}. All the piecies are mounted on one long rail and can be moved anytime above it. %TODO inserisci immagine schema_esperimento.png \begin{minipage}{\textwidth} \vspace{0.3cm} \hspace{-0.8cm} \begin{minipage}{0.45\textwidth} It is first needed an image, so a green laser beam with wavelenght $\lambda = 532$ \si{\nano\metre} followed by a beam expander with focal distance $f_1 = 12$ \si{\milli\metre} are first mounted on the rail to obtain a large beam of light. Secondly, a lens with focal distance $f_2=100$ \si{\milli\metre} is positionned at distance $11.5$ \si{\centi\metre} from the laser beam, wich is approximatively the sum $f_1+f_2$ where there would be the combined focals of the two lenses. \end{minipage} %\hspace{0.005\textwidth} \begin{minipage}{0.65\textwidth} \includegraphics[width=\textwidth]{schema_esperimento.png} -\captionof{figure}{4 focals configuration for creating a Fourier plane and an inverse image.} +\captionof{figure}{4 focals configuration used to project a Fourier plane and an inverse image.} \label{4f} \end{minipage} \vspace{0.3cm} \end{minipage} As it follows by the first principle in section \ref{prince}, as the beam passes through the focal point, it then becames parallel to the rail. \begin{minipage}{\textwidth} \vspace{0.3cm} \hspace{-1cm} \begin{minipage}{0.5\textwidth} \includegraphics[width=\textwidth]{foto/house_reverse.png} -\captionof{figure}{Reverse image of Fourier's home} +\captionof{figure}{Reverse image of Fourier's home.} \label{rev_img} \end{minipage} \hspace{0.01\textwidth} \begin{minipage}{0.55\textwidth} Thirdly, a filter with the image of a house is positionned in front of the first lens in orther to get an image. This filter consist now in the object with transmittance $f(x,y)$ as it is studied in the therotical section and it will be called {\it Fourier's home} from now on. Then, following the schema in Figure \ref{4f}, one lens with focal distance $f_3 = 125$ \si{\milli\metre} is positionned at $f_3$ from the Fourier's home. At a distance of $f_3$ further can be found the Fourier plane. Finally $1,25$ \si{\centi\metre} further another lens with focal distance $f_3$ is mounted on the rail so that an inverse image of the house is projcted on a screen as it follows by relation (\ref{g}) (see Figure \ref{rev_img}). \end{minipage} \vspace{0.3cm} \end{minipage} \section{Results of the experiment and Discussions} \subsection{Lens law (\ref{d}) verification} To verify the law that appears on the leftside in (\ref{d}), it is first measured the leght of the real image in the focal plane of the second lens in the 4f configuration (see Figure \ref{4f}) which is $l = 0.5$ \si{\centi\metre}. Then for a fixed distance $d$ from a lens with $f=50$ \si{\milli\metre}, $d^{'}$ is measured by the rightside equation in (\ref{d}) and the screen is positionned at this distance from the lens. Measuring the lenght $l_{img}$ of the projected image, the experimental extension coefficient can be measured as $\gamma_{exp} = \frac{l_{img}}{l}$ and compared to the value of $\gamma_{ref}$ measured by using equation (\ref{gamma}). The error on all distance measurements is $\Delta x = 0.1$ \si{\centi\metre}. Table \ref{verification} shows several measurements for this experiment. It can be noticed that $\gamma_{ref}$ values do not differ from $\gamma_{ref}$ values which fall in the error interval. This shows that the leftside law in (\ref{d}) holds. \begin{table}[H] \centering \begin{tabular}{|c|c|c|c|c|} \hline $d$ [\si{\centi\metre}] & $d'$ [\si{\centi\metre}] & $l_{img}$ [\si{\centi\metre}] & $\gamma_{ref}$ & $\gamma_{exp}$ \\ \hline $7.0 $ & $17.5$ & $1.3$ & $2.5 \pm 0.1$ & $2.6 \pm 0.4$ \\ \hline $8.0 $ & $13.3$ & $0.9$ & $1.7 \pm 0.1$ & $1.8 \pm 0.5$ \\ \hline $9.0 $ & $11.3$ & $0.7$ & $1.3 \pm 0.1$ & $1.3 \pm 0.7$ \\ \hline $10.0$ & $10.0$ & $0.6$ & $1.0 \pm 0.1$ & $1.2 \pm 0.4$ \\ \hline $11.3$ & $9.0 $ & $0.5$ & $0.8 \pm 0.1$ & $1.0 \pm 0.4$ \\ \hline $13.3$ & $8.0 $ & $0.3$ & $0.6 \pm 0.1$ & $0.6 \pm 0.3$ \\ \hline $17.5$ & $7.0 $ & $0.2$ & $0.4 \pm 0.1$ & $0.4 \pm 0.3$ \\ \hline \end{tabular} \caption{Measurments of the extension coefficient $\gamma$ for different values of $d$} \label{verification} \end{table} \subsection{Fourier plane} The Fourier's home is an image composed by four quadrants, each one of them composed by parallel -lines following one direction per quadrant (see Figure \ref{fourier_house}). %TODO inserisci foto dove si vede la casa intera con le linee +lines following one direction per quadrant (see Figure \ref{fourier_house}). These lines form a pattern of different spatial frequencies in different direction: $\nu_x$, $\nu_y$ or a composition of both $\nu_x + \nu_y$. Between the two convergent lenses in configuration 4f, there is the Fourier Plane where each of these frequencies is ``projected on a spot'' on the plane. The image of the Fourier transform can be extended mounting a convergent lens in front of in order to obtain a bigger image as it is shown on Figure \ref{fourier_plane}. \begin{figure}[h] \centering %\hspace{-0.6cm} \begin{subfigure}{0.45\textwidth} \includegraphics[width=\textwidth]{foto/fourier_plane.pdf} \caption{Extended image of the Fourier plane.} %TODO \label{fourier_plane} \end{subfigure} \hspace{0.08\textwidth} \begin{subfigure}{0.45\textwidth} \includegraphics[width=0.9\textwidth]{foto/house_border.pdf} \caption{Extended image of Fourier's home.} %TODO \label{fourier_house} \end{subfigure} \caption{Fourier plane and Fourier's home sectors, highlighted per direction of frequencies in Fourier plane.} \end{figure} On Figure \ref{fourier_plane} can be observed that several points are allinged concentrically on 4 different directions. In fact each direction derives from one sector of Fourier's home, for exemple the vertical lines on Figure \ref{fourier_house} are separated in frequencies only on direction $\vec{e}_x$ so this sector ``projects'' the horizontal spots on the Fourier plane. \subsubsection{Verifing relation (\ref{nu}).} \begin{minipage}{\textwidth} \vspace{0.3cm} \hspace{-0.8cm} \begin{minipage}{0.45\textwidth} By counting the number of all visible lines for each sector of Fourier's home, the minimum non-zero frequency can be measured as $\nu_{i} = N/l$, where $N$ is the number of lines and $l$ the lenght of the sector when the image isn't extended. Then, by extending the image of the Fourier Plane, these frequencies can be measured by relation (\ref{nu}) measuring the distance $x^{'}$ -between one point in the first ``circle'' and the center. Just one measurement is needed because it can be verified that this distance -is equal for each point in this circle. -% TODO commenta poco qui +between one point in the first ``circle'' and the center. Table \ref{frequency} shows that these measurement don't coincide, but the order +of the two measurments is the same. This divergence can be explained by the errors in lenght measurements, that occurs even in +the determination of $\gamma$, affected also by the precision of the one that takes the measurments. \end{minipage} \hspace{0.05\textwidth} \begin{minipage}{0.4\textwidth} \includegraphics[width=\textwidth]{foto/house_elab.pdf} -\captionof{figure}{House image elaborated in order to clearly count lines} +\captionof{figure}{House image elaborated in order to clearly count lines.} \label{elab} \end{minipage} \vspace{0.3cm} \end{minipage} % TODO continua a commentare qui se ti prolunghi \begin{table}[H] \centering \begin{tabular}{cc|c|c|c|c|c|c|} \cline{2-8} \multicolumn{1}{ c| }{} & $N$ & $l$ [\si{\milli\metre}] & $x^{'}$ [\si{\centi\metre}] & $\gamma$ & $f$ [\si{\milli\metre}] & $\lambda$ [\si{\nano\metre}] & $\nu_x$ [\si{1/\kilo\metre}] \\ \hline % \multicolumn{1}{|c|}{Line counting} & \multirow{2}{*}{$24 \pm 1$} & \multirow{2}{*}{$2.5 \pm 0.5$} & \multirow{2}{*}{$0.5 \pm 0.1$} & \multirow{2}{*}{$7.25$} & \multirow{2}{*}{$50 \pm 1$} & \multirow{2}{*}{$532 \pm 1$} & $9.6 \pm 4.2$ \\ \cline{1-1} \cline{8-8} -\multicolumn{1}{|c|}{Circle radius} & & & & & & & $23.3 \pm 3.1$ \\ +\multicolumn{1}{|c|}{Relation \ref{nu}} & & & & & & & $23.3 \pm 3.1$ \\ %\hline \Xhline{4\arrayrulewidth} % \end{tabular} \caption{Frequency computed with two different methods} \label{frequency} \end{table} \subsubsection{Filtering the Fourier Plane.} Once obtained the Fourier Plane, one can try to filter light directly on this one and see -what appens to the real image reconstructed by the second lens in configuration 4f (Figure \ref{}) and then +what appens to the real image reconstructed by the second lens in configuration 4f and then extendend using a convergent lens with $f = 100$ \si{\milli\metre}. \paragraph{Filtering frequencies in one direction.} By using two filters composed by a diagonal or horizontal/vertical line, one can filter a whole set of frequencies -on one of the 4 directions illustrated in Figure \ref{}. This result in a real image where all of the spatial frequencies in +on one of the 4 directions illustrated in Figure \ref{fourier_plane}. This results in a real image where all of the spatial frequencies in that direction have been filtered. As in Fourier's home each sector is composed by lines following only one direction, this -means that this kind of filter hides one of the sectors letting pass light from the three others. -For exemple by filtering the points in $\vec{e}_{y}$ direction where $\nu_x = 0$ (see Figure \ref{}), -the sector with horizontal lines in the Fourier's home get hidden (see Figure \ref{}). +means that this kind of filter hides one of the sectors and let pass light from the three others. +For exemple by filtering the points in $\vec{e}_{x}$ direction where $\nu_y = 0$ (see Figure \ref{filter_dir}), +the sector with horizontal lines in the Fourier's home get hidden (see Figure \ref{filter_res}). \begin{figure}[h] \centering %\hspace{-0.6cm} \begin{subfigure}{0.45\textwidth} \includegraphics[width=\textwidth]{foto/miss_bottom_right.png} -\caption{Result on the house image} +\caption{Result on the house image.} \label{filter_res} \end{subfigure} \hspace{0.08\textwidth} \begin{subfigure}{0.45\textwidth} \includegraphics[width=\textwidth]{foto/fourier_horiz.png} -\caption{Filtering the horizontal direction} +\caption{Filtering the horizontal direction.} \label{filter_dir} \end{subfigure} -\caption{Application of a filter on the fourier plane and resulting visualization} +\caption{Application of a filter on the fourier plane and visualization of the result.} \label{filtering} \end{figure} \paragraph{Filtering a few or all except one frequencies.} With a filter consinsting on a flat black surface with one small hole, it is possible to chose only a small set of frequencies to let pass trough in the Fourier plane. Doing this it is possible to visualize only one spatial frequency that composes the signal coming from Fourier's Home. \begin{minipage}{\textwidth} \vspace{0.3cm} \hspace{-0.8cm} \begin{minipage}{0.35\textwidth} -An example can be observed in Figure \ref{}, or a more interesting one in Figure \ref{top_left_img} where only flat surfaces can be observed. This image was in fact obtained +An example can be observed in Figure \ref{top_left_img}, or a more interesting one in Figure \ref{lowpass_img} + where only quite flat surfaces can be observed. This image was in fact obtained by letting through the central point in the Fourier plane where $\nu_x = \nu_y = 0$, and the filter surface acts like a low-pass filter. In a similar way a large black point acts like a high-pass filter if it hides the central points -on Fourier plane. The result of the latter can be observed on Figure \ref{} where the lines that compose Fourier's home +on Fourier plane. The result of the latter can be observed on Figure \ref{elab} where the lines that compose Fourier's home are much neat than in other pictures. \end{minipage} \hspace{0.03\textwidth} \begin{minipage}{0.3\textwidth} \includegraphics[width=\textwidth]{foto/house_passebas.png} -\captionof{figure}{House image when a low-pass filter is applied} +\captionof{figure}{House image when a low-pass filter is applied.} \label{lowpass_img} \end{minipage} \hspace{0.03\textwidth} \begin{minipage}{0.3\textwidth} \includegraphics[width=\textwidth]{foto/top_left.pdf} -\captionof{figure}{Image shown letting pass a frequency only} +\captionof{figure}{What appears of Fourier's home letting pass one frequency only.} \label{top_left_img} \end{minipage} \vspace{0.3cm} \end{minipage} -% TODO aggiungere passe bas e top-left \paragraph{A brief discussion.} We've observed how it is possible to hide frequencies from a light signal. In fact it has been prooven that working on the fourier plane it is possible to create every type of frequency filter that is needed. This kind of manipulations on a signal are really usefull and common in the field of signal treatement. For example it directly follows from the report that this tecnique can be developped to create light filters for spatial frequencies in order to hide noise in images, or in radio techniques where superposition of frequencies and smart filtering permits to transmit audible signal over radio wave \cite{ref3}. \section{Conclusion} Throughout the different experiments and measurements, the geometrical relations for convergent lenses have been verified as the results were quite consistent with theory. The most interesting part is the projection of the real image on the frequency plane, also called here {\it Fourier plane}. This kind of manpulation is in fact at the basis of image processing as it permits to easily treat the light signal by filtering precise frequencies. As mentionned in the introduction, Fourier optics is in fact an usefull tool for studying light waves and treating light signals, and much more applications can be found in every kind of field as spectroscopy, holography, pattern recognition and so on \cite{ref4}. \section{Annexes} \subsection{About the error} \paragraph{Incertitude on distances} All measured distances taken with the ruler have the incertitude $\Delta x = 0.1$ \si{\centi\metre}. \begin{itemize} \item \textbf{Nota bene}: this incertitude is never shown on values written on tables because it's the same for all of them. \end{itemize} \paragraph{Incertitude on $\gamma$} \begin{itemize} \setlength\itemsep{0.7em} \item $\Delta \gamma_{ref} = \frac{\Delta d'}{d} + \frac{d'}{d^2} \Delta d$ \item $\Delta \gamma_{exp} = \frac{\Delta L}{l} + \frac{L}{l^2} \Delta l$ \end{itemize} where $L$ is the length of the projected image on the screen and $l$ is the real size of the house. \paragraph{Incertitude on $\nu_x$} \begin{itemize} \item $\Delta \nu_x(N, l) = \frac{\Delta N}{l} + \frac{N}{l^2} \Delta l$ \item $\Delta \nu_x(x^{'}, \lambda, f, \gamma) = \frac{\Delta x^{'}}{\lambda f \gamma} + \frac{x^{'}}{(\lambda f \gamma)^2} \cdot (\lambda \gamma \; \Delta f + f \gamma \; \Delta \lambda + f \lambda \; \Delta \gamma)$ \end{itemize} \section{Literature References} % Bibliographie \begin{thebibliography}{99} \bibitem{ref1} Dr. D. MARI, Dr. I. TKALCEC, Métrologie Lessons TP I, Optique I : {\it Optique géométrique, Réfraction, Lentilles} Faculté des sciences de base, Section de physique. \bibitem{ref2} Notice des TP de physique F10 : Optique de Fourier. \bibitem{ref3} HARALD, Brune, Cours de Physique IV, Faculté des Sciences de Base, Section de Physique, Ecole Polytechnique Fédérale de Lausanne. \bibitem{ref4} TYSON, Robert K. Principles and applications of Fourier Optics, Chapter 7: Practical Applications. University of North Carolina at Charlotte, USA, August 2014. \end{thebibliography} \end{document} diff --git a/physnum/rap6/rapport.synctex.gz b/physnum/rap6/rapport.synctex.gz new file mode 100644 index 0000000..b29f9f9 Binary files /dev/null and b/physnum/rap6/rapport.synctex.gz differ diff --git a/physnum/rap6/rapport.tex b/physnum/rap6/rapport.tex index f51966b..8e4d8d4 100644 --- a/physnum/rap6/rapport.tex +++ b/physnum/rap6/rapport.tex @@ -1,674 +1,674 @@ % debut d'un fichier latex standard \documentclass[a4paper,12pt,twoside]{article} % pour l'inclusion de figures en eps,pdf,jpg \usepackage{graphicx} % quelques symboles mathematiques en plus \usepackage{amsmath} % le tout en langue francaise \usepackage[francais]{babel} % on peut ecrire directement les caracteres avec l'accent % a utiliser sur Linux/Windows \usepackage[utf8]{inputenc} % a utiliser sur le Mac %\usepackage[applemac]{inputenc} % pour l'inclusion de links dans le document \usepackage[colorlinks,bookmarks=false,linkcolor=blue,urlcolor=blue]{hyperref} \usepackage{wrapfig} \usepackage{mathtools} \DeclarePairedDelimiter\abs{\lvert}{\rvert} % module pour les graphiques \usepackage{tikz} \usepackage{pgfplots} \usepackage{pgfplotstable} \usepackage{verbatim} % module pour positionner les images \usepackage{float} \usepackage{subcaption} \usepackage{siunitx} \usepackage{xcolor} %\usepackage[bottom = 0.5in]{geometry} \usepgfplotslibrary{external} \pgfplotsset{compat=newest} \tikzexternalize[prefix=tikzext/] \paperheight=297mm \paperwidth=210mm \setlength{\textheight}{235mm} \setlength{\topmargin}{-1.2cm} % pour centrer la page verticalement %\setlength{\footskip}{5mm} \setlength{\textwidth}{15cm} \setlength{\oddsidemargin}{0.56cm} \setlength{\evensidemargin}{0.56cm} \pagestyle{plain} % quelques abreviations utiles \def \be {\begin{equation}} \def \ee {\end{equation}} \def \dd {{\rm d}} \newcommand{\mail}[1]{{\href{mailto:#1}{#1}}} \newcommand{\ftplink}[1]{{\href{ftp://#1}{#1}}} %pgfplots setup \makeatletter \pgfplotsset{ /pgfplots/flexible xticklabels from table/.code n args={3}{% \pgfplotstableread[#3]{#1}\coordinate@table \pgfplotstablegetcolumn{#2}\of{\coordinate@table}\to\pgfplots@xticklabels \let\pgfplots@xticklabel=\pgfplots@user@ticklabel@list@x }, % layer definition layers/my layer set/.define layer set={ background, main, foreground }{ % you could state styles here which should be moved to % corresponding layers, but that is not necessary here. % That is why wo don't state anything here }, % activate the newly created layer set set layers=my layer set } \makeatother % % latex SqueletteRapport.tex % compile la source LaTeX % xdvi SqueletteRapport.dvi & % visualise le resultat % dvips -t a4 -o SqueletteRapport.ps SqueletteRapport % produit un PostScript % ps2pdf SqueletteRapport.ps % convertit en pdf % pdflatex SqueletteRapport.pdf % compile et produit un pdf % ======= Le document commence ici ====== \begin{document} % Le titre, l'auteur et la date \title{Physique Numérique I - II - Exercice 6} \date{\today} \author{Ancarola Raffaele, Cincotti Armando\\{\small \mail{raffaele.ancarola@epfl.ch}}, {\small \mail{armando.cincotti@epfl.ch}}} \maketitle \tableofcontents % Table des matieres % Quelques options pour les espacements entre lignes, l'identation % des nouveaux paragraphes, et l'espacement entre paragraphes \baselineskip=16pt \parindent=15pt \parskip=5pt \section{Introduction} %----------------------------------------- En physique, les systèmes étudiés les plus intéressant, souvent ne possèdent pas de solution détérminable analityquement. Pour ceci il est nécéssaire d'utiliser des schémas numériques et de discrétiser le problème pour pouvoir l'étudier correctement. Dans le rapport qui suit, un système éléctrodynamique sera étudié en utilisant la méthode des éléments finis pour la discrétisation du problème. Cette méthode sera donc expliquée et illustrée, ensuite testée et en fin employée sur un problème non trivial pour pouvoir en déduire différentes caractéristiques intéressantes. \section{Analyse du problème} \paragraph{Potentiel éléctrique.} Soit un diéléctrique dans un volume $\Omega$ avec bord $\partial\Omega$, il est cherché une fonction décrivant le potentiel eléctrique $\phi(\vec{x})$ tel qu'il est à $V_0$ sur $\partial\Omega$. Le milieu est de constante diéléctrique $\epsilon_r(\vec{x})$ et avec densité de charges libres $\rho_{lib}(\vec{x})$. Une des équations de Maxwell pour l'éléctrodynamique veut : \begin{equation}\label{maxwell} \nabla \cdot \vec{D}(\vec{x}) = \rho_{lib}(\vec{x}) \end{equation} où $\vec{D}(\vec{x}) =- \epsilon_0\epsilon_r(\vec{x})\nabla\phi(\vec{x})$ est le champs d'induction éléctrique, aussi appelé champs de déplacement. \subsection{Forme variationelle du problème.} En choisissant une fonction $\eta(\vec{x})$ quelconque de régularité $C^1$ et 0 sur $\partial\Omega$, et en utilisant les propritée véctorielle de $\nabla \cdot$, de l'équation (\ref{maxwell}) il est possible de déduire les rélations suivantes : \begin{align*} \nabla \cdot (\eta(\vec{x})\epsilon_r(\vec{x})\nabla\phi(\vec{x})) -\nabla\eta(\vec{x}) \cdot \epsilon_r\nabla \phi(\vec{x}) = \eta(\vec{x}) \nabla \cdot (\epsilon_r(\vec{x})\nabla\phi(\vec{x}))\\ \Rightarrow\quad \nabla\cdot( \eta(\vec{x})\epsilon_r(\vec{x}) \nabla\phi(\vec{x})) -\nabla\eta(\vec{x})\cdot \epsilon_r(\vec{x})\nabla\phi(\vec{x}) = - \eta(\vec{x})\frac{\rho_{lib}(\vec{x})}{\epsilon_0} \Leftrightarrow \\ \iiint \limits_{\Omega} \nabla\cdot(\eta(\vec{x})\epsilon_r(\vec{x}) \nabla\phi(\vec{x})) -\nabla\eta(\vec{x})\cdot \epsilon_r(\vec{x})\nabla\phi(\vec{x}))dV = -\iiint \limits_{\Omega} \eta(\vec{x})\frac{\rho_{lib}(\vec{x})}{\epsilon_0} dV\\ \Leftrightarrow \iint \limits_{\partial\Omega} \eta(\vec{x})\epsilon_r(\vec{x}) \nabla\phi(\vec{x}) \cdot d\vec{\sigma} - \iiint \limits_{\Omega} \nabla\eta(\vec{x})\cdot \epsilon_r(\vec{x})\nabla\phi(\vec{x})) dV = -\iiint \limits_{\Omega} \eta(\vec{x})\frac{\rho_{lib}(\vec{x})}{\epsilon_0} dV\\ \end{align*} Pour la dérnière équivalence il a été employé le théorème de la divergence. Comme par définition $\eta(\vec{x})$ est 0 sur $\delta\Omega$ on en déduit en fin : \begin{equation}\label{general} \iiint \limits_{\Omega} \nabla\eta(\vec{x})\cdot \epsilon_r(\vec{x})\nabla\phi(\vec{x})) dV = \iiint \limits_{\Omega} \eta(\vec{x})\frac{\rho_{lib}(\vec{x})}{\epsilon_0} dV \end{equation} \paragraph{Description du système.} Le système étudié ici voit $\Omega$ comme étant un domaine cylindrique de rayon $R$ et longeur $L_z$ , et par les symmétries cylindriques du système, les fonctions définies sur ce domaine dépendent uniquement de $r$, et on choisit $\eta = \eta(r)$. Donc en passant par un changement en cordonnées cylindriques, en sachant que la jacobiènne d'un tel changement a pour détérminant $r$, l'équation (\ref{general}) devient : \begin{equation*} \iiint \limits_{\Omega(r,\theta,z)} \epsilon_r(r) \frac{\partial\eta(r)}{\partial r}\frac{\partial\phi(r)}{\delta r}r drd\theta dz = \iiint \limits_{\Omega(r,\theta,z)} \eta(r)\frac{\rho_{lib}(r)}{\epsilon_0}r drd\theta dz \end{equation*} Comme les fonctions à l'intérieur de l'intégrale dépendent uniquement de r, cette équation peut être simplifiée pour obtenir : \begin{equation}\label{cylindriques} \int \limits_{0}^{R} \epsilon_r(r) \frac{\partial\eta(r)}{\partial r}\frac{\partial\phi(r)}{\delta r}r dr = \int \limits_{0}^{R} \eta(r)\frac{\rho_{lib}(r)}{\epsilon_0}r dr \end{equation} \subsection{Discretisation par la méthode des éléments finis.} Pour discrétiser le problème, l'on commence en partageant l'intervalle $[0,R]$ par des plus petits intavalles $[r_i,r_{i+1}]$, pour $i = 1,2,...n$. En suite les fonctions $\eta(r)$ et $\phi(r)$ sont approximées par une somme discrète de fonctions chapeau $\Lambda_i$ ayant valeur 1 sur $r_i$, linéaires dans l'intervalle $[r_{i-1}, r_{i+1}]$ et 0 ailleur. \begin{equation}\label{discr} \phi(r) = \sum_{k=1}^{n}\phi_k\Lambda_k(r) \; \quad \quad \; \eta(r) = \sum_{k=1}^n \eta_k\Lambda_k(r) \end{equation} Ainsi faisant, l'équation (\ref{cylindriques}) peut s'écrire comme suit : \begin{equation*} \sum_{i=1}^n \eta_i \sum_{j=1}^n (\int \limits_{0}^R \epsilon_r(r)r\frac{\partial \Lambda_i(r)}{\partial r}\frac{\partial \Lambda_j(r)}{\partial r} dr)\phi_j = \sum_{i=1}^n \eta_i (\int \limits_{0}^{R} \Lambda_i(r)\frac{\rho(r)}{\epsilon_0}r dr) \quad \quad \forall\;\; \eta(r) = \sum_{k=1}^n \eta_k\Lambda_k(r) \end{equation*} Comme cette équation est valide pour toute forme de $\eta(r)$, elle peut être simplifiée pour obtenir : \begin{equation}\label{eqdiscrete} \sum_{j=1}^n (\int \limits_{0}^R \epsilon_r(r)r\frac{\partial \Lambda_i(r)}{\partial r}\frac{\partial \Lambda_j(r)}{\partial r} dr)\phi_j = \int \limits_{0}^{R} \Lambda_i(r)\frac{\rho(r)}{\epsilon_0}r dr \quad \quad \forall i \end{equation} Il est intéressant en fin de rémarquer que l'équation (\ref{eqdiscrete}), en considérant i de 1 à $n+1$, est équivalente à l'équation vectorielle $A\vec{\phi}=\vec{b}$ où $A$ est une matrice carrée symmétrique dont les coéfficients sont : \begin{equation}\label{Aij} A_{i,j} = \int \limits_{0}^R \epsilon_r(r)r\frac{\partial \Lambda_i(r)}{\partial r}\frac{\partial \Lambda_j(r)}{\partial r} dr \end{equation} où par convetion l'indice i est indice de ligne et j de colonne. Les coéfficient de $\vec{b}$ sont donnés par : \begin{equation}\label{bi} b_i = \int \limits_{0}^{R} \Lambda_i(r)\frac{\rho(r)}{\epsilon_0}r dr \end{equation} \paragraph{Détermination de $A_{i.j}$ et $b_i$.} \begin{equation*} \end{equation*} L'équation (\ref{Aij}) peut s'écrire comme $A_{i,j} = \sum_{k=1}^n \int \limits_{r_k}^{r_{k+1}} \epsilon_r(r)r\frac{\partial \Lambda_i(r)}{\partial r}\frac{\partial \Lambda_j(r)}{\partial r} dr$. De cette forme, par la définition des fonctions $\Lambda_j$ il peut être facilement déduit que les seules coéfficients non nuls de la matrice A sont ceux pour lequels $j=i=k$, $j=i=k+1$, $j=k, i=k+1$ et $j=k+1 ,i=k$, car les dérivées de ces fonctions ne sont jamais non-nulles au même temps dans d'autres cas. De ça il est possible de déduire que la matrice $A$ est tridiagonale. Pour faciliter le calcul des intègrales, la fonctions $f(r) = \epsilon_r(r)r$ est approximée par une ``valeur moyenne'' $\epsilon_{r,k}^{*}$ sur l'intérval $[r_k, r_{k+1}]$, qui corréspond à la valeur par laquelle la fonction $f(r)$ serait approximée à l'intérieur de l'intégrale lorsque l'on applique la formule du trapèze mixte à celle du point milieu (voir équation (\ref{trapezio2})). Il ne reste plus que détérminer la valeur de l'intégrale pour les quatres cas non nuls. Soit $h_k = r_{k+1} - r_k$ la longuer de l'intérval $[r_k, r_{k+1}]$, alors $\Lambda_k = -\frac{1}{h_k}$ et $\Lambda_{k+1} =\frac{1}{h_k} $ sur cet interval, il en suit directement que $\int \limits_{r_k}^{r_{k+1}} \frac{\partial \Lambda_k(r)}{\partial r}\frac{\partial \Lambda_k(r)}{\partial r} dr = \frac{1}{h_k}$ et $\int \limits_{r_k}^{r_{k+1}} \frac{\partial \Lambda_k(r)}{\partial r}\frac{\partial \Lambda_{k+1}(r)}{\partial r} dr = -\frac{1}{h_k}$. Ceci implique donc les équations suivantes: \begin{align*} & A_{i,i} = \frac{\epsilon_{r,i-1}^{*}}{h_{i-1}} + \frac{\epsilon_{r,i}^{*}}{h_{i}} \;\; i=2,... n; & & A_{1,1} = \frac{\epsilon_{r,1}^{*}}{h_{1}} & A_{n+1,n+1} = \frac{\epsilon_{r,n}^{*}}{h_n} \\ & A_{i,i+1}= -\frac{\epsilon_{r,i}^{*}}{h_{i}}\; \; i = 1,... n; & & A_{i,i-1} = -\frac{\epsilon_{r,i-1}^{*}}{h_{i-1}}\;\; i = 2,... n+1; \end{align*} \begin{equation}\label{coefficienti} A_{i,j} = 0 \quad \forall j \notin \{i-1,i,i+1\} \end{equation} Pour la détérmintation des coéfficient $b_i$ le même raisonnement est fait, et l'on rémarque que l'intégrale à l'équation (\ref{bi}) est non-nulle seulement sur les interavalles $[r_{i-1},r_{i}]$ et $[r_i,r_{i+1}]$ pour $i = 2,... n$ , de même sur l'intérval $[r_1,r_2]$ pour $i=1$ et $[r_{n},r_{n+1}]$ pour $i=n+1$. L'intégrale dans ce cas est simplement approximée par la formule du trapèze mixte à celle du point milieu (\ref{trapezio}) donnant les coéfficient qui suivent : \begin{align*} b_i = h_{i-1}[p\frac{\rho_{lib}(r_i)r_i}{2\epsilon_0} + (1-p)\frac{\rho_{lib}(r_{i-1/2})r_{i-1/2}}{2\epsilon_0}] \\ +\;\;h_{i}[p\frac{\rho_{lib}(r_i)r_i}{2\epsilon_0} + (1-p)\frac{\rho_{lib}(r_{i+1/2})r_{i+1/2}}{2\epsilon_0}] \quad \; i = 2,... n. \end{align*} \begin{equation*} b_1 = h_1[p\frac{\rho_{lib}(r_1)r_1}{2\epsilon_0} + (1-p)\frac{\rho_{lib}(r_{1+1/2})r_{1+1/2}}{2\epsilon_0}] \end{equation*} \begin{equation}\label{bcoeff} b_{n+1} = h_{n}[p\frac{\rho_{lib}(r_{n+1})r_{n+1}}{2\epsilon_0} + (1-p)\frac{\rho_{lib}(r_{n+1/2})r_{n+1/2}}{2\epsilon_0}] \end{equation} où $r_{i +1/2} = \frac{r_{i+1} + r_i}{2}$ et $r_{i -1/2} = \frac{r_{i} + r_{i-1}}{2}$, et où $p \in [0,1]$ est choisi en configuration avec les autres paramètres constants du système. En fin, pour imposer la condition limite $V(R) = V_0$ dans l'équation véctoriel on pose en dérnière ligne $A_{n+1,n} = 0$, $A_{n+1,n+1} = 1$ et $b_{n+1} = V_0$. L'algoritme permettant en fin de détérminer la valeur de $\vec{\phi}$ dans l'équation $A \vec{\phi} = \vec{b}$ est décrit à page 87 du polycopier de cours \cite{ref2}. Par l'application de cet algorithme, il est possible d'obtenir le vecteur de composantes $\phi_i = \phi(r_i)$, pour avoir une {\it mappe} des valeurs de $\phi$ sur les points discrétisées du système. \subsection{Champs éléctrique} Lorsque le champs scalaire $\phi(r)$ à été détérminé par discrétisation en un vecteur $\vec{\phi}$, il est intéréssant de détérminer aussi les champs véctoriel dérivant de ce potentiel. Le but ici est de détérminer la méthode numérique permettant d'obtenir le champs de déplacement $D(r) = \epsilon_0\epsilon_r(r)E(r)$ avec $E(r) = -\nabla\phi = \frac{\partial \phi}{\partial r}(r)$. Par la méthode des éléments finis, comme décrit à la ligne (\ref{discr}) on obtient $E(r) =- \sum_{j} \phi_j \frac{\partial\Lambda_j}{\partial r}(r)$. Lorsque l'on cherche d'évaluer $E(r^{'})$ pour $r^{'} \in ]r_k,r_{k+1}[$ l'on observe que seulement deux fonctions de base $\Lambda_j$ contribuent à la somme ( $\Lambda_k$ et $\Lambda_{k+1}$), ce qui implique : \begin{equation}\label{D} D(r^{'})= \epsilon_0\epsilon_r(r^{'})E(r^{'})\;\;\;\; \text{ où } \;E(r^{'}) = \frac{\phi_{k+1} - \phi_k}{h_k} \quad\quad \forall \; \; r^{'} \in {]r_k, r_{k+1}[} \end{equation} Il est possible d'observer donc que le champs éléctrique est approximé par une formules des différences finies entre $r_k$ et $r_{k+1}$. \paragraph{Verification de la loi de Maxwell (\ref{maxwell}).} Maintenant que le champs de déplacement $D(r)$ peut être détérminé, il est intéressant de tester l'algorithme en vérifiant la loi de maxwell pour l'éléctrodynamique du système. Par l'équation (\ref{maxwell}) et par la forme cylindrique de l'omperateur de divergence ($\nabla \cdot$) il est possible d'obtenir la rélation suivante : \begin{equation}\label{maxdiscr} \frac{1}{r}\frac{\partial}{\partial r}(r D(r))=\rho_{lib}(r) \end{equation} Soit $\alpha(r) = r D(r)$, pour vérifier numériquement l'équation (\ref{maxdiscr}) il est nécessaire d'évaluer numériquement $\frac{\partial \alpha(r)}{\partial r}$. Pour faire ainsi, sont d'abord évaluée les valeur $\alpha(r_{k+1/2})$ pour $k=1,... n$, car la foction est discontinue sur tout $r_k$ par construction de $\frac{\partial \phi}{\partial r}$ (Rappel: $r_{k+1/2} = \frac{r_{k+1} +r_k}{2}$). Ainsi, on approxime la dérivée partielle de $\alpha$ évaluée dans le point $r_{mid,k} = \frac{r_{k+1/2} + r_{k+3/2}}{2}$, lequel corréspond au point intémédaire entre les points pour lesquels on a calculé $\alpha(r)$, par $\frac{\partial\alpha(r_{mid,k})}{\partial r} = \frac{\alpha(r_{k+3/2}) - \alpha(r_{k+1/2})}{r_{k+3/2} - r_{k+1/2}}$. Ainsi faisant, numériquement il ne reste qu'à vérifier la rélation suivante : \begin{equation}\label{alphamid} \frac{1}{r_{mid,k}} \frac{\partial\alpha(r_{mid,k})}{\partial r} = \rho_{lib}(r_{mid,k}) \end{equation} \section{Simulations Numériques} Pour discrétiser le problème on partage l'intérval $[0,R]$ en $N_1$ intervalles entre $r = 0$ et $r = b$, et en $N_2$ intervalles entre $r = b$ et $r = R$. La condition de bord est $\phi(R) = V_0$ avec $V_0$ donné. \subsection{Cas trivial} Pour tester l'algorithme est simulé le cas trivial où $\rho_{lib}(r)= \epsilon_0$ et $\epsilon_{r}(r) = 1$, avec $V_0 = 0$, $b = 0.06$ \si{\metre}, $R = 0.12$ \si{\metre} et $N_1 = N_2$. \paragraph{Solution analytique.} De l'équation (\ref{maxdiscr}) l'on déduit : \begin{equation}\label{solving} -\frac{1}{r}\frac{\partial}{\partial r}(r \epsilon_0\epsilon_r(r)\frac{\partial \phi(r)}{\partial r})=\rho_{lib}(r) \end{equation} En appliquant les hypothèse pour le cas trivial plus la condition de bord, il est possible de résoudre l'équation différentielle pour obtenir la solution analytique suivante : \begin{equation}\label{solution} \phi(r) = V_0 + \frac{1}{4}(R^2 - r^2) \end{equation} L'on observe de l'équation (\ref{solution}) que pour ce cas trivial $\phi(0) = \frac{R^2}{4}$. \begin{minipage}{\textwidth} \hspace{-0.05\textwidth} \begin{minipage}{0.45\textwidth} \paragraph{Étude de convergence} Ayant obtenu la valeur réelle du potentiel en $r=0$, il est facile d'étudier la convergence du schéma numérique en ce point, en évaluant l'erreur $|\phi_0 - \phi(0)|$ en fonction de la discretisation faite, c'est à dire en fonction du nombre de points $N$ par lequel l'interval $[0,R]$ est discrétisé. Le graphe en Figure \ref{conv_triv_fig} montre que la valeur simulée $\phi_0$ converge effectivement vers $\phi(0)$. L'algorithme donc converge à l'ordre 2, car étant le graphe en log log cet ordre correspond à la pente de la courbe. \end{minipage} \hspace{0.02\textwidth} \begin{minipage}{0.6\textwidth} \resizebox{\textwidth}{!}{ \input{graphs/conv_triv.tex} } \captionof{figure}{Test de convérgence en $\phi(r = 0)$} \label{conv_triv_fig} \end{minipage} \end{minipage} % TODO templates for layout %\begin{minipage}{\textwidth} %\hspace{-0.05\textwidth} %\begin{minipage}{0.45\textwidth} % % HERE %\end{minipage} %\hspace{0.02\textwidth} %\begin{minipage}{0.55\textwidth} %\resizebox{\textwidth}{!}{ % % HERE %} %\end{minipage} %\end{minipage} % % % %\begin{figure}[h] %\hspace{-0.1\textwidth} %\begin{subfigure}{0.55\textwidth} %\resizebox{\textwidth}{!}{ % \input{graphs/phi_nontriv.tex} %} %\end{subfigure} %\hspace{0.02\textwidth} %\begin{subfigure}{0.55\textwidth} %\resizebox{\textwidth}{!}{ % \input{graphs/E_nontriv.tex} %} %\end{subfigure} %\end{figure} \subsection{Cas non trivial} Maintenant l'on a $V_0 = 0$, $b = 0.02$ \si{\metre}, $R = 0.12$ \si{\metre}, et \begin{equation}\label{epsir} \epsilon_r(r) = \begin{cases} 1 & (0\leq r <b)\\ 8 - 6\frac{r - b}{R - b} & (b \leq r \leq R) \end{cases} \end{equation} \begin{equation}\label{rhoo} \rho_{lib}(r) = \begin{cases} \epsilon_0 a_0 (1 - (\frac{r}{b})^2) & (0\leq r <b)\\ 0 & (b \leq r \leq R) \end{cases} \end{equation} avec $a_0 = 10^4$ \si{\volt\per\metre^2} \paragraph{Solutions numérique.} Sur les graphes en Figure \ref{phi_subfig} et \ref{E_subfig} il est possible d'observer respectivement les soultions numériques pour le potentiel $\phi(r)$ et le champs éléctriques $E(r)$ dans le cylindre. Le graphe \ref{phi_subfig} illustre le potentiel en fonction de r, et il est possible d'y observer que ceci prends deux formes distinctes. Ceci est du naturellement à la discontinuité du milieu entre le vide et le diéléctrique et cette discontinuité peut s'observer plus clairement sur le graphe \ref{E_subfig} du champs éléctrique. Dans le premier interval $\phi(r)$ subit une descente approximativement parabolique. En passant par une estimation de l'équation (\ref{solving}) lorsque l'on y injecte l'équation (\ref{rhoo}), il est en effet possible de determiner que $\phi(r)$ descend à l'ordre 4 sur $r$. En effet la situation est similaire au cas trivial si l'on considère le domaine représentant le cylindre vide de rayon $b$. \begin{figure}[h] \hspace{-0.15\textwidth} \begin{subfigure}{0.65\textwidth} \resizebox{\textwidth}{!}{ \input{graphs/phi_nontriv.tex} } \caption{$\phi(r)$} \label{phi_subfig} \end{subfigure} \hspace{0.02\textwidth} \begin{subfigure}{0.65\textwidth} \resizebox{\textwidth}{!}{ \input{graphs/E_nontriv.tex} } \caption{$E(r)$} \label{E_subfig} \end{subfigure} \caption{Graphe de $\phi(r)$ et $E(r)$ pour différents $N_1$} \label{phiE_fig} \end{figure} Sur l'intérval $[b,R]$, $\phi(r)$ descend plus lentement et sa dérivée $E(r)$ tend donc à se minimiser et elle est discontinue en $r=b$ comme anticipé. Ce comportement peut s'expliquer par le fait que le diéléctrique s'oppose à l'action du champs éléctrique par l'éffet d'une polarisation qui en contient donc les effets. \begin{minipage}{\textwidth} \hspace{-0.05\textwidth} \begin{minipage}{0.45\textwidth} En confrontant les deux graphes il est possible de vérifier en fin le bon comportement et cohérence physique de $E(r)$ qui approxime bien le gradient de $\phi(r)$. Pour finir, pour une discretisation assez grande une valeur de $\phi(r=b)$ est cherchée pour ensuite étudier la convergence en $N_2 \propto N_1$ du potentiel en cette valeur pour $r=b$. Le graphe en Figure \ref{conv_nontriv_fig} illustre qu'en effet l'algorithme converge aussi et de nouveau à l'ordre 2 comme dans le cas trivial. \end{minipage} \hspace{0.02\textwidth} \begin{minipage}{0.6\textwidth} \resizebox{\textwidth}{!}{ \input{graphs/conv_nontriv.tex} } \captionof{figure}{Test de convérgence de $\phi(r = b)$} \label{conv_nontriv_fig} \end{minipage} \end{minipage} \subsubsection{Verifier (\ref{maxwell}).} En utilisant les différences finies pour l'operateur $d/dr$ il est possible de calculer la divergence du champ de déplacement $\nabla \cdot D$ et vérifier ainsi numériquement l'équation (\ref{maxwell}). Le graphe en Figure \ref{divD_subfig} illustre en effet la solution numérique de $\nabla \cdot D$ et le graphe \ref{divDrho_subfig} montre l'ordre de l'erreur $|\nabla \cdot D - \rho_{lib}|(r)$ où $\rho_{lib}(r)$ est défini en (\ref{rhoo}). %TODO, inserire grafico Div D e rholib \begin{figure}[h] \hspace{-0.15\textwidth} \begin{subfigure}{0.65\textwidth} \resizebox{\textwidth}{!}{ \input{graphs/divD.tex} } \caption{$\nabla \cdot D$} \label{divD_subfig} \end{subfigure} \hspace{0.02\textwidth} \begin{subfigure}{0.65\textwidth} \resizebox{\textwidth}{!}{ \input{graphs/divD_diff_rho.tex} } \caption{Différence entre $\nabla \cdot D$ et $\rho_{lib}$} \label{divDrho_subfig} \end{subfigure} \caption{Verification de l'équation (\ref{maxwell})} \label{divDrho_fig} \end{figure} %TODO Correggere Sur le graphe \ref{divDrho_subfig} il est possible d'observer que la rélation (\ref{div}) est vérifiée car l'ordre de grandeur de l'erreur est assez petit pour être négligeable. \subsubsection{Charges de polarisation.} Lorsque un diélèctrique est traversé par un champs éléctrique $E(r)$, il se polarise pour contraster l'effet de ce dernier. Elles aparaissent donc des charges de polarisation sur l'intèrface vide-diélèctrique, dont la charge totale est $Q_{pol}$ et il est possible d'en évaluer numériquement la densitée $\rho_{pol}(r)$ par la rélation $\rho_{pol}(r) = \epsilon_0 \nabla \cdot E (r) - \nabla \cdot D (r)$. \paragraph{Valeurs de référence pour les charges de polarisation.} À partir de l'équation différentielle (\ref{maxdiscr}), il est possible de trouver une éxpréssion analytique da $\alpha$ en sachant que $\alpha(0) = 0$ et en intégrant sur le domaine $[0, r]$: \begin{equation}\label{alpha} \alpha(r) = \int \limits_{0}^{r} s \; \rho_{lib}(s) ds = \begin{cases} \frac{1}{2} a_0 \; \epsilon_0 \; r^2 \; (1 - \frac{r^2}{2b^2}) & \text{si } 0 \leq r < b \\ \frac{1}{4} a_0 \; \epsilon_0 \; b^2 & \text{si } b \leq r \leq R \end{cases} \end{equation} En considerant maintenant la grandeur $\rho_{pol}(r)$, en intègrant sur un domaine volumique $\Omega$ à forme de tube mince de rayon moyen $b$ et d'épaisseur $h_- + h_+$ (trés petit), on obtien en utilisant le \textit{théorème de la divérgence}: \begin{equation*} Q_{pol} := \iiint \limits_\Omega \nabla \cdot (\epsilon_0 \vec{E} - \vec{D}) dV = \epsilon_0 \iint \limits_{\partial \Omega} \vec{E} (1 - \epsilon_r) \; d\vec{\sigma} \end{equation*} Puisque $\vec{E}$ est projecté selon la diréction radiale par rapport au centre du cilindre, alors il est toujours pérpendiculaire à la surface $\partial \Omega$. La notation véctorielle n'est donc plus nécessaire, ce qui correspond en coordonnées cilindriques à: \begin{equation*} Q_{pol} = L_z \epsilon_0 \int \limits_{0}^{2\pi} \left. r \; E(r) (1 - \epsilon_r(r)) \right\rvert_{r = b - h_-}^{r = b + h_+} \; d\theta \end{equation*} où $L_z$ est l'hauteur du cilindre. -Puisque $\epsilon_r(r) = 1$ si $r < b$, la valeur évaluée à $b-h_-$ s'annulle et +Puisque $\epsilon_r(r) = 1$ si $r < b$, la valeur évaluée en $r = b-h_-$ s'annulle et si l'on considère que $E(r) = \frac{\alpha(r)}{\epsilon_0 \epsilon_r(r) r}$, on obtien l'éxpréssion suivante en faisant tendre $h_+$ vers $0$: \begin{equation}\label{qpol_ref} -Q_{pol} = 2 \pi \; L \; \alpha(b+h_+) \left(\frac{1}{\epsilon_r(b+h_+)} - 1\right) -= - \frac{7}{16} \pi \; L \; \epsilon_0 \; a_0 \; b^2 +Q_{pol} = 2 \pi \; L_z \; \alpha(b+h_+) \left(\frac{1}{\epsilon_r(b+h_+)} - 1\right) += - \frac{7}{16} \pi \; L_z\; \epsilon_0 \; a_0 \; b^2 \end{equation} -D'autre côté, au niveau de la discrétisation numérique, en applicant la règle du -trapèze pour l'intégrant $r \rho_{pol}$ évalué sur le même intérval qu'en avance +D'un autre côté, au niveau de la discrétisation numérique, en applicant la règle du +trapèze pour l'intégrant $r \rho_{pol}(r)$ évalué sur le même intérval qu'en avance et en supposant que $\rho_{pol}(r) \approx 0 \quad \forall r \neq b$, l'approximation numérique qui en dérive est la suivante: \begin{equation}\label{qpol_num} -Q_{pol} = \iiint \limits_{\Omega} \rho_{pol} \;dV = 2\pi \; L \int \limits_{b-h_-}^{b+h_+} r \; \rho_{pol} \;dr \approx -\pi \; L \; r_{N_b} \; \rho_{pol}(r_{N_b}) \cdot (r_{N_b+1} - r_{N_b-1}) +Q_{pol} = \iiint \limits_{\Omega} \rho_{pol} \;dV = 2\pi \; L_z \int \limits_{b-h_-}^{b+h_+} r \; \rho_{pol} \;dr \approx +\pi \; L_z \; r_{N_b} \; \rho_{pol}(r_{N_b}) \cdot (r_{N_b+1} - r_{N_b-1}) \end{equation} où $N_b$ est tel que $\rho_{pol}(r_{N_b}) \neq 0$ donne la singularité charactèristique des charges de polarisation. \paragraph{Solution numérique} Le graphe en Figure \ref{rhopol_pic_subfig} illustre $\rho_{pol}(r)$ pour différentes discretisations. Il est possible d'observer que la densitée est non nulle seleument en proximité de l'intérface, et -que le pic de $\rho_{pol}$ est négatif en impliquant que les charges de polarisation ont un signe négatif. +que le pic de $\rho_{pol}(r)$ est négatif en impliquant que les charges de polarisation ont un signe négatif. En effet $E(r)$ a un sens positif (voir graphe \ref{E_subfig}), donc par opposition à cela le diéléctrique cumule des charges négatives sur son intérface. En suite sur le graphe \ref{rhopol_pic_subfig} il est possible d'observer que le pic augmente lineairement en valeur absolu lorsque des $N$ plus grands sont pris pour la discretisation du système. Ceci est du au fait que physiquement les charges de polarisation se cumulent uniquement sur l'intérface du diéléctrique en $r=b$, mais par discretisation il est impossible de representer cela. Pourtant la charge totale de polarisation $Q_{pol}$, laquelle est donnée par l'intégrale sur $r$ de la courbe $(r, \; \rho_{pol}(r))$, se conserve numériquement. \begin{figure}[h] \hspace{-0.15\textwidth} \begin{subfigure}{0.65\textwidth} \resizebox{\textwidth}{!}{ \input{graphs/rhopol.tex} } \caption{$\rho_{pol}(r)$} \label{rhopol_subfig} \end{subfigure} \hspace{0.02\textwidth} \begin{subfigure}{0.65\textwidth} \resizebox{\textwidth}{!}{ \input{graphs/rhopol_augm.tex} } \caption{Pic de dénsité de charges en $r = b$ en fonction de $N_1$} \label{rhopol_pic_subfig} \end{subfigure} \caption{Observations sur les pics de $\rho_{pol}(r)$.} \label{rhopol_fig} \end{figure} En effet, pour des discretisation plus {\it fines} les charges de polarisation restent dans l'intervalle approchant $r=b$, et donc pour un interval plus mince le pic de densité en ce point augmente en valeur absolue en conservant $Q_{pol}$. \begin{minipage}{\textwidth} \hspace{-0.1\textwidth} \begin{minipage}{0.4\textwidth} En fin, en intégrant donc sous la courbe $\rho_{pol}(r)$ pour différentes discretisations de $N_1$ et $N_2$, il est possible de chercher une convergence sur la valeur totale de $Q_{pol}$. Puisque la charge dépend aussi de la hauteur $L_z$ du cylindre, cette valeur à été fixé à $1$ \si{\metre}. Ce qu'on observe par la figure \ref{q_pol_fig} c'est que les valeurs numériques convergent vers la valeur de référence $Q_{pol_{ref}}$ calculée précedemment par l'éxpréssion (\ref{qpol_ref}) et qui vaut: \begin{equation*} Q_{pol} = -4.8677 \cdot 10^{-11} \; \text{\si{\coulomb}} \end{equation*} \end{minipage} \hspace{0.02\textwidth} \begin{minipage}{0.7\textwidth} \resizebox{\textwidth}{!}{ \input{graphs/qpol_conv.tex} } \captionof{figure}{Test de convérgence sur la valeur de $Q_{pol}$.} \label{q_pol_fig} \end{minipage} \end{minipage} \section{Conclusions} Grace à la méthode des éléments finis, il a été possible de construire un algorithme pour la détérmination numérique d'une solution à un problème à condition de bords non trivial. La façon la plus immédiate de tester un tel algorithme est celle de tester sa cohérence physique par exemple en vérifiant une des lois de Maxwell de l'éléctrodynamique, ou même de le tester sur un problème dont la solution analytique est connue explicitement. De même une fois verifiée la convergence de l'algorithme, il est possible de l'utiliser pour obtenir numériquement la solution à un problème non trivial pour pouvoir en déduire différentes observations, et par exemple en tirer des infos intéresantes comme la charge de polarisation totale du dièlèctrique dans le système étudié. Celle ci est une procédure d'étude très récurrente est typique en physique, lorsque le système étudié ne possède pas de solutions pouvant être détérminé analytiquement. \section{Annexes} \subsection{Formule du trapèze et du point du milieu} Pour approximer l'intégrale unidimensionnelle d'une fonction $f(x)$ dans l'intérval $[x_k,x_{k+1}]$ de longuer $h_k = x_{k+1} - x_k$ il est possible d'employer la formule du trapèze mixte à celle du point du milieu, pour un facteur $p \in [0,1]$: \begin{equation}\label{trapezio} \int \limits_{x_k}^{x_{k+1}} f(x) dx \approx h_k [p\frac{f(x_k)+f(x_{k+1})}{2} + (1-p)f(\frac{x_k+x_{k+1}}{2})] \end{equation} Pour obtenir cette formule la fonction $f(x)$ est approximée à une {\it valeur moyenne} qui prend la fonction dans l'intérval d'intégration : \begin{equation}\label{trapezio2} f_{moy}(x) = p\frac{f(x_k)+f(x_{k+1})}{2} + (1-p)f(\frac{x_k+x_{k+1}}{2}) \end{equation} Ainsi faisant en intégrant sur cette valeur constante il en résulte le facteur $h_k$ apparaissant dans l'équation (\ref{trapezio}). \begin{thebibliography}{99} \bibitem{ref2} L. VILLARD, PHYSIQUE NUMERIQUE I - II ,Swiss Plasma Center, Faculté des Sciences de Base, Section de Physique, Ecole Polytecnique Fédérale de Lausanne. \bibitem{ref1} ECOLE POLYTECNIQUE FEDERALE DE LAUSANNE, Semestre de printemps 2019, Physique Numérique I - II - Exercice 5. \bibitem{ref3} L. VILLARD, Cours de Physique Numérique pour physiciens, Faculté des Sciences de Base, Section de Physique, Ecole Polytecnique Fédérale de Lausanne. \end{thebibliography} \end{document} %%%% THE END %%%%