Corda Vibrante

Grupo : Gustavo H. Guesser, Joshua L. Kipper, Marcos Pasa.

A equação da onda

No estudo das oscilações, é comum entre os físicos o emprego de modelos simples como representações "prototípicas" de certos padrões básicos observados na natureza. A eficácia desses modelos simples decorre de duas características fundamentais: a capacidade de serem compreendidos em detalhes minuciosos e a habilidade de reproduzir comportamentos semelhantes aos de situações reais e mais complexas, auxiliando na compreensão destas, pelo menos em termos qualitativos.

Quando se trata de comportamentos oscilatórios, o modelo mais comumente utilizado é o do sistema massa-mola. Por outro lado, para descrever comportamentos ondulatórios unidimensionais, o modelo simples mais difundido é o da corda vibrante.

Dedução

Consideremos uma corda esticada, como a corda de um violão, por exemplo. Suponhamos que esta corda tenha um comprimento $L$ e que suas extremidades estejam fixas nos pontos $x=0$ e $x=L$. Além disso, vamos supor que a corda tenha uma densidade linear uniforme representada por $\mu ={\frac {\Delta {m}}{\Delta {x}}}$, e que esteja sob uma tensão constante $T$.

Assumindo que a corda realize vibrações transversais apenas na direção $y$ (embora possa ter vibrações transversais na direção $z$, as quais vamos ignorar aqui), podemos representar a configuração da corda em qualquer instante de tempo no plano $xy$ por uma função $y(x,t)$.

Vamos considerar mais algumas suposições:

I) A tensão $T$ que estica a corda é tão alta que podemos negligenciar a força gravitacional sobre ela;

II) A corda é perfeitamente elástica, ou seja, não oferece resistência a dobras;

III) Os deslocamentos da corda, que ocorrem apenas na direção $y$, são de pequena amplitude.

Com essas suposições em mente, estamos prontos para abordar o problema da corda vibrante. Em qualquer instante de tempo, um segmento arbitrário da corda estará na posição geral indicada pela figura abaixo:

Posição genérica de um pedaço qualquer da corda

A massa do pequeno segmento de corda de comprimento $\Delta {x}$ é

$\Delta {m}=\mu \Delta {x}.$

As componentes horizontal e vertical da força resultante atuando sobre esse segmento são:

$F_{x}=T\cos {(\theta +\Delta {\theta })}-T\cos {(\theta )}$

e

$F_{y}=T\sin {(\theta +\Delta {\theta })}-T\sin {(\theta )}.$

Estamos assumindo que a corda não se move na direção $x$ (apenas na direção $y$). Isso implica que a força resultante na direção $x$ é nula ($F_{x}=0$). Substituindo essa condição na equação, obtemos:

$\cos {(\theta +\Delta {\theta })}=\cos {(\theta )}.$

E a força resultante na direção $y$, $F_{y}$, é dada de acordo com a segunda lei de Newton:

$F_{y}=(\mu \Delta {x})a_{y}=(\mu \Delta {x}){\frac {\partial ^{2}{y}}{\partial {t^{2}}}},$

então:

$T\sin {(\theta +\Delta {\theta })}-T\sin {(\theta )}=(\mu \Delta {x}){\frac {\partial ^{2}{y}}{\partial {t^{2}}}},$

ou

$\sin {(\theta +\Delta {\theta })}-\sin {(\theta )}=\Delta {x}{\frac {\mu }{T}}{\frac {\partial ^{2}{y}}{\partial {t^{2}}}}.$

Vamos prosseguir dividindo os dois lados da equação acima pelo termo $\cos {(\theta )}$, sabendo que $\cos {(\theta )}=\cos {(\theta +\Delta {\theta })}$:

${\frac {\sin {(\theta +\Delta {\theta })}}{\cos {(\theta +\Delta {\theta })}}}-{\frac {\sin {(\theta )}}{\cos {(\theta )}}}={\frac {\Delta {x}}{\cos {(\theta )}}}{\frac {\mu }{T}}{\frac {\partial ^{2}{y}}{\partial {t^{2}}}}\to \tan {(\theta +\Delta {\theta })}-\tan {(\theta )}={\frac {\Delta {x}}{\cos {(\theta )}}}{\frac {\mu }{T}}{\frac {\partial ^{2}{y}}{\partial {t^{2}}}}.$

Considerando que o coeficiente angular da reta tangente a uma função em um dado ponto do seu domínio é igual à derivada da função nesse ponto, temos:

${\frac {\partial {y}}{\partial {x}}}(x+\Delta {x},t)-{\frac {\partial {y}}{\partial {x}}}(x,t)={\frac {\Delta {x}}{\cos {(\theta )}}}{\frac {\mu }{T}}{\frac {\partial ^{2}{y}}{\partial {t^{2}}}}.$

Dividindo ambos os lados por $\Delta {x}$, teremos do lado esquerdo:

${\frac {{\frac {\partial {y}}{\partial {x}}}(x+\Delta {x},t)-{\frac {\partial {y}}{\partial {x}}}(x,t)}{\Delta {x}}}$

e no limite que $\Delta {x}\to 0$, a expressão se torna uma derivada parcial segunda ${\frac {\partial ^{2}{y}}{\partial {x^{2}}}}$, nos permitindo reescrever a equação como:

${\frac {\partial ^{2}{y}}{\partial {x^{2}}}}(x,t)={\frac {1}{\cos {(\theta )}}}{\frac {\mu }{T}}{\frac {\partial ^{2}{y}}{\partial {t^{2}}}}(x,t).$

E, como última manipulação, vamos supor agora que os deslocamentos da corda são pequenos, o que implica que os ângulos associados a esses deslocamentos também são pequenos: $\theta <<1$. Com essa condição, $\cos {(\theta )}\approx 1$ e chegamos em:

${\frac {\partial ^{2}{y}}{\partial {x^{2}}}}(x,t)={\frac {\mu }{T}}{\frac {\partial ^{2}{y}}{\partial {t^{2}}}}(x,t)$

onde o termo ${\frac {\mu }{T}}$ tem dimensão de ${\frac {1}{c^{2}}}$ e $c$ é a velocidade de propagação de ondas na corda esticada:

${\frac {\partial ^{2}{y}}{\partial {x^{2}}}}(x,t)={\frac {1}{c^{2}}}{\frac {\partial ^{2}{y}}{\partial {t^{2}}}}(x,t)\Leftrightarrow {\frac {\partial ^{2}{y}}{\partial {t^{2}}}}(x,t)=c^{2}{\frac {\partial ^{2}{y}}{\partial {x^{2}}}}(x,t)\Leftrightarrow y_{tt}=c^{2}y_{xx}.$

Chegando a conlusão de que essa é a expressão para a corda vibrante e que a velocidade de propagação aumenta com a tensão na corda e diminui com a sua inércia (massa por unidade de comprimento).

Solução Analítica

Começaremos tomando as condições de contorno da corda:

$y(0,t)=0$

e

$y(L,t)=0$,

indicando que as extremidades da corda permanecem fixas.

Além disso, tomaremos também as seguintes condições iniciais para a posição e velocidade da corda:

$y(x,0)=f(x)$, $0\leq x\leq L$

e

${\frac {\partial {y}}{\partial {t}}}(x,0)=y_{t}=g(x)$, $0\leq x\leq L$,

onde $f(x)$ e $g(x)$ são funções tais que

$f(0)=f(L)=0$

e

$g(0)=g(L)=0.$

Iniciaremos o estudo das vibrações em uma corda elástica admitindo que a velocidade inicial da corda é nula, ou seja:

$g(x)=0$, $\forall 0\leq x\leq L.$

Em outras palavras, estamos tratando do problema:

${\begin{aligned}{\begin{cases}y_{tt}=c^{2}y_{xx},\\y(0,t)=0,\\y(L,t)=0,\\y(x,0)=f(x),&0\leq x\leq L,\\y_{t}=g(x),&0\leq x\leq L.\\\end{cases}}\end{aligned}}$

Em que $f(0)=f(L)=0$ descreve a configuração da corda no início.

Usaremos o método da separação de variáveis para resolver o problema admitindo que $y$ pode ser escrita como

$y(x,t)=X(x)T(t)$

na qual $X$ depende apenas de $x$ e $T$ depende apenas de $t$.

Aplicando na equação diferencial parcial, temos:

${\frac {X''}{X}}={\frac {1}{c^{2}}}{\frac {T''}{T}}=-\lambda ,$

onde $\lambda$ é uma constante de separação.

Logo, obtemos

$X''+\lambda X=0$

e

$T''+c^{2}\lambda T=0$

Aplicando as condições de contorno, encontramos o problema

$X''+\lambda X=0$, $X(0)=X(L)=0,$

cuja solução é dada por

$X_{n}(x)=\sin {\left({\frac {n\pi x}{L}}\right)},\ {\lambda }_{n}=\left({\frac {n\pi }{L}}\right)^{2},\ n=1,2,3,...$

Com as constantes de separação acima, podemos obter a EDO

$T''+{\omega }^{2}T=0,\ \omega ={\frac {n\pi c}{L}},$

cujas soluções são

$T(t)=k_{1}\cos {(\omega t)}+k_{2}\sin {(\omega t)}.$

Pelo fato da velocidade inicial da corda ser nula, deduzimos que

$y_{t}(x,0)=X(x)T'(0)=0,\ \forall 0\leq x\leq L,\ \Rightarrow T'(0)=0.$

Como

$T'(t)=-\omega k_{1}\sin {(\omega t)}+\omega k_{2}\cos {(\omega t)},$

temos

$T'(0)=0\ \Rightarrow k_{2}=0.$

Assim, as soluções fundamentais da equação da onda, envolvendo as condições de contorno e a segunda condição inicial são dadas por:

