Exibir código-fonte para Métodos Computacionais para Estudo de Fenômenos Ondulatórios

''' Filssen Schereiber, João Roth e Lucas Oliveira '''

Este trabalho tem como objetivo realizar um estudo introdutório da equação da onda, adotando uma abordagem computacional que proporciona uma base teórica sólida e intuitiva. Inicialmente, é realizada a implementação de uma solução numérica (FTCS) para a equação da advecção, visando o estudo da propagação de uma perturbação e explorando aspectos computacionais fundamentais. Em seguida, a equação da onda é aplicada em uma dimensão, permitindo uma investigação detalhada de fenômenos ondulatórios, tais como reflexão, interferência e refração. A abordagem é estendida para duas dimensões, incluindo a análise dos efeitos de contorno nesse contexto. Por fim, o trabalho aprofunda-se no estudo do fenômeno de difração, contribuindo para uma compreensão mais enriquecedora dos aspectos complexos dos fenômenos ondulatórios em múltiplas dimensões.


= Introdução =

O que é uma onda exatamente? De acordo com 'A Student’s Guide to Waves', a literatura apresenta diversas definições. Contudo, a característica mais comum que define uma onda é sua natureza como uma perturbação, uma alteração no estado de equilíbrio que, inicialmente, permanece inalterado. Quando uma fonte externa, como um objeto vibrante ou uma força inicial, perturba a condição de equilíbrio de um meio, essa perturbação é transmitida de partícula para partícula ao longo do meio. Esse processo de transmissão de energia ocorre sem que as partículas individuais do meio se desloquem significativamente de suas posições de equilíbrio. Em outras palavras, a energia da perturbação é transferida, mas as partículas do meio não se movem em conjunto com a onda.

O estudo da física das ondas desempenha um papel crucial, permeando inúmeras áreas da ciência e da tecnologia ao proporcionar uma compreensão profunda dos fenômenos ondulatórios e suas aplicações práticas. A física das ondas desempenha um papel essencial em nossa vida cotidiana. Seu impacto é evidente na tecnologia de comunicação, onde as ondas de rádio e micro-ondas são a espinha dorsal de redes globais de telefonia móvel e internet sem fio.

Com base no exposto, propõe-se a realização de um estudo abrangente sobre a modelagem de fenômenos ondulatórios. Inicialmente, serão apresentados pontos importantes acerca da equação da onda. Em seguida, dedicaremos atenção ao estudo da propagação de uma perturbação a partir da equação de advecção. Nesse contexto, abriremos a discussão sobre a abordagem computacional, e, tangencialmente, comentaremos sobre outros métodos numéricos para resolver a advecção. Posteriormente, faremos um estudo sobre a solução numérica da equação da onda, permitindo a modelagem de fenômenos ondulatórios, como reflexão, interferência, refração e difração.


= Equação da Onda =

A função de uma onda é aquela que especifica o valor da perturbação em cada ponto e instante ao longo de seu percurso. Para iniciar nossa análise, abordaremos o exemplo mais simples, em que a propagação ocorre unicamente em uma direção, como é o caso das ondas transversais em uma corda.

== Dedução ==

