Corda Vibrante

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_{i}\right\}_{i=0}^{N}$ e $\left\{t_{i}\right\}_{j=0}^{M}$ . Assim obtemos uma grade de malhas de pontos $(t_{i},x_{i})$ onde a função $y(t,x)$ tem seus valores aproximados por $y(t_{i},x_{i})$ 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(t+\Delta t,x)-y((t+\Delta t)-\Delta t,x)}{\Delta t}}-{\frac {y(t,x)-y(t-\Delta t,x)}{\Delta t}}\right)={\frac {y(t+\Delta t,x)+y(t-\Delta t,x)-2y(t,x)}{(\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 do calor, 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 )$ e a identidade trigonométrica fundamental, 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 esse termo é real. 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 $\Delta x/\Delta t$.

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

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 :

${\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)-b{\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$;

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 (também chamadas de 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

Agora, considerando uma corda vibrante, 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 esse 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 vamos definir 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 exercer 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, então

$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$

As linhas verticais pretas e tracejadas representam as possíveis frequências da solução do caso ideal, 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$

Podemos observar a influência da rigidez no som produzido rodando a simulação com $\epsilon =10^{-}8$ (10 vezes maior do que o valor anterior) e olhar o espectro produzido

Espectro para math>\epsilon=10^-8</math>

O efeito mais notável é o deslocamento para esquerda das frequência 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 aqui0 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.