$y_{n}(x,t)=\sin {\left({\frac {n\pi x}{L}}\right)}\cos {(\omega t)}=\sin {\left({\frac {n\pi x}{L}}\right)}\cos {\left({\frac {n\pi ct}{L}}\right)},$

para $n=1,2,3,...$ .

Por fim, a superposição das soluções fundamentais nos fornece

$y(x,t)=\sum _{n=1}^{\infty }c_{n}\sin {\left({\frac {n\pi x}{L}}\right)}\cos {\left({\frac {n\pi ct}{L}}\right)}$

e a condição inicial $y(x,0)=f(x)$ fornece

$y(x,0)=\sum _{n=1}^{\infty }c_{n}\sin {\left({\frac {n\pi x}{L}}\right)}=f(x).$

Então, admitindo que $f(x)$ é uma função ímpar com período $T=2L$, concluímos que os coeficientes satisfazem

$c_{n}={\frac {2}{L}}\int _{0}^{L}f(x)\sin {\left({\frac {n\pi x}{L}}\right)}dx,\ n=1,2,3,...$

Juntando as informações, temos que a solução do problema

${\begin{aligned}{\begin{cases}y_{tt}=c^{2}y_{xx},\\y(0,t)=0,\\y(L,t)=0,\\y(x,0)=f(x),&0\leq x\leq L,\\y_{t}=g(x),&0\leq x\leq L.\\\end{cases}}\end{aligned}}$

é dada por

$y(x,t)=\sum _{n=1}^{\infty }c_{n}\sin {\left({\frac {n\pi x}{L}}\right)}\cos {\left({\frac {n\pi ct}{L}}\right)},$

em que

$c_{n}={\frac {2}{L}}\int _{0}^{L}f(x)\sin {\left({\frac {n\pi x}{L}}\right)}dx,\ n=1,2,3,...$

Adaptação da equação da onda para uma corda real

Até o presente momento, tudo que foi apresentado diz respeito a uma corda ideal. Como é de se esperar, cordas reais terão algumas perdas por atrito (amortecimento) e também uma rigidez. Para termos ideia de como se dá esse comportamento real, adicionaremos esses novos dois termos na composição da nossa equação da onda:

${\frac {\partial ^{2}y}{\partial t^{2}}}=c^{2}\left({\frac {\partial ^{2}y}{\partial x^{2}}}-\epsilon L^{2}{\frac {\partial ^{4}y}{\partial x^{4}}}\right)-2b{\frac {\partial y}{\partial t}}\Leftrightarrow y_{tt}=c^{2}(y_{xx}-\epsilon L^{2}y_{xxxx})-2by_{t},$

onde o termo da derivada de 4ª ordem vem da rigidez da corda, o qual é controlado por $\epsilon$, $L$ é o comprimento da corda e $-2by_{t}$ é o amortecimento sofrido pela corda ao longo do tempo de propagação da onda.

Utilizaremos essa equação futuramente para compreensão visual do comportamento de uma corda mais próxima da realidade.

Método FTCS

O método FTCS (Forward Time Central Space) é uma abordagem progressiva no tempo e centrada no espaço. Em outras palavras, ao lidarmos com uma função de duas variáveis, buscamos uma solução futura em termos do tempo, centrada em uma vizinhança espacial. Para a resolução numérica da equação diferencial unidimensional da onda, optamos por trabalhar exclusivamente com o método explícito.

Corda Ideal

Método Explícito

Para aplicar o método, é necessário inicialmente discretizar tanto as variáveis espaciais quanto as variáveis temporais. Seja a função $y(t,x):\mathbb {R} ^{2}\mapsto \mathbb {R}$, sejam os intervalos $\left[0,X\right],\left[0,T\right]\subset \mathbb {R}$ e $N,M\in \mathbb {N}$, discretizamos os intervalos em espaçamentos iguais $\Delta x=X/N$, $\Delta t=T/M$ de tal maneira que obtemos as sequências crescentes monótonas $\left\{x_{j}\right\}_{j=0}^{N}$ e $\left\{t_{i}\right\}_{i=0}^{M}$ . Assim obtemos uma grade de malhas de pontos $(t_{n},x_{j})$ onde a função $y(t,x)$ tem seus valores aproximados por $y(t_{n},x_{j})$ denotado $y_{j}^{n}$. Uma abordagem comum envolve a discretização da primeira derivada no tempo utilizando uma diferença finita progressiva baseada na expansão em séries de Taylor até a primeira ordem :

${\frac {\partial y}{\partial t}}\approx {\frac {y(t+\Delta t,x)-y(t,x)}{\Delta t}}={\frac {y_{j}^{n+1}-y_{j}^{n}}{\Delta t}}$

Discretizando a derivada temporal segunda utilizando um diferença finita regressiva :

