Um exercício interessante sobre EDP’s

Primeiramente, eu gostaria de disponibilizar o material da minha aula de 18/10, sobre  EDP’s lineares de segunda ordem e também alguns exercícios resolvidos sobre EDP’s.

Na aula de exercícios do dia 20/10 surgiu uma discussão muito interessante durante a resolução de um exercício, mais especificamente o item (k) do exercício (2) do capítulo 6 do livro-texto: resolver a EDP

y u_x + x u_y = xy ,

com as condições de contorno

u(0,y) = e^{-y^2} para y>0 e u(x,0) = e^{-x^2} para x>0.

A solução geral desta equação terá a forma (façam as contas!):

u(x,y) = \frac{x^2}{2} + f \left( x^2 - y^2 \right) ,

com a função f(u) a ser determinada pelas condições do problema. Usando a condição no contorno x = 0, obtemos a solução

u_1(x,y) = \frac{x^2}{2} + e^{-(y^2 - x^2)} .

Usando a condição em y = 0, obtemos uma outra solução, dada por

u_2(x,y) = \frac{y^2}{2} + e^{-(x^2 - y^2)} .

E agora? O que fazer? Bem, precisaremos “colar” essas duas soluções para formar a solução correta do problema.

Na solução que eu levei para a aula (que está no PDF de exercícios que disponibilizarei ao final da postagem), eu simplesmente assumi que a solução u_1 seria válida na região y > x, enquanto a solução u_2 seria válida na região x > y. Isto me pareceu razoável, visto a maneira como as condições de contorno foram definidas e a simetria entre as variáveis x e y na equação.

Entretanto, surgiu entre os alunos uma dúvida muito pertinente: por quê? Por que escolher justamente esta curva, a reta x=y para acoplar as soluções?

Notemos que a curva característica x^2-y^2 = a gera uma família de hipérboles. Poderíamos então escolher uma hipérbole qualquer para unir as soluções? O problema então deixaria de ter uma única solução?

hiperboles

A figura acima mostra uma família de hipérboles. Vamos escolher uma delas e dividir o plano x-y em duas regiões, de modo que as soluções u_1 e u_2 valham, respectivamente, na região acima e abaixo da curva escolhida. No caso da figura acima, a região 0 < x < a está englobada na parte de cima da curva, que corresponderia à solução gerada pela condição de contorno em y = 0; no entanto, a condição em x = 0 gera a solução u_2, que é diferente de u_1. Ou seja: a escolha de um a \neq 0 geraria uma solução inconsistente!!! Deste modo, a única escolha posível para a constante seria a = 0, de onde obteremos a curva x=y.

Uma terceira maneira de justificar esta escolha seria impor que a solução seja contínua na origem. Daí, obteríamos de maneira direta que f(0) = 1 em ambas as soluções. O que resultaria novamente na curva x^2 - y^2 = 0, isto é, x = y, uma vez que o domínio de nossa EDP é o quadrante x \ge 0, y \ge 0 de \mathbb{R} ^2 .

Deste modo, a solução deste problema é

u(x,y) = \begin{cases} \frac{x^2}{2} + e^{-(y^2 - x^2)} & x < y \\ e^{-(x^2 - y^2)} & x > y \end{cases} .

Abaixo, a solução u(x,y) vista a partir de dois pontos de vista. Note a simetria apresentada pela solução.

graf2graf1

Bom teste para todos! 🙂

Advertisements

A KdV linear e o caso dos “5 pontos extras”

Voltemos à questão 3 da prova, a KdV linear

\displaystyle \partial_t u + \partial_{xxx}u = 0 , \quad x\in \mathbb{R},\quad t\ge 0,

com condição inicial u(x,0)=f(x), sendo a transformada de Fourier F(k) de f(x) conhecida. A solução deste problema será (ver aqui)

\displaystyle u(x,t) = \frac{1}{\sqrt{2\pi}}\int_{-\infty}^\infty F(k) e^{-ik(x+k^2t)}dk

