Esquema Explícito Semi-Analítico para a Solução da Equação da Onda Unidimensional com Condições de Contorno Naturais

RESUMO

Este artigo aborda o desenvolvimento de um esquema explícito semi-analítico destinado à solução da equação da onda unidimensional não homogênea. Uma malha espaço-tempo foi construída através da relação incremental $c\mathrm{\Delta }t=\mathrm{\Delta }x$, de tal modo que volumes de controle foram formados a partir das retas características da referida equação. O desenvolvimento se deu a partir da forma integral da lei de conservação sobre estes volumes de controle. O esquema deduzido possui a propriedade de erro de truncamento local nulo, mesmo nos casos não homogêneo ou de solução generalizada (visto que satisfaz a lei integral). O método proporciona facilidade na inclusão das condições iniciais e de contorno, inseridas também sem necessidade de técnicas de aproximação. Diante dos desenvolvimentos e dos experimentos numéricos realizados, concluímos que o esquema proposto é uma excelente técnica numérica, com ótima acurácia e robustez para a resolução do problema de onda linear unidimensional.

Palavras-chave:
Equação da onda; Erro de truncamento nulo; Condições de Neumann

ABSTRACT

This paper discusses the development of an semi-analytical explicit scheme for solving the non-homogeneous one-dimensional wave equation. A space-time mesh was constructed through the incremental relation $c\mathrm{\Delta }t=\mathrm{\Delta }x$, such that control volume were formed from the characteristic lines of said equation. Schemes were developed from the integral form of the conservation law on these control volumes. They have the property of null local truncation error, even in cases not homogeneous or generalized solution. The method facilitates the inclusion of the initial and boundary conditions, also inserted without the need for approximation techniques. Considering the developments and numerical experiments performed, we conclude that the proposed scheme is an excellent numerical technique, with great accuracy and robustness to solve the one-dimensional linear wave problem.

Keywords:
Wave Equation; Null Truncation Error; Neumann Boundary Condition

1 INTRODUÇÃO

Seja $\Omega =\left(a,b\right)\times {\mathrm{ℝ}}_{+}^{*}\subset {\mathrm{ℝ}}^2$ e consideremos a equação hiperbólica unidimensional

$\frac{{\partial }^{2}u}{\partial t^2}-c^2\frac{{\partial }^{2}u}{\partial x^2}=f,\text{em}\mathrm{\Omega },$ (1.1)

onde a constante $c>0$ é a velocidade de propagação da onda; o termo f representa uma fonte externa; x é a direção no espaço em coordenadas cartesianas; t o tempo e a variável u um campo de pressão ⁴ ou um campo de deslocamento ⁶.

A motivação para este trabalho é o fato de que, no caso homogêneo ($f=0$), a equação (1.1) admite um esquema de diferenças finitas explícito com erro de truncamento localmente nulo ¹¹:

$\frac{u_i^{j+1}-2u_i^j+u_i^{j-1}}{\Delta t^2}=c^2\left(\frac{u_{i+1}^j-2u_i^j+u_{i-1}^j}{\Delta x^2}\right),1\le i\le n-1,j\ge 1,$ (1.2)

bastando, para isso, impor que^* $c\mathrm{\Delta }t=\mathrm{\Delta }x$. De fato, supondo que $u\in C^{\infty }\left(\Omega \right)$ e usando série de Taylor, podemos substituir a solução verdadeira u(x, t) em (1.2), para obtermos

${\overline{)\left\{\frac{{\partial }^{2}u}{\partial t^2}+2\sum _{k=2}^{\infty }\frac{\Delta t^{2\left(k-1\right)}}{\left(2k\right){\displaystyle !}}\frac{{\partial }^{2k}u}{\partial t^{2k}}\right\}}}_{\left(x_i,t_j\right)}={\overline{)\left\{c^2\frac{{\displaystyle {\partial }^{2}u}}{{\displaystyle \partial x^2}}+2c^2\sum _{k=2}^{\infty }\frac{{\displaystyle \Delta x^{2\left(k-1\right)}}}{{\displaystyle \left(2k\right)!}}\frac{{\partial }^{2k}u}{\partial x^{2k}}\right\}}}_{\left(x_i,t_j\right)}$ (1.3)

Mas, como u satisfaz (1.1) e na hipótese de que $u\in C^{\infty }\left(\Omega \right)$, segue-se então, pelo teorema de Schwarz, que $\frac{{\partial }^{2k}u}{\partial t^{2k}}=c^{2k}\frac{{\displaystyle {\partial }^{2k}u}}{{\displaystyle \partial x^{2k}}},k\ge 1$. Assim, os dois lados de (1.3) diferem por

${\mathcal{T}}_{\Psi }^2{\overline{)\left(x_i,t_j\right)=2\sum _{k=2}^{\infty }\frac{{\gamma }_k}{\left(2k\right){\displaystyle !}}\frac{{\partial }^{2k}u}{\partial x^{2k}}}}_{\left(x_i,t_j\right)},\text{em que}{\gamma }_k=\Delta t^{2\left(k-1\right)}{c}^{2k}-c^2\Delta x^{2\left(k-1\right)}.$ (1.4)

Portanto, se $c\mathrm{\Delta }t=\mathrm{\Delta }x$, então a equação (1.4), que descreve o erro de truncamento local para a aproximação (1.2), será igual a zero, isto é, ${\mathcal{T}}_{\Psi }^2\left(x_i,t_j\right)=0$. Em outras palavras, uma função u de classe C^∞ (Ω) que satisfaça a equação (1.1), satisfará também o esquema (1.2), quando $c\mathrm{\Delta }t=\mathrm{\Delta }x$. Por este motivo, esta relação é chamada de razão de ouro ¹¹, e espera-se uma acurácia perfeita, visto que o esquema é perfeitamente consistente. Observe que ao se impor $c\mathrm{\Delta }t=\mathrm{\Delta }x$, a equação (1.2) se torna:

$u_i^{j+1}=u_{i-1}^j+u_{i+1}^j-u_i^{j-1},1\le i\le n-1,j\ge 1,$ (1.5)