${\frac {\partial ^{2}y}{\partial t^{2}}}\approx {\frac {\partial }{\partial t}}\left({\frac {y_{j}^{n+1}-y_{j}^{n}}{\Delta t}}\right)\approx {\frac {1}{\Delta t}}\left[{\frac {y\left(t+\Delta t,x\right)-y\left(\left(t+\Delta t\right)-\Delta t,x\right)}{\Delta t}}-{\frac {y\left(t,x\right)-y\left(t-\Delta t,x\right)}{\Delta t}}\right]={\frac {y\left(t+\Delta t,x\right)+y\left(t-\Delta t,x\right)-2y\left(t,x\right)}{(\Delta t)^{2}}}={\frac {y_{j}^{n+1}+y_{j}^{n-1}-2y_{j}^{n}}{(\Delta t)^{2}}}$

note que o processo retorna uma derivada de diferenças finitas centrada, além disso obtemos um aproximação de segunda ordem. O mecanismo de discretização das derivadas espaciais é completamente análogo, a derivada segunda centrada :

${\frac {\partial ^{2}y}{\partial x^{2}}}\approx {\frac {y_{j+1}^{n}+y_{j-1}^{n}-2y_{j}^{n}}{(\Delta x)^{2}}}$

Aplicando as expressões das diferenciais na equação da onda :

${\frac {\partial ^{2}y}{\partial t^{2}}}=c^{2}{\frac {\partial ^{2}y}{\partial x^{2}}}\mapsto {\frac {y_{j}^{n+1}+y_{j}^{n-1}-2y_{j}^{n}}{(\Delta t)^{2}}}=c^{2}{\frac {y_{j+1}^{n}+y_{j-1}^{n}-2y_{j}^{n}}{(\Delta x)^{2}}}$

defina $a:=c\Delta t/\Delta x$ e explicite a variável no estado temporal futuro, obtemos :

$y_{j}^{n+1}=2\left(1-a^{2}\right)y_{j}^{n}-y_{j}^{n-1}+a^{2}\left(y_{j+1}^{n}+y_{j-1}^{n}\right)$

uma equação explicita no tempo $n+1$ dependendo somente de $n$ e $n-1$.

Estabilidade

Para analisar a estabilidade da solução numérica para a equação da onda, utilizaremos o método de Von Neumann. O procedimento consiste em utilizar a transformada de Fourier para determinar a estabilidade do esquema, a solução problema tem a forma :

$y_{j}^{n+1}(\xi )=A(\xi \Delta x,\Delta x,\Delta t)y_{j}^{n}(\xi )$

onde $\xi \in \left[-\pi /\Delta x,\pi /\Delta x\right]$ e $A$ é o fator de amplificação, observe que podemos estudar a estabilidade da aproximação somente analisando o fator de amplificação.

• Teorema : Se o fator de amplificação não depender de $\Delta x$ nem $\Delta t$, então a aproximação é estável se $\left|A(\xi \Delta x)\right|\leq 1$.

Vejamos a estabilidade da aproximação de diferenças finitas para a equação da onda, utilizando as componentes de Fourier $y_{j}^{n}\mapsto B^{n}e^{ikj\Delta x}$, obtemos :

$B^{n+1}e^{ikj\Delta x}=2\left(1-a^{2}\right)B^{n}e^{ikj\Delta x}-B^{n-1}e^{ikj\Delta x}+a^{2}\left(B^{n}e^{ik(j+1)\Delta x}+B^{n}e^{ik(j-1)\Delta x}\right)$

simplificando os termos em comum :

$B^{n+1}+B^{n-1}=B^{n}\left[2\left(1-a^{2}\right)+a^{2}\left(e^{ik\Delta x}+e^{-ik\Delta x}\right)\right]$

note que $e^{ik\Delta x}-e^{-ik\Delta x}=2cos(k\Delta x)$,

${\frac {B^{n+1}+B^{n-1}}{B^{n}}}=2\left(1-a^{2}\right)+2a^{2}cos(k\Delta x)$

utilizando a identidade trigonométrica $cos(\theta )=cos^{2}(\theta /2)-sen^{2}(\theta /2)$ e a identidade trigonométrica fundamental $1=cos^{2}(\theta )+sen^{2}(\theta )$, obtemos :

$A+{\frac {1}{A}}=-2\left[2a^{2}sen^{2}\left({\frac {k\Delta x}{2}}\right)-1\right]$

onde $A=B^{n+1}/B^{n}$ é o fator de amplificação. Defina $d:=2a^{2}sen^{2}\left(k\Delta x/2\right)-1$. Portanto obtemos a equação quadrática com soluções :

$A^{2}+2dA+1=0\Rightarrow A_{1,2}=-d\pm {\sqrt {d^{2}-1}}$

se $\left|d\right|>1$, então :

$-1>d>1\Rightarrow -1>2a^{2}sen^{2}\left(k\Delta x/2\right)-1>1\Rightarrow 0>2a^{2}sen^{2}\left(k\Delta x/2\right)>2\Rightarrow 0>\left(asen\left(k\Delta x/2\right)\right)^{2}>1$