Os “5 pontos extras” comentados na prova correspondem a “resolver” essa integral para F(k)  arbitrário. Obviamente, a ideia aqui é usar o teorema de convolução, isto é, escrever a solução como

\displaystyle u(x,t) =\int_{-\infty}^\infty f(s) g_t(x-s)\, ds

sendo

\displaystyle g_t(x) =  {\cal F}^{-1}\left[e^{-ik^3t}\right]= \frac{1}{\sqrt{2\pi}}\int_{-\infty}^\infty   e^{-i(k^3t+kx)}dk.

Nossa equação é real e linear, portanto podemos sempre pegar as combinações lineares pertinentes de soluções a fim de se obter uma solução real. Isto nos permite escrever, sem perda de generalidade

\displaystyle g_t(x) =   \sqrt{\frac{2}{\pi}}\int_{0}^\infty   \cos (k^3t+kx)\, dk.

Façamos agora a mudança de variável k = w/\sqrt[3]{3t}, o que resulta (para t>0)

\displaystyle g_t(x) =    \frac{\sqrt{2}}{\sqrt{\pi}\sqrt[3]{3t}} \int_{0}^\infty   \cos \left(\frac{w^3}{3}+\frac{wx}{\sqrt[3]{3t}}\right)\, dw.

Essa integral é conhecida. Ela nos remete a chamada Função de Airy

\displaystyle {\rm Ai}(y) = \frac{1}{\pi }\int_{0}^\infty   \cos \left(\frac{w^3}{3}+ yw\right)\, dw,

cujas propriedades são muito bem conhecidas. É evidente que nossa solução será

\displaystyle g_t(x) =    \frac{\sqrt{2\pi}}{\sqrt[3]{3}}  t^{-\frac{1}{3}}{\rm Ai}\left(\frac{x}{\sqrt[3]{3t}} \right) ,

válido para todo t>0. O limite t\to 0 corresponde a uma delta, veja a convolução acima.

Ora, quer dizer que para ganhar os 5 pontos tínhamos que saber de cor a integral que dava a função de Airy. Sim. Sinto muito. C’est la vie. Eu sabia. Se eu sei, porque vocês não podem saber? Na verdade, eu não sabia de cor, mas tinha uma vaga lembrança que as funções de Airy  tinham transformadas de Fourier que iam com k^3. Tendo uma vaga lembrança, consulta-se as fontes pertinentes! Quando eu era estudante de Métodos Matemáticos (disciplinas chamada  Física Matemática em lugares com menos complexos), a referência fundamental eram as tabelonas russas.  Gradshteyn and Ryzhik e Abramowitz and Stegun marcaram indelevelmente todas as gerações anteriores a de vocês. Ainda há exemplares na biblioteca (não sei por quanto tempo ainda), consultem, é como fazer “arquelogia” dos métodos matemáticos… Mas vocês não precisam mais disso. Cliquem aqui, e vocês verão como eu me certifiquei que deveria aparecer uma função de Airy no problema. Essas ferramentas são muito uteis, não as menosprezem.

Equações diferenciais e de diferenças

Não vou conseguir até a prova fazer o resumo que havia pensado sobre aplicações da Transformada de Laplace  às equações de diferenças, sinto muito. Porém, vou fazer aqui o problema que me pediram hoje em sala. Farei o 13.a, tenho certeza que com a solução vocês conseguirão fazer o b. O 12 é mais fácil. Qualquer dúvida, escrevam aqui que respondo até a prova. Notem que o problema 8 está resolvido nas notas da Stefania. O problema 13.a é resolver a equação diferencial de diferenças:

y'(t) - y(t-1) = t^2, \quad y(t)=0, t\le 0.

Calculando-se a transformada de Laplace de ambos os lados, tem-se

\displaystyle  sY(s) - e^{-s}Y(s) = \frac{2}{s^3}

onde usou-se a propriedade do deslocamento. A transformada de Laplace será portanto

\displaystyle F(s) = \frac{2}{s^3}\times \frac{1}{s-e^{-s}}

Vamos calcular a inversa de duas maneiras: integrando-se no plano complexo e usando-se as propriedades elementares da TL.