Quando $f\ne 0$, no entanto, não é possível obter ${\mathcal{T}}_{\Psi }^2\left(x_i,t_j\right)=0$, mesmo que se imponha $c\mathrm{\Delta }t=\mathrm{\Delta }x$. O objetivo deste trabalho consiste, por tanto, em desenvolver um esquema de diferenças finitas explícito com erro de truncamento localmente nulo para a equação não homogênea, bem como desenvolver esquemas para a inclusão das condições iniciais e de contorno (Neumann) que também contemplem esta propriedade de erro local nulo.

Para este fim, organizamos este trabalho a partir da seção 2, onde apresentamos o desenvolvimento de um esquema explícito de passo múltiplo para a solução da equação (1.1). Uma propriedade interessante decorre dessa dedução: soluções de classe C²(Ω), para a referida equação, também satisfazem o esquema resultante. Na subseção 2.1 realizamos tratamentos sobre as condições iniciais, com o intuito de se calcular o primeiro passo de tempo com a melhor acurácia possível. O fim das deduções ocorre com a seção 2.2, onde condições de contorno do tipo Neumann são consideradas. Experimentos numéricos que ilustram a precisão e a robustez das formulações são apresentados junto à seção 3. Encerramos com as conclusões na seção 4.

2 DESENVOLVIMENTO

O ponto chave para o desenvolvimento deste trabalho é a forma na qual discretizamos o domínio espaço-tempo. Considere uma malha Ψ com $n+1$ pontos de discretização (n elementos) no espaço, tal que $x_i=a+i\Delta x,0\le i\le n,t_j=j\Delta t,0\le j\le q$, com $2\le n,q\in \mathrm{\mathbb{N}}$ e $\Delta x=\left(b,a\right)/n$, sendo q o número de avanços no tempo. É conhecido que a equação (1.1) possui duas retas características: $x-ct=\xi$ e $x+ct=\eta$, em que ξ e η são duas constantes. Faremos o incremento espacial depender do incremento temporal por meio da relação $c\mathrm{\Delta }t=\mathrm{\Delta }x$, fornecida pelas retas características, conforme Figura 1a, ao mesmo tempo em que formamos os volumes de controle característicos, conforme Figura 1b.

Figura 1:
(a) Discretização do domínio $\overline{\Omega }$ por meio das características; (b) Formação dos volumes de controle característicos.

Integrando a equação (1.1), ao longo do volume de controle V representado pela Figura 1b e supondo que u seja uma função de classe $C^2\left(\overline{\Omega }\right)$, podemos aplicar o teorema da divergência de Gauss ⁹, de modo a obtermos

${\int }_{V}div\mathbf{F}dV={\oint }_S\mathbf{F}\mathbf{·}\mathbf{N}ds=\sum _{i=1}^4{\int }_{L_i}\mathbf{F}\mathbf{·}{\mathbf{N}}_{i}ds={\int }_{V}fdV,$ (2.1)

em que $\mathbf{F}=\left(-c^2{u}_x,u_t\right)$. Note que F é de classe ${\left(C^1\left(\overline{\Omega }\right)\right)}^2=C^1\left(\overline{\Omega }\right)\times C^1\left(\overline{\Omega }\right)$.

A fim de obtermos o fluxo ao longo da face L ₁, observamos que ${\mathbf{C}}_1\left(s\right)=\left(x_i+c\left(s-t_{j-1}\right),s\right)$, com $t_{j-1}\le s\le t_j$, é uma representação paramétrica desse segmento, de modo que seu respectivo vetor normal é ${\mathbf{N}}_1=\left(1,-c\right)$. Com isso, o fluxo ℱ₁ será dado por