o que é um aburdo, pois $asen\left(k\Delta x/2\right)\in \mathbb {R}$. Se $\left|d\right|\leq 1$ , então :

$-1\leq d\leq 1\Rightarrow -1\leq 2a^{2}sen^{2}\left(k\Delta x/2\right)-1\leq 1\Rightarrow 0\leq 2a^{2}sen^{2}\left(k\Delta x/2\right)\leq 2\Rightarrow 0\leq \left(asen\left(k\Delta x/2\right)\right)^{2}\leq 1$

note que $sen\left(k\Delta x/2\right)\leq 1$, logo $a\leq 1$. Portanto para o método ser condicionalmente estável :

${\frac {c\Delta t}{\Delta x}}\leq 1$

Embora $a\leq 1$ garante a estabilidade do processo, $a=1$ é a melhor escolha possível para a estabilidade do método, pelos seguintes motivos :

• Quando $a=1$, os termos de ordem superior que são descartados na discretização da equação da onda são amplamente cancelados. Essa compensação implica que o método numérico se torna mais estável e menos suscetível a erros numéricos significativos, ou seja, não ocorrerão variações abruptas devido aos termos de ordem maiores;
• Quando $a=1$, a perturbação na corda se propaga exatamente um passo espacial a cada passo de tempo. Esse comportamento assegura que a velocidade da perturbação seja adequadamente representada pelo algoritmo numérico. Isso ocorre devido a relação entre a velocidade física da corda $c$ e a velocidade da perturbação $v=\Delta x/\Delta t$. A animação abaixo contém soluções da equação da onda ideal, obtidas utilizando o método FTCS com diferentes a's. Nela fica bem claro que a velocidade de propagação das ondas na solução para $a=0,5$ é metade da solução para $a=1$

Soluções obtidas para a equação da onda ideal com diferentes a's.

${\frac {c\Delta t}{\Delta x}}={\frac {c}{v}}$

Observe que além de obtermos o intervalo de estabilidade $a\in \left[0,1\right]$, também obtemos uma restrição na discretização das malhas temporais e espaciais, limitando assim $N,M\in \mathbb {N}$ a um subconjunto que satisfaça a condição de estabilidade.

Corda Não Ideal

Método Explícito

Optaremos novamente por realizar a análise utilizando um método explícito. Em geral, o desenvolvimento teórico para a discretização da equação diferencial da onda com amortecimento e rigidez é análogo ao caso ideal, portanto, não serão apresentados com detalhes os cálculos necessários para a implementação do método. A derivada quarta centrada é aproximada utilizando diferenças finitas :

${\frac {\partial ^{4}y}{\partial x^{4}}}\approx {\frac {y_{j+2}^{n}+y_{j+1}^{n}+y_{j}^{n}+y_{j-1}^{n}+y_{j-2}^{n}}{\left(\Delta x\right)^{4}}}$

transformando a EDP não homogênea para um espaço discretizado usando a aproximação de diferenças finitas centrada :

${\frac {\partial ^{2}y}{\partial t^{2}}}=c^{2}\left({\frac {\partial ^{2}y}{\partial x^{2}}}-\epsilon L^{2}{\frac {\partial ^{4}y}{\partial x^{4}}}\right)-2b{\frac {\partial y}{\partial t}}\mapsto y_{j}^{n+1}=\left[\left(2-2a^{2}-6\epsilon a^{2}E^{2}-F\right)y_{j}^{n}-y_{j}^{n-1}+a^{2}\left(1+4\epsilon E^{2}\right)\left(y_{j+1}^{n}+y_{j-1}^{n}\right)-\epsilon a^{2}E^{2}\left(y_{j+2}^{n}+y_{j-2}^{n}\right)\right]\left({\frac {1}{1-F}}\right)$

onde $E:=L/\Delta x$ e $F:=b/\Delta t$. Novamente obtemos um equação explicita no tempo $n+1$ que depende somente de $n$ e $n-1$;

Estabilidade

A abordagem que empregamos para demonstrar a estabilidade da equação da onda não homogênea é completamente análoga ao caso homogêneo. Em outras palavras, utilizamos o método de Von Neumann em conjunto com o teorema da estabilidade. Ao aplicar novamente as componentes de Fourier, representadas por $y_{j}^{n}\mapsto B^{n}e^{ikj\Delta x}$, obtemos uma equação quadrática que impõe a seguinte restrição :

$a={\frac {c\Delta t}{\Delta x}}\leq {\frac {1}{4}}$

logo, o método é condicionalmente estável no intervalo $a\in \left[0,1/4\right]$.

Convergência, Consistência e Estabilidade

Até este ponto, ainda não discutimos o quão bem o esquema de diferenças finitas descreve a solução da equação diferencial. Vamos abordar isso agora. A propriedade mais essencial que um esquema deve possuir é que, à medida que os espaçamentos na malha discreta diminuem, a aproximação para a solução deve se aproximar da solução real do sistema. Portanto, um esquema é considerado uma boa aproximação se

