Estudo de um sistema dinâmico.
\documentclass[brazilian, 16pt, a4]{article}
\usepackage[brazil]{babel}
\usepackage[utf8]{inputenc}
\usepackage{verbatim}
\usepackage{float}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{graphicx}
\usepackage{color}
\usepackage{url}
\pagestyle{empty}
\hoffset 0 cm
\evensidemargin 0 cm
\oddsidemargin 0 cm
\setlength{\textwidth}{16 cm}
\title{Trabalho 5 \\ Métodos Computacionais da Física B}
\author{Esteban Gerling - 00231223}
\begin{document}
\maketitle
\section*{Introdução}
\paragraph{}Este trabalho é sobre o estudo de um mapa, em específico sobre o mapa de \textit{Arnold tongue}\cite{arnold-wiki} descrito na equação (1) apresentada abaixo:
\begin{equation}\label{eq1}
\theta_{n+1}= \theta_n + \omega - \frac{k}{2 \pi}sen(2 \pi \theta_n ),
\end{equation}onde o parâmetro k será estudado no intervalo de $0$ até $4\pi$, os valores de $\theta_n$ serão estudados no intervalo $[0,1]$ e o parâmetro $\omega$, a princípio, será um valor fixo $\omega = 1/3$.
\\
\paragraph{}Uma aplicação do problema descrito acima é no sincronismo de diodos túnel (\textit{resonant-tunneling diode}), que é um tipo de semicondutor extremamente rápido capaz de operar em frequências da ordem de $GHz$ através da utilização de um efeito da mecânica quântica chamado de tunelamento\cite{arnold-wiki}\cite{artigo}.
\paragraph{}O parâmetro $\omega$ utilizado neste estudo será fixado, conforme informado acima, para obter-se um \textit{mode-locking}, parâmetro este que está relacionado com o ruído de sinal\cite{arnold-wiki}, e estudar a resposta do sistema à variação do parâmetro $k$.
\paragraph{}
\section{Encontrando os pontos fixos de primeira ordem}
\paragraph{}Para encontrar os pontos fixos precisamos que $\theta_{n+1}=\theta_n$, ou seja:
\begin{equation}\label{eq2}
\omega = \frac{k}{2 \pi}sen(2 \pi \theta_n ).
\end{equation}
\paragraph{}Resolvendo a equação para $\theta_n$ temos:
\begin{equation}\label{eq3}
\theta_n = \frac{1}{2\pi}arcsen\left(\frac{2 \pi \omega}{k}\right)=\theta^*.
\end{equation}
\paragraph{}Para estes valores de $\theta_n=\theta^*$ encontramos pontos fixos de primeira ordem dependentes de $k$.
\subsection{Determinando o intervalo de estabilidade}
\paragraph{}Para encontrarmos o intervalo de estabilidade dos pontos fixos temos que o módulo da derivada da função $\theta_{n+1}=f(\theta_n)$ deve ser menor do que $1$, ou seja, $f'(\theta_{n})$deve estar entre $-1$ e $1$. Isto deve-se ao fato de que analisando a função em um ponto fixo, quando ocorre uma pequena perturbação ($\theta^*$ vai para um ponto $\theta_n + \epsilon_n$) a função retorna para um ponto fixo $\theta^*$, onde $\epsilon_n$ é uma pequena perturbação. Assim seguindo as iterações para o próximo valor de $\theta_n$, o próximo ponto seria $\theta_{n+1}+\epsilon_{n+1}$
\paragraph{}Pode-se verificar, através de uma expansão de $f(\theta_n+\epsilon_n)$ em torno de $\theta_n=\theta^*$, que para que exista estabilidade deve ocorrer que $|\epsilon_n|>|\epsilon_{n+1}|$, ou seja, quando $n$ aumenta a função se aproxima do ponto fixo $\theta^*$.
\paragraph{}Temos então que a derivada da função \ref{eq1} a ser estudada é:
\begin{equation}\label{eq4}
f'(\theta_n)=1-k(cos(2\pi\theta_n)).
\end{equation}
\paragraph{} Substituindo a equação \ref{eq3} em \ref{eq4} e fazendo o módulo menor do que $1$ temos o intervalo de estabilidade para pontos fixos de primeira ordem:
\begin{equation}\label{eq5}
|1-k[cos\left(arcsen\left( \frac{2\pi \omega}{k} \right)\right)]|<1.
\end{equation}
\paragraph{}Para encontrar valores de $k$ que satisfaçam esta inequação, foi feito um programa que varia o parâmetro $k$ no intervalo $[0, 4\pi]$ e testa a equação \ref{eq5} para cada $k$ e também testa se para estes valores a variável $\theta_n$ está dentro do domínio a ser estudado que é $[0,1]$, e guarda os valores de $k$ e $\theta_n$ para os quais estes testes são verdadeiros.
\paragraph{}A base do programa em $C$ é a seguinte:
\begin{verbatim}
for (k=(0);k<=(4*M_PI);k=k+0.01)
{
estab=modulo(x,k,omega);
dom=intervalo(x,k,omega);
if((estab<1)&&(dom>0)&&(dom<1))
{
fprintf(arq1,"%lf %lf\n",dom,k);
n++;
}
printf("%f\n",k);
}
fclose(arq1);
\end{verbatim}
\paragraph{}Onde a variável $x$ é o $\theta_n$ e as funções $estab$ e $dom$ são funções das equações \ref{eq5} e \ref{eq3} respectivamente, são as seguintes:
\begin{verbatim}
double modulo (double a, double b, double c) //a=x, b=k, c=omega;
{
return (fabs(1-((b)*(cos(asin((2*M_PI*c)/(b)))))));
}
double intervalo (double d, double e, double f) //d=x, e=k, f=omega;
{
return ((1.0/(2.0*M_PI))*(asin((2*M_PI*f)/(e))));
}
\end{verbatim}
\paragraph{}Para valores que o módulo do termo $arcsen\left( \frac{2\pi \omega}{k} \right)$ resultaria em um valor maior do que $1$ não existe resultado, pois o valor do $seno$ deve estar compreendido no intervalo $[-1,1]$, isto ocorre para valores de $k$ entre $-2.086371$ e $2.093629$.
\paragraph{}Com a execução do programa verificou-se que o intervalo de estabilidade é entre $k > \frac{2\pi}{3}$ e $k \approx 0.92\pi$.
\begin{figure}[H]
\begin{center}
\includegraphics[scale=0.5, angle = -90]{figura2.eps}
\caption{Gráficos da variação dos valores de $k$ para verificar $\theta$ dentro do intervalo do domínio $[0,1]$. Verificação de pontos fixos}\label{fig2}
\end{center}
\end{figure}
\paragraph{}Na figura abaixo podemos verificar, para valores iniciais de $\theta_n=0.1$, o que ocorre conforme se varia $k$ dentro do intervalo de estabilidade quando $n$ aumenta.
\begin{figure}[H]
\begin{center}
\includegraphics[scale=0.5, angle = -90]{figura1.eps}
\caption{Gráficos da variação de $k$ com $n$ para valores de $k$ dentro do intervalo de estabilidade.}\label{fig1}
\end{center}
\end{figure}
\paragraph{}Conforme verificado na Figura \ref{fig1}, para um número $n$ grande de iterações, os valores de $\theta_n$ convergem para uma solução assimptótica. Foi executado o programa para $n=1000$, porém só foi plotado o intervalo de interesse $n$ entre $0$ e $50$, pois para valores maiores de $n$ podia-se ver que a solução era assimptótica, mas não se enchergava a parte interessante que é os valores de $\theta_n$ oscilando e convergindo.
\paragraph{}Na Figura \ref{fig4} abaixo pode ser visto o que ocorre quando $k=0.92\pi$ (muito próximo ao limite do intervalo de estabilidade):
\begin{figure}[H]
\begin{center}
\includegraphics[scale=0.5, angle = -90]{figura4.eps}
\caption{Gráficos da evolução de $\theta_n$ com $n$ para valor de $k=0.92\pi$ - próximo ao limite do intervalo de estabilidade.}\label{fig4}
\end{center}
\end{figure}
\paragraph{}Como esperado, a função converge para uma solução assimptótica, mesmo que demorando mais tempo.
\paragraph{}Na Figura \ref{fig3} abaixo, pode-se ver o que acontece quando $k$ está fora do intervalo de estabilidade determinado anteriormente ($k<\frac{2\pi}{3}$):
\begin{figure}[H]
\begin{center}
\includegraphics[scale=0.5, angle = -90]{figura3.eps}
\caption{Gráfico da variação de $k$ versus $n$ para valores de $k$ fora do intervalo de estabilidade. $k=2.0$, $ k<\frac{2\pi}{3}$}\label{fig3}
\end{center}
\end{figure}
\paragraph{}Como era de se esperar, para valores de $k$ fora do intervalo determinado de estabilidade, a função tende à situação caótica quando $n$ aumenta. Neste caso até $n=10$ a função estava dentro do domínio, na próxima iteração os valores começaram a crescer de forma instável (situação caótica).
\section{Encontrando Pontos fixos de $2^a$ ordem}
\paragraph{}Para encontrar pontos fixos de segunda ordem, fazemos $f(f(\theta_n)) = \theta_n$, que significa $\theta_{n+2}=\theta_n$, e para esta igualdade encontramos regimes cíclicos e sua estabilidade conforme foi feito na seção anterior.
\paragraph{}Assim fazemos a substituição $f(f(\theta_n))$:
\begin{equation}\label{eq6}
f(f(\theta_n)) = \theta_n +2\omega - \frac{k}{2\pi} \left[sen(2\pi\theta_n)+sen(2\pi(\theta_n+\frac{k}{2\pi}sen(2\pi\theta_n)))\right],
\end{equation}
\paragraph{}Igualando $f(f(\theta_n)) = \theta_n$:
\begin{equation}\label{eq7}
\frac{4\pi\omega}{k} = sen(2\pi\theta_n)+sen(2\pi(\theta_n+\frac{k}{2\pi}sen(2\pi\theta_n))).
\end{equation}
\paragraph{}Equação que se tornaria um tanto difícil, senão impossível, de isolar somente os termos $\theta_n$ em um lado da igualdade. Porém podemos ter uma ideia do que acontece com a função \ref{eq1} para valores de $k$ maiores do que o intervalo de estabilidade para os pontos fixos de primeira ordem determinado anteriormente se executarmos um programa que roda a função \ref{eq1} com valores de $k$ a partir de $0.92\pi$.
\paragraph{}O programa utilizado para tal tarefa é o mesmo identificado anteriormente na seção para verificação da estabilidade de pontos fixos, porém agora com o parâmetro $k$ variando de $0.92\pi$ até $4\pi$.
\paragraph{}Obteve-se os seguintes resultados para um dos primeiros valores de $k>0.92\pi$:
\begin{figure}[H]
\begin{center}
\includegraphics[scale=0.5, angle = -90]{figura5.eps}
\caption{Gráfico de $\Theta_n$ por $n$ com $k=2.99$, $k>0.92\pi$.}\label{fig5}
\end{center}
\end{figure}
\paragraph{}No gráfico da Figura \ref{fig5} pode-se verificar que $\theta_n$ oscila entre dois valores fixos quando $n$ aumenta, isto significa uma órbita de período $2$, que são os pontos fixos de segunda ordem. Com isso pode-se concluir que já para $k>0.92\pi$ obtém-se pontos fixos de segunda ordem.
\paragraph{}Vamos agora verificar até que valor de $k>0.92\pi$ tem-se soluções de órbita de segunda ordem. Para isto, uma forma mais prática de se encontrar o intervalo de estabilidade dos pontos fixos é fazendo o diagrama de bifurcações, no qual varia-se o parâmetro $k$ de $0$ até $4\pi$ e para cada valor individual de $k$ se faz um número grande de $n$ iterações, salvando-se os últimos passos de $\theta_n$, os quais já chegaram na situação assimptótica. Então \textit{plota}-se o gráfico dos valores de $k$ por $\theta_n$.
\section{Diagrama de Bifurcações}
\paragraph{} Com o diagrama de bifurcações podemos verificar para quais valores do parâmetro $k$ ocorrem bifurcações, que são aumento do período da órbita.
\paragraph{}A base do programa utilizado para a elaboração do diagrama de bifurcações foi o seguinte:
\begin{verbatim}
neq=700; //numero de iteracoes para cada valor de k
nprod=300; //quantidade de primeiros valores de theta que serao descartados
x0=0.1; //theta inicial
x=x0;
li=(0); // k inicial
lf=(4*M_PI); // k final
dl=(lf-li)/400.0; // divisao de intervalos de k
l=li;
while(l<lf) // while para a variacao de k
{
x=x0; // reinicia o valor do theta zero
n=1; // reinicia o valor do tempo n
while(n<neq) // inicia os calculos de theta
{
x=(x+(1.0/3.0)-((l/(2*M_PI))*(sin(2*M_PI*x))));
if(n>nprod) // so salva os 400 utlimos valores de theta para cada k
{
fprintf(arq1,"%d %.12lf %.12lf\n",n,x,l);
}
n=n+1;
}
l=l+dl;
}
\end{verbatim}
\paragraph{}A partir deste programa, foi gerado o seguintes gráfico de $\theta_n$ versus $k$:
\begin{figure}[H]
\begin{center}
\includegraphics[scale=0.5, angle = -90]{figura6.eps}
\caption{Diagrama de bifurcações de $\Theta_n$ por $k$.}\label{fig6}
\end{center}
\end{figure}
\paragraph{}Na Figura \ref{fig6} pode-se verificar o intervalo de convergência que foi encontrado na primeira seção, para $ \frac{2\pi}{3}<k< 0.92\pi$, e também que a partir de $k$ maior que este intervalo os valores de $\theta_n$ começam a oscilar entre duas soluções até $k \approx 3.2$, que é onde ocorre novamente bifurcação das soluções para quatro valores. Após esta bifurcação que ocorre no intervalo de $k$ $[3.2,3.3]$, o sistema entra em estado de caos.
\paragraph{}Na figura \ref{fig7} abaixo, podemos verificar o diagrama de bifurcações completo, para o intervalo de $k$ $[0,4\pi]$.
\begin{figure}[H]
\begin{center}
\includegraphics[scale=0.5, angle = -90]{figura7.eps}
\caption{Diagrama de bifurcações de $\Theta_n$ por $k$. Intervalo de $k$ $[0;4\pi]$}\label{fig7}
\end{center}
\end{figure}
\paragraph{}É possível notar que existem vários intervalos em que o sistema está numa situação caótica, porém existem ``ilhas" de estabilidade em alguns pontos em meio ao caos. Isto será especificado e aprofundado na próxima seção. Também pode-se notar que para o intervalo de $\theta_n$ entre $[-1,1]$ a região onde existem mais soluções assimptóticas são aquelas para as quais os valores de $k$ estão entre $2$ e $4$, pois é onde tem maior densidade de pontos.
\section{Expoente de Lyapunov}
\paragraph{}A taxa com que as distâncias entre duas trajetórias aumenta ou diminui com o tempo está relacionada com uma quantidade chamada de \textit{expoente de Lyapunov}.
\paragraph{}O expoente de Lyapunov é definido como:
\begin{equation}\label{lyapunov}
\lambda_L=\lim_{n \rightarrow \infty}\frac{1}{n}\sum\limits_{i=1}^{n}ln|f'(\theta_i)|.
\end{equation}
\paragraph{}A análise do valor de $\lambda_L$ é a seguinte:
se $\lambda_L>0$, as trajetórias vizinhas se distanciam umas das outras conforme o tempo $n$ avança, e caracterizam um comportamento caótico;
se $\lambda_L<0$, as trajetórias convergem para um valor fixo ou um limite cíclico, caracterizam estabilidade do sistema, elas se aproximam.
\paragraph{}Foi feito um programa para o cálculo do expoente de Lyapunov para a sistema estudado para o intervalo de $k$ de $[0,4\pi]$. O programa realiza o cálculo de 400 valores de $\theta_i$ e salva os últimos $300$, que é onde o sistema já alcançou algum ponto fixo ou período cíclico. O programa calcula o expoente de Lyapunov para cada valor de $k$. A base do programa é a seguinte:
\begin{verbatim}
while(l<lf) // while para variacao de k
{
x=x0; // retorna a posicao inicial
n=1;
while(n<neq) // inicia os neq=400 calculos de theta e lyapunov para um dado k
{
x=(x+(1.0/3.0)-((l/(2*M_PI))*(sin(2*M_PI*x)))); // calcula theta
lyap = lyap + log(fabs(derivada(x,l,omega))); // soma os ln da derivada de theta_i
if(n>nprod) // guarda somente os ultimos pontos de theta
{
if((fabs(x))<=1) fprintf(arq1,"%d %.12lf %.12lf\n",n,x,l);
}
n=n+1;
}
lyap = (lyap)/neq; // finaliza o calculo do expoente de lyapunov para este k
fprintf(arq2,"%.12lf %.12lf\n",l,lyap);
l=l+dl;
lyap=0.0; //zera o valor de lyapunov, para calcular o do proximo k
}
\end{verbatim}
\paragraph{}Onde a função \textit{derivada} corresponte à equação \ref{eq4} e é a seguinte:
\begin{verbatim}
double derivada (double d, double e, double f )
{
return (1-((e)*(cos(2*M_PI*d)))); //d=theta, e=k, f=omega;
}
\end{verbatim}
\paragraph{}Na figura \ref{fig8} podemos verificar a comparação do diagrama de bifurcações com o expoente de Lyapunov para cada valor de $k$, o que condiz com o que foi verificado anteriormente, que para os valores de $k$ em que há pontos fixos e fases cíclicas de $\theta_n$ o expoente de Lyapunov é negativo, e em alguns casos quando ele está na iminência de passar de $0$ e ficar positivo o sistema entra em bifurcação das soluções de $\theta_n$ (fase cíclica) e o expoente de Lyapunov fica negativo novamente.
\begin{figure}[H]
\begin{center}
\includegraphics[scale=0.5, angle = -90]{figura8.eps}
\includegraphics[scale=0.53, angle = -90]{figura7.eps}
\caption{Expoente de Lyapunov e o Diagrama de bifurcações de $\Theta_n$ por $k$. Intervalo de $k$ $[0;4\pi]$}\label{fig8}
\end{center}
\end{figure}
\paragraph{}Agora vamos mudar o \textit{range} dos gráficos para o intervalo estudado anteriormente, de $k$ entre $\frac{2 \pi}{3}$ e $0.92 \pi$:
\begin{figure}[H]
\begin{center}
\includegraphics[scale=0.48, angle = -90]{figura10.eps}
\includegraphics[scale=0.45, angle = -90]{figura9.eps}
\caption{Expoente de Lyapunov e o Diagrama de bifurcações de $\Theta_n$ por $k$. Intervalo de $k$ $[2;3.8]$}\label{fig9}
\end{center}
\end{figure}
\paragraph{}Conforme esperado, podemos verificar que as bifurcações condizem com expoentes de Lyapunov que estão muito próximos de zero e que para valores negativos do expoente de Lyapunov encontramos pontos fixos e fases cíclicas para $\theta_n$, e para valores positivos do expoente de Lyapunov encontramos o sistema em uma fase caótica.
\section{Conclusões Gerais}
\paragraph{}Após as análises realizadas, pode-se concluir que é possível prever fases cíclicas, pontos fixos e fases caóticas de um sistema dinâmico através de cálculos analíticos e verificar com cálculos numéricos através do diagrama de bifurcações e do cálculo do expoente de Lyapunov comprovando os resultados esperados. Mesmo que o sistema seja difícil, senão impossível, de ser estudado analíticamente, pode-se utilizar métodos de cálculo numérico para prever situações deste sistema.
\begin{thebibliography}{9}
\bibitem{arnold-wiki}
``Arnold tongue" Disponível em ``\url{http://en.wikipedia.org/wiki/Arnold_tongue}"
\bibitem{artigo}
Romeira, Bruno ``Chaotic Dynamics in Resonant Tunneling Optoelectronic Voltage Controlled Oscillators" Departmento de Fis., Univ. do Algarve, Faro, Portugal
\end{thebibliography}
\end{document}