assignment 1 2016, Ordinära differentialekvationer, MMG511 / MVE162
Author
Elias Kamyab
Last Updated
6 years ago
License
Creative Commons CC BY 4.0
Abstract
Chalmers, Göteborg universitet
Chalmers, Göteborg universitet
\documentclass{article}
\usepackage[T1]{fontenc} % Svenska bokstäver
\usepackage[utf8]{inputenc} % Teckenkodning i filen
\usepackage[swedish]{babel} % Svenska rubriker i innehållsförteckning etc
\usepackage{amsfonts} % Innehåller t.ex. fonterna \mathbb
\usepackage{listings}
\usepackage{color} %red, green, blue, yellow, cyan, magenta, black, white
\definecolor{mygreen}{RGB}{28,172,0} % color values Red, Green, Blue
\definecolor{mylilas}{RGB}{170,55,241}
\usepackage{fancyvrb}
\usepackage{amsmath} % Innehåller mycket användara kommandon för matte
\usepackage{amsthm} % Ger omgivningar för att skapa satser och bevis
\usepackage{graphicx} % för att infoga bilder
% Ofta vill man definiera egna kommandon för sådana saker som man skriver ofta
\usepackage{sidecap}
\usepackage{listings}
\newcommand{\tb}{\textbackslash} % Gör att den late kan skriva bara \tb
\newcommand{\R}{\mathbb{R}} % Ger reella talen
\newcommand{\Z}{\mathbb{Z}} % Ger heltalssymbolen
% Skapa olika typer av satser som man vill använda i sin fil (kräver amsthm)
% Det första är namnet på omgivningen som använder när man skapar den
% och det andra är det som skrivs ut.
\newtheorem{sats}{Sats}
\newtheorem{prop}{Proposition}
\newtheorem{lemma}{Lemma}
% Kommentarer inleds med procenttecken och visas inte
% Här kommer en titel och författare att visa i början av dokumnetet
% och i sidhuvudena (om man vill)
\title{assignment 1, Ordinära differentialekvationer, MMG511 / MVE162}
\author{Elias Kamyab Orvar}
\date{27 april 2016} % \today ger dagens datum
\begin{document}
%Skapa titelinfo
\maketitle
%Lägga till rubriker för sektioner etc.
\section{Inledning till rapporten}
Projekten utgår ifrån att man ska använda två klassiska metoder: en fri pendulum och en pendulum med periodisk oscillating pivot, den sistnämnda heter Kapitza pendulum, döpt efter den rysska fysikern $ Pyotre Kapitza $. Båda modellerna skall testa med friktion och utan friktion, Syftet med projekten är att i alla fyra fallen undersöka hur $ {\theta}$ och $\dot{\theta} $ kan förändras med tiden med hjälp av ordinära differential ekvationer med initailvärdena $ {\theta} \approx 0 \ eller\ \approx {\pi} $ och en starthastigheten $\dot{\theta} \approx 0 $. \newline Målet på projekten är att hitta stabiliteten och instabiliteten efter en viss tid på respektive jämvikstpunkter dvs. stabiliseringen vid start på ${\theta} \approx {\pi}$ och destabiliseringen vid start på ${\theta} \approx 0 $. Det är tänkt att varje fall analyseras numeriskt, analytisk och grafiskt. Som hjälpmedel har vi Matlab ODE45.
\section{Analys av problemen}
\subsection{beteende av fria pendulum utan friction}
En fri pendulum utan friktion har ekvationen
\begin{equation}
l\ddot{\theta} = g sin{\theta} \end{equation}
där g är tyngdaccelerationen i ${m}/s{2}$ och l är längden.\newline
vi ska först icekdimensionalisera ekvationen då gör vi först och främst en variabelbyte ${\tau} = {\sqrt[]{g/l}}.t$ vilket leder till att $\dot{\theta} = \frac{d{\theta}}{dt} = \frac{d{\theta}d{\tau}}{d{\tau}dt} = \frac{d{\theta}}{d{\tau}}\ {\sqrt[]{g/l}}\ t.$\newline
$\ddot{\theta} = \frac{d^{2}{\theta}d^{2}{\tau}}{d{\tau}^{2}dt^{2}} = \frac{d^{2}{\theta}}{dt^{2}}\ {g/l} $\newline
$l\ddot{\theta}$ = g sin${\theta} \iff \ddot{\theta}$ = $\frac{g}{l}\ sin{\theta} \iff $ (efter variabelbyte) $\frac{d^{2}{\theta}}{dt^{2}} = -sin{\theta}. $ Nu har vi inga enheter (l och g i det fallet) kvar i ekvationen; alltså är ekvationen dimensionslös både på HL och VL. \newline
Vi kan sätta $\dot{x1} = x2$ och $\dot{x2} = -sin(x1). $ Då får vi en system som \newline $ \begin{bmatrix} x1' \\ x2' \end{bmatrix} $ = $ \begin{bmatrix} 0 & 1 \\ -sin({\theta}) & 0 \end{bmatrix} \begin{bmatrix} x1 \\ x2 \end{bmatrix}$
\bigbreak
Med hjälp av Matlab ode45 plottar vi phase portrait av den fria pendulum utan friktion med start positionerna (vinkeln x-led) som $[0, \frac{\pi}{2}, \frac{3\pi}{4}$ och ${\pi}]$ och med starthastigheterna $[0, \frac{\pi}{8}].$ Som grafiskt resultat får vi figur 1.
Bilden ser mer eller mindre ut som en cirkel och säger oss om vi startat i noll eller ${\pi}$ med hastigheten noll stannar vi på respektive punkt (alltså stabila stationära punkter) men däremot om vi startar i vinklar mellan $(0, {\pi})$ får vi cirkelliknande öga där ju närmare noll vi startar desto finare cirklar vi får.
Ändrar vi starthastigheten till $\frac{\pi}{8}$ får vi andra typer av banor tx. vid startpositioner på $ {\pi} $ flyger banan iväg.\newline
Vi tar slutsatsen att för att denna punkt ${\theta} = 0$ är stabil (med eller utan små hastighet) alltså pendulum med start på 0 (utan friktion) är stabil däremot vid startvilnek på ${\pi}$ med små hastighet är en pendulum instabil (OBS: om hastigheten $\dot{\theta} = 0$ då är även ${\pi}$ en stabil punkt).
\begin{figure}[!ht]
\centering
\includegraphics[width=1\textwidth]{FreePendWithoutFriction.jpg}
\caption{\label{fig:frog}Olika kvalitativa types of banor (icke-linjäriserad).}
\end{figure}
En godtycklig ekvation som ger denna reslutat är $\overline{\dot{x}} = A\overline{x} $ Där A är en kvadratisk matris med trace(A) = 0 och Det(A) >0. ex. $ \begin{bmatrix} 0 & {\alpha} \\ {\beta} & 0 \end{bmatrix} $ där
$ {\alpha} , {\beta} \neq 0$
och har olika tecken. En godtyckligt ekvations system ges av
\newline
$ \begin{cases}
\dot{x}= {\alpha}y \\
\dot{y}= -{\beta}x \\
\end{cases} \ $ (svar till fråga 1)
\bigbreak
\bigbreak
vi kan linjäriserar den dimesionslösa ekvationen $\ddot{\theta} = -sin{\theta} $ genom att sätta $sin{\theta} = {\theta}$ när vinkeln ${\theta} = 0$ och $sin{\theta} = -{\theta}$ när vinkeln $ {\theta} = {\pi} $ då får vi $\ddot{\theta} = \pm {\theta} $
syftet med linjärisering är att kunna spproximera stabiliteten på ett bättre och enklare sätt.
matrisen blir A = $ \begin{bmatrix} 0 & x(2) \\ -x(1) & 0 \end{bmatrix} $ och systemets utseende $\newline \begin{bmatrix} x1' \\ x2' \end{bmatrix} $ = $ \begin{bmatrix} 0 & 1 \\ -1 & 0 \end{bmatrix} \begin{bmatrix} x1 \\ x2 \end{bmatrix}. Det(A-{\lambda}I) = 0$
$ \implies \lambda_1,_2 = \pm i $ där $Re({\lambda})=0 $ \newline
$ \begin{bmatrix} 0 & 1 \\ -1 & 0 \end{bmatrix} \begin{bmatrix} 1 \\ 0 \end{bmatrix} = \begin{bmatrix} 0 \\ -1 \end{bmatrix} $ och $ \begin{bmatrix} 0 & 1 \\ -1 & 0 \end{bmatrix} \begin{bmatrix} 0 \\ 1 \end{bmatrix} = \begin{bmatrix} 1 \\ 0 \end{bmatrix} $ vilket betyder med pilarna neråt på punkten $ \begin{bmatrix} 1 \\ 0 \end{bmatrix} $ och pilarna är höger på punkten $ \begin{bmatrix} 1 \\ 0 \end{bmatrix} $
\newline Resultat från analytiskt beräkning tyder på en center.
Här stödjer även numeriska beräkningen av linjäriserade ekvationen min analytiska resultat. Bilden man får från Matlab visar en center oberoende på vilka startpositionen och hastigheter man tar. På figure 2 användas;
$ \theta $ (startpositioenerna)
$ =[0, \frac{pi}{4}] $ och hastigheter
$ \dot{\theta} = {[0, \frac{pi}{8}]} $\newline
\begin{figure}[!ht]
\centering
\includegraphics[width=0.7\textwidth, height=0.3\textwidth]{LinjPendWithoutFrik00.jpg}
\caption{\label{fig:frog}phase portrait på linjäriserad pendulum utan friktion.}
\end{figure} Tille skillnad från icke-linjäriserade fallet får vi här bara en typ av banor dvs alla banor symboliserar cirklar (center).\newline
Vid start på ${\theta} = {\pi}$ får vi $ \begin{bmatrix} x1' \\ x2' \end{bmatrix} $ = $ \begin{bmatrix} 0 & 1 \\ 1 & 0 \end{bmatrix} \begin{bmatrix} x1 \\ x2 \end{bmatrix}. Det(A-{\lambda}I) = 0$
$ \implies \lambda_1,_2 = \pm 1 $ med egenvektorerna $ \begin{bmatrix} 1 \\ 1 \end{bmatrix} $ och $ \begin{bmatrix} -1 \\ 1 \end{bmatrix} $
alltså en sadelpunkt, se figur 3!
\begin{figure}[!ht]
\centering
\includegraphics[width=0.7\textwidth, height=0.3\textwidth]{sadell.jpg}
\caption{Fas portrait, linjäriserad pendulum utan friktion med start på ${\theta}={\pi}.$}
\end{figure}
\subsection{Beteende av fria pendulum med friktion}
En fri pendulum med friktion har ekvationen \begin{equation}
l\ddot{\theta} = -l{\gamma}\dot{\theta} - g sin{\theta}. \end{equation}
där ${\gamma}$ är friktionen. \newline
Vi kör samma metod (variabelbyte) här som vi gjorde i fallet uan friktion för att icke-diemsionalisera och linjärisera ekvationen. \newline
Som resultat får vi $\ddot{\theta} = -\frac{\gamma}{\sqrt[]{g/l}} - sin{\theta} $. Här gör vi ytterligare en variabelbyte och sätter $ {\gamma^{*}} = \frac{\gamma}{\sqrt[]{g/l}}$ då våran friktion-force beror enbart på en variabel dvs. ${\gamma^{*}}$ och den dimensionslösa ekvationens utseende blir $\ddot{\theta} = -{\gamma^{*}{\theta}} - sin{\theta}. $ \newline
Systemet ser ut \newline $ \begin{bmatrix} x1' \\ x2' \end{bmatrix} $ = $ \begin{bmatrix} 0 & 1 \\ -sin({\theta}) & {\gamma^{*}} \end{bmatrix} \begin{bmatrix} x1 \\ x2 \end{bmatrix}$ \newline
vi kan linjäriserar den dimesionslösa ekvationen $\ddot{\theta} = -{\gamma^{*}} -sin{\theta} $ genom att sätta $sin{\theta} = {\theta}$ när vinkeln ${\theta} = 0$ och $sin{\theta} = -{\theta}$ när vinkeln $ {\theta} = {\pi} $ då får vi \newline $\ddot{\theta} = -{\gamma^{*}\dot{\theta}} \pm {\theta} $ \newline
Tillbaka till systemet efter lijarisering $ \begin{bmatrix} x1' \\ x2' \end{bmatrix} $ = $ \begin{bmatrix} 0 & 1 \\ -1 & -{\gamma^{*}} \end{bmatrix} \begin{bmatrix} x1 \\ x2 \end{bmatrix}$
$Det(A - {\lambda}I) = 0$ och om vi tar ${\gamma}^{*} = 0.1$ (litet)
$ \implies \lambda_1,_2 \approx -0.05 \pm i $ där $Re({\lambda})< 0 $ alltså en stabil spiral. \newline
$ \begin{bmatrix} 0 & 1 \\ -1 & -0.1 \end{bmatrix} \begin{bmatrix} 1 \\ 0 \end{bmatrix} = \begin{bmatrix} 0 \\ -1 \end{bmatrix} $ och $ \begin{bmatrix} 0 & 1 \\ -1 & -0.1 \end{bmatrix} \begin{bmatrix} 0 \\ 1 \end{bmatrix} = \begin{bmatrix} 1 \\ -0.1 \end{bmatrix} $ vilket betyder med pilarna neråt på punkten $ \begin{bmatrix} 1 \\ 0 \end{bmatrix} $ och pilarna är höger på punkten $ \begin{bmatrix} 1 \\ 0 \end{bmatrix} $
bilderna från matlab styrker min analys, se figur 4!
\begin{figure}[!ht]
\centering
\includegraphics[width=1\textwidth,height=0.4 \textwidth]{FrePenWithFrinonLin.jpg}
\caption{\label{fig:frog}Olika kvalitativa types of banor (icke-linjäriserad).}
\end{figure}
\begin{figure}[!ht]
\centering
\includegraphics[width=0.8\textwidth,height=0.4 \textwidth]{feiPenWithFriLin.jpg}
\caption{\label{fig:frog} banorna (linjäriserad).}
\end{figure}
Startar vi i $0 \ eller \ {\pi}$ (på jämviktspunkterna) i det icke-linjäriserade fallet med starthastigehten exakt 0 är punkterna stabila (banorna sitter stille) och det händer ingenting, på startvinklarna mellan $0\ och\ {\pi}$ får vi en stabilt spiral som mer eller mindre är cirkellikt. Däremot knuffar vi med lite hastighet $\dot{\theta} = {\pi}/8 $ får vi inte länge stabilitet på $ {\theta} = {\pi}\ $(flyger iväg) men fortfarande stabilt på ${\theta} = 0$.\newline
På den linjäriserade fallet får vi en stabil spiral som har samma form oberoende av startvinklar och starthastigheter, se figur 5! Därför kan vi alltid approximera bättre beteendet om stabilitet eller instabilitet med linjäriserade ekvationer.
\subsection{pendulum med vertically oscillating pivot}
\begin{equation}
l\ddot{\theta} = -l{\gamma}\dot{\theta}-[g + \ddot{\epsilon}(t)] sin{\theta} \iff
\end{equation} ($ {\epsilon}(t)= Acos({\omega}t) $ är en periodisk funktion med A och ${\omega}$ som konstanter)
$ \ddot{\theta} = -{\gamma}\dot{\theta}-[\frac{g}{l} - \frac{A\omega}{l}^{2} cos(wt)] sin{\theta}$. Variabelsubstitution (annat än tidigare) $t' = {\omega}t \Rightarrow dt' = {\omega}dt$ vilket ger oss $ {\omega}^{2} \ddot{\theta} = -{\gamma}{\omega}\dot{\theta}-[\frac{g}{l} - \frac{A\omega}{l}^{2} cos(t')] sin{\theta} \iff \newline \ddot{\theta} = -\frac{\gamma}{\omega}\dot{\theta}-[\frac{g}{l{\omega}^{2}} - \frac{A}{l} cos(t')] sin{\theta}$.
Sätter vi $ x1 = {\theta}, x2 = \dot{\theta} $ får en ett system av första ordnings differentialekvationer: \newline
$\begin{cases}
x1'= x(2) \\
x2' = [\frac{g}{l{\omega}^{2}} \mp \frac{A}{l} cos(t')] sin(x(1)) -\frac{\gamma}{\omega} x(2) \end{cases}$ \newline
för fallet utan friktion är ekvationen:
\begin{equation}
l\ddot{\theta} = -[g + \ddot{\epsilon}(t)] sin{\theta} \iff
\end{equation} Efter variabelsubstiturion får vi
$ \ddot{\theta} = -[\frac{g}{l{\omega}^{2}} \mp \frac{A}{l} cos(t')] sin{\theta} $ \newline
OBS: vi kan även ha både plustecken och minustecken framför oscillationstermen eftersom funktionen är periodisk så kommer vi principiellt få samma beteende oavsett tecknet. \newline
Tar vi $ {\alpha} = \frac{g}{l{\omega}^{2}} $ och $ {\epsilon} = \frac{A}{2l} $ och $ {\beta} = \frac{\gamma}{\omega} $ får vi tabellen nedan som visar dimensionslösa differential ekvationer
\begin{center}
\begin{tabular}{| l | l | l |}
\hline
& utan friktion & med friktion \\ \hline
Vertikal oscillating pivot & $ \ddot{\theta} = -[{\alpha} + 2{\epsilon}\cos(t')]sin({\theta}) $ & $ \ddot{\theta} = -{\beta}\dot{\theta}-[{\alpha} + 2{\epsilon}\cos(t')]sin({\theta}) $ \\ \hline
\end{tabular}
\end{center}
För linjärisering noterar vi först att för små värden på ${\theta}, \ $ för ${\theta}$ kring 0, har vi då $sin({\theta}) \approx 0$ och för ${\theta}$ kring $ {\pi}$ har vi då $ sin({\theta}) \approx {\pi}-{\theta} $ därför ändrar vi minus tecknet framför [ ] till plus. då får vi de linjäriserade diffrenetialekvationerna nedan:
\begin{center}
\begin{tabular}{| l | l | l |}
\hline
Vertikal oscillating pivot & utan friktion & med friktion \\ \hline
$ {\theta} \approx 0 $ & $ \ddot{\theta} = -[{\alpha} + 2{\epsilon}\cos(t')]{\theta} $ & $ \ddot{\theta} = -{\beta}\dot{\theta}-[{\alpha} + 2{\epsilon}\cos(t')]{\theta} $ \\ \hline
$ {\theta} \approx {\pi} $ & $ \ddot{\theta} = +[{\alpha} + 2{\epsilon}\cos(t')]({\theta}-{\pi}) $ & $ \ddot{\theta} = -{\beta}\dot{\theta}+[{\alpha} + 2{\epsilon}\cos(t')]({\theta}-{\pi}) $ \\ \hline
\end{tabular}
\end{center}
För ${\theta} \approx 0$ får vi systemet: \begin{equation}
\overline{\dot{x}} = \begin{bmatrix} x1' \\ x2' \end{bmatrix} = \begin{bmatrix} 0 & 1 \\ -{\alpha} - 2{\epsilon}\cos(t') & -{\beta} \end{bmatrix} \begin{bmatrix} x1 \\ x2 \end{bmatrix} = A(t')\overline{x}
\end{equation}
För ${\theta} \approx {\pi}$ får vi systemet:
\begin{equation}
\overline{\dot{x}} = \begin{bmatrix} x1' \\ x2' \end{bmatrix} = \begin{bmatrix} 0 & 1 \\ {\alpha} + 2{\epsilon}\cos(t') & -{\beta} \end{bmatrix} \begin{bmatrix} x1 \\ x2 \end{bmatrix} = A(t')\overline{x}
\end{equation}
där i fallet utan friktion tar vi $-{\beta} = 0$
Först tittar vi på den oscillerande pendulum med friktion med ${\theta} \approx {\pi}$ och en starthastighet $\approx 0$. Parametrarna här har samma värde givna värde som tidigare dvs. $l =$ längden $=2$ och amplituden $ 10\% $ av $l; A = 0.2 $. Vi testa med liten friktion ${\gamma} = -0.1$, samtidigt tar vi ett längre tidsintervall; $t'$ eller ${\tau} = [0, 300]$, för då blir lättare och säkrare att finna stabiliteten.\newline
Med hjälp av en for-loop på matlab som loopar genom ${\omega} \in [0, 150]$ får vi $\forall \ {\omega} \geq 31.6$ har vi stabilitet och för $ {\omega} < 31.6$ får vi instabila. alltså ${\omega} \geq 31.6$ är gränsen.
För fallet utan friktion när vi loopar som ovan för ${\omega} \in [0, 150]$ får vi samma svar dvs $ \forall {\omega} \geq 32 $ får vi stabilitet, se figur 6!
\begin{figure}[!ht]
\centering
\caption{Oscillerande pendulum med friktion kring ${\theta} \approx {\pi}$ och ${\omega} \approx 31.6 $}
\includegraphics[width=1\textwidth]{PivotPendWithFrikStartiPI.jpg}
\end{figure}
Nu ska vi starta med jämviktspunkten ${\theta} \approx 0$ och $\dot{\theta} \approx 0$ på oscillerande pendulum utan friktion för att hitta ett/flera intervall för ${\omega} = {\omega}_{res}$ där pendulum blir instabilt. Vi kan anta först att instabilitets gränsen vid start nära noll är ${\pi/4}$ och $-{\pi/4}$ och samtidigt tar vi tidsintervallet stort och låta Matlab loopar genom ${\omega} = [0, 20]$ med små steglängd. Som resultat får vi oilka intervall: \newline i) ${\omega}_{res} \in [4, 4.85]$, \ ii) ${\omega}_{res} \in [9.05, 9.10]$ och \ \ iii) ${\omega}_{res} \in [16.05, 20]$, se figur 7! \newline
Observera att vi kan hitta fler intervall för ${\omega}_{res}$ om vi loopar genom större intervall än ${\omega} = [0, 20]$
\begin{figure}[!ht]
\centering
\includegraphics[width=0.8\textwidth,height=0.27 \textwidth]{ResfrPivPenutanFrikN_rNoll.jpg}
\caption{Pivot pendulum utan friktion kring ${\theta}, \dot{\theta} \approx 0 $ och ${\omega} \approx 2 \sqrt[]{g/l}$}
\end{figure}
Med fallet med ett litet och specifik friktion ${\gamma} = -0.1$ fast med samma tidsintervallm och parametrar som tidigare får vi
${\omega}_{res} \in [4.05, 4.75]$ vilket kan tolkas att friktionen kan på sätt och viss kan försvåra instabiliteten, se figur 8!
\begin{figure}[!ht]
\centering
\includegraphics[width=0.8\textwidth,height=0.27 \textwidth]{ResPendWithFrik.jpg}
\caption{Pivot pendulum med friktion kring ${\theta}, \dot{\theta} \approx 0 $ och ${\omega} \approx 2 \sqrt[]{g/l}$}
\end{figure}
Vi kan välja ${\omega}_{res} = 2{\sqrt[]{g/l}}=4.4294$: som en specifik resonans frekvens (alltså 2 gånger den naturliga frekvensen) och försöka testa det med fallet utan pivot (både med friktion och utan friktion) med ${\theta}, \dot{\theta} \approx 0 $. Som relutat får jag ingen instabilitet (både med eller utan friktion) med given ${\omega}_{res} = 2{\sqrt[]{g/l}} $. Alltså det skilljer sig mellan pendulum med pivot och utan pivot i det här fallet. Pendulum utan pivot är alltid stabilt, det kan beror på att det saknar pivot.
\begin{figure}[!ht]
\centering
\caption{Utan pivot och friktion med ${\omega} \approx 2 \sqrt[]{g/l}$ kring ${\theta}, \dot{\theta} \approx 0$}
\includegraphics[width=0.8\textwidth,height=0.2 \textwidth]{Upp5UtanPivotUtanFrik.jpg}
\end{figure}
\begin{figure}[!ht]
\centering
\caption{Utan pivot och men med friktion med ${\omega} \approx 2 \sqrt[]{g/l}$ kring ${\theta}, \dot{\theta} \approx 0$}
\includegraphics[width=0.8\textwidth,height=0.25 \textwidth]{Upp5CompUtanPivotmedFrik.jpg}
\end{figure}
\subsection{Floquet teori}
Tillbaka till funktion (3) och funktion (4) och ekvationssystemen (5) och (6)
\begin{center}
\begin{tabular}{| l | l | l |}
\hline
Vertikal oscillating pivot & utan friktion & med friktion \\ \hline
$ {\theta} \approx 0 $ & $ \ddot{\theta} = -[{\alpha} + 2{\epsilon}\cos(t')]{\theta} $ & $ \ddot{\theta} = -{\beta}\dot{\theta}-[{\alpha} + 2{\epsilon}\cos(t')]{\theta} $ \\ \hline
$ {\theta} \approx {\pi} $ & $ \ddot{\theta} = +[{\alpha} + 2{\epsilon}\cos(t')]({\theta}-{\pi}) $ & $ \ddot{\theta} = -{\beta}\dot{\theta}+[{\alpha} + 2{\epsilon}\cos(t')]({\theta}-{\pi}) $ \\ \hline
\end{tabular}
\end{center}
De ekvationerna ovan är redan linjäriserade kring ${\theta}$.
Vi gör en variabelsubstitution ${\delta}={\theta}-{\pi}$ då ekvationen kan skrivas som
$\ddot{\delta}=-{\beta}\dot{\delta} + [{\alpha} + 2 {\epsilon}cos({\tau})]{\delta} $ där i fallet utan friktion sätter vi ${\beta} = 0. \newline $ För enkelhetens skull tar vi ${\tau} = t'$ dvs. tiden.
Vi observera att $A({\tau})$ på (5) och (6) är periodeiska med perioden $2{\pi}.\newline$
Enligt Floquet sats.352 [s.64 Sze-Bi Hsu]; om ${\Phi}({\tau})$ är en fundamental matris till ett peridisk linjärt system så är även ${\Phi}({\tau} + T)$ är en fundamental matris till samma system. Dessutom $\exists$ en konstant matris $R \in \mathbb{C}^{2*2}$ (samma storlek som A) och en periodisk matris $P({\tau}) = P({\tau} + T)$, $P({\tau}) \in \mathbb{C}^{2*2}$ sådant att $ {\Phi}({\tau}) =P(\tau)\mathrm{e}^{TR}$ där T är perioden $2{\pi}$ (samma period)
\newline
Då gäller i så fall att $\exists$ en icke-singular matris $C \in \mathbb{C}^{2*2}$ sådant att ${\Phi}({\tau} + T) = {\Phi}({\tau})C. \newline$
Enligt Sats.3.5.1 [s.62 Sze-Bi Hsu] kan C skrivas som som $\mathrm{e}^{RT}$ där $R \in \mathbb{C}^{2*2}$ är icke-singulär ty $C \in \mathbb{C}^{2*2}$ och icke-singulär. \newline
Vi låter ${\Phi}(0) = I_{2*2}$ där $I_{2*2}$ är identitet matrisen av samma storlek som A, då även $P(0) = I_2 = P(T)$ dvs. ${\Phi}(T) = \mathrm{e} ^{RT}$ vilket $ \implies {\Phi}(0)^{-1} {\Phi}(T) = {\Phi}(T) = \mathrm{e}^{RT} $, Detta kallas Monodromy matris. I Det allmänna fallet har monodromy matris formen
${\Phi}({\tau})^{-1} {\Phi}({\tau}+T) = \mathrm{e}^{RT}. \newline $
Egenvärdena av Monodromy matrisen ${\mu}_1$, ${\mu}_2$ kallas karakteristiska multiplar och egenvärdena till vår matris $R,\ {\rho}_1$, ${\rho}_2$ kallas karakteritiska exponenterna. Vi använder nu Abels sats som säger oss om ${\Phi}({\tau})$ är en lösning till $\bar{x'} = A({\tau}) \bar{x} $ till (5) och (6). så är
$ \newline det{\Phi}(t + s) = det{\Phi}(s)\mathrm{e}^{\int_s^t trace(A(u))du}.\newline$
I fallet med friktion på numeriks beräknigen tar vi tiden $T \in [0, 2{\pi}]$ med 100 steglängd där egentligen plockar vi ut den sista dvs. matrisen vi får när tiden når $2{\pi}$. Detta är en $2*2$ matris. \newline
Låter vi ${\alpha} = \frac{g}{l{\omega}^{2}}$ och ${\epsilon} = \frac{A}{2l}$ loopar genom $I_{\alpha} \in [0, 5]$ och $I_{\epsilon} \in [-3, 3]$ med 100 steg på varje får vi 10 000 olika monodromy matriser.\newline
Vi vill nu hitta de karektaristiska multiplarna som är de
två lösningarna till $\newline det({\Phi}(T)-{\mu}I^{*})=0 \newline $
$\newline det({\Phi}(T)-{\mu}I^{*})= \ \
\begin{vmatrix}
{\phi}_{11} - {\mu} & {\phi}_{12} \\ {\phi}_{21} & {\phi}_{22} - {\mu}
\end{vmatrix} = {\mu}^{2} - {\mu}({\phi}_{11} + {\phi}_{22}) + {\phi}_{11}{\phi}_{22} - {\phi}_{12}{\phi}_{21} = {\mu}^{2} - {\mu}* trace({\Phi}(T)) + det({\Phi}(T)) = 0. \newline \newline$
Vi löser andragradsekvationen: \newline
$ {\mu}_{1,2} = \frac{1}{2}[trace({\Phi}(T)) \pm \sqrt[]{trace({\Phi}(T))^{2} - 4det({\theta}(T))}].\newline$
Vi låter nu ${\varphi} = trace({\Phi}(T))$ och eftersom $det({\Phi}(T)) = \mathrm{e}^{-{\beta}T}$ så har vi
$ \newline {\mu}_{1,2} = \frac{1}{2}[{\varphi} \pm \sqrt[]{{\varphi}^{2} - 4 \mathrm{e}^{-{\beta}T}}].\newline$
Vi får sammanlagt 10000 st ${\mu}_1, {\mu}_2$ och utav dem plockar vi endast sådana som uppfyller kravet på nedanstående sats.
Enligt Sats 3.5.5 [s.65 Sze-Bi Hsu] så är lösningen till (5) och (6) asymptotiskt stabil om och endast om $ |{\mu}_i| \leq 1$
$ i = 1,2. \newline $
Startar vi i ${\theta} = 0$ utgår vi från systemet (5).\newline
Startar vi i ${\theta} = {\pi}$ utgår vi från ekvations systemet (6). alltså skillnaden blir bara tecknet. \newline Vi tar ${\beta}=\frac{\gamma}{\omega} = 0.01$; alltså en fix friktion för enkelhetens skull på båda fallen. \newline
Som numeriska resultat får vi med start på ${\theta} = {\pi}$ stabilitet på $\newline {\alpha} = (0,\ 0.0505,\ 0.101,\ 0.1515,\ 0.2525,\ 0.3030,\ 0.4545)$, där den minsta ${\alpha}< 0.0505$; alltså nära noll.
${\alpha}=\frac{g}{l{\omega}^{2}} \implies$ ${\omega} = \sqrt[]{\frac{g}{l{\alpha}}}$. Tar vi längden som 2 motsvarar små ${\alpha}$ (nära noll) mot ${\omega} > 32$. Alltså får vi samma svar som uppgift 4.\newline
Stabilitet på $ {\epsilon} = (\pm 2.212,\pm 0.5757,\pm 0.4545,\pm 0.3939,\pm 0.3333, \pm 0.2727,\newline \pm 0.2121,\pm 0.1515,\pm 0.0909,\pm 0.0303)$ Vi får ett symmetrisk intervall på ${\epsilon}-leden. \newline$ ${\epsilon} = \frac{A}{2l} \implies amplituden A = {\epsilon}*2l$.
Vi plockar den minsta ${\epsilon} = \pm 0.0303$ med samma längd = 2 får vi minsta amplituden för stabiliteten som $A = \pm 0.1212. \newline$
\bigbreak
Som numeriska med start på $ {\theta} = 0$ (samma intervall för ${\alpha}$ och ${\epsilon}$ och samma friktion) får vi nedanstående stabilitets diagrammet. \newline
Tittar vi noggrant på bilden ser vi att vi har inga prickor mellan intervall
${\alpha}\in (0.2, 0.3) $. \newline
Enligt uppgift 5 (fallet med friktion) fick vi en instabilitets intervall (vid start ${\theta} = 0$) ${\omega}_{res} \in [4.05; 4.75]. \newline$
Här har vi ${\alpha}=\frac{g}{l{\omega}^{2}}$. Stoppa vi de givna $[4.05; 4.75]$ istället för ${\omega}$ får vi ett intervall för ${\alpha} = [0.2174; 0.299]$ vilket matchar väl med figur12 ty på bilden mellan intervall ${\alpha}\in (0.2, 0.3) $ finns inga prickor. \newline
Plockar vi den minsta ${\epsilon}$ dvs. $|{\epsilon}| = 0.0303$. ${\epsilon} = \frac{A}{2l} \implies$ amplituden $A = {\epsilon}*2l = 0.1212. \newline$
Dock är jag inte $100\%$ säker på att $A = 0.1212$ är den minimala amplituden som är nödvändligt för att observera instabilitets egenskapen i det här fallet.
\begin{figure}[!ht]
\centering
\caption{De markerade områden är kombination av ${\alpha}$ och ${\epsilon}$ där vi har stabilitet kring ${\theta} = {\pi}$ }
\includegraphics[width=1\textwidth,height=0.5 \textwidth]{up7smallFrikUpward.jpg}
\end{figure}
\begin{figure}[!ht]
\centering
\caption{De markerade områden är kombination av ${\alpha}$ och ${\epsilon}$ där vi har stabilitet kring ${\theta} = 0$}
\includegraphics[width=1.1\textwidth,height=0.5 \textwidth]{up8smallFrikDownward.jpg}
\end{figure}
\bigbreak
\bigbreak
\bigbreak
\section{Matlab codes}
\begin{verbatim}
%% 2 Free pendulum med friktion
clc; clear all;
% for system parametets
%g = 9.81;
%L = 5.5; % perioden beror på längd och friktion
Gm = 0.1; % friktion
f = @(t, x) [x(2); -Gm*x(2) - sin(x(1))];
y1 = linspace(-2, 8, 50);
y2 = linspace(-2, 2, 50);
[x, y] = meshgrid(y1, y2);
u = zeros(size(x));
v = zeros(size(x));
t = 0;
for i = 1: numel(x)
Yprim = f(t, [x(i); y(i)]);
u(i) = Yprim(1);
v(i) = Yprim(2);
end
quiver(x, y, u, v, 'black');
figure(gcf)
title('Free pendulum med friktion med startvinklar [0 pi/2 (3/4)*pi pi]
och starthastigheter [0 pi/8]')
xlabel('positioner i radien')
ylabel('hastigheter')
axis tight equal;
hold on
for y21 = [0 pi/8]; % för olika hastigheter från 0 till pi
for y20 = [0 pi/2 (3/4)*pi pi] % för olika vinklar/positioner från 0 till pi
[t, X] = ode45(f, linspace(0, 30), [y20; y21]);
plot(X(:, 1), X(:, 2))
plot(X(1, 1), X(1, 2), 'bo') %starting point
plot(X(end, 1), X(end, 2), 'ks') % ending point
end
end
hold off
% 2 olika typer
%% 1 free pendulum utan friktion
clc; clear all;
% for system parametets
%g = 9.81;
%L =3;
f = @(t, Y) [Y(2); -sin(Y(1))]; % def. of diff. ekv.
y1 = linspace(0, 25, 20);
y2 = linspace(-2, 2, 20);
[x, y] = meshgrid(y1, y2);
u = zeros(size(x));
v = zeros(size(x));
t = 0;
for i = 1: numel(x)
Yprim = f(t, [x(i); y(i)]);
u(i) = Yprim(1);
v(i) = Yprim(2);
end
quiver(x, y, u, v, 'black');
figure(gcf)
title('Free pendulum utan friktion med olika vinklar [0 (3/4)*pi pi]
och hastigheter [0 pi/8]')
xlabel('position')
ylabel('hastighet')
axis tight equal;
hold on
for y21 = [0 pi/8] % hastigheten thetha prim
for y20 = [0 (3/4)*pi pi] % vinkel theta eller positionen
[t, X] = ode45(f, linspace(0, 20), [y20; y21]); % tiden tar slut,
vi kan testa med längre tider med.
plot(X(:, 1), X(:, 2))
plot(X(1, 1), X(1, 2), 'bo') %starting point
plot(X(end, 1), X(end, 2), 'ks') % ending point
end
end
hold off
% 4 olika typer ??
figure(2)
plot(t, X(:, 1))
%% 3 linearization eq. utan friction
clc; clear all;
Y(1) = 0;
f = @(t, Y) [Y(2); -Y(1)]; % def. of diff. ekv.
y1 = linspace(-1, 1, 20);
y2 = linspace(-1, 1, 20);
[x, y] = meshgrid(y1, y2);
u = zeros(size(x));
v = zeros(size(x));
t = 0;
for i = 1: numel(x)
Yprim = f(t, [x(i); y(i)]);
u(i) = Yprim(1);
v(i) = Yprim(2);
end
quiver(x, y, u, v, 'black');
figure(gcf)
title('Free pendulum utan friktion med start vinklar [0 pi/4]
och start hastigheter [0 pi/8]')
xlabel('positioner mätt i radien')
ylabel('hastigheter')
axis tight equal;
hold on
for y21 = [0 pi/8]
for y20 = [0 pi/4]
[t, X] = ode45(f, linspace(0, 20), [y20; y21]); % tiden tar slut,
vi kan testa med längre tider med.
plot(X(:, 1), X(:, 2))
plot(X(1, 1), X(1, 2), 'bo') %starting point
plot(X(end, 1), X(end, 2), 'ks') % ending point
end
end
hold off
% bara en typ av lösning
%% 3. sadellpunkt
dxdt = @(t, X) [X(2); X(1)];
% Producing vectorfield: -------------------------------------
y1 = linspace(-2,2,20);
y2 = linspace(-2,2,20);
[x,y] = meshgrid(y1,y2);
u = zeros(size(x));
z = zeros(size(x));
t=0;
for i = 1:numel(x)
fx = dxdt(t,[x(i); y(i)]); % x' = f(x, t)
u(i) = fx(1);
z(i) = fx(2);
end
quiver(x,y,u,z,'r');
hold on
axis tight equal;
axis([-2, 2, -2, 2])
H = linspace(0, 10);
for x10 = linspace(0, pi, 5)
for x20 = linspace(0, pi, 5)
initial = [x10, x20];
[ts,x] = ode45(dxdt,H,initial); % Lek runt med initialvärden
och H för att rita ut lösningar.
plot(x(:,1),x(:,2), 'b')
end
end
hold off
%% 3. Linearization eq. with friction
clc; clear all;
Gm = 0.1; % friktion
f = @(t, x) [x(2); -Gm*x(2) - x(1)];
y1 = linspace(-4, 4, 20);
y2 = linspace(-4, 4, 20);
[x, y] = meshgrid(y1, y2);
u = zeros(size(x));
v = zeros(size(x));
t = 0;
for i = 1: numel(x)
Yprim = f(t, [x(i); y(i)]);
u(i) = Yprim(1);
v(i) = Yprim(2);
end
quiver(x, y, u, v, 'black');
figure(gcf)
title('Free pendulum med friktion med start vinklar [0 pi/2 (3/4)*pi pi]
och start hastigheter [0 pi/8]')
xlabel('positioner mätt i radien')
ylabel('hastigheter')
axis tight equal;
hold on
for y21 = [0 pi/8]
for y20 = [0 pi/2 (3/4)*pi pi]
[t, X] = ode45(f, linspace(0, 30), [y20; y21]);
plot(X(:, 1), X(:, 2))
plot(X(1, 1), X(1, 2), 'bo') %starting point
plot(X(end, 1), X(end, 2), 'ks') % ending point
end
end
hold off
% bara en typ av lösning
%stabilt spiral, kolla räkningarna där du fick egenvärdena och
%egenvektorerna
%% 4. pendulum with priodically oscillating pivot utan friktion close
%% to upwards, we check the stability For-Loop
clc; clear all;
g = 9.81;
A = 0.2;
L = 2; % perioden beror på längd och friktion
tht = (15/16)*pi;
thtP = 0;
t = linspace(0, 300, 1000);
w1 = [];
% w = 31.5999; % for att frekvens lika eller större än 32 får vi stabilisitet om
% we börjar uppwards med hastogheten=0 och positionen nära pi
for w = 0:150
f = @(t, x) [x(2); - sin(x(1)) + ((A*w^2)/g)*cos(t*w*sqrt(L/g))*sin(x(1))];
[T, x] = ode45(f, t, [tht; thtP]);
condition = (3/4)*pi < x(:, 1) & x(:, 1) < (5/4)*pi;
if all(condition)
w1 = [w1, w];
end
end
subplot(2,1,1)
plot(x(:, 1), x(:, 2))
title('Pivot pendulum med friktion=-0.1 och startpunkten nära pi och start hastigheter 0')
xlabel('Tiden')
ylabel('Positioner')
%% uppgift 5; pendulum with vertically oscillation (utan friction) close
%% to downward. we check the instability
clc; clear all;
g = 9.81;
A = 0.2;
L = 2;
tht = 0.1; % start-positionen/vinkel
thtP = 0.1; % theta prim/hastigheten
w1 = [];
t = linspace(0, 3000, 1000);
%w = sqrt(g/L)*2;
for w = 0:0.05:20 % värdena 4. 4.1 4.2 4.3 4.4
f = @(t, x) [x(2); -g/(L*w^2)*sin(x(1)) + (A/L)*cos(t)*sin(x(1))];
[T, x] = ode45(f, t, [tht; thtP]);
condition = pi/4 < x(:, 1) | x(:, 1) < -pi/4;
if (any(condition))
w1 = [w1, w];
end
end
subplot(2,1,1)
plot(T, x(:, 1))
title('Resonansfrekvensen=4.4294 för Pivot pendulum, utan friktion med
startpunkten och start hastigheter nära noll')
xlabel('Tiden')
ylabel('Positioner')
%% 5 plot för instability resonans (vertically oscillation med friction)
%% close to downward.
clc; clear all;
g = 9.81;
A = 0.2;
L = 2;
Gm = -0.1;
tht = 0.1; % start-positionen/vinkel
thtP = 0.1; % theta prim/hastigheten
t = linspace(0, 3000, 1000);
w1 = [];
w = sqrt(g/L)*2;
%for w = 0:0.05:20
f = @(t, x) [x(2); (Gm/w)*x(2)-g/(L*w^2)*sin(x(1)) + (A/L)*cos(t)*sin(x(1))];
[T, x] = ode45(f, t, [tht; thtP]);
% condition = pi/4 < x(:, 1) | x(:, 1) < -pi/4;
% if (any(condition))
% w1 = [w1, w];
% end
%end
subplot(2,1,1)
plot(T, x(:, 1))
title('Pivot pendulum med friktion=-0.1 och startpunkten nära noll och
start hastigheter nära noll')
xlabel('Tiden')
ylabel('Positioner')
%% 5. Jämförelse med free pendulum utan friktion och med friktion lägg till
%% resonans friktionen Wres = 2sqrt(g/L)
clc; clear all;
% for system parametets
g = 9.81;
L =2;
force = 0.1;
Wres = 2*sqrt(g/L);
f = @(t, Y) [Y(2); -force*Y(2)-Wres*sin(Y(1))]; % def. of diff. ekv.
y1 = linspace(0, 25, 20);
y2 = linspace(-2, 2, 20);
[x, y] = meshgrid(y1, y2);
u = zeros(size(x));
v = zeros(size(x));
t = 0;
for i = 1: numel(x)
Yprim = f(t, [x(i); y(i)]);
u(i) = Yprim(1);
v(i) = Yprim(2);
end
quiver(x, y, u, v, 'black');
figure(gcf)
title('Free pendulum när vi lägger till resonans frekvensen (med friktion) med
olika start vinkel och starthastigheten nära noll')
xlabel('position')
ylabel('hastighet')
axis tight equal;
hold on
for y21 = [0.1] % hastigheten thetha prim
for y20 = [0.1] % vinkel theta eller positionen
[t, X] = ode45(f, linspace(0, 3000, 1000), [y20; y21]); % tiden tar slut,
vi kan testa med längre tider med.
plot(X(:, 1), X(:, 2))
plot(X(1, 1), X(1, 2), 'bo') %starting point
plot(X(end, 1), X(end, 2), 'ks') % ending point
end
end
hold off
% 4 olika typer ??
figure(2)
plot(t, X(:, 1))
title('Free pendulum när vi lägger till resonans frekvensen (med friktion) med
olika start vinkel och starthastigheten nära noll')
xlabel('tid')
ylabel('position')
%% 6 pendulum with vertically oscillation (without friction) close
%% to downward. we check the instability
clc; clear all;
g = 9.81;
A = 0.2;
L = 2;
tht = 0.1; %positionen/vinkel
thtP = 0.1; %theta prim/hastigheten
t = linspace(0, 5000,5000);
w = 100; %w = sqrt(g/L)*2; %för w >= 151 börjar det hända nånting
f = @(t, x) [x(2); -g/(L*w^2)*x(1) + (A/L)*cos(t)*x(1)];
[T, x] = ode45(f, t, [tht; thtP]);
subplot(2,1,1)
plot(T, x(:, 1))
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function DY = upp7(t, y, a, eps)
%g = 9.81;
%Gm = -0.01;
%L = 2;
%A = 0.2;
%w = 155;
%w = sqrt(g/(a*L));
%A = eps * L;
%a = g/(L*w^2);
%eps = A/L;
%force = Gm/w;
force = -.025;
DY = zeros(2,1);
DY(1) = y(2);
DY(2) = force*y(2)+(a+2*eps*cos(t))*y(1); %på uppgift 8 tag -(g/(L*w^2) istället
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% uppgift7 linjäriserade vertically oscillation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% with friction around the upward, theta lika med pi
tic
T = linspace(0, 2*pi, 100);
a = linspace(0, 5, 100);
eps = linspace(-3, 3, 100);
a1 = [];
eps1 = [];
for k = 1:length(a)
for j = 1:length(eps)
[t,y1] = ode45(@(t,y) upp7(t, y, a(k), eps(j)),T, [1, 0]);
[t,y2] = ode45(@(t,y) upp7(t, y, a(k), eps(j)),T, [0, 1]);
Phi = [y1(end,:); y2(end,:)];
% Here we calculate the characteristic multipliers,
% i.e. the eigenvalues
mu1(k, j) = 1/2*(trace(Phi))+(sqrt((((1/2)*trace(Phi))^2-det(Phi)))); %100x100
mu2(k, j) = 1/2*(trace(Phi))-(sqrt((((1/2)*trace(Phi))^2-det(Phi)))); %100x100
% Here we check the stability conditions
if (abs(mu1(k, j)) <= 1 && abs(mu2(k, j)) <=1 )
plot(a(k), eps(j), '*');
a1 = [a1, a(k)];
eps1 = [eps1, eps(j)];
hold on
%elseif abs(trace(Phi)) <= 2
%plot(a(k), eps(j), '*');
%a1 = [a1, a(k)];
%eps1 = [eps1, eps(j)];
% skit(j, k) = eps(j);
% hold on
%else
end
end
end
xlabel('a');
ylabel('eps');
axis([-2,5,-3,3])
toc
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function DY = upp8(t, y, a, eps)
%g = 9.81;
%Gm = -0.01;
%L = 2;
%A = 0.2;
%w = 155;
%w = sqrt(g/(a*L));
%A = eps * L;
%a = g/(L*w^2);
%eps = A/L;
%force = Gm/w;
force = -.001;
DY = zeros(2,1);
DY(1) = y(2);
DY(2) = force*y(2)-(a+2*eps*cos(t))*y(1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 8. linearized with vertically oscillation pivot with small
%% friction around the downward theta = 0
tic
T = linspace(0, 2*pi, 100);
a = linspace(0, 5, 100);
eps = linspace(-3, 3, 100);
a2 = [];
eps2 = [];
for k = 1:length(a)
for j = 1:length(eps)
[t,y1] = ode45(@(t,y) upp8(t, y, a(k), eps(j)),T, [1, 0]); % y1 = 100x2
[t,y2] = ode45(@(t,y) upp8(t, y, a(k), eps(j)),T, [0, 1]); % y2 = 100x2
Phi = [y1(end,:); y2(end,:)]; % den är monodromy matris
% Here we calculate the characteristic multipliers,
% i.e. the eigenvalues
mu1(k, j) = 1/2*(trace(Phi))+(sqrt((((1/2)*trace(Phi))^2-det(Phi)))); % 100x100
mu2(k, j) = 1/2*(trace(Phi))-(sqrt((((1/2)*trace(Phi))^2-det(Phi)))); % 100x100
% Here we check the stability conditions
if (abs(mu1(k, j)) <= 1 && abs(mu2(k, j)) <=1 )
plot(a(k), eps(j), 'o');
a2 = [a2, a(k)];
eps2 = [eps2, eps(j)];
hold on
% elseif abs(trace(Phi)) <= 2
%plot(a(k), eps(j), 'o');
%a2 = [a2, a(k)];
%eps2 = [eps2, eps(j)];
% hold on
% else
end
end
end
xlabel('a');
ylabel('eps');
axis([-2,5,-3,3])
toc
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function dy = mathieu(t, y, a, eps) %ft ,f,gt,g)
dy = [y(2); -a - eps * cos(t)];
%function dydt = myode(t,y,ft,f,gt,g)
%[T y] = ode45(@(t,y) mathieu(t,y,ft,f,gt,g),Tspan,IC); % Solve ODE
%[T,y1]=ode45(@(t2,v)mathieu(t2,v,a(k),eps(j)),t,v1);
%dy = [y(2); -a - eps * cos(t)];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\end{verbatim}
%Avsluta med en innehållsförteckning
\tableofcontents
\end{document}