Trabalho 5: Métodos Computacionais da Física B
Author:
Esteban Gerling
Last Updated:
8년 전
License:
Creative Commons CC BY 4.0
Abstract:
Estudo de um sistema dinâmico.
\begin
Discover why 18 million people worldwide trust Overleaf with their work.
Estudo de um sistema dinâmico.
\begin
Discover why 18 million people worldwide trust Overleaf with their work.
\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}