Física Matemática I: aula do dia 22/11/2016

Davi C. Rodrigues

EDP separável (sessão 6.5 do Cahill) (continuação)

Equação de onda bidimensional e o problema do tambor: funções de Bessel e Neumann

Apresentação do problema

Voltemos à equação abaixo:

(∇^2)ϕ (x,t) - (1/(v^2))(Overscript[ϕ, ..])(x,t) =0

Note que a constante real v precisa ter dimensão de velocidade.

Note ainda que esta equação é invariante por inversão temporal e por inversão espacial. Ademais, o sinal relativo é muito importante, pois caso fosse o contrário, não seria possível encontrar soluções oscilatórias simultanemanete para o tempo e para o espaço (requisito essencial para descrever, por exemplo, uma corda que oscila).

Não falei antes, mas essa mesma equação, para v = 1, pode ser escrita na forma

□ϕ (x,t) =0,

Em que o operador d’Alembertiano é definido por □ = (∇^2) - (∂_t)(∂_t).

Contudo, como antes dito, v tem dimensão de velocidade, e como pode algo medido, por exemplo, em metros por segundo valer 1? Note que 1 é diferente de 1 m/s.  Se bem me lembro, já comentei isso com vocês em sala de aula, mas se não souberem entrem em contato.

Aqui trataremos de uma abordagem razoavelmente geral à equação de onda em duas dimensões, mas é bom ter um problema físico para nos nortear. Para isso consideraremos o chamado problema do tambor, ou seja, o problema de uma película bidimensional capaz de vibrar mas que encontra-se presa em suas bordas. Em outras palavras, é o análogo bidimensional do problema da corda que pode vibrar e está presa em suas extremidades.

Seja

ϕ (x,t) =ϕ (r,φ,t) = R(r) Θ(θ) T(t).

A eq.(1) pode ser expressa por

(1/r)(∂_r) (r (∂_r)ϕ ) + (1/(r^2))((∂_θ)^2) ϕ - (1/(v^2))((∂_t)^2)ϕ = 0,

logo

