Homec4science
Diffusion yearII_reports (master) Browse

rSTICAZZI/physnum/rap2/RapportExercice2_Ancarola_Cincotti/RapportExercice2_Ancarola_Cincotti.tex
fe2a2725821dmaster

Last Change
Raw File

RapportExercice2_Ancarola_Cincotti.tex
View Options

% debut d'un fichier latex standard
\documentclass[a4paper,9pt,twoside]{article}
\usepackage{comment}
% 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}
\usepackage[T1]{fontenc}
% 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}
% module pour les graphiques
\usepackage{tikz}
\usepackage{pgfplots}
\usepackage{pgfplotstable}
% module pour positionner les images
\usepackage{float}
\usepackage{subcaption}
%\usepackage[bottom = 0.5in]{geometry}
\usepgfplotslibrary{external}
\pgfplotsset{compat=1.3}
%\tikzexternalize
\paperheight=300mm
\paperwidth=210mm
\setlength{\textheight}{220mm}
\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}}}
%
% 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{Particule dans un champ élécromagnétique}
\date{\today}
\author{Raffaele Ancarola, Cincotti Armando\\{\small \mail{raffaele.ancarola@epfl.ch} \small \mail{cincotti.armando@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} %------------------------------------------
Le mouvement d'une particule chargée dans un champs éléctromagnétique est régit par l'action de la force de Lorentz sur cette particule.
Ce système physique est analytiquement simple à détérminer, et pourtant, il est facile d'obtenir des valeurs analytiques à travers lesquels est possible tester
la stabilité d'un schéma d'intégration numérique.
Ce rapport illustréra donc la détérmination analityque du système, et ensuite étudiera la simulation numérique du système à travers trois schéma d'intégration numérique différents : Euler, Euler-Cromer, Runge-Kutta d'ordre 2.
Elles seront étudiées ensuite, la convergence et la stailité des différents schémas numérique en les confrontant avec les valeurs analytiques.
Le but de l'étude est en effet celui de tester les schéma pour pouvoir les confronter, et savoir ensuite lequel est plus susceptible d'être utile pour l'étude d'un système physique plus complexe.
\section{Schémas numériques}
Ici sont illustrés les algorithmes pour les trois schémas d'intégration numérique cités en introduction.
\paragraph{Euler}
\begin{align} \label{algoeuler}
\vec{r}(t_{i+1}) &= \vec{r}(t_i) + \vec{V}(t_i) \cdot dt \\
\vec{V}(t_{i+1}) &= \vec{V}(t_i) + \vec{a}(\vec{V}(t_i), t_i) \cdot dt
\end{align}
\paragraph{Euler-Cromer}
\begin{align} \label{algoeulercromer}
\vec{V}(t_{i+1}) &= \vec{V}(t_i) + \vec{a}(\vec{V}(t_i), t_i) \cdot dt \\
\vec{r}(t_{i+1}) &= \vec{r}(t_i) + \vec{V}(t_{i+1}) \cdot dt \\
\end{align}
\paragraph{Runge-Kutta d'ordre 2}
\begin{align} \label{algorungekutta}
\vec{k_r} &= (\vec{V}(t_i) + \frac{1}{2} \vec{a}(\vec{V}(t_i), t_i) \cdot dt) \cdot dt \\
\vec{r}(t_{i+1}) &= \vec{r}(t_i) + \vec{k_r} \\
\vec{k_1} &= \vec{a}(\vec{V}(t_i), t_i) \\
\vec{k_2} &= \vec{a}(\vec{V}(t_i) + \frac{1}{2}\vec{k_1}, t_{i+1/2}) \cdot dt \\
\vec{V}(t_{i+1}) &= \vec{V}(t_i) + \vec{k_2}
\end{align}
\section{Détermination Analytique du système}
On étudie le mouvement d'un proton chargé dans le plan {\it xOy} plongé dans un champ éléctromagnétique uniforme.
\subsection{Donnés}
\begin{itemize}
\item Masse du proton $m = 1.6726 \cdot 10^{-27} kg$.
\item Charge du proton $q = 1.6022 \cdot 10^{-19} C$.
\item La vitesse initiale ${\bf v_0} = (v_{x0},v_{y0})$.
\item Le champ éléctrique uniforme ${\bf E} = {\it E}{\bf y}$.
\item Le champs magnétique uniforme ${\bf B} = {\it B}{\bf z}$.
\end{itemize}
\subsubsection{Rélations}
La seule force appliquée au proton dans le système, est la force de {\it Lorentz}:
\begin{equation}\label{lorentz}
\vec{F} = q(\vec{E} + \vec{v} \times \vec{B})
\end{equation}
l'action de {\it B} sur la particule change en fonction de la position {\it x} de celle ci, en suivant la rélation:
\begin{equation}\label{B}
B(x) = B_0\cdot(1 + \kappa x)
\end{equation}
Le facteur $\kappa $ apparaitra dans une des séction quant il faudra l'éxpliciter
Le moment magnétique d'une particule chargée soumise à l'action d'un champ magnétique est donné par la rélation:
\begin{equation}\label{mu}
\mu = \frac{mv^2}{2B}
\end{equation}
où {\it v} est le module de la vitesse prise par la particule.
\subsection{Calculs Analytique}
Par le deuxième principe de Newton, on obtien l'équation différentielle suivante en
considerant l'action de la force de {\it Lorentz} sur la particule:
\begin{align} \label{equadiff}
\vec{F_{res}} &= m \cdot \frac{d\vec{V}}{dt} \Rightarrow m \frac{d\vec{V}}{dt} = q (\vec{E}
+ \cdot \vec{V} \times \vec{B}) \\
\Leftrightarrow &
\begin{cases}
\frac{dV_x}{dt} &= q E_x + \omega V_y \\
\frac{dV_y}{dt} &= q E_y - \omega V_x \\
V_x &= \frac{dx}{dt} \\
V_y &= \frac{dy}{dt}
\end{cases} \quad \omega = \frac{qB}{m}
\end{align}
Pour résoudre analytiquement le système d'équation différentielles on suppose $E = 0$ et $B = B_0$ et on obtient, pour les deux prémières équations du système (\ref{equadiff}), l'équation vectorièlle suivante.
\begin{align}
\Leftrightarrow
\frac{d\vec{V}}{dt} &=
\begin{pmatrix}
0 & \omega \\
-\omega & 0
\end{pmatrix} \cdot \vec{V}
\end{align}
Pour résoudre l'équation on cherche les valeurs propres de la matrice:
\begin{equation} \label{eigenvalues}
\begin{vmatrix}
\lambda & \omega \\
-\omega & \lambda
\end{vmatrix} = \lambda^2 + \omega^2 = 0 \Rightarrow \lambda = \pm i\omega
\end{equation}
Puisque la matrice se laisse diagonaliser par des valeur complexes, par un théorème
analytique la solution de l'équation différentielle est de la forme:
\begin{equation} \label{solutionform}
\begin{cases}
V_x(t) &= A \cos(\omega t) + C \sin(\omega t) \\
V_y(t) &= B \cos(\omega t) + D \sin(\omega t)
\end{cases}
\end{equation}
avec $A$, $B$, $C$, $D$ à déterminer selon les condition initielles.
En posant, en fait $\vec{V}(t = 0)$ on obtien:
\begin{equation} \label{initspeed}
\vec{V}(0) =
\begin{pmatrix}
A \\
B
\end{pmatrix} :=
\begin{pmatrix}
V_{x_0} \\
V_{y_0}
\end{pmatrix}
\end{equation}
Pour trouver $C$ et $D$ il faut dériver l'expression de la solution (\ref{solutionform}) et puis
appliquer la condition $\frac{d\vec{V}}{dt}(t = 0)$:
\begin{equation} \label{initacc}
\frac{d\vec{V}}{dt}(0) = \omega \cdot
\begin{pmatrix}
C \\
D
\end{pmatrix} = \omega \cdot
\begin{pmatrix}
V_{y_0} \\
-V_{x_0}
\end{pmatrix}
\end{equation}
On obtient, enfin, l'expréssion de la vitesse de la particule:
\begin{equation} \label{speedexpr}
\vec{V}(t) =
\begin{pmatrix}
V_{x_0} \cos(\omega t) + V_{y_0} \sin(\omega t) \\
V_{y_0} \cos(\omega t) - V_{x_0} \sin(\omega t)
\end{pmatrix}
\end{equation}
Si on applique enfin la trasformation $V_0 = \sqrt{V_{x_0}^2 + V_{y_0}^2}$ et
$\phi = \arcsin(\frac{V_{x_0}}{V_0}) = \arccos(\frac{V_{y_0}}{V_0})$ on a
une forme plus agréable:
\begin{equation} \label{speedbestexpr}
\vec{V}(t) = V_0 \cdot
\begin{pmatrix}
\sin(\omega t + \phi) \\
\cos(\omega t + \phi)
\end{pmatrix}
\end{equation}
Une fois q'on a la fonction de la vitesse, il reste seulement à intégrer l'éxpression
pour obtenir la position en fonction du temps et, dans un deuxième moment,
déterminer la constante d'intégration $\vec{C_0}$ selon les conditions initiales.
\begin{align} \label{positionexpr}
\vec{r}(t) &= \int \vec{V(t)}dt + \vec{C_0} = \frac{V_0}{\omega}
\begin{pmatrix}
- \cos(\omega t + \phi) \\
\sin(\omega t + \phi)
\end{pmatrix} + \vec{C_0} \\
\vec{r}(t = 0) &= \frac{V_0}{\omega}
\begin{pmatrix}
-\cos(\phi) \\
\sin(\phi)
\end{pmatrix} + \vec{C}_0 :=
\begin{pmatrix}
x_0 \\
y_0
\end{pmatrix} \\
\Rightarrow \vec{C}_0 &=
\begin{pmatrix}
x_0 + \frac{V_0}{\omega} \cos(\phi) \\
y_0 - \frac{V_0}{\omega} \sin(\phi)
\end{pmatrix} =
\vec{r}(0) + \frac{1}{\omega}
\begin{pmatrix}
0 & 1 \\
-1 & 0
\end{pmatrix} \cdot \vec{V}(0)
\end{align}
Ce qui montre bien que la trajectoire est circulaire de rayon
$R = \frac{V_0}{\omega}$ et de centre $\vec{C}_0$.
\paragraph{Cas centré à l'origine (éxercice 2.1.3)}
C'est suffisant de poser $\vec{C_0} = 0$ et le resultat suivant est obtenu:
\begin{equation} \label{centered}
\begin{pmatrix}
x_0 \\
y_0
\end{pmatrix} = - \frac{1}{\omega}
\begin{pmatrix}
0 & 1 \\
-1 & 0
\end{pmatrix}
\begin{pmatrix}
0 \\
v_0
\end{pmatrix} = - \frac{m v_0}{qB} \hat{x} \approx \begin{pmatrix} -0.0014 \\ 0 \end{pmatrix}
\end{equation}
On conclut en observant que la tajectoire de la particule est décrite par le système suivant:
\begin{align}
\begin{cases}\label{traj}
x(t) &= x_0 cos(\omega t) \\
y(t) &= \frac{v_0}{\omega} sin(\omega t)
\end{cases}
\end{align}
\section{Simulations}
\subsection{Schéma de Euler, Euler Cromer et Runge-Kutta d'ordre 2}
En implémentant trois schémas numériques différents (\ref{algoeuler}) (\ref{algoeulercromer}) et (\ref{algorungekutta}),
et en les faisant travailler sur le système d'équation (\ref{equadiff}),
elle peut être lancée la simulation d'un proton dans un champ magnétique
pour les trois schémas d'intégration et en comparer les resultat.
Dans ce cas le valeurs numériques prises sont: $v_0 = 4 \cdot 10^5 \: m/s $ et $ B = B_0 = 3T $ et $E = 0$
Comme démandé dans la consigne, il est choisit un temps de simulation $t_{fin}$ corréspondant à
5 périodes de revolution du proton dans sa trajectoire circulaire, c'est à dire $ T = 5 \cdot \frac{2\pi}{\omega} \approx 1.0932 \cdot 10^{-7} s$ .
Les résultats des simulations pour chaque schéma figurent dans les graphes suivants.
\begin{figure}[H]
\centering
\caption{Trajectorie parcourue sur $n$ pas de temps}
\label{trajNsteps}
\begin{subfigure}{0.35\textwidth}
\includegraphics[width=\textwidth]{ScreenMatlab/22/Euler/Graphe1x.png}
\caption{Euler n-steps=$1000$}
\label{mileuler}
\end{subfigure}
\hspace{0.1\textwidth}
\begin{subfigure}{0.35\textwidth}
\includegraphics[width=\textwidth]{ScreenMatlab/22/Euler/Graphe3x.png}
\caption{Euler n-steps=$16000$}
\label{seizeuler}
\end{subfigure}
\begin{subfigure}{0.35\textwidth}
\includegraphics[width=\textwidth]{ScreenMatlab/22/EulerCromer/Graphe1x.png}
\caption{Euler-Cromer n-steps=$500$}
\label{cineulc}
\end{subfigure}
\hspace{0.1\textwidth}
\begin{subfigure}{0.35\textwidth}
\includegraphics[width=\textwidth]{ScreenMatlab/22/RungeKutta2/Graphe1x.png}
\caption{Runge-Kutta n-steps=$500$}
\label{cinrun}
\end{subfigure}
\end{figure}
En observant la les graphes (\ref{mileuler}) et (\ref{seizeuler}), il est évisent que pour Euler la trajectoire du proton reste une spirale malgré le pas de temps.
En effet, même si le {\it dt} est assez petit (par exemple $dt= t_{fin}/16000$), à chaque pas de temps Euler cumule une certaine fraction d'erreur, en faisant poursuivre la trajectoire sur la tangente de la trajectoire circulaire sans ajouster le pas du proton.
On voit, pour les graphes (\ref{cineulc}) et (\ref{cinrun}), que soit le schéma Euler-Cromer, que pour Runge-Kutta d'ordre 2, la trajectoire simulée est circulaire pour un dt rélativement grand.
En effet ces deux schémas ajustent le calcul en fonction du pas précedent, dans le cas d'Euler-Cromer, et en fonction du pas qui suit dans le cas de Runge-Kutta 2. Pour Euler-Cromer la trajectoire n'est pas parfaitement circulaire, mais on peut imaginer que, pour un pas de temps plus petit, la trajectoire sera ajustée encore mieux.
\subsubsection{Étude de convérgence} \label{Etconv}
Comme on choisit un temps de simulation $t_{fin}$ corréspondant
à 5 périodes de revolution, la particule devrait se trouver,
áprès ce temps, à la position initiale $(x_{fin} , y_{fin}) = (x_0 , y_0 )$ .
En prénant, pour des pas de temps différents la valeur finale
simulée pour les trois schémas, on obtient les deux graphes suivants,
qui illustrent la convergence des schéma en fonction du pas de temps
$dt$ pour les valeurs $x_{fin}$ et $y_{fin}$.
En fin, pour un étude de convergence plus détaillé, on prend l'erreur maximale de la
trajectoire simulée par rapport à celle analytique décrite par la rélation (\ref{traj})
en fonction du pas de temps $\triangle t$ et on le représente sur un graphe loglog.
% Primo grafico
%
\begin{figure}[H]
\centering
\begin{subfigure}{0.48\textwidth}
\input{graphs/convx.tex}
\caption{Étude sur $x_{fin}$}
\label{x}
\end{subfigure}
\begin{subfigure}{0.48\textwidth}
\input{graphs/convy.tex}
\caption{Étude sur $y_{fin}$}
\label{y}
\end{subfigure}
\caption{Étude de convérgence des algorithmes d'intégration}
\label{convall}
\end{figure}
\begin{figure}[H]
\includegraphics[width = 0.9\textwidth]{etudesup.eps}
\caption{Étude sur $e_{max}$}
\label{emax}
\end{figure}
Les graphes en figure (\ref{convall}) illustrent la convergence des trois schéma numériques en fonction de dt.
Une rémarque évidente est que, comme il peut être observé par les graphes (\ref{x}) et (\ref{y})
, soit pour la valeur de $x_{fin}$ que $y_{fin}$, le schéma d'Euler est celui qui converge le moins rapidement.
Ceci était prévisible en observant que le schéma d'Euler ne permet pas d'obtenir
la bonne trajectoire numériquement, comme vu dans la séction précedente.
Par les graphe (\ref{x}) et (\ref{y}) il peut être constaté donc que les trois schémas convergent vers les valeurs analytiques.
En observant la pente des courbes dans le graphe (\ref{emax}), l'ordre de convergence de chaque schéma peut être déterminé.
Ceci à cause des propriété du logarithme. En effet, une pente de valeur n correspond à un degré de convergence d'ordre n, car $log(x^n) = n\cdot log(x)$.
Ceci dit, il peut être constaté que Euler et Euler-Cromer convergent à l'ordre 1, et Runge-Kutta d'ordre 2, comme son nom l'indique, à l'ordre 2.
Il est intéressant de noter aussi, que pour des dt assez grands, le schéma d'Euler cumule son erreur de façons exponentielle. Ce caractère sera observé, plus
expiclitement, dans la section suivante.
Pour conclure, le schéma qui convege le plus rapidement, et qui présent un moindre $e_{max}$, est celui de Runge-Kutta d'ordre 2, suivit par Euler-Cromer et enfin par Euler.
\subsection{Étude de stabilité}
Comme les notions de convergence et de stabilitée sont différentes, il reste à tester cette dérnière pour les trois schémas.
Pour le faire, faut rémarquer que, dans le système étudié les forces appliquée à la particule (voir rélation (\ref{lorentz}))
sont toutes consérvative. Or l'énérgie est conservée.
Un simple étude de stabilité consiste donc, à observer comment varie l'erreur sur la conservation
de l'énérgie pour chaque schéma à un nombre de pas de temps
donné.
% Graph Euler 1
%
\begin{figure}[H]
\caption{Étude de stabilité de {\it Euler}}
\label{stabeuler}
\centering
\begin{subfigure}{0.48\textwidth}
\input{graphs/stabeuler.tex}
\caption{Energie pour Euler à n-steps=500}
\label{stabe}
\end{subfigure}
\begin{subfigure}{0.48\textwidth} % no spacing here
\input{graphs/stabeulerbetter.tex}
\caption{Energie pour Euler à n-steps=16000}
\label{stabem}
\end{subfigure}
\end{figure}
% Graph Euler Cromer
%
\begin{figure}[H]
\caption{Étude de stabilité de {\it Euler cromer} et {\it Runge Kutta}}
\label{stabother}
\centering
\begin{subfigure}{0.45\textwidth}
\input{graphs/stabeulercromer.tex}
\caption{Energie pour Euler-Cromer à n-steps=16000}
\label{stabec}
\end{subfigure}
\hspace{0.05\textwidth}
\begin{subfigure}{0.45\textwidth}
\includegraphics[width = \textwidth]{ScreenMatlab/22/RungeKutta2/Energie.PNG}
\caption{Energie pour Runge-Kutta à n-steps=16000}
\label{stabrk}
\end{subfigure}
\end{figure}
Les graphes (\ref{stabe}) et (\ref{stabem}) correspondent à l'enérgie calculée par l'algorithme d'Euler, pour un n-steps = 500 et 16000.
On observe que le schéma n'est pas stable, c'est à dire que l'erreur est croissant, même strictement: $e(t_{i+1}) > e(t_i)$ où $t_{i+1} > t_i$.
Il peut être par contre observé que l'erreur cumulée, dépend du choix du pas de temps {\it dt}:
en effet pour un dt rélativement grand l'erreur se cumule rapidement, de façons exponentielle (voir graphe (\ref{stabe}),
et pour un dt rélativement petit la croissance de l'erreur semble être linéaire (voir graphe (\ref{stabem})).
Dans le cas d'Euler Cromer (voir graphe (\ref{stabec})), l'erreur maximale sur l'énérgie est constant, donc le schéma est stable.
Dans le cas de Runge-Kutta (voir graphe (\ref{stabrk})), l'erreur augmente linéairement, donc le schéma n'est pas stable sous les conditions du système
, mais pourtant l'erreur augmente très lentement.
En effet sur le graphe (\ref{stabrk}) on observe que l'ordre de grandeur de l'erreur est $10^{-16}$ ce qui est très petit, donc l'erreur reste négligeable.
\section{Application 1, $\vec{E} \times \vec{B}$}
Dans cette séction, au traver le schéma de Runge-Kutta d'ordre 2, sera étudié le comportement d'un proton dans un champs éléctromagnétique.
Les valeurs numériques prises sont : $E = 6 \cdot 10^4 V/m , B = B_0 = 3T , v_{x0} = 0, v_{y0} = 4 \cdot 10^5$ et $t_{fin} \approx 1.0932 \cdot 10^{-7} s$ donc, comme décrit avant, la particule conplétera 5 révolution.
\begin{figure}[H]
\caption{Trajectorie de la particule chargée dans un champ éléctro-magnétique}
\label{trajEplusB}
\centering
\begin{subfigure}{0.4\textwidth}
\includegraphics[width=\textwidth]{ScreenMatlab/23/GrapheRinert.PNG}
\caption{Position}
\end{subfigure}
\hspace{0.05\textwidth}
\begin{subfigure}{0.4\textwidth}
\includegraphics[width=\textwidth]{ScreenMatlab/23/GrapheVVinert.PNG}
\caption{Vitesse}
\end{subfigure}
\begin{subfigure}{0.4\textwidth}
\includegraphics[width=\textwidth]{ScreenMatlab/23/Graphexinert.PNG}
\caption{Trajectorie sur un plan cartésien}
\end{subfigure}
\hspace{0.05\textwidth} % nessuno spazio qui, altrimenti va a capo
\begin{minipage}{0.5\textwidth}
Comme il peut être observé dans les graphes en figure (\ref{trajEplusB}), notamment dans le graphe c),
la particule décrit un cercle en translation horizontale.
Ceci permet de déduir que l'action de la force du champ éléctrique, combinée à celle du champ magnétique,
résulte en une translation de vitesse constante en direction $\hat{x}$.
La vitesse de translation horizontale du système est donc décrite par la rélation $\vec{v_e} = \frac{\vec{E} \times \vec{B}}{B^2}$.
En effet, si l'on reprèsente la position de la particule dans le référentiel non inértiel à vitesse constante $\vec{v_e}$,
il peut être observé que la trajectoire est en effet circulaire dans ce nouveau référentiel. Ceci est illustré par les graphes qui suivent.
\end{minipage}
\end{figure}
\begin{figure}[H]
\centering
\begin{subfigure}{0.4\textwidth}
\includegraphics[width=\textwidth]{ScreenMatlab/23/GrapheRNinert.PNG}
\caption{En fonction du temps}
\end{subfigure}
\hspace{0.08\textwidth}
\begin{subfigure}{0.4\textwidth}
\includegraphics[width=\textwidth]{ScreenMatlab/23/GraphexNinert.PNG}
\caption{Sur un plan cartésien}
\end{subfigure}
\caption{Trajéctoire dans un référentiel non-inértiel}
\label{trajEplusBrel}
\end{figure}
\subsection{Étude de convérgence}
\begin{figure}[H]
\input{graphs/convrk.tex}
\caption{Étude de convergence de l'algorithme {\it Runge-Kutta}}
\label{convrkEplusB}
\end{figure}
Dans la figure (\ref{convrkEplusB}) on observe que le schéma utilisé, Runge-Kutta d'ordre 2, converge.
L'étude de convergence ici est fait de la même façons que dans la section précedente (\ref{Etconv}), en étudiant donc la convergence
sur les valeurs de $x_{fin}$ et $y_{fin}$. La convergence ne dépend donc pas du système étudié, seulement du schéma.
Ceci est en effet vrai pour tout schéma d'intégration numérique.
\subsection{Étude de stabilité}\label{stabr}
Aussi dans ce cas, de façon générale, les forces aplliquées au protons sont consérvatives, donc l'énergie reste constante.
Un étude de stabilité du schéma peut être donc fait sur l'erreur de l'énergie mécanique calculée par le schéma.
\begin{figure}[H]
\input{graphs/stabrkapp.tex}
\caption{Étude de stabilité de {\it Runge-Kutta}}
\label{stabrkEplusB}
\end{figure}
Dans le graphe en figure (\ref{stabrkEplusB}), il peut être observé que l'erreur maximale sur le calcul de lénérgie mécanique n'augmente pas.
On en déduit que le schéma est stable pour ces condition de système.
De plus, on rémarque que l'ordre de grandeur de l'erreur est très petit: $10^{-16}$, c'est à dire le même observé dans la figure (\ref{stabother}) b)).
L'erreur est donc négligéable.
\section{Application 2, $\nabla B$}
Maintenant, il sera considéré un champ magnétique d'intesité variable $\vec{B} = B(x)\hat{z}$, où $B(x)$ est détérminé par la rélation (\ref{B}).
Les valeurs numériques prises sont: $\kappa = 100 m^{-1} , B_0 = 3T , E = 0 , v_{x0} , v_{y0 = 4} \cdot 10^{5} m/s$ et le même $t_{fin}$ pris dans les sections précedentes.
Les simulations sont faites avec l'algorithme de Runge-Kutta d'ordre 2 illustré par la rélation (\ref{algorungekutta}),
et seront simulés le comportement d'un proton ($q = e$) et antiproton ($q = -e$) sujets au système décrit ici.
{\bf N.B.} en simulation, l'antiproton est positionnné à une position initiale $x_{0antiprot} = 0.0014$, l'inverse additif de $x_{0prot}$.
On peut observer les résultats suivants:
\begin{figure}[H]
\centering
\begin{subfigure}{0.4\textwidth}
\includegraphics[width=\textwidth]{ScreenMatlab/24/Graphex2prot.PNG}
\caption{Trajectoire proton}
\label{trajprot}
\end{subfigure}
\hspace{0.05\textwidth}
\begin{subfigure}{0.4\textwidth}
\includegraphics[width=\textwidth]{ScreenMatlab/24/Graphex2antiprot.PNG}
\caption{Trajectoire antiproton}
\label{trajaprot}
\end{subfigure}
\begin{subfigure}{0.4\textwidth}
\includegraphics[width=\textwidth]{ScreenMatlab/24/GrapheR2prot.PNG}
\caption{$x(t)$ et $y(t)$ proton}
\label{xyprot}
\end{subfigure}
\hspace{0.05\textwidth}
\begin{subfigure}{0.4\textwidth}
\includegraphics[width=\textwidth]{ScreenMatlab/24/GrapheR2antiprot.PNG}
\caption{$x(t)$ et $y(t)$ antiproton}
\label{xyaprot}
\end{subfigure}
\begin{subfigure}{0.4\textwidth}
\includegraphics[width=\textwidth]{ScreenMatlab/24/GrapheV2prot.PNG}
\caption{$v_x(t)$ et $v_y(t)$ proton}
\label{vprot}
\end{subfigure}
\hspace{0.05\textwidth}
\begin{subfigure}{0.4\textwidth}
\includegraphics[width=\textwidth]{ScreenMatlab/24/GrapheV2antiprot.PNG}
\caption{$v_x(t)$ et $v_y(t)$ antiproton}
\label{vaprot}
\end{subfigure}
\caption{Trajectoires proton et antiproton sous un champ magnétique à intensité variable}
\label{trajBvariable}
\end{figure}
Par l'étude des graphes dans la figure (\ref{trajBvariable}) il peut être observé, que l'action d'un champ magnétique d'intensité variable se traduit dans une translation verticale du cercle (voir graphes (\ref{trajprot}) et (\ref{trajaprot})). En étudiant les graphes (\ref{xyprot}) et (\ref{xyaprot}) on peut observer que, par le fait que le proton et antiprotons ont des charges de signe opposé, l'action du champ magnétique variable est aussi opposée pour les deux particules. En effet, la trajectoire du prontons subit une translation verticale vers le haut (voir graphe (\ref{xyprot})),
et celle de l'antiproton vers le bas (voir graphe (\ref{xyaprot})).
Pour finir, en observant les graphes (\ref{vprot}) et (\ref{vaprot}), on rémarque que les vitesses prises par les deux particules sont essentiellement les même, car la rotation des particule suis le sens et l'intensité du champ magnétique. La seule différence est donnée par une légère translation verticale des courbes $v_y(t)$ vers le haut pour le proton (voir graphe (\ref{vprot})) et vers le bas pour l'antiproton (voir graphe (\ref{vaprot})). Ceci implique la direction et le sens de la translation subits par la trajectoire.
Pour la courbe $v_x(t)$ on observe un déphasage de $\pi \omega$, mais ceci est du au fait que les deux particules ont une position initiale opposée par simmétrie avec l'axe $x = 0$. $v_{xprot}(t) = - v_{xantiprot}$ car les particules sont contraintes à suivre une trajectoire centrée en $x = 0$ par les conditions imposées.
\subsection{Moment magnétique}\label{magn}
Le moment magnétique d'une particule chargée est donné par la rélation (\ref{mu}).
Comme $B$ varie en fonction de $x$, et que $x$ est une fonction sinusoidale, alors aussi $\mu$ est une fonction sinusoidale.
Comme la trajectoire des particules est une spirale, la valeur du module de la vitesse peut être approximée à une constante $v(t) \approx v_0$.
La valeur moyenne de $\mu$ est prise quand $B(x) = B(0) = B_0$, donc $\mu_{moy} = \frac{mv_0^2}{2B_0} \approx = 4.4603 \cdot 10^{-17} J/T$.
La valeur analytique du moment magnétique est donnée par $\mu(t) = \frac{mv_0^2}{2B(x(t))}$, où $B(x(t)) = B_0(1+\kappa x_0cos(\omega t))$.
\subsubsection{Fonction analytique}
\begin{figure}[H]
\begin{subfigure}{0.5\textwidth}
\input{graphs/pantipanalytique.tex}
\caption{Moment magnétique $\mu$ des deux particules}
\label{pantipanalytique}
\end{subfigure}
\hspace{0.05\textwidth}
\begin{minipage}{0.4\textwidth}
En figure (\ref{pantipanalytique}) on observe une portion d'un sinus, pour le proton et l'antiproton, qui a en effet comme valeur moyenne $\mu_{moy}$
mais une oscillation très ample par rapport au temps $t$, ceci est donné par le fait que l'ordre de grandeur de $t , 10^{-7}$
est très petit par rapport à celui de $\omega , 10^8$ . Il faut par contre tenir compte du fait que la solution analytique est approximée pour
$\omega = qB_0/m$ alors qu'une meilleure solution analytique devrait tenir compte de $\omega(B(x))$.
\end{minipage}
\caption{Moment magnétique d'un proton et un antiproton}
\label{muanalytique}
\end{figure}
\subsubsection{Fonction donnée par la simulation}
\begin{figure}[H]
\begin{subfigure}{0.45\textwidth}
\input{graphs/proton.tex}
\caption{Proton}
\label{muproton}
\end{subfigure}
\begin{subfigure}{0.45\textwidth}
\input{graphs/antiproton.tex}
\caption{Antiproton}
\label{muantiproton}
\end{subfigure}
\caption{Moment magnétique simulé par l'algorithme Rung-Kutta}
\label{moment}
\end{figure}
Les courbes décrites dans les graphes (\ref{muproton}) et (\ref{muantiproton}), représentent des courbe sinusoidales de même amplitude et valeur moyenne
que celles du graphe précedent, mais avec une oscillation beaucoup moins ample. Ceci peut être du à l'erreur d'instabilité de l'algorithme, ou à l'approximation
de $\omega$ faite avant. Pourtant la constatation liée à la valeur moyenne $\mu_{moy}$ et le fait que l'ordre de grandeurs des oscillations est très petit,
on constate que l'on peut quand même faire confiance au schéma en faisant les bonne approximation. Le schéma va très bien si nous intéresse juste la valeur moyenne du moment magnétique des particules étudiées.
\section{Conclusions}
Le but de ce rapport étant d'etudier le comportement de différents schémas d'intégration numérique.
Pour cela un système physique rélativement simple, dont les solutions analytiques pouvaient être détérminée facilement, à été choisi pour
confroter les calcul algorithmique aux solutions analytiques.
Par les observations faites dans la totalité du rapport, on peut constater que, les trois schéma choisis, c'est à dire Euler Euler-Cromer et Runge-Kutta, convergent.
On sait d'entre autre que Runge-Kutta converge plus rapidement des deux autres, et qui est plus précis. De plus on sait que le schéma d'Euler n'est pas assez précis même pour des rélativement petits pas de temps.
L'ordre de convergence pour les trois schémas est réspéctivement, 1 pour Euler et Euler-Cromer, 2 pour Runge-Kutta.
Grace à l'étude de ce système, on peut étudier aussi la stabilité des schéma, et on a pu constater que:
Euler n'est jamais stable, Euler Cromer est stable dans les même conditions, et Runge-Kutta est stable seulement sous certaine conditions, comme on a pu voir dans la section \ref{stabr}.
Pourtant, en ce concentrant sur ce dérnier schéma, il a été observé, que même quand il n'est pas stable, l'erreur est la plupart du temps négligéable, en effet on obtient des résultat presque analytiquement parfait, comme par exemple dans la section \ref{magn}.
Pour conclur, par l'étude de ce système physique, on a pu tester différents schémas d'intégration numérique. Les résultats de ces test seront utile dans la suite
lorsque il faudra choisir le bon schéma pour l'étude d'un certain système physique.
\end{document} %%%% THE END %%%%
c4science · Help