1) Propriedades elementares da TL (estratégia recomendada para a prova)

Vamos escrever a transformada como

\displaystyle F(s) = \frac{2}{s^4}\times \frac{1}{1-\frac{e^{-s}}{s}}

Supondo-se \Re s> 1 e usando-se a série geométrica, teremos

\displaystyle F(s) = 2\sum_{k=0}^\infty \frac{e^{-ks}}{s^{4+k}}.

Vamos olhar cada um desses termos:

\displaystyle F_k(s) = \frac{2}{s^{3+k}} \times \frac{e^{-ks}}{s}.

Notem que

\displaystyle {\cal L}^{-1}\left[\frac{2}{s^{3+k}}  \right] = \frac{2t^{k+2}}{(k+2)!}

\displaystyle {\cal L}^{-1}\left[\frac{e^{-ks}}{s}  \right] =  \theta(t-k)

e da propriedade da convolução, tem-se

\displaystyle {\cal L}^{-1}\left[F_k(s)\right] = 2\int_0^t \frac{(t-\tau)^{k+2}}{(k+2)!}\theta(\tau -k)  d\tau = 2\frac{(t-k)^{k+3}}{(k+3)!},

para t>k, caso contrário é zero. A resposta do livro segue naturalmente daqui.

Um ponto interessante, e de certa forma sutil, é que é suficiente conhecer a transformada de Laplace para \Re s> 1 para se resolver o problema. (Por que? Vejam a fórmula da inversa e relembre as propriedades das funções analíticas e suas continuações no plano complexo.)

2) Integração no plano complexo

Da fórmula de Mellin, temos

\displaystyle f(t) = \frac{1}{2\pi i} \int_{\gamma -i\infty}^{\gamma -i\infty} e^{ts}F(s)\, ds =\frac{1}{\pi i} \sum_{k=0}^\infty \int_{\gamma -i\infty}^{\gamma -i\infty} \frac{e^{(t-k)s}}{s^{4+k}} \, ds  

para um certo  \gamma >0. Vamos calcular a integral para um dado k fixo. Há duas possibilidades. Para t<k, a integral será zero, pois ao fecharmos o contorno pela direita, não encontramos nenhum polo, e portanto a integral se anula pelo argumento usual. Para t>k, por outro lado, a integral engloba um polo de ordem 4+k, e pelo cálculo usual do teorema dos resíduos, obtemos exatamente o resultado do item anterior. (Confiram!)

Notem que se quiséssemos calcular a integral sem utilizar a série geométrica, teríamos dificuldades de cara com a localização dos zeros de s-e^{-s}. De longe, seria a menos problemática… Notem que estes zeros correspondem aos pontos do plano complexo tais que ze^{z}=1 e, portanto são os pontos z=W(1) sendo W(z) a chamada função W de Lambert, que já foi motivo de diversão em cursos passados de variáveis complexas e, pasmem, Cálculo Numérico.  Esta função é um tanto complicada.

Sobre o trabalho para quarta-feira

Quarta-feira é a nossa P1 e o dia combinado para a entrega do trabalho. Porém, se necessário, podemos deixar para outra data, por exemplo, o dia 4/10, quarta-feira. Acho que ninguém aqui tem dúvidas acerca da finalidade deste tipo de atividade extra.

Duas perguntas me foram feitas sobre como, a partir do roteiro proposto, encontrar a forma explícita da raiz quadrada da transformada de Fourier (e também do caso fracionário mais geral). As perguntas foram:

  1. Como encontrar a transformada integral a partir das suas autofunções?
  2. Como, dada uma certa transformada, encontrar suas autofunções?

Ambas questões envolvem equações integrais, algo poco discutido nos cursos que vocês tiveram até agora, mas extremamente importante pelas mais diversas aplicações que estas equações têm. Não faz parte da nossa ementa, é um tópico que normalmente se discute em Métodos III, uma optativa que nunca está nas listas das preferidas da turma e, portanto, raramente é oferecida.  😦