$\left(\Delta t,\Delta x\right)\rightarrow 0\Rightarrow \left|y\left(t_{n},x_{j}\right)-y_{j}^{n}\right|\rightarrow 0$

onde $y\left(t_{n},x_{j}\right)$ é a solução exata e $y_{j}^{n}$ a solução aproximada. Se tais condições são satisfeitas, dizemos que o esquema de equações diferenciais finitas é convergente. Entretanto, geralmente, provar a convergência não é uma tarefa trivial. No entanto, podemos demonstrá-la de maneira indireta, utilizando estabilidade e consistência. Se um esquema é estável e consistente, então a convergência do esquema é garantida. Dada uma equação diferencial $Dy=f$, dizemos que o esquema de diferenças finitas $D_{\Delta t\Delta x}y=f$ é consistente com a equação diferencial se, para qualquer função suficientemente diferenciável $\phi (t,x)$,

$\left(\Delta t,\Delta x\right)\rightarrow 0\Rightarrow D\phi -D_{\Delta t\Delta x}\phi \rightarrow 0$

Mostremos a consistência da aproximação de diferenças finitas para a equação da onda homogênea. Seja a EDP da onda

$Dy\left(t,x\right)={\frac {\partial ^{2}y}{\partial t^{2}}}-c^{2}{\frac {\partial ^{2}y}{\partial x^{2}}}$

e seja a equação de diferenças finitas

$D_{\left(\Delta t\Delta x\right)}y\left(t,x\right)={\frac {y_{j}^{n+1}+y_{j}^{n-1}-2y_{j}^{n}}{(\Delta t)^{2}}}-c^{2}{\frac {y_{j+1}^{n}+y_{j-1}^{n}-2y_{j}^{n}}{(\Delta x)^{2}}}$

pela aproximação que demonstramos anteriormente, temos :

${\begin{aligned}&{\frac {\partial ^{2}y}{\partial t^{2}}}-c^{2}{\frac {\partial ^{2}y}{\partial x^{2}}}={\frac {y_{j}^{n+1}+y_{j}^{n-1}-2y_{j}^{n}}{(\Delta t)^{2}}}+O\left(\Delta t^{2}\right)-c^{2}{\frac {y_{j+1}^{n}+y_{j-1}^{n}-2y_{j}^{n}}{(\Delta x)^{2}}}+O\left(\Delta x^{2}\right)\\&{\frac {\partial ^{2}y}{\partial t^{2}}}-c^{2}{\frac {\partial ^{2}y}{\partial x^{2}}}-\left({\frac {y_{j}^{n+1}+y_{j}^{n-1}-2y_{j}^{n}}{(\Delta t)^{2}}}-c^{2}{\frac {y_{j+1}^{n}+y_{j-1}^{n}-2y_{j}^{n}}{(\Delta x)^{2}}}\right)=O\left(\Delta t^{2}\right)+O\left(\Delta x^{2}\right)\\&Dy\left(t,x\right)-D_{\left(\Delta t\Delta x\right)}y\left(t,x\right)=O\left(\Delta t^{2}\right)+O\left(\Delta x^{2}\right)\end{aligned}}$

como estamos supondo que os termos de ordem 2 são desprezíveis, $O\left(\Delta t^{2}\right)\rightarrow 0$ e $O\left(\Delta x^{2}\right)\rightarrow 0$. Logo :

$Dy\left(t,x\right)-D_{\left(\Delta t\Delta x\right)}y\left(t,x\right)\rightarrow 0$

portanto o esquema de diferenças finitas é consistente. Como na seção anterior já havíamos provado a estabilidade condicional do método, podemos concluir que o esquema é condicionalmente convergente (é consistente e condicionalmente estável) e portanto representa uma boa aproximação para a solução exata da equação diferencial.

Análise espectral

Uma possível forma para quantitativamente analisar o som gerado por uma corda vibrante é estudar as frequências que compõem o seu movimento, técnica essa chamada de análise espectral. Antes de prosseguirmos vamos recapitular alguns resultados da álgebra linear

Supremacia da álgebra linear

O seguinte conjunto $\mathbb {R} ^{\mathbb {R} }=\{f~|~f:\mathbb {R} \to \mathbb {R} \}$ é o espaço de funções reais de uma variável. Esse conjunto é um espaço vetorial, logo podemos utilizar toda a artilharia da álgebra linear, em especial, estamos interessados no sub-espaço gerado pela base $B=\{sen(\omega t)/{\sqrt {\pi }},cos(\omega t)/{\sqrt {\pi }}\}_{\omega \in \mathbb {R^{+}} }$ , pois elementos de $B$, interpretados como sinais sonoros, representam um frequência pura de valor $f=\omega /(2\pi )$. Dessa forma, um sinal arbitrário $s(t)$ pode ser escrito em termos das frequências puras que o formam