$$\begin{gathered} {\mathcal{F}}_1={\int }_{L_1}\mathbf{F}\mathbf{·}{\mathbf{N}}_{1}ds={\int }_{L_1}\left(-c^2{u}_x,u_t\right)·\left(1,-c\right)ds=-c{\int }_{L_1}\left(c{u}_x+u_t\right)ds, \\ =-c{\int }_{L_1}\left(u_x,u_t\right)·\left(c,1\right)ds=-c{\int }_{t_{j-1}}^{t_j}\nabla u\left({\mathbf{C}}_1\left(s\right)\right)·{\mathbf{C}}_1^{\text{'}}\left(s\right)ds, \\ =-c\left[u\left(x_{i+1},t_j\right)-u\left(x_i,t_{j-1}\right)\right]. \end{gathered}$$ (2.2)

Os demais fluxos podem ser encontrados de modo análogo:

${\mathcal{F}}_2=c\left[u\left(x_i,t_{j+1}\right)-u\left(x_{i+1},t_j\right)\right],{\mathcal{F}}_3=-c\left[u\left(x_{i-1},t_j\right)-u\left(x_i,t_{j+1}\right)\right]$ (2.3)

e

${\mathcal{F}}_4=c\left[u\left(x_i,t_{j-1}\right)-u\left(x_{i-1},t_j\right)\right].$ (2.4)

A integral sobre o termo fonte f(x, t) será

${\int }_{V}fdV={\int }_{t_j}^{t_{j+1}}{\int }_{x_i-c\left(t_{j+1}-\tau \right)}^{x_i+c\left(t_{j+1}-\tau \right)}f\left(s,\tau \right)dsd\tau +{\int }_{t_{j-1}}^{t_j}{\int }_{x_i-c\left(\tau -t_{j-1}\right)}^{x_i+c\left(\tau -t_{j-1}\right)}f\left(s,\tau \right)dsd\tau .$ (2.5)

Substituindo as equações (2.2)-(2.5) na equação (2.1), teremos

$u\left(x_i,t_{j+1}\right)=u\left(x_{i-1},t_j\right)+u\left(x_{i+1},t_j\right)-u\left(x_i,t_{j-1}\right)+F\left(x_i,t_j\right),$ (2.6)

para $1\le i\le n-1,j\ge 1$, em que

$F\left(x_i,t_j\right)=\frac{1}{2c}{\int }_{t_j}^{t_{j+1}}{\int }_{x_i-c\left(t_{j+1}-\tau \right)}^{x_i+c\left(t_{j+1}-\tau \right)}f\left(s,\tau \right)dsd\tau +\frac{1}{2c}{\int }_{t_{j-1}}^{t_j}{\int }_{x_i-c\left(\tau -t_{j-1}\right)}^{x_i+c\left(\tau -t_{j-1}\right)}f\left(s,\tau \right)dsd\tau ,$ (2.7)

ou, na notação padrão de diferenças finitas:

$u_i^{j+1}=u_{i-1}^j+u_{i+1}^j-u_i^{j-1}+F_i^j,\mathrm{para}1\le i\le n-1,j\ge 1,$ (2.8)

em que $F_i^j$ é dado pela equação (2.7).

Observe que aproximação alguma foi utilizada na dedução da equação (2.8). Dessa forma, podemos concluir o seguinte

Teorema 2.1. Seja$u\in C^2\left(\overline{\Omega }\right)$uma solução daequação (1.1)e considere um conjunto de pontos$\left(x_i,t_j\right)\in \overline{\Omega }$tais que$x_i=a+i\Delta x,0\le i\le n,t_j=j\Delta t,0\le j\le q$, com$2\le n,q\in \mathrm{\mathbb{N}}$sendo q o número de avanços no tempo,$\Delta x=\left(b-a\right)/n$e$c\mathrm{\Delta }t=\mathrm{\Delta }x$. Nestas condições, u também é solução daequação (2.8).

Observe que se $f=0$, então a equação (2.8) se torna idêntica a (1.5). Por isso, é possível compreender (2.8) como um esquema explícito de diferenças finitas. Decorre do resultado anterior o seguinte

Corolário 2.1.1. O esquema dado pelaequação (2.8)possui erro de truncamento local nulo para funções u de classe$C^2\left(\overline{\Omega }\right)$.

Proof. [Demonstração] Seja $u\in C^2\left(\overline{\Omega }\right)$ uma solução em sentido clássico para a equação (1.1). Pelo teorema 2.1, u também é solução da equação (2.8). Sendo assim:

$$\begin{gathered} {\mathcal{T}}_{\psi }\left(x_i,t_j\right)=\left[u_{tt}\left(x_i,t_j\right)-c^2{u}_{xx}\left(x_i,t_j\right)-f\left(x_i,t_j\right)\right]-\left[u\left(x_i,t_{j+1}\right)-u\left(x_{i-1},t_j\right)-u\left(x_{i+1},t_j\right)+u\left(x_i{t}_{j-1}\right)-F\left(x_i,t_j\right)\right], \\ =0-0=0. \end{gathered}$$

◻

Quando a integral sobre f for de difícil avaliação, pode-se utilizar uma aproximação por quadratura. O esquema produzido é, em todo caso, de passo múltiplo, pois para obter a solução no tempo $j+1$ é necessário conhecê-la nos tempos j e $j-1$. Isso implica, que em termos globais, ainda existe um erro inerente que é devido às condições de contorno e iniciais que devem ser consideradas. Observemos, entretanto, que se $u_i^1$ fosse conhecido para $0\le i\le n$, então $u_i^j,2\le j\le q$ seria determinado de modo exato pela equação (2.8). Em outras palavras, para obtermos um esquema com erro globalmente nulo precisamos determinar $u_i^1,u_0^j$ e $u_n^j$ (isto é, as soluções no primeiro passo de tempo e nos contornos) de modo exato.

2.1 Cálculo do primeiro passo de tempo: uso das condições iniciais

A equação (2.8) não pode ser utilizada para $j=0$, tendo em vista que para este caso $u_i^{j-1}$ estaria indefinido. O passo inicial precisa ser calculado por meio de alguma expressão auxiliar. A fim de encontrar esta expressão auxiliar, seja o volume de controle triangular definido conforme Figura 2. Suponha ainda que as condições iniciais prescritas sejam

$u\left(x,0\right)=\varphi \left(x\right)\text{e}{u}_t\left(x,0\right)=\Psi \left(x\right).$ (2.9)

Figura 2:
Volume de controle auxiliar triangular para o cálculo da solução numérica no primeiro passo de tempo.

Integrando a equação (1.1) sobre este volume de controle e aplicando novamente o teorema da divergência de Gauss:

${\int }_{V}fdv={\int }_V\mathrm{div}\mathbf{F}dV={\oint }_S\mathbf{F}\mathbf{·}\mathbf{N}ds=\sum _{i=1}^3{\int }_{L_i}\mathbf{F}\mathbf{·}{\mathbf{N}}_{i}ds.$ (2.10)

Então, de forma análoga ao caso anterior, temos que os fluxos sobre as faces do triângulo serão:

• • Sobre a face L ₁:

${\mathcal{F}}_1^0={\int }_{L_i}\mathbf{F}\mathbf{·}{\mathbf{N}}_{1}ds=-{\int }_{x_{i-1}}^{x_{i+1}}{u}_t\left(x,t_0\right)dx=-{\int }_{x_{i-1}}^{x_{i+1}}\Psi \left(x\right)dx.$ (2.11)

• • Sobre a face L ₂:

${\mathcal{F}}_2^0={\int }_{L_2}\mathbf{F}\mathbf{·}{\mathbf{N}}_{2}ds=c{\int }_{L_2}{u}_{x}dx+u_{t}dt=c\left[u\left(x_i,t_1\right)-\varphi \left(x_{i+1}\right)\right].$ (2.12)

• • Sobre a face L ₃:

${\mathcal{F}}_3^0={\int }_{L_3}\mathbf{F}\mathbf{·}{\mathbf{N}}_{3}ds=-c{\int }_{L_3}{u}_{x}dx+u_{t}dt=-c\left[\varphi \left(x_{i-1}\right)-u\left(x_i,t_1\right)\right].$ (2.13)

• • A integral de domínio sobre o termo fonte é dada por:

${\int }_{V}f\left(x,t\right)dV={\int }_{t_0}^{t_1}{\int }_{x_i-c\left(t_1-\tau \right)}^{x_i+c\left(t_1-\tau \right)}f\left(s,\tau \right)dsd\tau .$ (2.14)

Então, substituindo as equações (2.11)-(2.14) em (2.10), teremos

$u\left(x_i,t_1\right)=\frac{\varphi \left(x_{i+1}\right)+\varphi \left(x_{i-1}\right)}{2}+\frac{1}{2c}{\int }_{x_{i-1}}^{x_{i+1}}\Psi \left(x\right)dx+\frac{1}{2c}{\int }_{t_0}^{t_1}{\int }_{x_i-c\left(t_1-\tau \right)}^{x_i+c\left(t_1-\tau \right)}f\left(s,\tau \right)dsd\tau ,$ (2.15)

com $1\le i\le n-1$.

A equação anterior é conhecida como a solução de D’Alembert para o problema de Cauchy e determina os valores de u no primeiro passo de tempo de modo exato. Se as condições de contorno são essenciais (ou de Dirichlet), então, após a determinação de u no primeiro passo de tempo com a equação (2.15), basta utilizarmos a equação (2.8) para determinarmos a variável nos demais passos de tempo.

O problema de valor inicial formado pelas equações (1.1) e (2.9) com condições de contorno de Dirichlet $u\left(a,t\right)=g\left(t\right)$ e $u\left(b,t\right)=h\left(t\right)$ é completamente determinado pelo esquema formado pelas equações (2.8) e (2.15). Em outras palavras, tal esquema representa a solução clássica avaliada nos pontos da malha.

Abordamos, na seção a seguir, a situação em que as condições de contorno dadas são do tipo natural (ou de Neumann).

2.2 Condições de contorno naturais

2.2.1 Contorno esquerdo ( $x=a$ )

Se as condições de contorno forem naturais (ou de Neumann), então o valor da variável u será desconhecido em um ou em ambos os contornos. Supondo que $u_x\left(a,t\right)=g\left(t\right)$ seja conhecido, então no contorno esquerdo, consideremos o volume de controle triangular representado pela Figura 3a. Ao integrar a equação (1.1) sobre este volume, encontramos, por meio do mesmo raciocínio, a seguinte equação

$$\begin{gathered} u\left(x_0,t_{j+1}\right)=2u\left(x_1,t_j\right)-u\left(x_0,t_{j-1}\right)+c{\int }_{t_{j-1}}^{t_{j+1}}g\left(\tau \right)d\tau + \\ \frac{1}{c}{\int }_{t_j}^{t_{j+1}}{\int }_{x_0}^{x_0+c\left(t_{j+1}-\tau \right)}f\left(s,\tau \right)dsd\tau + \\ \frac{1}{c}{\int }_{t_{j-1}}^{t_j}{\int }_{x_0}^{x_0+c\left(\tau -t_{j-1}\right)}f\left(s,\tau \right)dsd\tau ,\text{para}j\ge 1. \end{gathered}$$ (2.16)

Figura 3:
Volumes de controle auxiliares utilizados no contorno esquerdo para (a) $j\ge 1$ e (b) $j=0$.

Observamos que a equação anterior não pode ser utilizada quando $j=0$, tendo em vista que $u\left(x_0,t_{j-1}\right)$ não está definido. Para solucionar este problema, repetimos o raciocínio sobre a célula dada pela Figura 3b. A expressão resultante após as simplificações é

$u\left(x_0,t_1\right)=\varphi \left(x_1\right)+\frac{1}{c}{\int }_{x_0}^{x_1}\Psi \left(x\right)dx+c{\int }_{t_0}^{t_1}g\left(\tau \right)d\tau +\frac{1}{c}{\int }_{t_0}^{t_1}{\int }_{x_0}^{x_0+c\left(t_1-\tau \right)}f\left(s,\tau \right)dsd\tau ,$ (2.17)

para $i=j=0$.

2.2.2 Contorno direito ( $x=b$ )

Se $u_x\left(b,t\right)=h\left(t\right)$ é conhecido, então no contorno direito, utilizamos o volume de controle auxiliar triangular designado pela Figura 4a que, de modo análogo aos casos anteriores, geram após todas as simplificações:

$$\begin{gathered} u\left(x_n,t_{j+1}\right)=2u\left(x_{n-1},t_j\right)-u\left(x_n,t_{j-1}\right)+c{\int }_{t_{j-1}}^{t_{j+1}}h\left(\tau \right)d\tau + \\ \frac{1}{c}{\int }_{t_j}^{t_{j+1}}{\int }_{x_n-c\left(t_{j+1}-\tau \right)}^{x_n}f\left(s,\tau \right)dsd\tau + \\ \frac{1}{c}{\int }_{t_{j-1}}^{t_j}{\int }_{x_n-c\left(\tau -t_{j-1}\right)}^{x_n}f\left(s,\tau \right)dsd\tau ,\text{para}j\ge 1. \end{gathered}$$ (2.18)

Figura 4:
Volumes de controle auxiliares utilizados no contorno direito para (a) $j\ge 1$ e (b) $j=0$.

Observamos, novamente, que a equação (2.18) anterior não pode ser utilizada quando $j=0$, pois $u\left(x_n,t_{j-1}\right)$ também não está definido. Repetimos novamente o raciocínio sobre a célula representada pela Figura 4b de modo à obtermos

$u\left(x_n,t_1\right)=\varphi \left(x_{n-1}\right)+\frac{1}{c}{\int }_{x_{n-1}}^{x_n}\Psi \left(x\right)dx+c{\int }_{t_0}^{t_1}h\left(\tau \right)d\tau +\frac{1}{c}{\int }_{t_0}^{t_1}{\int }_{x_n-c\left(t_1-\tau \right)}^{x_n}f\left(s,\tau \right)dsd\tau ,$ (2.19)

para $i=n,j=0$.

Os pares de equações (2.16)-(2.17) e (2.18)-(2.19) são os esquemas de avanço no tempo utilizados nos contornos esquerdo e direito, respectivamente, e são deduzidos sem inserção de erro de truncamento. Uma possível fonte de erro ocorre quando as funções ψ, g, h ou f são difíceis de se integrar, neste caso, os esquemas se tornam sujeitos aos erros das fórmulas de quadratura.

3 EXPERIMENTOS NUMÉRICOS

Nesta seção consideramos alguns exemplos numéricos típicos para a equação da onda (1.1). A menos que se diga o contrário, fixamos os parâmetros relativos ao domínio em $a=0$, $b=2$ e $L=b-a=2$. A velocidade da onda c foi variada conforme o problema. O incremento temporal Δt obedece, conforme salientado nas deduções anteriores, a relação $\Delta t=\Delta x/c$, sendo $\Delta x=L/n$. Os resultados numéricos foram comparados com as soluções analíticas ou espectrais dos respectivos problemas, que podem ser todos encontradas em ¹⁰.

Exemplo 3.1 (Barra com condições de contorno de Dirichlet variável) O primeiro experimento considerado tem por objetivo avaliar a técnica proposta de inclusão das condições iniciais. Para isso, seja a equação (1.1) na forma homogênea com as prescrições $u\left(a,t\right)=u\left(x,0\right)=u_t\left(x,0\right)=0$ e $u\left(b,t\right)=\mathrm{sen}\left(\omega t\right)$, sendo $\omega =\pi {\left(3L\right)}^{-1}$. Substituindo estes dados nas equações (2.8) e (2.15), obtemos, respectivamente,

$u_i^{j+1}=u_{i-1}^j+u_{i+1}^j-u_i^{j-1},1\le i\le n-1,1\le j\le q,$ (3.1)

$u_i^1=\frac{u_{i-1}^0+u_{i+1}^0}{2},1\le i\le n-1.$ (3.2)

A equação (3.1) é utilizada no interior da malha, enquanto que (3.2) é utilizada no primeiro passo de tempo. É importante salientar o número reduzido de operações envolvidas: 2 operações (1 soma e 1 divisão) para os cálculos do primeiro passo de tempo e 2 operações (1 soma e 1 subtração) nos demais passos. Os resultados numéricos foram obtidos utilizando $n+1=3$ pontos ao longo do espaço, incrementos $\Delta x=1$ e $c=1$. Foram tomados incrementos grosseiros para comprovar que o esquema é indiferente a magnitude do espaçamento. Na Figura 5a está apresentado o gráfico da solução numérica (pontos) e espectral (traço contínuo), avaliada ao longo do tempo sobre o único ponto interior da malha $x_1=1$. A Figura 5b contém uma ampliação da mesma.

Figura 5:
(a) Solução espectral versus solução numérica para a equação linear da onda unidimensional homogênea com condições de contorno do tipo Dirichlet avaliada no ponto $x=1$; (b) Ampliação da Figura 5a.

A solução espectral para este problema pode ser determinada por série de Fourier ¹⁰. A Tabela 1 computa os desvios relativo e absoluto avaliados no tempo $t=195$ s (observamos na Figura 5b que este é um dos picos do gráfico), calculados nas normas L ¹ e L ², dadas por

${||\mathbf{x}||}_1=\Delta x\sum _{i=0}^n\left|x_i\right|e{||\mathbf{x}||}_2={\left(\Delta x\sum _{i=0}^n{\left|x_i\right|}^2\right)}^{1/2},$ (3.3)

respectivamente. Inicialmente truncamos a solução espectral em $N=29$ termos, quando o critétio ${\hat{u}}_k-{\hat{u}}_{k-1}|<{10}^{-20}$ foi atingido (em que ${\hat{u}}_k$ é coeficiente de Fourier). A diferença absoluta e relativa concentrou-se em torno de 10^-3 (ver Tabela 1), não corroborando com a propriedade de erro localmente nulo.

Tabela 1:
Comparação entre as soluções numérica e espectral (truncada em diversos termos N), avaliadas no tempo $t=195$ s, para o problema da barra com condições de contorno de Dirichlet.

N	Relativo (L 2 )	Relativo (L 1 )	Absoluto (L 2 )	Absoluto (L 1 )
16	4,46032 × 10-3	3,16282 × 10-3	5,88330 × 10-3	5,88330 × 10-3
29	2,36138 × 10-3	1,67426 × 10-3	3,11899 × 10-3	3,11899 × 10-3
128	6,17200 × 10-4	4,37564 × 10-4	8,16149 × 10-4	8,16149 × 10-4
1.024	7,81780 × 10-5	5,54227 × 10-5	1,03414 × 10-4	1,03414 × 10-4
8.192	9,78851 × 10-6	6,93934 × 10-6	1,29489 × 10-5	1,29489 × 10-5
131.072	6,11918 × 10-7	4,33805 × 10-7	8,09492 × 10-7	8,09492 × 10-7
524.288	1,52981 × 10-7	1,08453 × 10-7	2,02375 × 10-7	2,02375 × 10-7

Com o objetivo de explorar o erro numérico, fixamos todos os parâmetros bem como a solução numérica obtida pelo sistema (3.1)-(3.2) e variamos (apenas) a quantidade de termos no truncamento da solução espectral. As diferenças relativas e absolutas também estão computadas na Tabela 1. Observamos que quanto mais termos são considerados na solução em série, menor se torna a diferença relativa e absoluta entre a solução numérica e espectral (salientamos que a solução numérica está fixa). Neste sentido, os erros apresentados na referida tabela são constituídos por erros de truncamento da série e erros computacionais de arredondamento.

O próximo exemplo retrata a aplicação do método proposto na equação não homogênea em conjunto com condições de contorno do tipo Neumann.

Exemplo 3.2 (Problema não homogêneo com condições de Neumann) Seja agora o problema com condições iniciais: $u\left(x,0\right)=0$ e $u_t\left(x,0\right)=\omega \left[\cos \left(\lambda x\right)-x^2\right]$; e condições de contorno de Neumann: $u_t\left(a,t\right)=0$ e $u_t\left(b,t\right)=-4\mathrm{sen}\left(\omega t\right)$, para $\omega =\pi {\left(3L\right)}^{-1}$ e $\lambda =\pi L^{-1}$. Considere, ainda, que a equação (1.1) seja não homogênea com termo fonte

$f\left(x,t\right)=\mathrm{sen}\left(\omega t\right)\left[\cos \left(\lambda x\right)\left(c^2{\lambda }^2-{\omega }^2\right)+{\omega }^2\left(x^2-t\right)+2c^2\right]+2\omega \cos \left(\omega t\right).$ (3.4)

Este problema tem solução analítica^† $u\left(x,t\right)=\mathrm{sen}\left(\omega t\right)\left[\cos \left(\lambda x\right)-x^2-t\right]$. Para o cálculo da solução numérica utilizamos $n+1=3$ pontos: o contorno esquerdo ($x=a$), um ponto interior ($x=a+\Delta x$) e o contorno direito ($x=b$). Dessa forma os incrementos utilizados foram $\Delta x=1,0$ e $\Delta t=0,25$. O esquema utilizado envolve todas as equações deduzidas nos parágrafos anteriores (equações (2.8) e (2.15)-(2.19)). As integrais sobre o termo fonte foram avaliadas por quadratura de Gauss-Legendre (5 pontos-pesos). A Figura 6 contém os gráficos da solução numérica nos três pontos da malha, computados no intervalo $0\mathrm{s}\le t\le 200\mathrm{s}$. Embora a função cresça ao longo do tempo, a solução numérica se mantém estável, sem dissipação ou dispersão, conforme gráfico constante na Figura 6d, ampliação das Figuras 6a, 6b e 6c.

Figura 6:
Solução analítica versus solução numérica da equação linear da onda unidimensional não homogênea, com termo fonte dado pela equação (3.4) e condições de contorno do tipo Neumann, avaliada em pontos fixos do espaço e velocidade de onda $c=4\mathrm{m}/\mathrm{s}$.

Os erros relativo e absoluto calculados nas normas L ₁ e L ₂ no tempo fixo em $t=200$ s, constantes na Tabela 2, podem ser justificados devido ao acúmulo dos erros de arredondamento.

Tabela 2:
Erros relativo e absoluto da solução numérica avaliada no tempo $t=200$ s para o problema não homogêneo com condições de contorno de Neumann.

Relativo (L 2 )	Relativo (L 1 )	Absoluto (L 2 )	Absoluto (L 1 )
4,163428 × 10-13	3,405960 × 10-13	1,238718 × 10-10	1,755041 × 10-10

No exemplo a seguir exploraremos a presença do erro de arredondamento computacional através de um problema não homogêneo com condições de Dirichlet.

Exemplo 3.3 (Problema não homogêneo com condições de Dirichlet) Tomando a equação (1.1) com termo fonte $f\left(x,t\right)=2\left(\alpha -\beta c^2\right)$ e prescrita por:

$u\left(x,0\right)=\mathrm{sen}\left(\lambda x\right)+\beta x\left(x-L\right),u_t\left(x,0\right)=0\text{e}u\left(0,t\right)=u\left(L,t\right)=\alpha t^2,$ (3.5)

obtemos a solução

$u\left(x,t\right)=\mathrm{sen}\left(\lambda x\right)\cos \left(\omega t\right)+\beta x\left(x-L\right)+\alpha t^2,\text{com}\lambda =2\pi /L\text{e}\omega =2\pi c/L.$ (3.6)

Os erros absolutos e relativos computados na norma L ₂ estão dispostos nas Tabelas 3 e 4. Neste experimento, a solução numérica foi iterada q vezes até o tempo 20.000 s sobre variados comprimentos de domínio L e velocidades de propagação de onda c. Utilizamos uma malha espacial com $n+1=3$ pontos (apenas os valores do ponto central são desconhecidos). Na Tabela 3, dispomos os resultados obtidos quando utilizamos incrementos temporais Δx/c com representação binária exata. Por exemplo, para $\Delta t=5/20$ (quando $\Delta x=5$), a representação na base 2 será 0,01₂. Neste caso, o erro computado é nulo independentemente do α adotado nas simulações (ver Tabela 3 $\alpha =0$ e $\alpha =-5$).

Tabela 3:
Erros relativo e absoluto na norma L ₂ da solução numérica avaliada no tempo $t=20.000\mathrm{s}$, para o problema da barra com condições de contorno de Dirichlet, com $\beta =3$ e incrementos temporais na forma de dízimas binárias não periódicas.

α	L	Δx	q	c	Relativo (L 2 )	Absoluto (L 2 )
-5 ou 0	2	1	20.001	1	0,00000 × 100	0,00000 × 100
	10	5	80.001	20	0,00000 × 100	0,00000 × 100
	100	50	801	2	0,00000 × 100	0,00000 × 100
	1000	500	1.601	40	0,00000 × 100	0,00000 × 100

Tabela 4:
Erros relativo e absoluto na norma L ₂ da solução numérica avaliada no tempo $t=20.000\mathrm{s}$, para o problema da barra com condições de contorno de Dirichlet, com $\beta =3$ e incrementos temporais na forma de dízimas binárias periódicas.

α	L	Δx	q	c	Relativo (L 2 )	Absoluto (L 2 )
-5	2	1	6.001	0,3	2,75302 × 10-16	9,53674 × 10-7
	10	5	76.001	19	4,81779 × 10-16	3,73184 × 10-6
	100	50	441	1,1	4,12953 × 10-16	1,01152 × 10-5
	1000	500	1.421	35,5	1,37634 × 10-15	1,06624 × 10-4
0	2	1	6.001	0,3	8,16431 × 10-17	2,44929 × 10-16
	10	5	76.001	19	3,26572 × 10-18	5,47679 × 10-16
	100	50	441	1,1	8,57665 × 10-20	4,54846 × 10-15
	1000	500	1.421	35,5	3,26572 × 10-22	5,47679 × 10-15

Não é o que ocorre, por outro lado, quando os incrementos temporais formam dízimas periódicas binárias, por exemplo, para $c=0,3$ e $\Delta t=1/0,3=3,{\overline{3}}_{10}$ (decimal), cuja representação binária é $11,{\overline{01}}_2$ (ver Tabela 4). Persiste, nestes casos, a presença de erro numérico atribuído ao arredondamento computacional. Quando a solução (3.6) se torna ilimitada (tomando $\alpha =-5$, por exemplo), o erro absoluto se torna mais sensível aos erros de arredondamento.

Exemplo 3.4 (Problema com descontinuidade nos dados iniciais) Considere agora o problema, também homogêneo, que modela uma corda inicialmente em repouso percutida por um martelo plano de largura 2δ. Os dados de contorno (Dirichlet) para este problema são todos nulos e as condições iniciais

$u\left(x,0\right)=0\mathrm{e}{u}_t\left(x,0\right)=\left\{\begin{array}{l}\nu ,\mathrm{se}\left|x-\xi \right|\le \delta ,\\ 0,\mathrm{se}\left|x-\xi \right|>\delta ,\end{array}\right.$ (3.7)

onde $0<\xi <L$.

O objetivo deste problema é avaliar a presença de descontinuidades a partir dos dados iniciais. Os gráficos nas Figuras 7a e 7c apresentam o comportamento da solução numérica ao longo do tempo nos três pontos interiores da malha ($n+1=5$), com incremento espacial $\Delta x=0,5$ e velocidade de onda $c=3\mathrm{m}/\mathrm{s}$. Estas figuras dispõem, também, os gráficos referentes à solução espectral avaliadas nos respectivos pontos.

Figura 7:
(a)-(c) Solução espectral versus solução numérica do problema da corda percutida por um martelo plano, avaliadas nos 3 pontos interiores da malha, com parâmetros $\xi =1,0\mathrm{m}$, $\delta =0,5\mathrm{m}$ e $\nu =2,5\mathrm{m}/\mathrm{s}$; (b)-(d) Ampliações das Figuras 7a e 7c, respectivamente.

A solução espectral, obtida por série de Fourier ⁵, foi truncada em $N=4$ termos a partir do critério $\left|{\hat{u}}_k-{\hat{u}}_{k-1}\right|<{10}^{-50}$, em que ${\hat{u}}_k$ é o coeficiente de Fourier. As Figuras 7b e 7d representam as ampliações das Figuras 7a e 7c, respectivamente. É evidente, a partir destas ampliações, a fragilidade da solução em série nas regiões de inexistência de derivadas (as “pontas” ou “quinas” dos gráficos): a solução obtida pelo esquema proposto não sofre influência nestes pontos.

Consideremos o gráfico constante na Figura 8a que apresenta as soluções numérica e espectral avaliadas ao longo do domínio espacial, $a\le x\le b$, com t fixo em 1 s. Observamos que existem dois pontos com rápida mudança caracterizando a inexistência de derivadas de primeira ordem. Ampliando-se a vizinhança desses pontos, Figuras 8b e 8c, observa-se que a solução espectral se torna oscilatória. Importante também é notar a constância produzida pela solução numérica que não oscila na região desses pontos.

Figura 8:
(a) Solução espectral versus solução numérica do problema da corda percutida por um martelo plano, avaliada no domínio espacial fixado o tempo em $t=1$ s; (b)-(c) Ampliações da Figura 8 nos pontos com rápida mudança e presença de oscilação da solução espectral.

Na Tabela 5, dispomos os desvios relativos e absolutos da solução numérica em relação à série truncada, calculados ambos nas normas L ₁ e L ₂ no tempo^‡ $t=29,\overline{6}\mathrm{s}$. Enquanto a solução numérica foi mantida fixa, com $n+1=5$ pontos, a solução em série foi variada no truncamento (N). De modo análogo ao exemplo 3.1, ocorre que quanto mais termos são considerados na solução em série, menor se tornam os desvios. Salientamos, mais uma vez, que o único parâmetro que está sendo variado é a quantidade de termos na solução em série. Isto corrobora com a afirmação de que os desvios computados se tratam de erros de truncamento da série e de arredondamento computacional.

Tabela 5:
Comparação entre as soluções numérica (com parâmetros fixos) e espectral (truncada em diversos termos N), avaliadas no tempo $t=29,\overline{6}\mathrm{s}$, para o problema da corda percutida por um martelo.

N	Relativo (L 2 )	Relativo (L 1 )	Absoluto (L 2 )	Absoluto (L 1 )
4	8,69334 × 10-2	7,71871 × 10-2	4,17778 × 10-2	4,53509 × 10-2
32	9,82114 × 10-3	8,33600 × 10-3	4,97375 × 10-3	5,17032 × 10-3
512	6,44142 × 10-4	5,27383 × 10-4	3,28540 × 10-4	3,29442 × 10-4
4.096	8,07554 × 10-5	6,59592 × 10-5	4,12076 × 10-5	4,12218 × 10-5
65.536	5,04921 × 10-6	4,12275 × 10-6	2,57665 × 10-6	2,57671 × 10-6
262.144	1,26233 × 10-6	1,03069 × 10-6	6,44178 × 10-7	6,44181 × 10-7
524.288	6,31165 × 10-7	5,15346 × 10-7	3,22090 × 10-7	3,22091 × 10-7

4 CONCLUSÃO

Este artigo abordou o desenvolvimento de um esquema numérico semi-analítico destinado à solução da equação da onda unidimensional não homogênea. Uma malha espaço-tempo foi construída através da relação incremental $c\mathrm{\Delta }t=\mathrm{\Delta }x$, de tal modo que volumes de controle foram formados a partir das retas características da equação. O desenvolvimento se deu a partir da forma integral da lei de conservação sobre estes volumes de controle. O esquema deduzido possui a propriedade de erro de truncamento local nulo, mesmo nos casos não homogêneo ou de solução generalizada (visto que satisfaz a lei integral). O método proporciona facilidade na inclusão das condições iniciais e de contorno, inseridas também sem necessidade de técnicas de aproximação.

As equações de avanço no tempo deduzidas contemplam expressões integrais que podem exigir aproximações por quadratura (como o exemplo3.2). Neste caso, podem ocorrer erros devido à fórmula de quadratura e aos arredondamentos. Nos experimentos em que se dispunha da solução em série (exemplos 3.1 e 3.2), observamos que o erro calculado é atribuído, principalmente, ao truncamento da série. Erro nulo é possível de ser obtido, desde que as integrais presentes nos esquemas possam ser avaliadas analiticamente (neste caso o erro de truncamento é exatamente nulo), bem como erros de arredondamento possam ser totalmente evitados (ver Tabela 3 do exemplo 3.3).

O problema de incluir aspectos que envolvam as curvas características é a dificuldade de generalização à maiores dimensões, conforme já salientava ^{2), (1), (3}.

Diante dos desenvolvimentos e dos experimentos numéricos realizados, todos sobre malhas que podem ser consideradas grosseiras, concluímos que o esquema desenvolvido é uma excelente ferramenta numérica, com ótima acurácia e robustez para resolução do problema de onda linear unidimensional.

AGRADECIMENTOS

O presente trabalho foi realizado com o apoio da CAPES. Os autores agradecem ao Programa de Pós Graduação em Métodos Numéricos em Engenharia, da Universidade Federal do Paraná (PPGMNE-UFPR), e ao Instituto Federal de Educação, Ciência e Tecnologia Catarinense, Campus Araquari, pelo apoio à pesquisa.

REFERÊNCIAS

• ¹
  S.C. Chang. The Method of Space-Time Conservation Element and Solution Elemen - A New Approach for Solving the Navier-Stokes and Euler Equations. Journal of Computational Physics, 119(2) (1995), 295-324. doi:10.1006/jcph.1995.1137. URL http://www.sciencedirect.com/science/article/pii/S0021999185711370
  » https://doi.org/10.1006/jcph.1995.1137» http://www.sciencedirect.com/science/article/pii/S0021999185711370
• ²
  S.C. Chang & W.W. To. A New Numerical Framework for Solving Conservation Laws - The Method of Space-Time Conservation Element and Solution Element. Technical Memo TM 104495, NASA, Lewis Research Center, NASA (1991). URL http://www.grc.nasa.gov/WWW/microbus.TM 104495
  » http://www.grc.nasa.gov/WWW/microbus.TM 104495
• ³
  S.C. Chang, X.Y. Wang & W.M. To. Application of the Space-Time Conservation Element and Solution Element Method to One-Dimensional Convection-Diffusion Problems. Journal of Computational Physics, 165(1) (2000), 189-215. doi:10.1006/jcph.2000.6610. URL http://www.sciencedirect.com/science/article/pii/S0021999100966105
  » https://doi.org/10.1006/jcph.2000.6610» http://www.sciencedirect.com/science/article/pii/S0021999100966105
• ⁴
  L.L. Fernandes, J.a.C.R. Cruz, C.J.C. Blanco & A.R.B. Barp. Modelagem Sísmica via Métodos das Diferenças Finitas: Caso da Bacia do Amazonas. Acta Amazonica, 39 (2009), 155 - 163. URL http://www.scielo.br/scielo.php?script=sci_arttext&pid=S0044-59672009000100015&nrm=iso
  » http://www.scielo.br/scielo.php?script=sci_arttext&pid=S0044-59672009000100015&nrm=iso
• ⁵
  D.G.d. Figueiredo. “Análise de Fourier e equações diferenciais parciais”. IMPA, Rio de Janeiro (2012), 274 pp.
• ⁶
  M.D. Greenberg. “Advanced Engineering mathematics”. Prentice Hall, New Jersey (1998), 1324 pp.
• ⁷
  A. Iserles. “A First Course in the Numerical Analysis of Differential Equations”. Cambridge University Press, New York (1996), 378 pp.
• ⁸
  R.J. Leveque. “Numeriacal methods for conservation laws”. Birkhäauser, Basel; Boston; Berlin (1992), 214 pp.
• ⁹
  E.L. Lima. “Curso de Análise vol. 2”. IMPA, Rio de Janeiro (2011), 546 pp.
• ¹⁰
  A.D. Polyanin. “Handbook of linear partial differential equations for engineers and scientists”. Chapman & Hall/CRC (2002), 667 pp.
• ¹¹
  G. Strang. “Introduction to Applied Mathematics”. Wellesley-Cambridge Press (1986).
• ¹²
  J.C. Strikwerda. “Finite Difference Schemes and Partial Differential Equations”. Siam, Philadelphia (1947), 435 pp.

• *
  Observe que a condição de estabilidade para este esquema é $\left|c\Delta t/\Delta x\right|\le 1$ ^8),(12),(7.
• †
  Esta solução foi obtida pelo método da solução fabricada, em que o termo-fonte é construído de tal forma que a função-solução satisfaça a EDP. Assim, a função-solução determina o termo fonte (3.4). Observação análoga pode ser feita para a solução em (3.6).
• ‡
  A partir da Figura 7b é possível observar que este é um dos tempos em que a curva da solução em série é distinta da numérica, esta observação motivou a escolha do tempo $t=29,\overline{6}\mathrm{s}$ para análise do erro.

Disponibilidade de dados

Citações de dados

S.C. Chang & W.W. To. A New Numerical Framework for Solving Conservation Laws - The Method of Space-Time Conservation Element and Solution Element. Technical Memo TM 104495, NASA, Lewis Research Center, NASA (1991). URL http://www.grc.nasa.gov/WWW/microbus.TM 104495

Datas de Publicação

• Publicação nesta coleção
  03 Jun 2019
• Data do Fascículo
  Jan-Apr 2019

Histórico

• Recebido
  13 Nov 2017
• Aceito
  22 Out 2018