(1/(r R ))(r  R' )' + (1/((r^2) Θ)) Θ'' - (1/((v^2) T))T'' = 0.

Resolução da parte temporal

Como a última parcela da última equação só depende de t, e as duas primeiras de θ e r, a única forma dessas parcelas serem iguas para quaisquer valores dessas variáveis é que

(1/((v^2) T))T'' =  - (α^2),

em que α é certa constante real. A EDO acima é a mesma que encontramos para a parte temporal do problema unidimensional. Portanto já sabemos que, dado que buscamos soluções oscilatórias, a constante que aparece na parte direita é necessariamente negativa. Por isto, já usamos uma notação que força que o lado direito seja sempre negativo para quaquer α real. Como visto, a solução para T pode ser escrita como

T(t) = A Cos(α v t + δ),

em que A e δ são constantes de integração.

Resolução da parte angular

Combinando as eqs. (5) e (6) temos,

(r/( R ))(r  R' )'  +(r^2)(α^2) = -(1/Θ) Θ''.

Acima temos uma equação cujo lado esquerdo só depende de r e o direito só depende de θ, e portanto deve existir constante real β tal que

(1/Θ) Θ'' = β.

A função Θ precisa ser periódica, em particular Θ(0) = Θ(2π) (caso contrário, ϕ(r,0,t) ≠ ϕ(r,2π,t), o que é um absurdo para o que pretendemos descrever fisicamente). Se β for positivo, Θ terá crescimento exponencial, o que impede a periodicidade de 2π. Se β for negativo, a solução geral para Θ é

Θ = B Cos((Sqrt[-β]) θ + (OverTilde[δ])).

Logo, para obter a periodicidade de 2π,

(Sqrt[-β]) = m, em que m  ∈   ℕ. (trocar para inteiro)

Assim,

Θ = B Cos(m θ + (OverTilde[δ])).

Resolução da parte radial

A eq. (8) pode agora ser escrita da seguinte forma

(r/( R ))(r  R' )'  +(r^2)(α^2) = - β =  (m^2).

Logo,

r (r  R' )'  +(r^2)(α^2) R -  (m^2) R = 0,

R'' + (1/r)R'+ ((α^2) - ((m^2)/(r^2)))R = 0.

No caso da corda, ao resolvermos a parte radial, encontramos soluções que demandavam o uso de funções especiais, contudo eram funções especiais bem conhecidas: seno ou cosseno. A equação diferencial acima também tem solução em termos de funções especiais, mas veremos que essas são chamadas de funções de Bessel e de Neumann. A eq. diferencial acima é essencialmente (a menos da manipulação a seguir) a eq. diferencial de Bessel.

Tal como pode ser visto em diversas fontes (por exemplo, http://mathworld.wolfram.com/BesselDifferentialEquation.html), a eq. de Bessel é

Assim, para obter a forma acima, falta uma simples mudança de coordenadas, a saber, a eq. (14) pode também ser escrita como,

(1/(α^2))R'' (r)+ (1/((α^2)r))R'(r)+ (1 - ((m^2)/((α^2)(r^2))))R(r) = 0.

E portanto,

R''(x) + (1/x)R'(x)+ (1 - ((m^2)/(x^2)))R(x) = 0,

em que x=α r, e assim temos uma equação que é idêntica à eq. de Bessel.

É natural que, com o passar do tempo, caso vocês não usem a eq. de Bessel com frequência, venham a se esquecer da forma precisa dela. Entretanto, pode-se lembrar que essa equação é obtida da parte radial da equação da onda bidimensional.

As soluções da eq. (16) são em geral dadas pela combinação linear de duas funções LI (como seria natural de se esperar), e assim com duas constantes de integração (para m fixo). Essas funções são as funções de Bessel de Primeiro tipo (ou função de Bessel simplesmente), e as funções de Bessel de segundo tipo (ou funções de Neumann).

Para qualquer que seja o valor de m a solução é

DSolve[{f''[x] + (1/x) f'[x] + (1-(( (m^2))/(x^2)))f[x] == 0}, f[x], x]

Note que a solução acima é a mais geral para m fixo. Note que para obter essa solução com o Mathematica nada foi dito com respeito à natureza de m , realmente, a solução acima é válida para qualquer m, seja inteiro, real ou complexo.  Contudo, considerando o problema do tambor que temos em mente, m deve ser inteiro. Vimos que da análise da parte angular, descobre-se que m é inteiro, mas seu valor não é fixado. Ou seja, a solução de interesse para a resolução do problema da equação de onda bidimensional com simetria axial depende da soma de 2 coeficientes arbitrários para cada m.

Funções de Bessel e de Neumann

Acima, BesselJ[m,k r] representa a função de Bessel do primeiro tipo, que é comumente denotada por (Subscript[J, m])(k r). BesselY[m,k r] representa a função de Bessel do segundo tipo, também chamada de função de Neumann, e é comumente denotada por (Subscript[N, m])(k r) ou (Subscript[Y, m])(k r).

Sendo m inteiro, mostra-se que

Abaixo seguem exemplos gráficos das funções de Bessel.

Row[{Plot[{BesselJ[0,x], BesselJ[1,x], BesselJ[2,x], BesselJ[3,x]}, {x,-10,10}, PlotTheme->"Detailed", ImageSize->Medium],
Plot[{BesselJ[0,x], BesselJ[1,x], BesselJ[2,x],BesselJ[3,x]}, {x,-50,50}, PlotRange -> All,PlotTheme->"Detailed",ImageSize->Medium]}]

Exemplos gráficos de funções de Neumann:

Plot[{BesselY[0,x], BesselY[1,x], BesselY[2,x], BesselY[3,x]}, {x,-1,30},PlotTheme->"Detailed"]

Ambas são funções que oscilam em torno de zero, mas não são periódicas como senos ou cossenos, pode-se ver que a amplitude das oscilações decai com x.

Uma característica essencial que difere as funções de Bessel das de Neumann é a singularidade sempre presente nas últimas em x =0. 

Em x=0, todas as funções de Bessel são nulas, exceto pela (Subscript[J, 0])(x).

Esses limites são apresentados abaixo.

Limit[BesselJ[0,x], x-> 0, Direction -> -1]
Limit[BesselY[0,x], x-> 0, Direction -> -1]
Limit[BesselJ[n,x], x-> 0, Direction -> -1, Assumptions -> n>0]
Limit[BesselY[n,x], x-> 0, Direction -> -1, Assumptions -> n>  0]

No infinito, todas essas funções se anulam,

Limit[BesselJ[n,x], x-> ∞, Assumptions -> n≥ 0]
Limit[BesselY[n,x], x-> ∞, Assumptions -> n≥  0]

Devido ao comportamento singular das funções de Neumann, e dada nossa motivação física de descrever uma película que oscila, todos os coeficientes que multiplicam funções de Neumann tem de ser necessariamente nulos.

Assim, até o presente nomento nossa solução para um α fixo é

ϕ(r,ϕ, t)= (∑_(m=0)^∞)(Subscript[A, m α]) cos(α v t) cos(mθ) (Subscript[J, m])(α r),

em que fixamos as constantes δ e (OverTilde[δ]) como nulas (sem perda de generalidade -- sabe por quê?) e (Subscript[A, m α]) são constantes, que dependem de m e de α. Até o momento, α é um parâmetro real qualquer, mas falta impor mais uma condição. No caso do problema do tambor, não importa o comportamento das funções de Bessel no infinito, mas sim num raio fixo, digamos que o raio da película do tambor seja L, então temos a condição de contorno

ϕ(L,ϕ, t)= 0.

A fim de satisfazer a condição acima, teremos de buscar pelos valores de α tais que (Subscript[J, m])(α r) = 0.

No caso da corda, resolvemos problema análogo, mas lá tínhamos de resolver a equação sen(k L) = 0, que levou naturalmente a k = ((n π )/L). Como encontrar os zeros das funções de Bessel?

Este não é um problema passível de ser resolvido analiticamente, mas pode ser resolvido por métodos numéricos. Antigamente, para saber os zeros da função de Bessel buscava-se comumente por um livro que contivesse esses zeros tabelados. 


Com os softwares de hoje, é mais simples buscar diretamente pelas raízes numericamente, veja por exemplo,

FindRoot[BesselJ[0,x] == 0, {x, 1}] (* Procura por uma raiz usando x == 1 como ponto de partida*)

FindRoot[BesselJ[0,x] == 0, {x, 4} ] 

FindRoot[BesselJ[0,x] == 0, {x, 8}]

Compare os valores acima com o gráfico de (Subscript[J, 0])(x).

De forma mais automática e robusta, o Mathematica dispõe de uma função dedicada para os zeros da função de Bessel:

N[BesselJZero[0,1]]
N[BesselJZero[0,2]] (* São a primeira, a segunda e a terceira raizes de BesselJ[0,x]*)
N[BesselJZero[0,3]] (*Usa-se N aqui para forçar o Mathematica a calcular numericamente, do contrário ele simplesmente repete o símbolo BesselJZero[0,1]*)

É fácil fazer uma tabela com as 10 primeiras raízes de (Subscript[J, 0]), (Subscript[J, 1]), (Subscript[J, 2]), (Subscript[J, 3]), por exemplo. A saber:

N[Table[BesselJZero[i,j], {j,1,10}, {i,0,30}]] // TableForm

Sendo (Subscript[k, n,p]) o p-ésimo zero da função de Bessel (Subscript[J, n])(x), ou seja, (Subscript[J, n])(kn, p) = 0 para todo p, então vemos que os possíveis valores de α são

α = ((Subscript[k, n, p])/L).

E a solução geral do problema do tambor é

ϕ(r,θ, t)= (∑_(m=0)^∞)(∑_(p=1)^∞)(Subscript[A, m p]) cos(((Subscript[k, m,p])/L)v t) cos(mθ) (Subscript[J, m])(((Subscript[k, m,p])/L) r).

Tal como ocorreu com o caso da corda, só dizer que as oscilações são nulas nas extreminadades não é suficiente para determinar todas as constantes do problema. Isto é, note que todos os coeficientes (Subscript[A, m p]) são arbitrários. É bom lembrar que v é constante que depende do meio, mas é uma constante dada (aparece explicitamente na equação diferencial original). Fisicamente, é natural que os coeficientes (Subscript[A, m p]) continuem arbitrários, pois a película de um mesmo tambor pode vibrar de diferentes formas. Isto é, a solução acima impõe limites às soluções possíveis, mas não a fixa a ponto de fixar todas as constantes.

Caso particular: (Subscript[A, 01]) = 1 e demais coeficientes nulos. Ademais, sem perda de generalidade, podemos re-escalar a coordenada radial e o tempo t tais que possamos usar v = 1 e L =1. --- Isto está claro para você?

Assim, para esse caso particular,

Para uma questão de velocidade para gerar os gráficos, é bom calcular e salvar em variáveis os seguintes valores:

k01 = N[BesselJZero[0,1]];
k02 = N[BesselJZero[0,2]];
k11 = N[BesselJZero[1,1]];
k21 = N[BesselJZero[2,1]];
k22 = N[BesselJZero[2,2]];

Vamos chamar a solução na qual  (Subscript[A, 01])=1 e os demais são nulos de ϕ01, assim

ϕ01[r_,θ_,t_] = Cos[k01 t]  BesselJ[0,(k01 r)]

Ou seja,

(Subscript[ϕ, 01])(r,θ, t)= cos((Subscript[k, 0,1]) t)  (Subscript[J, m])((Subscript[k, 0,1]) r).

Num primeiro momento, vamos fixar o tempo, e seja t=0. Assim,

RevolutionPlot3D[ϕ01[r,θ,0],{r,0,1},{θ,0,2π}, AxesLabel-> "  ϕ(r,θ,0)" ]

Agora podemos considerar a evolução no tempo...

Animate[ RevolutionPlot3D[ϕ01[r,θ,t],{r,0,1},{θ,0,2π}, AxesLabel-> "  ϕ(r,θ,0)" ,Mesh->True, PlotRange->{{-1,1}, {-1,1}, {-1,1}}], {t,0,((2π)/k01)},SaveDefinitions->True,AnimationRunning->False]

Agora o caso  (Subscript[A, 02]) = 1 e os demais nulos

ϕ02[r_,θ_,t_] = Cos[k02 t]  BesselJ[0,(k02 r)]

Animate[ RevolutionPlot3D[ϕ02[r,θ,t],{r,0,1},{θ,0,2π}, AxesLabel-> "  ϕ(r,θ,0)" ,Mesh->True, PlotRange->{{-1,1}, {-1,1}, {-1,1}}], {t,0,((2π)/k02)},SaveDefinitions->True,AnimationRunning->False]

Agora o caso  (Subscript[A, 11]) = 1 e os demais nulos.

ϕ11[r_,θ_,t_] = Cos[k11 t]  Cos[θ] BesselJ[1,(k11 r)]

Animate[ RevolutionPlot3D[ϕ11[r,θ,t],{r,0,1},{θ,0,2π}, AxesLabel-> "  ϕ(r,θ,0)" ,Mesh->True, PlotRange->{{-1,1}, {-1,1}, {-1,1}}], {t,0,((2π)/k11)},SaveDefinitions->True,AnimationRunning->False]

Esta acima é uma oscilação com paridade ímpar. Note que  depende exclusicamente de uma função de Bessel de paridade ímpar (todos os (Subscript[J, n])(x) com n ímpar são funções ímpares)

ϕ21[r_,θ_,t_] = Cos[k21 t]  Cos[2 θ] BesselJ[2,(k21 r)];
Animate[ RevolutionPlot3D[ϕ21[r,θ,t],{r,0,1},{θ,0,2π}, AxesLabel-> "  ϕ(r,θ,0)" ,Mesh->True, PlotRange->{{-1,1}, {-1,1}, {-1,1}}], {t,0,((2π)/k21)},SaveDefinitions->True,AnimationRunning->False]

Agora o caso  (Subscript[A, 22]) = 1 e os demais nulos.

ϕ22[r_,θ_,t_] = Cos[k22 t]  Cos[2 θ] BesselJ[2,(k22 r)];
Animate[ RevolutionPlot3D[ϕ22[r,θ,t],{r,0,1},{θ,0,2π}, AxesLabel-> "  ϕ(r,θ,0)" ,Mesh->True, PlotRange->{{-1,1}, {-1,1}, {-1,1}}], {t,0,((2π)/k22)},SaveDefinitions->True,AnimationRunning->False]