O perfil da onda em uma corda em um determinado instante é equivalente à forma que a corda apresenta nesse momento, que é dada pela função <math>y(x, t)</math>. A perturbação assume a forma de uma onda progressiva, movendo-se como uma entidade coesa para a direita, mantendo sua configuração inalterada, com velocidade <math>v</math>. Dessa forma, ao acompanhar a onda em um referencial inercial diferente, em que <math>y(x,t)=y(x',t)</math>, onde <math>x'</math> representa o antigo <math>x</math> deslocado por <math>vt</math>, a relação entre os dois referenciais é estabelecida por

<center><math>
x' = x - vt
</math>.</center>

<center><math>
y' = y
</math>.</center>

de modo que, no referencial original,

<center><math>
y(x,t) = f(x - vt)
</math>.</center>
descreve uma onda progressiva, que se propaga para a direita, com velocidade <math>v</math>.

Para associar uma equação de movimento com a propagação da onda, vamos calcular a aceleração num dado ponto <math>x</math>. A velocidade e a aceleração em <math>x</math> se obtêm fixando <math>x</math> e derivando em relação ao tempo, o que corresponde a tomar derivadas parciais.

<center><math>
        V = \frac{\partial}{\partial t} y(x,t)
</math>.</center>

<center><math>
        A = \frac{\partial ^2}{\partial t^2} y(x,t)
</math>.</center>

pela regra da cadeia

<center><math>
\frac{\partial y}{\partial t} = \frac{\partial f}{\partial x'} \frac{\partial x'}{\partial t} = -v \frac{\partial f}{\partial x'}
</math>.</center>
Analogamente,

<center><math>
\frac{\partial ^2}{\partial t^2} = -v \frac{\partial}{\partial t} \frac{\partial f}{\partial x'} = -v \frac{\partial}{\partial x'} \frac{\partial f}{\partial x'} \frac{\partial x'}{\partial t}
</math>.</center>

ou seja,

<center><math>
\frac{\partial ^2 y}{\partial t^2} = v^2 \frac{\partial ^2 f }{\partial {x'}^2}
</math>.</center>

Derivando <math>y</math> em relação a <math>x</math> obtemos que:

<center><math>
        \frac{\partial y}{\partial x} = \frac{\partial f}{\partial x'}
</math>.</center>
<center><math>
        \frac{\partial ^2 y}{\partial x^2} = \frac{\partial ^2 f}{\partial {x'}^2}
</math>.</center>

Assim,

<center><math>
\frac{\partial ^2 y}{\partial t^2} - v^2 \frac{\partial^2 y }{\partial x^2} = 0  
</math>.</center>

Denominada equação da onda unidimensional, esta é uma das equações fundamentais da física.

== Solução Analítica ==

Agora, procederemos à análise do movimento de uma corda de comprimento <math>L</math> que vibra com extremidades fixas. Como essas extremidades são fixas, temos, nos pontos <math>x = 0</math> e <math>x = L</math>, os valores <math>y(0,t) = y(L,t) = 0</math>.

Além disso, supomos que a função <math>y(x,0) = f(x)</math> e <math>\frac{\partial y}{\partial t} = 0</math> para <math>t = 0</math>.

Com base nas observações realizadas, torna-se possível identificar um conjunto fundamental de soluções por meio da aplicação do método de separação de variáveis.

Suponhamos que existam duas soluções, uma dependendo exclusivamente de <math>x</math> e outra exclusivamente de <math>t</math>, de modo que possamos expressar a função <math>y(x,t)</math> como o produto dessas duas.

<center><math>
y(x,t) = X(x)T(t)
</math></center>

Assim, podemos escrever as derivadas parciais como

<center><math>
\begin{cases}
\frac{\partial^2 y}{\partial x^2} = \frac{X''}{X} \\
\frac{\partial^2 y}{\partial t^2} = \frac{T''}{T} \end{cases}
</math></center>

substituindo na equação da onda,

<center><math>
\frac{T''}{T} - v^2 \frac{X''}{X} = 0
</math></center>

Dessa forma, torna-se evidente que a razão entre as funções de <mathT</math> consiste em uma constante, assim como <math>v^2 \frac{X''}{X}</math>.

<center><math>
v^2 \frac{X''}{X} = \lambda ^2 \to X'' = \bigg(\frac{\lambda}{v}\bigg)^2 X
</math></center> '''(I)'''

<center><math>
\frac{T''}{T} = \lambda ^2  \to T''= \lambda^2 T
</math></center> '''(II)'''

Observa-se que as equações (I) e (II) são equações diferenciais ordinárias (EDOs) lineares de segunda ordem, cujos coeficientes podem depender da variável independente e, potencialmente, serem funções não constantes. Assim, é possível usarmos da teoria de Sturm-Liouville para resolver essas equações.

'''Resolvendo a Equação (I)'''

Se avaliarmos as condições de contorno para os possíveis valores de lambda, que devem ser maiores ou iguais a zero, obtemos apenas a solução trivial. Contudo, se lambda for negativo, surge uma solução para <math>X(x)</math> na forma

<center><math>
X(x) = a \cos \bigg( \frac{\lambda x}{v} \bigg) + b \sin \bigg( \frac{\lambda x}{v} \bigg)
</math></center>

Avaliando as condições de contorno, temos:

<center><math>
X(0) = a \cos (0) + b \sin (0) = a  \therefore a=0
</math></center>

<center><math>
X(L) = b \sin \bigg( \frac{\lambda L}{v} \bigg)
</math></center>

Certamente, observa-se que, para alcançarmos uma solução não trivial, é necessário que <math>\lambda = \bigg( \frac{nv\pi}{L} \bigg) </math>, onde <math>n \in \mathbb{N}</math>.

<center><math>
X_n (x) = \sum_{n=1}^{\infty} b_n \sin \bigg( \frac{n\pi x}{L} \bigg)
</math></center>

'''Resolvendo a Equação (II)'''

Analogamente, temos uma solução não trivial para lambda menor que zero.

<center><math>
T(t) = c \cos(\lambda t) + d \sin(\lambda t) \to T_n (t) = \sum_{n=1}^{\infty} c_n \cos \bigg(\frac{nv\pi t}{L} \bigg) + d_n \sin \bigg(\frac{nv\pi t}{L} \bigg)
</math></center>

Voltando para a solução geral:

<center><math>
y(x,t) = \sum_{n=1}^{\infty}  b_n \sin \bigg( \frac{n\pi x}{L} \bigg) \bigg(c_n \cos \bigg(\frac{nv\pi t}{L} \bigg) + d_n \sin \bigg(\frac{nv\pi t}{L} \bigg) \bigg)
</math></center>

Juntando as constantes

<center><math>
y(x,t) = \sum_{n=1}^{\infty} \sin \bigg( \frac{n\pi x}{L} \bigg) \bigg(A_n \cos \bigg(\frac{nv\pi t}{L} \bigg) + B_n \sin \bigg(\frac{nv\pi t}{L} \bigg) \bigg)
</math></center>

Avaliando as condições iniciais:

Para <math>t=0</math>, temos

<center><math>
y(x,0) =  \sum_{n=1}^{\infty} A_n \sin \bigg( \frac{n\pi x}{L} \bigg) = f(x)
</math></center>

É possível observar que <math>f(x)</math> representa uma série de Fourier de senos, na qual os coeficientes <math>A_n</math> são determinados por:

<center><math>
A_n = \frac{2}{L}\int_0^L f(x) \sin \bigg( \frac{n\pi x}{L} \bigg) dx
</math></center>

Avaliando <math>y'(x,0)</math>:

<center><math>
\frac{\partial }{\partial t} y(x,0) = \sum_{n=1}^{\infty} B_n \bigg( \frac{nv\pi}{L} \bigg) \sin \bigg( \frac{n\pi x}{L} \bigg) = g(x)
</math></center>

Analogamente a <math>f(x)</math>, <math>g(x)</math> é uma série de Fourier de senos. Portanto, podemos determinar <math>B_n</math> da seguinte maneira:

<center><math>
B_n = \frac{2}{nv \pi}\int_0^L g(x) \sin \bigg( \frac{n\pi x}{L} \bigg) dx
</math></center>

=== Implementação da Solução Analítica ===
[[Arquivo:eq onda analitica.png|thumb|upright=0.0|right|400px]]

Se considerarmos 10 como o número máximo de termos na série de Fourier de uma corda com comprimento <math>L=1</math> e velocidade <math>v=1</math>, é possível criar um mapa da amplitude em função das posições de cada ponto da corda e do tempo. Esse tipo de gráfico será recorrente em nossa análise em diferentes casos.

Neste caso, é possível observar que inicialmente a corda possui um máximo localizado no centro do eixo de posição, com mínimos iguais a zero nas extremidades do eixo. À medida que a função evolui no tempo, a amplitude diminui, passando por um instante onde em todos os pontos da corda a amplitude é zero, e prossegue para completar um meio período, onde apresenta um mínimo menor que zero centrado no eixo x e máximos iguais a zero nas extremidades da corda.

= A equação de advecção unidimensional =

Partimos da equação de advecção para estudarmos a propagação de uma perturbação.

A equação de onda está intimamente relacionada à chamada equação de advecção, que, em uma dimensão, assume a seguinte forma:

<center><math>
\frac{\partial  u}{\partial t} + v\frac{\partial u}{\partial x} = 0
</math></center>

A equação a seguir descreve a advecção passiva de um campo escalar <math>u(x,t)</math> transportado por um fluxo de velocidade constante <math>v</math>. Como a equação de advecção é um tanto mais simples do que a equação da onda, discutiremos aquela em primeiro lugar. A equação de advecção possui a solução formal:

<center><math>
u(x,t) = F(x-vt)
</math></center>

Essa solução descreve um pulso de forma arbitrária que é arrastado pelo fluxo a uma velocidade constante <math>v</math>, sem alterar sua forma.

== Esquema FTCS para equação da advecção == 

Para a abordagem numérica computacional do nosso estudo, utilizamos o esquema FTCS (Forward Time-Centered Space).

O método FTCS é amplamente empregado na resolução de equações diferenciais parciais. Trata-se de um método de diferenças finitas que adota uma abordagem explícita no tempo e central no espaço. Esse esquema apresenta uma acurácia de primeira ordem no tempo e de segunda ordem no espaço.

A discretização no tempo (Forward Time) pode ser feita utilizando a diferença finita para a derivada temporal:

<center><math>
\frac{\partial u}{\partial t} \approx \frac{u_i^{n+1} - u_i^n}{\Delta t}
</math></center>

onde <math>u_i^n</math> representa a amplitude da perturbação no ponto <math>i</math> do espaço e no instante <math>n</math>, e <math>\Delta t</math> é o passo de tempo.

A discretização no espaço (Central Space) pode ser feita utilizando a diferença finita para a derivada espacial:

<center><math>
\frac{\partial u}{\partial x} \approx \frac{u_{i+1}^n - u_{i-1}^n}{2\Delta x}
</math></center>

onde <math>\Delta x</math> é o espaçamento entre os pontos discretizados no espaço.

As equações acima podem ser reescritas como

<center><math>
u_i^{n+1}= u_i^n + \frac{r}{2} (u_{i+1}^n - u_{i-1}^n)
</math></center>

onde <math>r = v\frac{\Delta t}{\Delta x}</math>.

== Análise de estabilidade de von Neumann == 

A análise de estabilidade de von Neumann é uma abordagem teórica utilizada para avaliar a estabilidade numérica de esquemas de diferenças finitas, incluindo o esquema FTCS. Usando o método de von Neumann, supõe-se a seguinte solução tentativa para a EDP:

<center><math>
u(x, t) = A(t)e^{jkx}
</math></center>

onde <math>j = \sqrt{-1}</math>, <math>k</math> é o número de onda, e <math>A</math> é a amplitude da onda, uma função do tempo.

Assim, a solução no passo de tempo <math>n</math> e em <math>x = x_i = i\Delta x</math> é escrita como

<center><math>
A_ne^{jki\Delta x}
</math></center>

Substituindo esta solução tentativa na expressão da advecção discretizada temos:

<center><math>
A_{n+1}e^{jki\Delta x} = A_ne^{jki\Delta x} - \frac{u\Delta t}{2\Delta x}\left(A_ne^{jk(i+1)\Delta x} - A_ne^{jk(i-1)\Delta x}\right)
</math></center>(III)

Seja <math>\xi</math> a razão da amplitude da onda no passo de tempo <math>n + 1</math> em relação ao passo <math>n</math>. Portanto,

<center><math>
\xi = \frac{A_{n+1}}{A_n}
</math></center>

Dividir (III) por <math>A_n</math> resulta em

<center><math>
\xi e^{jki\Delta x} = e^{jki\Delta x} - \frac{u\Delta t}{2\Delta x}\left(e^{jk(i+1)\Delta x} - e^{jk(i-1)\Delta x}\right)
</math></center>

Se dividirmos esse resultado por <math>e^{jki\Delta x}</math> temos:

<center><math>
\xi = 1 - \frac{u\Delta t}{2\Delta x}\left(e^{jk\Delta x} - e^{-jk\Delta x}\right) = 1 - \frac{u\Delta t}{\Delta xj}\sin(k\Delta x)
</math></center>

Assim,

<center><math>
|\xi| = \sqrt{1 + \left(\frac{u\Delta t}{\Delta x}\sin(k\Delta x)\right)^2}
</math></center>

Isso implica que <math>|\xi| \geq 1</math> independentemente do passo de tempo <math>\Delta t</math> selecionado ou do passo espacial <math>\Delta x</math>, portanto, FTCS é incondicionalmente instável.

Apesar de ser incondicionalmente instável, o esquema funcionará bem em nossa modelagem, uma vez que o método apresenta também suas vantagens. O esquema FTCS é relativamente simples de implementar, o que facilita sua aplicação, principalmente para problemas bidimensionais ou tridimensionais, como pretendemos abordar posteriormente.

Em alguns casos específicos, o FTCS pode ser computacionalmente eficiente. Se a equação diferencial parcial for relativamente simples e o domínio de interesse for bem comportado, o FTCS pode oferecer uma solução razoável com custo computacional relativamente baixo.

É feita uma comparação entre outros métodos numéricos de diferenças finitas para solução da equação da advecção e comparado o erro relativo à solução analítica no Apêndice A.

== Implementação do esquema FTCS == 

[[Arquivo:ftcs eq advec.png|thumb|upright=0.0|right|400px]]

Vamos agora examinar a implementação do esquema FTCS na equação da advecção para visualizarmos a propagação de uma perturbação na forma de <math>1-\cos(x)</math> em um eixo <math>x</math> infinito. Em outras palavras:

* Condição inicial: <math>u(x,0) = 1-\cos(x)</math>;

* Condições de contorno periódicas.

Neste caso, podemos observar uma perturbação se propagando no sentido negativo do eixo <math>x</math> sem perda de amplitude ou mudança de forma. Isso acontece porque, com a equação da advecção, podemos analisar apenas a cinemática de uma perturbação, a qual, por sua vez, apenas relaciona os atributos como amplitude, comprimento de onda, frequência e velocidade.

Podemos observar, ainda, que a amplitude está relacionada à condição inicial em que o máximo de <matH>1-\cos(x)</math> é atingido quando <math>\cos(x) = -1</math>. Nesse contexto, considerando <math>0 < x < 2\pi</math>, percebemos que <math>\cos(x) = -1</math> ocorre exatamente no meio do eixo. Em consequência, a amplitude máxima, que é 2, é inicialmente alcançada no ponto central do eixo <math>x</math>.

= A equação da onda em uma dimensão =

A equação da onda, que em uma dimensão assume a forma

<center><math>
 \frac{\partial^2 u}{\partial t^2} = v^2 \frac{\partial^2  u}{\partial x^2}
</math></center>(15)

onde <math>u</math> geralmente é algum tipo de deslocamento ou perturbação, enquanto <math>v</math> é a velocidade da onda. A equação de onda possui a solução formal

<center><math>
 u(x,t) = F(x+vt) + G(x+vt),
</math></center>(16)

onde <math>F</math> e <math>G</math> são funções arbitrárias. A solução acima representa pulsos de onda de formato arbitrário que se propagam com velocidade <math>v</math> nas direções <math>+x</math> e <math>-x</math>, respectivamente, sem alterar a forma.

== Esquema FTCS para equação da onda unidimensional== 

Iniciamos o estudo da solução numérica da equação da onda a partir da implementação do esquema FTCS, pois, como comentado antes, este é o método de implementação mais prático para generalizarmos para mais dimensões espaciais.

As derivadas de segunda ordem podem ser substituídas por diferenças centrais. A aproximação de diferença mais amplamente utilizada para a derivada de segunda ordem é

<center><math>
\frac{\partial^2}{\partial t^2}u(x_i,t_n)\approx \frac{u_i^{n+1} - 2u_i^n + u_i^{n-1}}{\Delta t^2}
</math></center>

Uma aproximação semelhante para a derivada de segunda ordem na direção <math>x</math> é dada por

<center><math>
\frac{\partial^2}{\partial x^2}u(x_i,t_n)\approx \frac{u_{i+1}^{n} - 2u_i^n + u_{i-1}^{n}}{\Delta x^2}
</math></center>

Agora podemos substituir as derivadas e obter

<center><math>
\frac{u_i^{n+1} - 2u_i^n + u_i^{n-1}}{\Delta t^2} = v^2\frac{u_{i+1}^{n} - 2u_i^n + u_{i-1}^{n}}{\Delta x^2}.
</math></center>

Assumimos que <math>u^n_i</math> e <math>u^{n-1}_i</math> já foram calculados para <math>i=0,\ldots,N_x</math>. A única quantidade desconhecida é, portanto, <math>u^{n+1}_i</math>, que podemos resolver da seguinte forma:

<center><math>
u^{n+1}_i = -u^{n-1}_i + 2u^n_i + r^2 \left(u^{n}_{i+1}-2u^{n}_{i} + u^{n}_{i-1}\right),
</math></center>

onde introduzimos o parâmetro

<center><math>
r = v\frac{\Delta t}{\Delta x},
</math></center>

conhecido como número de Courant.

== Análise de estabilidade de von Neumann ==

Pela análise de Von Neumann, substituímos <math>u_n</math> por <math>g_n e^{i\theta}</math> na equação da onda. Após manipulações, a expressão resultante é:

<center><math>
g = 2 - g^{-1} + \frac{C}{2}(e^{i\theta} - 2 + e^{-i\theta})
</math></center>

Usando a identidade de Euler <math>e^{i\theta} + e^{-i\theta} = 2 \cos \theta</math>, chegamos a:

<center><math>
g^2 - 2[1 - \frac{C}{2}(1 - \cos \theta)]g + 1 = 0
</math></center>

Reescrevendo <math>\cos \theta = \cos^2\left(\frac{\theta}{2}\right) - \sin^2\left(\frac{\theta}{2}\right)</math>, a equação torna-se:

<center><math>
g^2 + 2(2C^2 \sin^2\left(\frac{\theta}{2}\right) - 1)g + 1 = 0
</math></center>

Expressamos isso como <math>g^2 + 2bg + 1 = 0</math>, onde <math>b = 2C^2 \sin^2\left(\frac{\theta}{2}\right) - 1</math>. As raízes são:

<center><math>
g = \begin{cases} -b \pm \sqrt{b^2 - 1}, & \text{se } |b| > 1 \\ -b \pm i\sqrt{1 - b^2}, & \text{se } |b| \leq 1 \end{cases}
</math></center>

Se <math>|b| > 1</math>, implica <math>C \leq 1</math>. Se <math>|b| \leq 1</math>, <math>C \leq 1</math> pela condição de Von Neumann, tornando o método condicionalmente estável pela condição CFL.

== Reflexão de onda em uma corda ==

Quando uma onda se propaga ao longo de uma corda e incide em uma de suas extremidades, ocorre a sua reflexão. Ao atingir um extremo fixo, que não pode oscilar, é possível observar que o pulso refletido apresenta uma orientação oposta à do pulso incidente. Podemos afirmar que o pulso reflete com inversão de fase, pois isso ocorre devido à reação do ponto fixo sobre a corda com uma força de mesmo módulo. Uma situação distinta ocorre quando a extremidade da corda é móvel, ou seja, não há troca de forças entre ela e a corda, resultando em uma reflexão sem inversão de fase.

Além disso, em um sistema fechado e sem dissipação de energia, podemos observar que as ondas refletidas, independentemente da fase, sofrem interferência construtiva e produzem picos de máximo ou mínimo, dependendo do contexto.

Podemos observar esse fenômeno ao implementar a solução da equação da onda unidimensional com condições de contorno fixas e periódicas, respectivamente, para uma corda com extremidades fixas e extremidades móveis.

=== Perturbação em uma corda com extremidades fixas ===

Neste caso, é possível observar que a corda iniciará a partir do seu ponto mais alto no centro do eixo. A perturbação se propagará igualmente para ambos os sentidos do eixo. Ao atingir as extremidades fixas, as duas perturbações, invertidas, deslocam-se em direção ao centro da corda, onde interferem de forma construtiva e produzem um mínimo de amplitude.

<center>
{|
|[[Arquivo:corda_ponta_fixa.gif|thumb|upright=0.0|center|Evolução temporal da amplitude pela posição no eixo x em uma corda com extremidades fixas|350px]]
|[[Arquivo:MA corda com pontas fixas.png|thumb|upright=0.0|center|310px]]
|}
</center>

=== Perturbação em uma corda com extremidades móveis ===

Com extremidades móveis, não se observa uma mudança de fase da perturbação inicial ao atingir os extremos da corda. No entanto, ao implementarmos condições de contorno periódicas, implicamos que, ao alcançar uma das extremidades da corda, o pulso retorna à outra extremidade, continuando a se propagar na mesma direção e com a mesma velocidade. Como a perturbação inicial gera dois pulsos se propagando em direções opostas a partir do centro do eixo, com velocidades iguais, percebemos a interação construtiva dos pulsos que alcançam os limites do eixo simultaneamente. Isso faz com que retornem a se propagar em direção ao centro, gerando máximos de amplitude em cada interação.

<center>
{|
|[[Arquivo:corda_ponta_solta.gif|thumb|upright=0.0|center|Evolução temporal da amplitude pela posição no eixo x em uma corda com extremidades interligadas por condições de contorno periódicas.|350px]]
|[[Arquivo:MA corda com pontas soltas.png|thumb|upright=0.0|center|310px]]
|}
</center>

=== Alteração no meio de propagação ===

A refração é um fenômeno que ocorre quando uma onda passa de um meio para outro com uma velocidade de propagação diferente. Esse fenômeno é mais comumente associado a ondas de luz, mas também pode ocorrer com ondas sonoras, ondas de rádio, ondas sísmicas e outras formas de ondas.

Se uma onda se propaga em uma corda e encontra uma mudança abrupta nas propriedades da corda, como a densidade ou a tensão, isso pode resultar em uma mudança na velocidade da onda. O grau de refração dependerá das diferenças nas propriedades do meio antes e depois da mudança.

<center>
{|
|[[Arquivo:muda_meio_de_propag.gif|thumb|upright=0.0|center|Evolução temporal da amplitude pela posição no eixo x em uma corda com extremidades fixas e diferentes densidades.|350px]]
|[[Arquivo:MA corda muda meio.png|thumb|upright=0.0|center|340px]]
|}
</center>


Em um modelo que considera a reflexão e a refração na junção das duas ondas [Elmore e Head], é demonstrado que a razão entre as amplitudes da onda refratada (<math>A_2</math>) e a amplitude da onda incidente (<math>A_1</math>) é dada por:

<center><math>
\frac{A_2}{A_1} = \frac{2\sqrt{\mu_1}}{\sqrt{\mu_1} +\sqrt{\mu_2}}
</math></center>

Usando a relação da velocidade de propagação com as variáveis de tensão e densidade podemos chegar em expressões para as amplitudes do pulso que irá voltar (<math>A_v</math>) e do pulso que irá passar para a outra corda (<math>A_p</math>).

<center><math>
A_p = \frac{2v_2}{v_1 + v_2}A_i 
</math></center>


<center><math>
A_v = \frac{|v_1 - v_2|}{v_1 + v_2}A_i
</math></center>

'''(IV)'''

onde <math>A_i</math> é a amplitude do pulso incidente.

No nosso caso, a perturbação inicial tem amplitude de 1, porém propaga-se isotopicamente nos dois sentidos do eixo x. O pulso que se desloca no sentido negativo do eixo interage imediatamente com o ponto fixo da corda e torna a propagar-se no sentido positivo, mas agora com uma fase. Assim, temos dois pulsos, um positivo e um negativo, ambos com amplitude igual a <math>0.5</math>, propagando-se em direção a outra corda de densidade diferente.

Com as velocidades de propagação nos dois meios, utilizando as fórmulas em (IV), podemos determinar que a amplitude do pulso que irá voltar é <math>0,167</math>, enquanto a do pulso que será refratado é <math>0,333</math>.

[[Arquivo:muda_meio_de_propag_analise.gif|thumb|upright=0.0|center|450px]]

=== Oscilador em uma corda ===

Agora, em vez de variarmos as condições de contorno para um pulso senoidal em uma corda, estamos avaliando um oscilador senoidal acoplado a corda com parâmetros físicos bem definidos: amplitude (A), frequência (f) e comprimento de onda ( λ ). Assim, ao invés de termos uma condição inicial que aplicada em um pequena porção da corda se propaga, iremos avaliar um oscilador da forma <math>A\sin(kx)</math> que define a amplitude inicial para cada ponto da corda, onde <math>k</math> é o número de onda definido como duas vezes <math>\pi</math> sobre o comprimento de onda. Vale notar ainda que λ é uma função da velocidade de propagação (propriedade da corda) e da frequência (propriedade do oscilado).

Seja <math>A = 1</math>, <math>v = 1</math> e <math>f =  0.13</math>, que resulta em um comprimento de onda igual a 7.7. Implementando um oscilador com esses parâmetros para uma corda com extremidades fixas, obtemos

[[Arquivo:oscilador com pontas fixas.png|thumb|upright=0.0|center|450px]]

Observamos que, ao ajustar arbitrariamente os parâmetros do oscilador, obtemos regiões mais escuras e mais claras, pouco organizadas. Nessas regiões, ocorrem interferências destrutivas e construtivas, respectivamente. As interferências ocorrem de forma pouco definida, uma vez que não há, nesse caso, um padrão bem definido para as interações que a onda sofrerá. No entanto, veremos que isso pode ser modificado.

=== Ondas estacionárias ===

A relação entre o comprimento de onda da onda incidente ( λ) e a distância disponível para a onda se refletir (<math>2L</math>, onde <math>L</math> é o comprimento da região confinada, da corda) deve permitir a formação de nós e ventres estacionários. Para ondas estacionárias em uma corda com extremidades fixas, os comprimentos de onda permitidos são dados por <math>\lambda_n = \frac{2L}{n}</math>, onde <math>n</math> é um número inteiro positivo.

Essa fórmula mostra que à medida que <math>n</math> aumenta, o comprimento de onda diminui. Cada valor de <math>n</math> corresponde a uma harmônica específica, que representa o número de meios comprimentos de onda que cabem na corda. A primeira harmônica (<math>n = 1</math>) tem um nó em cada extremidade da corda, a segunda harmônica (<math>n = 2</math>) tem dois nós, e assim por diante.

Dado que <math>\lambda = \frac{v}{f} </math> e <math>\lambda_n = \frac{2L}{n}</math>, é possível expressar um conjunto discreto de frequências para o nosso oscilador, as quais resultam em ondas estacionárias.

<center><math>
f_n = \frac{nv}{2L}
</math></center>

Dessa forma, podemos aplicar valores de <math>f</math> condicionados a gerar as condições necessárias para a formação de ondas estacionárias.

[[Arquivo:ondas_estacionarias.png|thumb|upright=0.0|center|650px]]

= A equação da onda em duas dimensões =

Agora, procederemos à avaliação da modelagem de ondas em duas dimensões. Nesse contexto, além dos fenômenos já observados, é possível perceber também a difração.

== Esquema FTCS para equação da onda bidimensional== 

Seja a equação da onda em duas dimensões dada por:

<center><math>
\frac{\partial^2 u}{\partial t^2} = v^2 \left( \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} \right)
</math></center>

A discretização desta equação utilizando o esquema FTCS envolve as seguintes aproximações:

<center><math>
\frac{\partial^2 u}{\partial t^2} \approx \frac{u^{n+1}_{i,j} - 2u^{n}_{i,j} + u^{n-1}_{i,j}}{\Delta t^2}
</math></center>


<center><math>
\frac{\partial^2 u}{\partial x^2} \approx \frac{u^{n}_{i+1,j} - 2u^{n}_{i,j} + u^{n}_{i-1,j}}{\Delta x^2}
</math></center>


<center><math>
\frac{\partial^2 u}{\partial y^2} \approx \frac{u^{n}_{i,j+1} - 2u^{n}_{i,j} + u^{n}_{i,j-1}}{\Delta y^2}
</math></center>

Assim, a equação da onda discretizada usando FTCS é dada por:

<center><math>
u^{n+1}_{i,j} = 2u^{n}_{i,j} - u^{n-1}_{i,j} + v^2 \left( \frac{\Delta t^2}{\Delta x^2} \cdot (u^{n}_{i+1,j} - 2u^{n}_{i,j} + u^{n}_{i-1,j}) + \frac{\Delta t^2}{\Delta y^2} \cdot (u^{n}_{i,j+1} - 2u^{n}_{i,j} + u^{n}_{i,j-1}) \right)
</math></center>

== Esquema FTCS para equação da onda bidimensional== 

Seja a equação da onda em duas dimensões dada por:

<center><math>
\frac{\partial^2 u}{\partial t^2} = v^2 \left( \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} \right)
</math></center>

A discretização desta equação utilizando o esquema FTCS envolve as seguintes aproximações:

<center><math>
\frac{\partial^2 u}{\partial t^2} \approx \frac{u^{n+1}_{i,j} - 2u^{n}_{i,j} + u^{n-1}_{i,j}}{\Delta t^2}
</math></center>


<center><math>
\frac{\partial^2 u}{\partial x^2} \approx \frac{u^{n}_{i+1,j} - 2u^{n}_{i,j} + u^{n}_{i-1,j}}{\Delta x^2}
</math></center>


<center><math>
\frac{\partial^2 u}{\partial y^2} \approx \frac{u^{n}_{i,j+1} - 2u^{n}_{i,j} + u^{n}_{i,j-1}}{\Delta y^2}
</math></center>

Assim, a equação da onda discretizada usando FTCS é dada por:

<center><math>
u^{n+1}_{i,j} = 2u^{n}_{i,j} - u^{n-1}_{i,j} + v^2 \left[ \frac{\Delta t^2}{\Delta x^2} \left(u^{n}_{i+1,j} - 2u^{n}_{i,j} + u^{n}_{i-1,j}\right) + \frac{\Delta t^2}{\Delta y^2} \left(u^{n}_{i,j+1} - 2u^{n}_{i,j} + u^{n}_{i,j-1}\right) \right]
</math></center>

Para <math>\Delta x = \Delta y</math>, segue que

<center><math>
u^{n+1}_{i,j} = \left(2-4a\right)u^{n}_{i,j} - u^{n-1}_{i,j} + a\left(u^n_{i-1,j}+u^n_{i+1,j}+u^n_{i,j-1}+u^n_{i,j+1}\right), \qquad a\equiv \left(v\frac{\Delta t}{\Delta x}\right)^2
</math></center>

ou, ainda,
<center><math>
u^{n}_{i,j} = \left(2-4a\right)u^{n-1}_{i,j} - u^{n-2}_{i,j} + a\left(u^{n-1}_{i-1,j}+u^{n-1}_{i+1,j}+u^{n-1}_{i,j-1}+u^{n-1}_{i,j+1}\right), \qquad a\equiv \left(v\frac{\Delta t}{\Delta x}\right)^2
</math></center>

== Reflexão de ondas em duas dimensões ==

Analogamente ao que analisamos para uma perturbação em uma corda, podemos agora avaliar como uma perturbação evolui ao se propagar em um plano em direção às suas extremidades, e como interage com os contornos, que podem ser fixos ou periódicos. Neste caso, temos uma perturbação senoidal centrada em <math>x=0,5</math> e <math>y=0,5</math>, propagando-se isotropicamente em direção aos limites <math>0</math> e <math>1</math> dos dois eixos.

=== Contorno fixo ===

=== Contorno periódico ===

= Apêndice A: Métodos numéricos para solução da equação de advecção =

[[Método Lax-Friedrich]]

[[Método Lax-Wendroff]]

[[Upwind Differencing]]

[[Método Leapfrog]]

[[Método Lax-Wendroff de dois passos]]

[[Comparação dos métodos numéricos]]