Grupo: Lucas Dória, Caetano Pires e Ânderson Rosa.
O objetivo deste trabalho é analisar três métodos diferentes de integração numérica para resolução de equações diferenciais parciais (EDP's). A equação aqui abordada é a equação da onda, a qual consiste em uma EDP hiperbólica de segunda ordem. Uma pequena noção dos métodos numéricos é dada, após a explicação do problema físico abordado. Os dados encontrados são comentados e comparados com os resultados esperados. Simulações também foram feitas para melhor visualização dos fenômenos. Por último, foi feita uma análise de erros e estabilidade para cada um dos métodos, assim como algumas conclusões foram inferidas.
Introdução
Equações diferenciais parciais hiperbólicas geralmente podem ser formuladas a partir de teoremas de conservação. Um exemplo é a equação do tipo:
onde
é o vetor de densidades da quantidade conservada, i.e.,
é o fluxo de densidade e
é um termo genérico representando fontes e/ou sumidouros.
Uma classe especial de equações hiperbólicas são as chamadas equações de adveção, na qual a derivada temporal da quantidade conservada
é proporcional à sua derivada espacial. Nesses casos,
é diagonal e dada por:
onde
é a matriz identidade.
Considerando apenas uma dimensão e com
, temos a equação de adveção:
onde
é a velocidade de propagação do pulso gerado. A equação admite uma solução analítica da forma
representando um pulso se movendo na direção
A equação da onda em uma dimensão é uma EDP hiperbólica de segunda ordem dada por
onde
é o deslocamento vertical de uma posição
no instante
e admite duas soluções, representadas por pulsos, na forma
e
Assumindo que
não depende da posição na equação da onda, nos restringimos a problemas lineares. Além disso, se escrevermos
então a equação da onda pode ser escrita como um sistema de três equações diferenciais de primeira ordem:
Em notação vetorial, o sistema acima pode ser reescrito na forma conservativa como:
onde
Condição CFL
Uma condição necessária mas nem sempre suficiente para a convergência de equações diferenciais parciais resolvidas a partir de métodos de diferença, conhecida como condição CFL [1], é formulada a partir do termo de “domínio de dependência”.
Considerando, por exemplo, a equação da advecção [1] em uma aproximação a partir de métodos de diferença em sua forma explícita
diz-se que o valor de
depende dos valores anteriores de
e
, e esses dois pontos dependem, novamente, de outros dois pontos em tempos anteriores. Todos esses pontos dependentes formam o domínio de dependência do valor
, representado abaixo (botar figura do domínio de dependência).
A condição CFL enuncia que a condição mínima para que haja estabilidade em métodos de diferenças o domínio de dependência da equação diferencia parcial, dado por sua equação característica, deve estar situado dentro do domínio de dependência do esquema numérico.
A partir dessa condição, define-se o número CFL como
onde
é um passo na malha temporal e
é um passo na malha espacial. Para casos unidimensionais, como o que será tratado aqui, a condição CFL é satisfeita se r
Um uso da condição CFL é determinar o tamanho do passo temporal, sabendo-se
e
:
onde o fator
.
O número CFL pode ser entendido como um fluxo numérico advectivo adimensionalisado pelas malhas espacial e temporal do problema. De um ponto de vista matemático, ele garante que o domínio numérico de dependência é sempre maior que o domínio físico. De um ponto de vista físico, é garantido que a velocidade de propagação de qualquer perturbação, como uma onda, seja menor ou ao menos igual que a velocidade de propagação numérica, fazendo com que a distância propagada pela perturbação não seja maior do que a divisão da malha espacial:
O Problema Físico
O Modelo de Corda Ideal
Para uma primeira abordagem da equação da onda, consideramos uma corda com suas extremidades fixas. Podemos dividir o comprimento
dessa corda em
intervalos de comprimentos iguais, dessa forma
. Cada intervalo é discretizado, representado por
,
. Também podemos dividir o tempo em intervalos iguais
e denotá-los como
,
. Abaixo temos uma esquematização das informações que precisamos para cada ponto da corda:
Tendo feita a discretização das variáveis, podemos aproximar a equação da onda por diferenciação finita, utilizando derivadas centradas da seguinte forma:
onde
representa o valor discretizado de
.
Assim, chegamos em uma equação discretizada:
.
Sabendo que essa discretização da equação da onda pode ser verificada como sendo o método Leapfrog (ver seção do método de Leapfrog), podemos resolver a equação para
para sabermos o deslocamento vertical de uma partição da corda no momento de tempo seguinte, assim obtendo
,
onde
Um Quadro Mais Realístico - O Modelo de Corda Rígida
Para nos aproximarmos de um modelo mais real, podemos adicionar um termo à equação original da onda que corresponde ao efeito de fricção em uma corda. A equação da onda mais geral com efeito de fricção pode ser escrita como [2]
onde
é a velocidade transversal de propagação do pulso na corda, dada pela relação
(sendo
a tensão na corda e
a densidade linear da mesma),
é um parâmetro adimensional de fricção que representa a rigidez da corda e
o comprimento da corda.
O parâmetro
é dado por
Falhou ao verificar gramática (erro de sintaxe): {\displaystyle \epsilon = \kappa² \frac{E S}{T L^2}}
,
onde
é o raio de giro da corda,
é o Módulo de Young e
a área da secção da corda.
Discretizamos a equação da seguinte maneira:
e resolvemos para
, obtendo:
O fato de essa discretização depender do deslocamento da corda em posições
e
implica em precisarmos simular "pontos fantasmas" quando integramos os extremos das cordas. Para fazermos isso, podemos ou utilizar a aproximação
ou podemos considerar esses "pontos fantasmas" como pontos presos e, portanto, sempre iguais a zero.
Os Métodos Utilizados
Foi realizada uma abordagem ao problema da corda real a partir de três métodos diferentes de integração numérica.
O método mais básico é chamado de FTCS (Forward-Time-Centered-Space) e consiste em duas expansões de Taylor ao redor do ponto
:
Subtraindo as duas expressões, encontramos a expressão
,
A qual podemos substituir na equação da onda, juntamente com a discretização da derivada parcial temporal. Temos então que, para um sistema linear de equações hiperbólicas:
Visto que essa última notação é mais genérica, ela será utilizada para a explicação dos métodos posteriores.
O Método de Lax-Friedrichs
O método de Lax-Friedrichs consiste em substituir o termo
com sua respectiva média espacial, i.e.,
. Logo, temos a seguinte equação de recorrência:
O Método de Leapfrog
Tanto o método FTCS quanto o método de Lax-Friedrichs são métodos de primeira ordem para a derivada temporal. Nessas circunstâncias,
deve ser significantemente menor do que
, muito abaixo do limite imposto pela condição de Courant.
Uma nova expressão para a derivada temporal, com precisão de segunda ordem é dada por
Substituindo a nova expressão acima no método de FTCS discutido anteriormente, encontramos o método de Leapfrog:
Como o método de Leapfrog foi o mais aplicado na resolução do problema da onda, é interessante um aprofundamento maior do método. Podemos adaptar o método de Leapfrog para o sistema de equações definido para a equação da onda ao fazermos a discretização ao redor de um ponto qualquer
no instante de tempo
:
Com essas discretizações, podemos utilizar as equações (1) e (2) para resolvermos para
apenas em termos de
,
e
. Fazendo as substituições dessas duas equações uma na outra, obtemos:
Dessa equação, chegamos a
Utilizando o mesmo raciocínio, podemos também resolver para
e obter
Com essas duas equações, podemos fazer uma integração utilizando o método de Euler para obter
, ou seja, o deslocamento de um determinado ponto
da corda no próximo instante de tempo:
Contudo, podemos fazer uma simples substituição das equações
e
nas equações
e
e, assim, obtemos que a representação de Leapfrog da equação da onda é dada pela discretização de segunda ordem da própria equação da onda, com
. Isso nos dá uma solução de "um passo", onde só precisamos efetuar o cálculo da equação discretizada pelo método da derivação finita, ou seja, para o caso da corda ideal, a equação
O Método de Lax-Wendroff
O método de Lax-Wendroff é a extensão do método de Lax-Friedrichs de segunda ordem. Calculamos o vetor
a partir de um passo médio de Lax-Friedrichs:
,
,
E encontramos os fluxos
a partir dos valores de
Logo, com um meio passo de Leapfrog, temos a expressão final do método:
Análise e Discussão dos Resultados
Escolhemos para simular quatro diferentes cordas: as cordas C2, C4 e C7 de um piano [3] e uma corda com dados pré-estabelecidos encontrados em [4]. Os dados das cordas C2, C4 e C7 estão na tabela abaixo.
|
C2
|
C4
|
C7
|
Comprimento ( )
|
|
|
|
Massa( )
|
|
|
|
Tensão ( )
|
|
|
|
Divisões
|
|
|
|
Amostragem de sinal
|
|
|
|
Parâmetro de fricção
|
|
|
|
Desses dados, temos que a densidade linear de massa das cordas é dada por
onde
é a massa e
o comprimento da corda,
onde
é a tensão na corda,
e
Para a corda com dados pré-estabelecidos, utilizamos uma corda com 2 metros de comprimento, com velocidade de propagação da onda sendo 300
,
de 0.01
,
e parâmetro de fricção de
.
Supondo a corda inicialmente em repouso, temos que em
a corda recebe em seu centro o equivalente à batida do martelo do piano. Supomos que esse estímulo possuía o formato aproximado de uma Gaussiana com amplitude de
do comprimento de cada corda. Então, com o estado inicial sendo um pulso com o formato de um pacote gaussiano e os dados da tabela, simulamos a propagação de ondas em cada uma das cordas.
Utilizando o método de Leapfrog, foi realizada uma primeira simulação para uma onda ideal em uma corda:
Simulação de uma onda em uma corda ideal. A corda simulada possuía comprimento de 2 metros e a amplitude inicial do pulso gaussiano tinha amplitude de 0.5 metros. Foi utilizado uma velocidade de propagação de 300 m/s,

e

.
Corda C2:
Simulação de uma corda C2 de um piano com os parâmetros da tabela acima.
Estado da corda C2 em diferentes instantes de tempo.
Corda C4:
Simulação de uma corda C4 de um piano com os parâmetros da tabela acima.
Estado da corda C4 em diferentes instantes de tempo.
Corda C7:
Simulação de uma corda C7 de um piano com os parâmetros da tabela acima.
Estado da corda C7 em diferentes instantes de tempo.
Corda com as definições do Giordano:
Simulação da corda exemplificada no livro Computational Physics [1]. Os parâmetros utilizados estão descritos no texto.
Estado da corda usada como exemplo pelo autor Giordano em diferentes instantes de tempo.
É possível observar que o tempo necessário para que, devido a fricção, os movimentos na corda diferenciem-se cada vez mais da solução para a corda ideal é menor para um tamanho de corda menor, ou seja, a fricção acaba sendo mais relevante mais rapidamente em cordas mais curtas.
Análise de Erro e Estabilidade dos Métodos
Estabilidade do método Leapfrog
Pela estabilidade de Von Neumann, podemos escrever que
Para
, a equação da discretização da onda pode ser reescrita como
Essa escolha com
nos dá a solução exata sem dispersão numérica, como pode ser verificado pela condição CFL. Contudo,
é válido somente no caso de uma corda ideal.
É conveniente escrever a condição acima em termos da amostragem de sinal
e a frequência fundamental da corda
, o que nos leva a
O Teorema de Nyquist diz que a frequência superior no espectro deve ser menor do que
para evitar serrilhamento e para garantir uma reconstrução única e contínua. Logo, no caso ideal quando as autofrequências da corda são igualmente espaçadas (
), a condição de Nyquist indica que o número máximo de frequências no espectro é
. Isso significa que essa condição pode ser usada para selecionar o número apropriado de pontos espaciais
para a corda. Entretanto, como
é um inteiro, apenas valores discretos da frequência natural podem ser obtidos sem erros de trucamento, ou seja, usando
. Como series discretas não costumam ser utilizadas, precisamos aceitar pequenos erros de truncamento para ajustar
, ou seja, utilizando
. No caso de uma onda com fricção, temos que
é um valor de boa estabilidade.
Referências Bibliográficas
- ↑ Luciano Rezzolla, "Numerical Methods for the Solution of Partial Differential Equations". Albert Einstein Institute, Max-Planck-Institute for Gravitational Physics, Potsdam, Germany.
- ↑ N. J. Giordano, "Computational Physics". Department of Physics, Purdue University. Upper Saddle River, New Jersey. Prentice-Hall, 1997.
- ↑ Antoine Chaign, Anders Askenfelt, "Numerical simulations of piano strings. I. A physical model for a struck string using finite difference methods". Signal Department, Paris, France. Department of Speech Communication and Music Acoustics; Royal Institute of Technology, Stockholm, Sweden.
- ↑ N. J. Giordano, "Computational Physics". Department of Physics, Purdue University. Upper Saddle River, New Jersey. Prentice-Hall, 1997.