Dado um operador integral com kernel G(k,x), o problema de auto-valores em questão corresponde a encontrar os pares autovalores/autofunções (\lambda_n,f_n) tais que

\displaystyle \int_a^b G(k,x)f_n(x)\, dx = \lambda_n f_n(k)

Este tipo de equação tem uma longa tradição da Física e, obviamente, também na Matemática. Chama-se equação integral de Fredholm (de segundo tipo) e sua teoria é riquíssima. Como uma leitura elementar e introdutória, sugiro o Capítulo 16 do Arfken. Por ora, adianto que este tipo de problema tem muitas semelhanças com o Problema de Sturm-Liuville (de fato, pode-se sempre escrever o problema de SL numa forma integral de Fredholm). Em particular, haverá um produto interno “canônico” associado ao operador integral, e as autofunções serão ortogonais em relação a esse produto.

Vamos começar com a questão 1. Conhecendo-se os autovalores/autofunções (\lambda_n,f_n) de um certo operador integral, podemos “reconstruir” o núcleo G(k,x)? A resposta é sim, se o conjunto de autofunções for completo (base para o espaço pertinente de funções). Eu os convido a testar esta solução:

\displaystyle   G(k,x) = \sum_n\lambda_n f_n(k)f_n(x)

Funciona, não? Além disso, é uma solução “óbvia” (depois que se é apresentado a ela, certo… 🙂 )? Isso não te remete as construções do seu professor de Mecânica Quântica de coisas do tipo \displaystyle \sum_n |x_n\rangle \langle x_n | ? Se estivéssemos em dimensão finita, como reconstruiríamos uma matriz a partir dos seus autovalores/autovetores? Pois estes exemplo são exatamente a mesma coisa…

No caso particular do nosso problema, as autofunções são as funções de Hermite, que, como vocês provavelmente sabem, são as autofunções de quem mesmo? Sim! Dele! O omnipresente! 😮  The harmonic oscillator rules!!!  O kernel da transformada de Fourier vem do oscilador harmônico! Fantástico! Para mostrar esse resultado, explorem a função geratriz dos polinômios de Hermite. Um pouco de álgebra e vocês terão o resultado!

A questão 2 é de certa maneira a complementar. Conhecendo-se o kernel G(k,x), como resolver a equação de autovalores/autofunções? Há várias maneiras. O kernel da transformada de Fourier é muito simples, G(k,x) = \frac{1}{\sqrt{2\pi}} e^{ikx}, não é difícil fazer algumas manipulações com derivas e encontrar uma EDO para as autofunções. O surpreendente e que esta construção também nos remeterá ao oscilador harmônico! (Óbvio, pois descobriremos exatamente as suas autofunções). Mais precisamente, aos famosos operadores “subir” e “descer” nível, que são a base dos operadores de criação/aniquilação da Teoria Quântica de Campos. Vocês não precisam fazer isso, mas se quiserem, go ahead, aposto que se divertirão.

Como falei na aula de sexta, meu interesse nesse assunto foi reavivado por este artigo de memórias. Nele, o autor, um célebre físico teórico que está se recuperando de uma grave doença, conta sua experiência de vida como físico. Na seção 2.4, ele conta suas lembranças de Feynman, segundo suas próprias palavras

I first met Feynman as an idol, not a person.

Conta que este problema foi proposto por Feynman para seus alunos no Caltech nos anos 70, e dentre eles,  estava o autor. Tenho certeza que vocês merecem a mesma deferência, portanto, ai está o problema para o deleite dos interessados… Só lamento não ser possível aqui ter uma disciplina como a Physics X…. 

A Transformada Inversa de Laplace

Como prometido, vou deixar aqui a construção da Transformada inversa de Laplace a partir da Transformada de Fourier. Lembrando, a transformada de Laplace de uma função f(x) é definida por

\displaystyle {\cal L}\left[ f(x)\right] = F(s) = \int_0^\infty e^{-sx}f(x)\, dx,

