| | | | |

6.3 Equação da Onda

Ajude a manter o site livre, gratuito e sem propagandas. Colabore!

Consideramos a equação da onda com condições iniciais dadas e condições de contorno de Dirichlet homogêneas

uttαuxx=0, 0<ttf,a<x<b, (6.39a)
u(0,x)=f(x),axb, (6.39b)
ut(0,x)=g(x),axb, (6.39c)
u(t,a)=u(t,b)=0, 0ttf, (6.39d)

onde u=u(t,x) é a incógnita com f, g e α>0 dadas.

Para a aplicação do Método das Diferenças Finitas (MDF), assumimos as discretizações: no tempo, t(k)=kht, j=0,1,,nt, ht=tf/nt; no espaço xi=a+(i1)hx, i=1,2,,nx+1, hx=(ba)/nx. Então, assumindo a notação ui(k)u(t(k),xi) usando a fórmula de diferenças finitas central D0,h22, obtemos a seguinte forma discreta da equação Eq. (6.39a)

ui(k1)2ui(k)+ui(k+1)ht2 (6.40)
αui(k)2ui(k)+ui+1(k)hx2=0,

para i=2,3,,nx, j=1,2,,nt1. Denotando λ:=αht2/hx2, rearranjando os termos e aplicando as condições de contorno, obtemos

u2(k+1)=2(1λ)u2(k)+λu3(k)u2(k1), (6.41a)
ui(k+1)=λui(k)+2(1λ)ui(k)+λui+1(k)ui(k1), (6.41b)
unx(k+1)=λunx1(k)+2(1λ)unx(k)unx(k1), (6.41c)

para i=2,3,,nx, j=1,2,,nt1. Ou, equivalentemente, na forma matricial

𝒖~(k+1)=A𝒖~(k)𝒖~(k1), (6.42)

para k=1,2,,nt1, onde 𝒖~(k)=(ui(k))i=2nx e A=[ai,j]i,j=1nx1,nx1 é a matriz tridiagonal de elementos