$s(t)={\frac {1}{\sqrt {\pi }}}\int _{0}^{\infty }[a(\omega )cos(\omega t)+b(\omega )sen(\omega t]d\omega$

E podemos extrair suas coordenadas (conhecidas como transformada de Fourier do sinal) ($a(\omega )$ e $b(\omega )$), fazendo o produto escalar com os elementos da base

${\begin{aligned}a(\omega )&=\int _{-\infty }^{\infty }s(t){\frac {1}{\sqrt {\pi }}}cos(\omega t)dt\\b(\omega )&=\int _{-\infty }^{\infty }s(t){\frac {1}{\sqrt {\pi }}}sen(\omega t)dt\end{aligned}}$

Se o domínio de $s(t)$ é limitado, digamos $t\in [0,T]$, então uma base infinita com cardinalidade enumerável (em contraste com a base anterior, que possui cardinalidade não enumerável) é suficiente para representar $s(t)$, uma possível base ortonormal é a seguinte: ${\bigg \{}{\sqrt {\frac {1}{T}}},{\sqrt {\frac {2}{T}}}sen(\omega _{n}t),{\sqrt {\frac {2}{T}}}cos(\omega _{n}t){\bigg \}}_{n\in \mathbb {N^{*}} }$, em que $\omega _{n}={\frac {n\pi }{T}}$. Dessa forma, a representação e coordenadas de $s(t)$ ficam

${\begin{aligned}s(t)&=a_{0}{\sqrt {\frac {1}{T}}}+{\sqrt {\frac {2}{T}}}\sum _{n=1}^{\infty }{\bigg [}a_{n}cos(\omega _{n}t)+b_{n}sen(\omega _{n}){\bigg ]}\\a_{0}&=\int _{0}^{T}s(t){\sqrt {\frac {1}{T}}}dt\\a_{n}&=\int _{0}^{T}s(t){\sqrt {\frac {2}{T}}}cos(\omega _{n}t)dt\\b_{n}&=\int _{0}^{T}s(t){\sqrt {\frac {2}{T}}}sen(\omega _{n}t)dt\\\end{aligned}}$

É impossível falar sobre bases enumeráveis de um sub-espaço de $\mathbb {R} ^{\mathbb {R} }$ sem representar esse canhão matemático com uma animação. Abaixo segue uma visualização que calcula as primeiras $N$ coordenadas ($a_{0},\dots ,a_{N}$ e $b_{1},\dots ,b_{N}$) de um sinal qualquer e sobrepõem a série obtida, incrementando $N$ até as duas curvas serem indistinguíveis a olho nu.

Animação de uma série de fourier. É interessante notar que a série converge muito bem com apenas uma dúzia de frequências.

Potência espectral

Adiante vamos ver que o sinal sonoro provém da vibração de um ponto específico da mesma, digamos em $x=x_{o}$, então a função que representa o sinal é $y(x_{o},t)$. Como estamos interessados nas frequências que compõem o sinal, será calculado a transformada de fourier de $y(x_{o},t)$ e definido que a potência da frequência $f=\omega /(2\pi )$ é $a_{\omega }^{2}+b_{\omega }^{2}$. A potência em função da frequência é o resultado da análise espectral.

Simulando uma corda de violão

Uma corda de violão geralmente é excitada por uma pancada dada por uma palheta ou pelo próprio dedo/unha do violonista. Essa pancada define uma condição inicial para a equação de onda. Uma suposição razoável da condição gerada é a seguinte

${\begin{aligned}y(x,t=0)&={\begin{cases}{\frac {h}{x_{0}}}x&,~0\leq x\leq x_{0}\\{\frac {h}{L-x_{0}}}(-x+L)&,~x_{0}\leq x\leq L\\\end{cases}}\\{\frac {\partial }{\partial t}}y(x,t=0)&=0\end{aligned}}$

supondo que a corda possui comprimento $L$ e a pancado ocorreu em $x_{o}$, causando uma deslocamento máximo $h$. A imagem a seguir ilustra o estado da corda logo após a excitação

Condição inicial de uma corda de violão.

Dado essa condição inicial, podemos evoluir temporalmente o estado da corda utilizando a equação da onda, mantendo as bordas fixas ($y(x=0,t)=y(x=L,t)=0$), mas como extrair som dessa simulação? Para responder essa pergunta precisamos saber como um violão gera som. Ao contrário do que inicialmente pareça, as ondas sonoras não são diretamente geradas pela vibração das cordas, mas sim da caixa do violão, que está diretamente conectada com as cordas em uma peça chamada ponte. A vibração das cordas gera uma força dependente do tempo que atua na caixa através da ponte, assim vibrando a caixa e gerando o som que escutamos. Portanto, para gerar som de forma realista, precisaríamos fisicamente simular a caixa, levando em consideração a sua geometria e as propriedades física do seu material, e então determinar as ondas de pressão que seriam geradas por essa vibração, o que está fora do escopo do presente trabalho. Felizmente, as seguintes simplificações vão nos permitir calcular as ondas de pressão:

• A força que a ponte exerce na caixa é aproximadamente proporcional a sua velocidade (Dinâmica aristotélica).
• A onda de pressão produzida pela caixa é aproximadamente proporcional a sua velocidade.

A força que a ponte exerce na caixa ($F(t)$) é a força que a corda exerce na ponte (pois a ponte está firmemente conectada na caixa), e essa força é proporcional a inclinação da corda, ou seja

$F(t)=T{\frac {\partial }{\partial x}}y(x=0,t)\approx T{\frac {y(\Delta x,t)-\overbrace {y(0,t)} ^{0}}{\Delta x}}={\frac {T}{\Delta x}}y(\Delta x,t)$

em que $T$ é a tensão na corda. Portando, o sinal da onda de pressão pode ser aproximado como o deslocamento de um ponto da corda próximo a ponte. Para gerar som com esse sinal, o tipo de arquivo WAV] é utilizado, pois o seu dado de entrada pode ser justante o sinal da onda de pressão.