quando a integral existe, sendo s\in \mathbb{C}. Façamos agora a substituição  s= u +iv, sendo uv reais. Lembrando as definições das Transformada de Fourier de uma função g(x)

\displaystyle {\cal F}\left[ g(x)\right] = \hat G(k) = \frac{1}{\sqrt{2\pi}}\int_{-\infty}^\infty g(x)e^{ikx}\, dx,

e de sua inversa

\displaystyle {\cal F}^{-1}\left[ \hat G(k)\right] = g(x) = \frac{1}{\sqrt{2\pi}}\int_{-\infty}^\infty \hat G(k)e^{-ikx}\, dk,

ficamos com

\displaystyle  F(u+iv) = \int_{-\infty}^\infty \underbrace{\left[e^{-st}\theta(t)f(t)\right]}_{\hat H_u(t)}e^{-ivt}\, dt = \sqrt{2\pi} h_u(v) ,

Calculando-se agora a transformada de Fourier de h_u(v) tem-se

\displaystyle \hat H_u(x)  =e^{-sx}\theta(x)f(x) =  \frac{1}{2\pi}\int_{-\infty}^\infty  F(u+iv) e^{ivx}dv ,

de onde temos

\displaystyle \theta(x)f(x) =  \frac{1}{2\pi}\int_{-\infty}^\infty  F(u+iv) e^{(u+iv)x}dv .

Introduzindo-se a variável z = u+iv, ficamos finalmente com

\displaystyle \theta(x)f(x) =  \frac{1}{2\pi i}\int_{u-i\infty }^{u+i\infty }  F(z) e^{xz}dz,

que é a fórmula usual. Notem que para garantir que a integral se anula para x<0, devemos escolher o valor de u de tal maneira que todos os polos de F(z) estejam a esquerda da reta vertical ao longo da qual a integral é calculada no plano complexo. Isto é fácil de ser verificado via teorema dos resíduos (façam!). Para x<0, a integral pode ser calculada a partir de um caminho fechado pela direita da reta. A única maneira de garantir que ela seja zero para qualquer x<0, é que não haja nenhum polo de F(z)  localizado a direita da reta vertical.  (Fechando-se o caminho da maneira usual, com um semicírculo, pode-se usar o Lema de Jordan para concluir que, no limite do raio grande, a contribuição do semicírculo vai a zero e, portanto, a integral que procuramos é igual a soma dos resíduos que estão a direita (para x0) da reta vertical. Façam os detalhes.).

Exercícios e uma dúvida

Atendendo aos pedidos que foram feitos na segunda-feira, a Stefânia organizou aqui todas suas notas com as soluções apresentadas nas aulas. Não deixem de agradecer-lhe  (clique aqui para mais detalhes sobre a regência do verbo agradecer… 🙂 )

Na mesma segunda-feira, o Matheus levantou um ponto correto: na solução do T3, que está aqui, a solução na forma (3) não existe para o caso do amortecimento crítico, para o qual \lambda_+ = \lambda_- e, portanto, não há duas soluções L.I. para comodar as constantes A e B. Notem que esse caso é precisamente o caso para o qual \omega(k) = 0,  no qual as equações (9) vão apresentar problema. Para o caso do amortecimento crítico, a solução geral será uma combinação linear das soluções e^{\lambda t} e t e^{\lambda t}. Rigorosamente, portanto, a solução apresentada, em particular as Eq. (9), só valem para o caso não crítico \omega(k) \ne  0. O caso crítico corresponde as valores \displaystyle k = \pm\frac{1}{c}\sqrt{\alpha^2 - 4\beta}, o que significa que esses valores devem ser retirados da integral (6). Porém, sabemos que se retirarmos um conjunto finito de pontos numa integral de Riemann, não alteraremos nada…

PS.: Notaram a importância de numerar as equações em um texto físico/matemático? Vejam que não só  as equações que vocês próprios fazem referência devem ser numeradas, de preferência TODAS devem ser numeradas, para que outras pessoas possam fazer referências a elas como fizemos no parágrafo acima.

Um erro muito instrutivo