ai,j={ai,i1=λ,1<inx1,ai,i=2(1λ),1inx1,ai,i+1=λ,1i<nx1. (6.43)

Para a inicialização, a Eq. (6.42) requer que conhecemos 𝒖~(0) e 𝒖~(1). A primeira, vem diretamente da condição inicial Eq. (6.39b), i.e.

𝒖~(0)=f(𝒙~), (6.44)

onde 𝒙~=(xi)i=2nx. Agora, aplicando a fórmula de diferenças finitas progressiva D+,h, temos da condição inicial Eq. (6.39c)

𝒖~(1)𝒖~(0)ht=g(𝒙~) (6.45)

ou, equivalentemente,

𝒖~(1)=𝒖~(0)+htg(𝒙~). (6.46)

De tudo isso, temos que a solução numérica da equação da onda pode ser computada com a seguinte iteração

𝒖~(0)=f(𝒙~), (6.47a)
𝒖~(1)=𝒖~(0)+htg(𝒙~), (6.47b)
𝒖~(k+1)=A𝒖~(k)𝒖~(k1), (6.47c)

para k=1,2,,nt1, com 𝒖(k)=(0,𝒖~,0).

Observação 6.3.1.

(Condição de Estabilidade e Erro de Truncamento.) Pode-se mostrar a seguinte condição de estabilidade [3, p. 487]

αhthx1. (6.48)

Com isso e para f e g suficientemente suaves, o esquema numérica (6.3) tem erro de truncamento O(ht2+hx2).

Exemplo 6.3.1.

Consideramos o seguinte problema

uttuxx=0, 0<t2, 0<x<1, (6.49a)
u(0,x)=0, 0x1, (6.49b)
ut(0,x)=πsen(πx), 0x1, (6.49c)
u(t,0)=u(t,1)=0, 0t2. (6.49d)

Sua solução exata é u(t,x)=sen(πt)sen(πx). A Figura 6.6 contém gráficos de comparação entra as soluções numérica e exata. Para a solução numérica, tomamos nt=40 (ht=0.05) e nx=10 (hx=0.1).

Refer to caption
Refer to caption
Figura 6.6: Gráficos comparativos das soluções numérica e exata do problema de onda do Exemplo 6.3.1.
1import numpy as np
2from numpy import pi, sin, cos
3
4# params
5nt = 40
6ht = 2./nt
7tt = np.linspace(0., 2., nt+1)
8
9nx = 10
10hx = 1./nx
11xx = np.linspace(0., 1., nx+1)
12
13# c.i.s
14def f(x):
15    return np.zeros_like(x)
16
17def g(x):
18    return pi*sin(pi*x)
19
20# axiliares
21lbda = ht**2/hx**2
22
23A = np.zeros(((nx-1), (nx-1)))
24A[0,0] = 2*(1. - lbda)
25A[0,1] = lbda
26for i in range(1,nx-2):
27    A[i,i-1] = lbda
28    A[i,i] = 2*(1 - lbda)
29    A[i,i+1] = lbda
30A[nx-2,nx-3] = lbda
31A[nx-2,nx-2] = 2*(1 - lbda)
32
33# laço no tempo
34## c.i.s
35u0 = f(xx)
36
37u1 = u0.copy()
38u1[1:-1] = u0[1:-1] + ht*g(xx[1:-1])
39
40u = u1.copy()
41for k in range(1,nt):
42
43    print(f'{k+1}: t = {tt[k+1]:f}')
44
45    u[1:-1] = A@u1[1:-1] - u0[1:-1]
46
47    u0 = u1.copy()
48    u1 = u.copy()

6.3.1 Exercício

E. 6.3.1.

Considere o problema

uttuxx=0, 0<t1.5, 0<x<1, (6.50a)
u(0,x)=0, 0x1, (6.50b)
ut(0,x)=πsen(πx), 0x1, (6.50c)
u(t,0)=u(t,1)=0, 0t1.5. (6.50d)

Sua solução exata é u(t,x)=sen(πt)sen(πx). Faça testes numéricos para determinar os passos ht e hx para os quais o esquema numérico (6.3) compute o valor de u(1.5,0.5) com 5 dígitos significativos corretos.

Resposta.

ht=2.5e3, hx=1.e2

E. 6.3.2.

Considere o problema

uttuxx=et(2+xx2), 0<t1, 0<x<1, (6.51a)
u(0,x)=xx2, 0x1, (6.51b)
ut(0,x)=x2x, 0x1, (6.51c)
u(t,0)=u(t,1)=0, 0t1. (6.51d)

Sua solução exata é u(t,x)=et(xx2). Implemente um esquema numérico semelhante ao (6.3) para computar soluções numéricas desse problema.

E. 6.3.3.

Considere o problema

uttuxx=et(2+xx2), 0<t1, 0<x<1, (6.52a)
u(0,x)=xx2, 0x1, (6.52b)
ut(0,x)=x2x, 0x1, (6.52c)
ux(t,0)=et, 0t1, (6.52d)
u(t,1)=0, 0t1. (6.52e)

Sua solução exata é u(t,x)=et(xx2). Implemente um esquema numérico semelhante ao (6.3) para computar soluções numéricas desse problema.

E. 6.3.4.

Considere o problema

uttuxx=0, 0<t2, 0<x<1, (6.53a)
u(0,x)=0, 0x1, (6.53b)
ut(0,x)=πsen(πx), 0x1, (6.53c)
u(t,0)=u(t,1)=0, 0t2. (6.53d)

Sua solução exata é u(t,x)=sen(πt)sen(πx). Baseado em (6.3), desenvolva um novo esquema numérico substituindo o passo (6.47b) por um esquema numérico de mais alta ordem.

Resposta.

Dica: use, por exemplo, um método de R-K-2.

E. 6.3.5.

Considere o problema

uttαuxx=0, 0<t1, 0<x<1, (6.54a)
u(0,x)=x(1x), 0x1, (6.54b)
ut(0,x)=0, 0x1, (6.54c)
u(t,0)=u(t,1)=0, 0t2. (6.54d)

Use o esquema numérico (6.3) para fazer testes numéricos para α=1.,0.5,0.1,0.01. É necessário ajustar os parâmetros ht e hx ao variar o parâmetro α? Justifique sua resposta.

Resposta.

Dica: consulte a Observação 6.3.1.


Envie seu comentário

As informações preenchidas são enviadas por e-mail para o desenvolvedor do site e tratadas de forma privada. Consulte a Política de Uso de Dados para mais informações. Aproveito para agradecer a todas/os que de forma assídua ou esporádica contribuem enviando correções, sugestões e críticas!