Agora tudo está pronto para fazermos uma simulação, vamos tentar reproduzir a nota lá ($f=440~Hz$). Primeiro precisamos descobrir qual deve ser a velocidade de propagação das ondas ($c$) para gerar a nota em questão. Resolvendo a equação da onda ideal, obtemos que as possíveis frequência que existem na solução são $f_{n}={\frac {nc}{2L}}$, então não é possível fazer a corda vibrar apenas com uma frequência, mas quando músicos se referem a uma nota, eles realmente estão se referindo a $f_{1}$, a primeira frequência que compõem o sinal, portanto

$f_{1}={\frac {c}{2L}}=440~Hz\Rightarrow c=L\cdot 880~Hz$

Assumindo $L=1~m$, temos que $c=880~m/s$. Resolvendo a equação da onda realista pelo método FTCS com os seguinte parâmetros

Parâmetros da condição inicial

Parâmetros	Valor
$x_{o}$	$0,5~L$
h	1 cm

Parâmetros da integração

Parâmetros	Valor
k	1/4
c	880 m/s
$\Delta x$	0,01 m
L	1 m
b	5 $s^{-1}$
$\epsilon$	$10^{-9}$

obtemos o seguinte resultado

Solução da equação de onda ideal e realista.

Na animação acima, por questão de comparação, também foi coloca a solução de uma simulação com os mesmos parâmetros, apenas com a modificação $b=\epsilon =0$, que seria o caso ideal (sem amortecimento e rigidez).

Realizando a análise espectral do deslocamento de um ponto próximo de $x=0$, obtemos o seguinte gráfico

Espectro do deslocamento de um ponto próximo de $x=0$. Linhas pretas tracejadas representam os harmônicos.

As linhas verticais pretas e tracejadas representam as possíveis frequências da solução do caso ideal (também chamadas de harmônicos), ou seja, a primeira linha está em $f=440Hz$ como era de se esperar, mas é notável que a segunda e quarta linhas possuem uma potência praticamente nula, isso ocorre porque a excitação inicial foi exatamente no meio, e as frequência correspondentes a essas linhas provém de ondas estacionário que possum um nodo em $x=L/2$, logo elas não foram excitadas pela condição inicial. Mudando o ponto de excitação inicial, podemos ver essas frequências aparecendo. Rodando a simulação novamente, mas com $x_{o}=0,2~L$, obtemos o seguinte espectro

Espectro para excitação inicial em $x_{0}=0,2~L$

Para observar a influência da rigidez, a simulação foi rodada com $\epsilon =10^{-8}$ (10 vezes maior do que o valor anterior) resultando no seguinte espectro

Espectro para $\epsilon =10^{-8}$

O efeito mais notável é o deslocamento para esquerda das frequência em relação as frequências que compõem o sinal no caso ideal.

Por fim, conforme descrito anteriormente, podemos gerar som com as dados da simulação, mas infelizmente essa wiki não nos permite upar arquivos de áudio. No entanto, áudio das simulações discutidas aqui podem ser encontrados no repositório do projeto ana-julia dentro da pasta "sound". Em especial, o efeito de deslocamento de frequências causado pela rigidez é bastante evidente escutando os áudios.

Bônus

Apenas por diversão, segue uma animação de um pacote gaussiano como condição inicial

Animação da solução da equação da onda com um pacote gaussiano como condição inicial.

Referências

1. Giordano, Nicholas, Nakanishin Hisao. Computacional Physics, Sencond Edition - 2006.
2. Strikwerda, John. Finite Diference Schemes and Partial Diferential Equations, Second Edition, SIAM - 2004.
3. Moraes, Alciney das Neves. Critério de Estabilidade de um Esquema Explícito em Diferenças Finitas para o Modelo de Placas de Mindlin-Timoshenko, Universidade Federal do Pará — 2019.
4. Grigoryan, Viktor. Finite differences for the wave equation. UC Santa Bárbara Mathematics - 2012.
5. Boyce, E., W. Equações Diferenciais Elementares e Problemas de Valores de Contorno, 11ª Edição - 2020.