A questão 2 do T2 é provavelmente o problema mais interessante que encontramos nos últimos tempos… Há um erro sutil na nossa discussão, um erro tão instrutivo, sob vários pontos de vista, que não irei corrigi-lo. Desnecessário dizer que não há nenhum impacto para a nota do T2. O erro passou desapercebido porque, provavelmente, ele some no limite do contínuo, onde temos uma boa intuição do problema. Estes são os erros mais perigosos, aqueles que “somem” nos limites nos quais conhecemos melhor o problema em questão. Eu só percebi o erro olhando as animações. Vocês verão que o salto na velocidade é meio estranho, parece dessincronizado. Sim, ele efetivamente está!

Relembrando, queremos resolver a EDO

\displaystyle \ddot x  = \alpha x(t)\sum_{k=-\infty}^\infty\Delta t \delta\left( t-k\Delta t\right),

que é equivalente a

\displaystyle \ddot x  = \alpha  \sum_{k=-\infty}^\infty x(k\Delta t)\Delta t \delta\left( t-k\Delta t\right),

cuja solução geral é

\displaystyle x(t) = x_k + \dot x_k (t-k\Delta t),\quad k\Delta t < t < (k+1)\Delta t

A continuidade da solução exige

\displaystyle x_{k+1} = x_k +  \Delta t \dot x_k,

e o salto da derivada é

\displaystyle \dot x_{k+1} = \dot x_k + \nu^2\Delta t x_{k+1}.

Aqui está o erro! Inicialmente, havia considerado

\displaystyle \dot x_{k+1} = \dot x_k + \nu^2\Delta t x_{k},

mas este caso corresponde à equação

\displaystyle \ddot x  = \alpha  \sum_{k=-\infty}^\infty x((k-1)\Delta t)\Delta t \delta\left( t-k\Delta t\right),

que é de fato uma equação com atraso, i.e., o lado direito e o esquerdo da equação envolvem a função em instantes de tempos diferentes! Vamos ver que era isto que estava por trás do comportamento amplificado que obtivemos. No limite do contínuo, ambos os casos coincidem.

A solução do problema correto, sem atraso, é

\displaystyle  \left(  \begin{array}{c}  x_{k+1}\\  \dot x_{k+1}  \end{array}  \right) =  \left(  \begin{array}{cc}  1& \Delta t\\  \alpha\Delta t & 1 + \alpha \Delta t^2  \end{array}  \right)  \left(  \begin{array}{c}  x_{k}\\  \dot x_{k}  \end{array}  \right) ,

e a partir daqui os desenvolvimentos são análogos. Porém, há agora uma diferença fundamental! Essa matriz tem determinante 1, o que implica que, como uma transformação linear, ela preserva áreas em \mathbb{R}^2. Além disso, seus dois autovalores são tais que \lambda_+\lambda_- = 1. Caso sejam complexos, temos que \lambda_- = \lambda_+^*  |\lambda_+|= 1,  o que exclui a possibilidade de termos oscilações amplificadas como no caso do problema com atraso. De fato, amplificação/supressão exponencial é um comportamento usual nas equações lineares com atraso. Os autovalores da nossa matriz são

\displaystyle \lambda_\pm = \frac{1}{2}\left( {2+\alpha\Delta t^2 \pm 2\Delta t \sqrt{\alpha + \frac{\alpha^2\Delta t^2}{4}}}\right),

de onde temos que o movimento oscilatório “conservativo” aparecerá sempre que 4\alpha + \alpha^2 \Delta t^2 < 0. O limite do contínuo \Delta t \to 0 é idêntico ao caso com atraso, e foi por isso basicamente que não notamos o erro.

Abaixo, estão as animações correspondentes aos casos sem atraso, os corretos. Espero que notem o grande avanço na qualidade das produções da PIXAAR. 🙂  Vejam que os saltos agora estão como esperávamos, comparem com o caso com atraso. O código python para este caso sem atraso está aqui.

 

Aqui, vocês têm o caso \Delta t = 1, com as duas soluções: a com e a sem atraso. Comparem e vejam as diferenças.