Introducao as equacoes diferenciais parciais

atraves da discretizacao

Armando G. M. NevesUniversidade Federal de Minas Gerais

Departamento de Matematica

2a. Bienal da SBM, Salvador, 2004

Sumario

1 Discretizacao 31.1 EDP’s e discretizacao . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31.2 Discretizacao das derivadas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5

2 A equacao da difusao (ou do calor) 82.1 Modelagem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82.2 Condicoes de contorno . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10

2.2.1 Condicoes de contorno de Dirichlet . . . . . . . . . . . . . . . . . . . . . . . . 112.2.2 Condicoes de contorno de Neumann . . . . . . . . . . . . . . . . . . . . . . . . 112.2.3 Outras condicoes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12

2.3 Removendo a discretizacao no tempo . . . . . . . . . . . . . . . . . . . . . . . . . . . 122.4 Solucao numerica de problemas de difusao . . . . . . . . . . . . . . . . . . . . . . . . 13

2.4.1 Problemas de Dirichlet . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132.4.2 Estabilidade e princıpio do maximo e do mınimo . . . . . . . . . . . . . . . . . 17

3 A equacao da onda 223.1 Modelagem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 223.2 Problemas tıpicos de equacao da onda . . . . . . . . . . . . . . . . . . . . . . . . . . 243.3 Solucao numerica de problemas de onda . . . . . . . . . . . . . . . . . . . . . . . . . 253.4 A condicao de estabilidade para a solucao numerica de problemas com a equacao da

onda . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29

4 A equacao de Laplace 334.1 Aplicacoes da equacao de Laplace . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33

4.1.1 Potencial eletrostatico . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 334.1.2 Solucoes estacionarias da equacao de difusao . . . . . . . . . . . . . . . . . . . 35

4.2 A equacao de Laplace discretizada . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 364.2.1 O laplaciano discreto . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 364.2.2 Problemas de Dirichlet . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 374.2.3 Problemas de Neumann . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38

1

Prefacio

Em geral, EDP’s sao estudadas a nıvel basico apos um primeiro curso sobre Equacoes DiferenciaisOrdinarias (EDO’s), utilizando-se para isto o instrumental das series de Fourier e o metodo deseparacao de variaveis, ver por exemplo [1, 2, 3, 6]. Nossa intencao nestas notas e a de apresentaruma introducao ao assunto, sem tocar nas series de Fourier e no metodo de separacao de variaveis.Como pre-requisito bastam algum conhecimento basico sobre EDO’s lineares, por exemplo ao nıvelde [1], e um pouco de Algebra Linear. Na verdade, para quase tudo nestas notas, basta somente oconhecimento do conceito de derivada parcial.

A ideia de escrever uma introducao “diferente”a um assunto ja tao amplamente tratado parte deinterrogacoes feitas pelo proprio autor a si mesmo quando era estudante de graduacao e estudavapela primeira vez as EDP’s. Por que os metodos para resolver EDP’s sao tao diferentes dos metodospara resolver EDO’s? Por que aparecem as series de Fourier em EDP’s e nao em EDO’s? De ondevieram estas tais condicoes de contorno que nao apareciam antes em EDO’s? Por que usualmenteas solucoes de EDP’s aparecem como series? Por que nao existe para EDP’s uma teoria geral deexistencia e unicidade tao poderosa quanto a correpondente para EDO’s?

Nao consegui responder a estas perguntas naquela epoca. Embora as respostas sejam simplese estejam espalhadas pelas presentes notas, a maior parte dos livros sobre EDP’s nao as fornece.Quando, como professor, comecei a ensinar EDP’s da forma tradicional a meus alunos, percebi queas tecnicalidades do metodo de separacao de variaveis obscureciam o entendimento do assunto.

Estas notas sao portanto a minha tentativa de fornecer respostas simples a questoes teoricas.Idealmente, o lugar para situar um curso baseado nestas notas seria imediatamente antes de umcurso tradicional sobre series de Fourier e o metodo de separacao de variaveis. Ou entao, intercalarestas notas ao curso tradicional.

Uma palavra e devida a respeito de metodos numericos para solucoes aproximadas de EDP’s.Estes metodos, assim como o presente texto, partem de uma abordagem discreta as EDP’s. Conformeacima explicado, meu objetivo principal nao e o de ensinar metodos numericos para EDP’s e simensinar EDP’s. Porem, ensinar EDP’s da maneira presente leva-nos tao perto de aprender algo sobremetodos numericos que seria um grande desperdıcio nao mencionar esta conexao.

Ouvi recentemente uma palestra de um eminente cientista da computacao, prof. Juris Hartmanis.Afirma este que, apesar de a Ciencia da Computacao ter-se originado na necessidade de cientistasde outras areas em realizar longos calculos, com o passar dos anos a Computacao adquiriu vidapropria no panorama das ciencias. Como ciencia autonoma, ela pode influenciar as outras ciencias.Entender as EDP’s atraves da discretizacao significa basicamente entende-las sob o ponto de vistade um computador. Isto leva certamente a um quadro incompleto, uma vez que o engenho humano ebastante superior ao das maquinas, mas de forma alguma a um conhecimento inutil. Concordo como prof. Hartmanis em que a Ciencia da Computacao pode influenciar a Matematica e este texto eum exemplo do que julgo uma influencia extremamente benefica deste tipo.

Agradeco aos alunos que foram “cobaias”de versoes preliminares deste texto e que me ajudarama melhora-lo. A comissao organizadora da 2a. Bienal da SBM pela maravilhosa oportunidade delecionar um curso no evento, o que docemente me obrigou a escrever estas notas. Aos professoresMichael O’Carroll e Giovanni Gallavotti, que em estagios diversos de minha formacao me ofereceramexemplos da utilidade de discretizar sistemas contınuos. Finalmente a minha mulher Marcia e a meufilho Francisco, que todos os dias me ensinam como e bom viver e conhecer.

Belo Horizonte, setembro de 2004.

2

Capıtulo 1

Discretizacao

1.1 EDP’s e discretizacao

Uma equacao diferencial parcial (doravante abreviada como EDP) e uma equacao que envolve asderivadas parciais de uma funcao incognita de varias variaveis. Se u e uma funcao de x e t, porexemplo, entao as equacoes ∂u

∂t= ∂2u

∂x2 ,∂2u∂t2

= u + ∂2u∂x2 ou cos u ∂u

∂t= ∂2u

∂x∂tsao exemplos de EDP’s.

As equacoes envolvendo funcoes incognitas de uma unica variavel e suas derivadas, equacoes estasprovavelmente ja estudadas pelo leitor em um primeiro curso de equacoes diferenciais, sao chamadasequacoes diferenciais ordinarias (porque envolvem derivadas ordinarias) e abreviadas doravante comoEDO’s.

Antes de introduzir as primeiras EDP’s que estudaremos, faremos uma pequena digressao em umassunto que e certamente conhecido dos leitores que ja fizeram um estudo de metodos aproximadosde solucao numerica de EDO’s: a discretizacao de funcoes.

Considere uma funcao real f definida em um intervalo [a, b]. Se f e contınua, entao para dar umadescricao bastante boa do comportamento de f nao e necessario especificar o valor de f em todosos pontos do seu domınio; basta especificar este valor em um conjunto X = x0, x1, x2, . . . , xN depontos do domınio, onde N e um inteiro positivo suficientemente grande e os elementos razoavelmentebem “espalhados”por todo o intervalo [a, b]. Da continuidade de f segue que o valor f(c) de f emum ponto c /∈ X sera proximo do valor de f no ponto xi ∈ X que esteja mais proximo de c. Se onumero n+1 de elementos em X for suficientemente grande e estes elementos forem bem espalhadosem X, entao pouco se perde em aproximar f(c) por f(xi). O processo de se trocar um conjuntocontınuo de informacoes (todos os valores de f) por um conjunto finito (o conjunto dos valoresf(x0), f(x1), f(x2), . . . , f(xn)) e dito uma discretizacao da funcao f . O conjunto X dos valoresutilizados para discretizar f e chamado a rede de discretizacao.

Uma maneira comum de se discretizar de maneira simples uma funcao sem ter que se preocu-par se os pontos x0, x1, x2, . . . , xN da rede estao ou nao bem espalhados no domınio e escolherestes pontos como sendo igualmente espacados. Divide-se o intervalo [a, b] em N sub-intervalos demesmo comprimento e escolhem-se os pontos da discretizacao como os extremos destes sub-intervalos.Chamando-se ∆x = b−a

Nao tamanho de cada sub-intervalo, os pontos da discretizacao sao dados por

xi = a + i∆x, i = 0, 1, 2, . . . , N . O tamanho ∆x dos sub-intervalos da discretizacao e chamadopasso de discretizacao. De agora em diante adotaremos sempre discretizacoes em redes de pontosigualmente espacados.

Por exemplo, suponha que f seja a funcao seno definida em [0, π]. Toda a informacao sobre osvalores de f no intervalo considerado esta contida em seu grafico. Podemos discretizar a funcao futilizando um passo ∆x = π

5. Isto significa trocar o grafico de f somente pelos 6 pontos indicados

na figura 1.1. Como se pode ver pelo grafico, o erro maximo nesta discretizacao e aproximadamente

3

π5

2 π5

3 π5

4 π5

π

0.2

0.4

0.6

0.8

1

Figura 1.1: Uma discretizacao do seno.

π5

2 π5

3 π5

4 π5

π

-2

-1.5

-1

-0.5

0.5

1

1.5

Figura 1.2: O seno e um polinomio podem ser identicos quando discretizados.

0,3.A discretizacao de uma funcao implica necessariamente em perda de informacao, uma vez que duas

funcoes muito diferentes podem tornar-se exatamente iguais em uma determinada discretizacao. Umexemplo e fornecido na figura 1.2, onde mostramos que o seno e uma outra funcao bastante distintacoincidem exatamente nos pontos determinados por uma discretizacao com passo ∆x = π

5.

E claro que se diminuirmos o valor do passo da discretizacao, duas funcoes que, discretizadas,coincidem com um determinado valor do passo, passarao provavelmente a ser distintas. E obvio poremque com um determinado passo de discretizacao havera sempre infinitas outras funcoes distintas quecoincidem com uma funcao dada.

E intuitivamente obvio que quanto menor for o passo de discretizacao, tanto menores serao oserros ao se aproximar o valor de f(c) em um ponto c /∈ X pelo valor f(xi), onde xi e o elemento deX mais proximo de c. Podemos dar uma expressao quantitativa para o erro de aproximacao se f for

4

de classe C1 em [a, b]. Neste caso, podemos usar o Teorema do Valor Medio para escrever

f(c) = f(xi) + f ′(ξ)(c− xi) ,

onde ξ e algum ponto de [a, b] entre c e xi. O erro |f(c)−f(xi)| cometido na discretizacao e portanto

|f(c)− f(xi)| = |f ′(ξ)| |c− xi| < M∆x , (1.1)

onde M e o maximo de |f ′(x)| para x ∈ [a, b]. Como o erro cometido na discretizacao e menorque alguma constante vezes ∆x, dizemos que o erro e da ordem de ∆x ou, simbolicamente, O(∆x).Escrevemos ainda,

f(c) = f(xi) + O(∆x) .

1.2 Discretizacao das derivadas

Estamos falando em discretizacao de funcoes com o intuito de aplicar este conceito no estudo deequacoes diferenciais. E portanto importante, tendo um conjunto discreto de dados representandoaproximadamente uma funcao, saber aproximar as derivadas desta funcao.

A maneira mais obvia de aproximar a derivada primeira de uma funcao discretizada e pelo proprioquociente de Newton que define esta derivada, ou seja, se conhecemos f(x0), f(x1), . . . , f(xN),aproximamos f ′(xi) como

f ′(xi) ≈ f(xi+1)− f(xi)

∆x. (1.2)

Esta aproximacao discreta para a derivada e chamada aproximacao para frente, uma vez que usamosum ponto da discretizacao de f a frente de xi para aproximar a derivada de f em xi. E igualmenterazoavel fazer uma aproximacao da derivada para tras, definida como

f ′(xi) ≈ f(xi)− f(xi−1)

∆x. (1.3)

O erro cometido nestas aproximacoes para as derivadas pode ser estimado de maneira semelhanteao resultado (1.1). No lugar do Teorema do Valor Medio, utilizaremos uma das suas generalizacoes,a formula de Taylor com resto na forma diferencial: se f e n + 1 vezes diferenciavel em [a, b] ec, d ∈ [a, b], entao existe ξ entre c e d tal que

f(d) = f(c) + f ′(c) (d− c) +1

2!f ′′(c) (d− c)2 + . . . +

1

n!f (n)(c) (d− c)n

+1

(n + 1)!f (n+1)(ξ) (d− c)n+1 . (1.4)

Supondo entao f de classe C2 e aplicando o resultado acima com n = 1, c = xi, d = xi+1,∆x = xi+1 − xi, temos que

f ′(xi) =f(xi+1)− f(xi)

∆x− 1

2f ′′(ξ) ∆x ,

donde o erro em utilizar a aproximacao (1.2) para se aproximar a derivada f ′ num ponto qualquer emenor que M ′∆x, onde M ′ e a metade do maximo de |f ′′(x)| para x ∈ [a, b]. Ou ainda, utilizandoa simbologia dos O’s, o erro cometido na aproximacao discreta para frente para a derivada primeirae O(∆x). Utilizando novamente (1.4), agora com n = 1, c = xi, d = xi−1, ∆x = xi − xi−1, obtem-seque o erro obtido com a discretizacao para tras (1.3) da derivada e tambem O(∆x).

5

Tendo obtido 2 maneiras distintas de aproximar a derivada discreta de uma funcao, podemos obteruma infinidade de outras aproximacoes para a derivada primeira tomando medias ponderadas compesos arbitrariamente escolhidos das aproximacoes (1.2) e (1.3). Como os erros nestas aproximacoessao O(∆x), entao em geral o erro em uma media da aproximacao para tras e da aproximacao parafrente sera tambem O(∆x), a menos que haja algum cancelamento entre os erros das 2 parcelas,situacao em que o erro podera ser menor. Este cancelamento de fato acontece se fizermos a mediacom pesos iguais das aproximacoes para tras e para frente, ou seja, se tomarmos

f ′(xi) ≈ 1

2

[f(xi+1)− f(xi)

∆x+

f(xi)− f(xi−1)

∆x

]

=f(xi+1)− f(xi−1)

2 ∆x. (1.5)

Esta aproximacao para a derivada primeira e, por razoes obvias, chamada de aproximacao cen-trada. Para calcular o erro cometido nesta aproximacao, supondo agora que f seja 3 vezes difer-enciavel, usamos a formula de Taylor (1.4) com n = 2 para obter

f(xi+1) = f(xi) + f ′(xi) ∆x +1

2f ′′(xi) (∆x)2 +

1

6f (3)(ξ1) (∆x)3

f(xi−1) = f(xi) − f ′(xi)∆x +1

2f ′′(xi) (∆x)2 − 1

6f (3)(ξ2) (∆x)3 .

Subtraindo as equacoes acima e dividindo por 2 ∆x tem-se que

f ′(xi) =f(xi+1)− f(xi−1)

2 ∆x− 1

6

[f (3)(ξ1) + f (3)(ξ2)

](∆x)2 , (1.6)

mostrando que o erro na aproximacao centrada para derivada primeira e O((∆x)2). Isto e muitovantajoso, pois quando ∆x e menor que 1, entao (∆x)2 e muito menor que ∆x e a aproximacaocentrada e melhor que ambas as aproximacoes precedentes para a derivada primeira.

Finalizamos esta secao mostrando como aproximar a derivada segunda da funcao f . Tomamos

f ′′(xi) ≈ f ′(xi+1)− f ′(xi)

∆x≈ 1

∆x

[f(xi+1)− f(xi)

∆x− f(xi)− f(xi−1)

∆x

]

=f(xi+1) + f(xi−1) − 2f(xi)

(∆x)2. (1.7)

Nos problemas a seguir pediremos ao leitor para provar que o erro em aproximar-se desta forma aderivada segunda de uma funcao e O((∆x)2).

Problemas do capıtulo 1

1.1 A funcao cujo grafico e mostrado na figura 1.2 juntamente com o grafico do seno e um polinomiop(x). Este polinomio foi escolhido como sendo o polinomio de menor grau possıvel que coincidisseexatamente com o seno nos pontos da discretizacao utilizada e tambem tal que p(π

2) = 1

2. Encontre

a formula para p(x).

1.2 (a) Faca tabelas de discretizacao da funcao seno no intervalo [0, π] utilizando varios passos dediscretizacao, por exemplo ∆x = π

4, π

5, π

6, π

7, π

8. Utilize estas tabelas para aproximar (para frente, para

tras e de maneira centrada) a derivada da funcao seno nos pontos da rede de discretizacao. Compare

6

o valor de cada uma das aproximacoes com o valor exato da derivada do seno em cada um dos pontosindicados e o erro cometido nas aproximacoes.

(b) Faca um grafico onde no eixo horizontal voce vai colocar o valor do passo da discretizacao eno vertical o erro maximo, dentre todos os pontos da rede, cometido na aproximacao para a frentepara a derivada do seno. Explique porque o grafico e aproximadamente uma reta de inclinacao 1

2

passando pela origem.(c) Repita a letra (a) para as aproximacoes para tras e centrada. Tire as conclusoes pertinentes.

1.3 Repita o problema 1.2 trocando a funcao seno pelo polinomio encontrado no problema 1.1. Porque os erros sao maiores neste caso que os correspondentes no caso da funcao seno?

1.4 (a) Utilize os dados do problema 1.2 para calcular, usando a equacao (1.7), aproximacoes paraa derivada segunda da funcao seno usando discretizacoes de passos diversos.

(b) Compare os resultados obtidos em (a) com os valores exatos da derivada segunda.(c) Calcule o erro maximo na aproximacao para a derivada segunda para cada valor de ∆x

utilizado na letra (a). Faca graficos do erro em funcao de ∆x e do erro em funcao de (∆x)2. Qualdestes graficos se aproxima mais de uma reta passando pela origem? O que pode ser concluıdo arespeito do erro da aproximacao (1.7) para a derivada segunda.

1.5 Deduza uma formula para o erro maximo na aproximacao discreta da derivada segunda dadapor (1.7).

Sugestao: Use (1.4) com d = xi+1, c = xi e n = 3 e depois repita para d = xi−1, c = xi e n = 3.Some os resultados.

7

Capıtulo 2

A equacao da difusao (ou do calor)

2.1 Modelagem

Quando colocamos uma quantidade de soluto em uma pequena regiao de uma solucao, por exemplouma colherada de acucar numa bacia com agua, sabemos da nossa experiencia quotidiana que o soluto,inicialmente concentrado proximo a regiao onde foi colocado, ira paulatinamente se espalhando aolongo de toda a solucao. Este fenomeno de “espalhamento”do soluto e chamado difusao.

Microscopicamente, a difusao e devida ao movimento termico desordenado das moleculas dosoluto. Gostarıamos de descrever a difusao, ou seja descrever a concentracao de soluto como funcaoda posicao e do tempo. Veremos que a concentracao sera dada como a solucao de uma EDP. EstaEDP, chamada equacao da difusao, sera deduzida a seguir a partir de hipoteses razoaveis construıdasa partir da explicacao microscopica qualitativa acima.

Descrevemos inicialmente a situacao fısica.

• A solucao esta contida em um tubo longo e fino de comprimento l e secao transversal de area A.Devido ao tubo ser longo e fino, pode-se desprezar qualquer movimento de soluto em direcoesperpendiculares ao tubo e portanto a posicao ao longo do tubo sera descrita por uma unicavariavel x, x ∈ [0, l].

• Ligados a ambas as extremidades do tubo estao reservatorios contendo um volume infinito desolucao a concentracao nula (solvente puro). Um reservatorio de volume infinito e impossıvel deser realizado na pratica, mas pode-se pensar em um cujo volume seja muito maior que o volumede solucao contido no tubo. A hipotese do reservatorio infinito com solucao a concentracao nulasera usada no nosso modelamento da difusao para afirmar que quando uma molecula de solutoalcanca uma das extremidades do tubo, caira imediatamente no reservatorio e, como este einfinito, nao mais retornara a solucao no tubo.

• A concentracao inicial de solucao e uma funcao dada f(x) da posicao x. A medida que passao tempo, espera-se que a concentracao ira mudar em cada ponto devido a difusao.

Se a concentracao no ponto x e no instante t e designada por u(x, t), entao o objetivo do estudoda difusao e exatamente encontrar esta funcao u(x, t) sendo dados a concentracao inicial f(x) e ainformacao sobre a presenca dos reservatorios com solucao a concentracao nula nos extremos do tubo.

Para simplificar o modelamento, em vez de estudarmos a situacao real em que as variaveis x e tda concentracao u(x, t) sao ambas contınuas, iremos introduzir discretizacoes em ambas as variaveiscom passos respectivamente iguais a ∆x e ∆t. Adotamos entao a seguinte notacao simplificada:usaremos ui,j no lugar de u(i∆x, j∆t) e fi em vez de f(i∆x), i = 0, 1, 2, . . . , N − 1, j = 0, 1, 2, . . .,

8

com ∆x = lN

. Para simplificar a linguagem, passaremos doravante a falar em posicao i e instante jem vez de posicao i∆x e instante j∆t.

A discretizacao do tempo significa que desprezaremos a passagem contınua do tempo e olharemospara a concentracao somente em um conjunto discreto de instantes separados uns dos outros porum intervalo temporal ∆t. Teremos uma situacao inicial em t = 0, teremos uma nova situacao emt = ∆t, uma nova situacao em t = 2∆t, 3∆t, . . ., ate onde for desejado resolver-se o problema.

A discretizacao da posicao significa que nao nos interessaremos pela posicao exata das moleculasdo soluto. Para isto, dividiremos toda a extensao do tubo em N sıtios de tamanho ∆x iniciandoem x = 0, x = ∆x, x = 2∆x, . . . , x = (N − 1)∆x = l − ∆x e especificaremos somente o ındicei = 0, 1, . . . , N − 1 do sıtio em que se localiza cada molecula.

A cada passagem do tempo discreto, havera uma redistribuicao das moleculas do soluto entreos sıtios. Dada uma regra de passagem do soluto de um sıtio para o outro, poder-se-a calculara concentracao da solucao em funcao da posicao e do tempo. Iremos admitir a regra fisicamenterazoavel exposta no paragrafo seguinte.

As moleculas do soluto somente podem saltar de um sıtio para um dos sıtios imediatamentevizinhos a direita ou a esquerda. Portanto, uma molecula que no instante j esteja no sıtio i, noproximo instante discreto, j + 1, ou permanecera onde estava, ou entao passara para um dentreos sıtios i − 1, ou i + 1. Como a passagem de soluto de um sıtio para o outro e consequencia domovimento termico desordenado, entao a passagem do sıtio i para o sıtio vizinho a direita i + 1 etao provavel quanto a passagem para o sıtio vizinho a esquerda i − 1. Alem do mais, a massa desoluto que passa do sıtio i para i± 1 e proporcional a massa de soluto no proprio sıtio i e portantoproporcional a concentracao em i. A massa de soluto que passa de i para i±1 sera ainda proporcionalao intervalo temporal ∆t, a area A da secao transversal do tubo e inversamente proporcional a ∆x.

Em resumo, se ∆mi→i±1,j representa a massa de soluto que no instante j passa do sıtio i parai+1 (e tambem para i− 1), entao a nossa hipotese fundamental sobre a qual esta baseada a equacaoda difusao e

∆mi→i±1,j =k A ∆t

∆xui,j , (2.1)

onde a constante de proporcionalidade positiva k e chamada constante de difusao e deve dependerdo soluto, do solvente e da temperatura.

A massa de soluto em i no instante j+1 e calculada por conservacao da massa: basta descontar damassa mij em i no instante j a parte dela que passou para os sıtios vizinhos, assim como acrescentara massa que veio dos sıtios vizinhos i± 1:

mi,j+1 = mij − (∆mi→i+1,j + ∆mi→i−1,j) + ∆mi+1→i,j + ∆mi−1→i,j .

Usando agora (2.1), temos

mi,j+1 = mij − k A ∆t

∆x2ui,j +

k A ∆t

∆x(ui+1,j + ui−1,j) .

Para fazermos aparecer concentracoes no lugar de massas na ultima equacao, dividimo-la pelo volumeA∆x de cada sıtio, obtendo

ui,j+1 = ui,j +k ∆t

(∆x)2(ui+1,j + ui−1,j − 2ui,j) (2.2)

= (1− 2s) ui,j + s (ui+1,j + ui−1,j) , (2.3)

onde o paramtero de estabilidade s > 0 e definido por

s = k∆t

(∆x)2(2.4)

9

e sera util mais adiante.Rearranjando os termos em (2.2), chegamos tambem a

ui,j+1 − ui,j

∆t= k

ui+1,j + ui−1,j − 2ui,j

(∆x)2, (2.5)

que chamaremos de a equacao da difusao discreta.Mostraremos na secao 2.4 como utilizar a equacao da difusao discreta, ou equivalentemente (2.3),

para resolver numericamente problemas de difusao. Por ora nos interessa terminar a deducao daEDP da difusao. Para isto, lembremo-nos que os passos de discretizacao ∆x e ∆t foram introduzidosde maneira artificial de maneira a simplificar a formulacao do problema da difusao, mas nao possuemnenhuma realidade fısica, uma vez que o tempo e o espaco sao contınuos. A continuidade do espacoe do tempo sao recuperadas ao fazermos o limite do contınuo ∆x → 0 e ∆t → 0. Neste limite, deacordo com os resultados da secao anterior sobre a discretizacao de derivadas, temos

ui,j+1 − ui,j

∆t−→ ∂u

∂t

eui+1,j + ui−1,j − 2ui,j

(∆x)2−→ ∂2u

∂x2.

No limite contınuo, (2.5) torna-se a EDP da difusao

∂u

∂t= k

∂2u

∂x2. (2.6)

Se em vez de um tubo contendo solucao a concentracao u(x, t) pensarmos em uma longa e finabarra condutora de calor e u(x, t) como sendo a temperatura no ponto x no instante t, entao amesma equacao (2.6) descrevera a temperatura u(x, t). Neste contexto, ela e chamada equacao docalor. Sabemos hoje que o calor e uma forma de energia associada ao movimento microscopicodas partıculas (moleculas). No inıcio do sec. XIX pensava-se que o calor fosse um fluido e quea temperatura seria sua concentracao. Usando esta hipotese, Fourier deduziu a equacao do calorusando um raciocınio contınuo completamente analogo ao que fizemos acima no caso discreto.

Apesar da falsidade de uma premissa, no caso do calor, e de hipoteses deliberadamente simplesusadas na deducao, a equacao do calor/difusao ja foi experimentalmente validada em inumerassituacoes, no sentido de que as suas solucoes concordam com observacoes experimentais.

Varios problemas simples e interessantes envolvendo a equacao da difusao podem ser exatamenteresolvidos usando-se series de Fourier e o chamado metodo de separacao de variaveis. Nao iremostratar estes importantes metodos no presente texto e remetemos o leitor interessado a algum dosvarios otimos livros sobre o assunto [1, 2, 3, 6].

2.2 Condicoes de contorno

Assim como os problemas de EDO’s, quase sempre problemas de EDP’s envolvem condicoes iniciais.No exemplo mencionado na secao anterior, os valores iniciais f(x) da concentracao u(x, t) constituemuma condicao inicial para o problema de difusao que descrevemos.

Alem da condicao inicial, de um ponto de vista fısico o problema descrito na secao anteriordepende fundamentalmente da informacao de que em cada extremidade do tubo existem reservatoriosde volume infinito com solucao a concentracao 0. Cada molecula de soluto que chega a uma dasextremidades do tubo e levada ao reservatorio ali presente e, por causa da condicao de volume infinito,

10

nunca mais retornara a solucao no tubo. Intuitivamente e de se esperar que, independentemente dacondicao inicial f(x), a concentracao de solucao em qualquer ponto do tubo tenda a 0 quando t →∞.

Podemos pensar em outros problemas realısticos de difusao com comportamento distinto. Por ex-emplo, em vez de estarem ligadas a reservatorios, as extremidades do tubo poderiam estar tampadas.Neste caso, nossa intuicao fısica diz que o soluto inicialmente concentrado de forma distinta ao longodo tubo deve se espalhar e a concentracao tender a um valor constante uniforme ao longo do tuboquando t → ∞. Pela condicao de conservacao da massa, o valor desta concentracao uniforme podeser calculado dividindo-se a massa total de soluto presente inicialmente, A

∫ l

0f(x) dx pelo volume

A l do tubo.Outras possibilidades ainda sao reservatorios de solucao com concentracao (possivelmente depen-

dente do tempo) especificada, ou ainda sistemas que bombeiem soluto para o interior do tubo ataxas dadas (tambem dependentes do tempo para maior generalidade). Ou mesmo um reservatorioem uma extremidade e uma bomba na outra.

O importante aqui e frisar que nossa intuicao espera que a solucao de um problema de difusao de-penda da especificacao de o que acontece em cada uma das extremidades do tubo. Esta especificacaoe feita pelas chamadas condicoes de contorno do problema.

2.2.1 Condicoes de contorno de Dirichlet

A condicao de reservatorio com volume infinito contendo solucao a concentracao g(t) em uma dasextremidades, digamos x = 0, forca que o valor da concentracao nesta extremidade do tubo seja g(t),ou seja, u(0, t) = g(t). Analogamente, um reservatorio com volume infinito ligado a x = l contendosolucao com concentracao h(t) forca com que se tenha u(l, t) = h(t). Condicoes de contorno destetipo, que especificam o valor em um ponto da funcao incognita u(x, t) de uma EDP sao chamadascondicoes de Dirichlet. Em particular, se g(t) ≡ 0, a condicao u(0, t) = g(t) e dita homogenea deDirichlet e analogamente para h(t) ≡ 0.

Um problema de Dirichlet para a equacao de difusao e portanto o de obter solucao para a equacaoda difusao respeitando uma condiccao inicial e condicoes de Dirichlet. Mais precisamente, o problemae o de encontrar uma funcao u : [0, l]× [0,∞) → R que obedeca

ut = k uxx, ∀x ∈ (0, l), ∀t > 0u(x, 0) = f(x), ∀x ∈ [0, l]u(0, t) = g(t), ∀t > 0u(l, t) = h(t), ∀t > 0

. (2.7)

No problema acima, a primeira equacao nos informa que se trata de um problema de difusao (oucalor), a segunda nos diz a condicao inicial e as duas ultimas representam a especificacao, para cadauma das extremidades da barra, da condicao de contorno correspondente.

O problema que descrevemos na secao anterior, com reservatorios infinitos de solvente nas duasextremidades do tubo, e portanto um problema do tipo acima com g(t) = h(t) ≡ 0, condicoeshomogeneas de Dirichlet em ambas as extremidades.

2.2.2 Condicoes de contorno de Neumann

O fluxo de massa φ(x, t) no ponto x e no instante t e definido como a massa de soluto por unidadede tempo atraves da secao transversal do tubo. Para descrever quantitativamente se o movimento desoluto em x e t e predominantemente da direita para a esquerda ou o contrario, moleculas atraves-sando a secao transversal da esquerda para a direita irao contar positivamente para φ(x, t), enquantoas que atravessarem a secao da direita para a esquerda irao contar negativamente.

11

No modelo discretizado, temos

φ(i∆x, j∆t) =∆mi→i+1,j − ∆mi+1→i,j

∆t

= −kAui+1,j − ui,j

∆x

limite contınuo−→ −kA∂u

∂x(x, t) , (2.8)

onde na passagem da primeira para a segunda linha utilizamos (2.1).Vemos assim que nos problemas de difusao, a derivada espacial ux da concentracao esta rela-

cionada como acima com o fluxo de massa. E comum em problemas realısticos que, em vez de seespecificar o valor da concentracao em uma extremidade do tubo, as condicoes especifiquem o valordo fluxo de massa nas extremidades. O exemplo mais comum e o de um tubo com extremidadesfechadas de modo a nao permitir passagem de soluto: o fluxo de massa e entao 0. Um problema dedifusao com um tubo com ambas as extremidades fechadas e entao modelado por um conjunto deequacoes como (2.7) em que as duas ultimas sao substituıdas pelas condicoes de contorno

ux(0, t) = 0, ∀t > 0ux(l, t) = 0, ∀t > 0

. (2.9)

Condicoes de contorno em que se especifica o valor da derivada da funcao incognita na direcao nor-mal a fronteira sao conhecidas como condicoes de contorno de Neumann. Em particular, as condicoesde tubo com extremidades fechadas (2.9) sao condicoes homogeneas de Neumann. Condicoes deNeumann nao-homogeneas (possivelmente dependentes do tempo) aparecem em problemas onde seespecifica a taxa de entrada ou saıda de soluto em uma extremidade.

2.2.3 Outras condicoes

Obviamente, podemos pensar em um tubo com solucao em que uma das extremidades e fechada,enquanto a outra esta ligada a um reservatorio. Uma das extremidades e entao modelada por umacondicao de Dirichlet e a outra por uma de Neumann.

Existem ainda problemas em que aparecem condicoes de fronteira especificando por exemplo ovalor da funcao αu + βux, onde α e β sao constantes, ou mesmo condicoes mais complicadas.

2.3 Removendo a discretizacao no tempo

Consideremos inicialmente o problema (2.7) de difusao com condicoes de Dirichlet homogeneas evoltemos atras na deducao das equacao da difusao. Se na passagem que levou de (2.5) a (2.6)removermos a discretizacao temporal sem remover a espacial, em vez de (2.7) chegaremos ao seguinteproblema:

dui

dt= k

(∆x)2(ui+1(t) + ui−1(t) − 2ui(t)), i = 1, 2, . . . , N − 1

ui(0) = f(i∆x), i = 1, 2, . . . , N − 1u0(t) = 0, t > 0uN(t) = 0, t > 0

, (2.10)

onde ui(t) = u(i∆x, t), i = 0, 1, 2, . . . , N corresponde a versao discretizada no espaco da concentracao.Observe que o problema de valor inicial com condicoes de contorno para a EDP da difusao torna-

se assim um problema para um sistema de EDO’s de 1a. ordem. Cada uma das N − 1 EDO’sescritas na primeira linha de (2.10) envolve alem da concentracao ui que comparece no lado direito

12

com sua derivada, as concentracoes nos sıtios vizinhos ui−1 e ui+1. As N − 1 equacoes na 2a. linhasao as condicoes iniciais para cada uma das EDO’s da 1a. linha. Observe ainda que as equacoesda 1a. linha com ındices i = 0 e i = N − 1 envolvem respectivamente funcoes u0(t) e uN(t) quenao comparecem em nenhuma outra das EDO’s. A 3a. e 4a. linhas de (2.10), originadas pelascondicoes de contorno de Dirichlet homogeneas em (2.7), sao especificacoes para estas funcoes. Deum ponto de vista matematico, temos um problema em que se pode aplicar com tranquilidade oconhecido teorema de existencia e unicidade para problemas de valor inicial em EDO’s. Observeporem que o numero de EDO’s do sistema e N − 1, onde N = l

∆xe um valor artificial introduzido

pela discretizacao espacial. Ao removermos esta, ∆x → 0 e portanto N → ∞. Ou seja, em algumsentido, uma EDP e como se fosse um sistema de infinitas EDO’s.

Alem do mais, as EDO’s em (2.10) sao todas lineares e, por causa da homogeneidade das condicoesde Dirichlet, homogeneas. Combinacoes lineares de solucoes deste sistema de EDO’s sao solucoes dosistema. Portanto o conjunto das solucoes do sistema e um espaco vetorial. Da teoria de existenciae unicidade, resulta que a dimensao deste espaco e N − 1. A solucao geral do sistema contemportanto N − 1 constantes a serem determinadas usando-se as condicoes iniciais ui(0) = f(i∆x),i = 1, 2, . . . , N − 1. Quando removermos a discretizacao espacial estaremos entao lidando com umespaco de solucoes de dimensao infinita. Para uma prova rigorosa disto, o leitor deve fazer oproblema2.4. Teremos ainda infinitas condicoes iniciais u(x, 0) = f(x), x ∈ (0, l) e , para determinar as infinitasconstantes que aparecem na solucao geral de modo a satisfazer as condicoes iniciais, o instrumentoa ser usado sao as series de Fourier.

Como ja dissemos antes, nao e de nosso interesse entrar neste texto nos problemas de comoresolver problemas de EDP’s. Paramos aqui com a compreensao de que estes problemas possuemanalogia com sistemas de infinitas EDO’s e devem envolver espacos de dimensao infinita.

2.4 Solucao numerica de problemas de difusao

Conforme dito anteriormente, podemos usar a versao (2.3) da equacao da difusao discreta como ummetodo numerico para aproximar a solucao de problemas de difusao. O leitor que tenha feito umprimeiro curso de EDO’s deve ter ouvido que nem todas as EDO’s podem ser resolvidas por metodosexatos. Quando isto acontece, uma das abordagens e aproximar numericamente as solucoes atravesde varios metodos, por exemplo os de Euler e Runge-Kutta, ambos tratados de forma introdutoriapor exemplo em [1].

Do ponto de vista de possibilidade de solucao exata, as EDP’s possivelmente se encontram emsituacao pior que a das EDO’s. Para resolver problemas que nao possam ser resolvidos de maneiraexata existem tambem metodos numericos para EDP’s. Iremos nesta secao descrever o metodo maissimples para se resolver problemas envolvendo a equacao da difusao. Este metodo, conhecido comometodo direto ou metodo das diferencas finitas pode-se considerar como uma versao para EDP’s dometodo de Euler para EDO’s.

2.4.1 Problemas de Dirichlet

Vamos explicar atraves de exemplos o procedimento para se resolver problemas de difusao comcondicoes de Dirichlet. A solucao de problemas de Neumann e um pouco mais trabalhosa e seratratada nos problemas 2.8 e 2.9.

Exemplo 1 Consideremos o problema de Dirichlet (2.7) homogeneo (g(t) = h(t) ≡ 0) com k = 1,l = 1 e f(x) = 100x(1 − x)2(1

3− x)2. O grafico da funcao f acima e mostrado na figura 2.1. Por

nossa intuicao fısica e de se esperar que a maior quantidade de soluto concentrada nas vizinhancas

13

0.2 0.4 0.6 0.8 1x

0.2

0.4

0.6

0.8

f HxL

Figura 2.1: Grafico da condicao inicial f(x) do exemplo 1. Observe que a concentracao inicial temdois maximos locais: um com valor maior em x ≈ 0, 706 e um com valor menor em x ≈ 0, 094.

do maximo absoluto de f em x ≈ 0, 706 tenda a se espalhar pelo resto do tubo, ao mesmo tempoque a quantidade total de soluto no tubo decresce por causa da perda de soluto para os reservatoriosnas duas extremidades do tubo.

Vamos usar (2.3) para enxergar o acontecimento deste fenomeno e tentarmos descobrir a escalade tempos em que ocorre. Primeiramente, precisamos discretizar f , escolhendo um valor adequadode ∆x. Embora percamos bastantes detalhes da condicao inicial, comecemos escolhendo por simpli-cidade ∆x = 1

4.

Precisamos ainda escolher o passo ∆t da discretizacao temporal. Para isto, lembre-se da deducaoda equacao da difusao (2.6), que (2.3) vem de uma discretizacao desta em que para a derivadatemporal e adotada uma aproximacao para a frente, enquanto que para a derivada segunda comrelacao a posicao usa-se a aproximacao (1.7). Lembrando que o erro na aproximacao da derivadatemporal e O(∆t) e que o erro na derivada espacial e O((∆x)2), vamos convencer o leitor de que einutil escolher ∆t muito maior ou muito menor que (∆x)2.

De fato, suponhamos que se queira acompanhar a solucao do problema de difusao ate um instantet = T . Para isto, o numero de passagens do tempo discreto e T

∆t. Se escolhermos um ∆t pequeno,

presumivelmente teremos um erro menor na aproximacao numerica, mas porem teremos que calcularum numero grande de passagens do tempo discreto, o que implica em grande trabalho computacional,como veremos. Se, ao contrario, escolhermos um valor grande para ∆t, podemos economizar trabalhocomputacional com um numero menor de passagens do tempo discreto, porem as custas de um maiorerro na discretizacao da derivada temporal.

O compromisso razoavel entre trabalho computacional e erro a ser adotado e escolher ∆t nemmuito maior, nem muito menor que (∆x)2. Caso ∆t ¿ (∆x)2, o trabalho computacional grandeem se calcular a derivada temporal com precisao e desperdicado porque o erro cometido na derivadaespacial sobrepuja o erro no calculo da derivada temporal. Podemos dizer que o erro no calculoda derivada espacial contamina o calculo da derivada temporal. Caso ∆t À (∆x)2, acontece ocontrario: o erro no calculo da derivada temporal contamina o calculo da derivada espacial.

Em nosso exemplo com ∆x = 14

iremos escolher, com base em nossa experiencia, ∆t = 140

. Oparametro s definido em (2.4) neste caso vale 2/5. Conforme argumentamos acima, para nao cometererros de truncamento muito grandes e nem desperdicar trabalho computacional, s deve ser da ordemda constante de difusao k.

14

i 0 1 2 3 4j0 0 0,0977 0,3472 0,8138 01 0 0,1584 0,4340 0,3016 02 0 0,2053 0,2708 0,2339 03 0 0,1494 0,2299 0,1551 04 0 0,1218 0,1678 0,1230 0

Tabela 2.1: Tabela de valores numericos ui,j da solucao do exemplo 1

A primeira coisa a fazer e discretizar f com passo ∆x = 14. Obtemos os valores ui,0 = f(i∆x), i =

0, 1, 2, 3, 4, que colocamos como a linha j = 0 da tabela 2.1. A proxima linha, j = 1, correspondentea t = ∆t pode ser obtida a partir da linha j = 0 e das condicoes de contorno. As condicoes decontorno homogeneas de Dirichlet especificam o primeiro (i = 0) e o ultimo (i = 4) valores da linhacomo sendo nulos. Os demais valores sao obtidos aplicando-se para i = 1, 2, 3 a equacao (2.3) comj = 0 e s = 2/5. Por exemplo,

u1,1 = (1− 2s) u1,0 + s (u0,0 + u2,0)

= 0, 2 · 0, 0977 + 0, 4(0 + 0, 3472) = 0, 1584 .

A linha j + 1 da tabela e obtida a partir da linha j usando o mesmo procedimento que utilizamospara obter a linha j = 1 a partir da condicao inicial na linha j = 0.

Analisemos o resultado. A concentracao no ponto x = 1/4 (ou seja, i = 1), inicialmente pequena,cresce quando j passa de 0 a 2 (t passa de 0 a t = 2∆t = 0, 050) e depois passa a decrescer.O crescimento inicial corresponde a chegada em x = 1/4 de soluto inicialmente localizado maisfortemente na metade direita do tubo. Em x = 1/2 (i = 2), a concentracao, pelo mesmo motivo,cresce quando j passa de 0 a 1 e depois decresce. Em x = 3/4 (i = 3), onde o soluto estava muitoconcentrado inicialmente, a concentracao decresce o tempo todo. A partir de j = 2 (t ≥ 0, 050) aconcentracao decresce em todos os pontos. Isto corresponde a perda de soluto para os reservatoriosde solvente puro em x = 0 e x = 1. E de se esperar que se continuarmos o calculo para valoresmaiores de j, este fenomeno deve continuar e a concentracao tender a zero em todos os pontos dotubo. O leitor e convidado a calcular algumas linhas a mais na tabela 2.1 e verificar isto.

De qualquer forma, o resultado acima analisado parece concordar bastante razoavelmente comnossa intuicao fısica a respeito do que deveria acontecer.

Exemplo 2 O fato de termos tomado ∆x = 1/4 no exemplo anterior nos fez perder muitos detalhesespaciais do fenomeno fısico. Para enxergarmos melhor estes detalhes, tomemos a mesma condicaoinicial do exemplo 1, os mesmos valores de k e ∆t e o mesmo intervalo de tempo, ou seja t ∈ [0; 0, 100].Tomemos porem ∆x = 1/8, a metade do valor do exemplo anterior. Com estes valores, s = 8/5. Osvalores da solucao numerica sao mostrados na tabela 2.2, analoga a tabela 2.1.

Os resultados sao pessimos! Ja na linha j = 1, em i = 1 aparece uma concentracao negativa,algo completamente sem significado em nosso modelo. E nas linhas seguintes, alem de mais valoresnegativos, as concentracoes parecem estar aumentando em valor absoluto. Esta solucao, que deveriadar em mais detalhes o que acontecia no exemplo 1, nao se parece em nada com este.

Antes de passar a um proximo exemplo, precisamos explicar o que aconteceu de errado no exemplo2. Afinal de contas, o exemplo 1, apesar de calculado com valores muito grandes de ∆x e ∆t, pareciaindicar que o metodo numerico funciona. No exemplo 2, simplesmente tentamos estudar o mesmoproblema com um valor menor de ∆x e o mesmo de ∆t. Quando o metodo das diferencas finitas

15

i 0 1 2 3 4 5 6 7 8j0 0 0,4154 0,0977 0,0254 0,3472 0,7477 0,8138 0,4011 01 0 -0,7576 0,4905 0,6559 0,4731 0,2127 0,0477 0,4196 02 0 2,4514 -1,2418 0,0988 0,3490 0,3653 0,9067 -0,8467 03 0 -7,3799 6,8121 -1,6458 -0,0252 1,2054 -2,7649 3,3134 04 0 27,1352 -29,4278 14,4799 -0,6493 -7,1160 13,3129 -11,7134 0

Tabela 2.2: Tabela de valores numericos ui,j da solucao do exemplo 2.

1

1t=0.06

1

1t=0.07

1

1t= 0.10

1

1t=0.03

1

1t=0.04

1

1t=0.05

1

1t=0.

1

1t=0.01

1

1t=0.02

Figura 2.2: Graficos de u em funcao de x relativos a varios valores de t no exemplo 3. Usamos∆x = 1

50e ∆t = 1

5000.

falha, como no exemplo 2, dizemos que e instavel. A explicacao para a falha no exemplo 2 e dada noteorema abaixo, demonstrado em [5]:

Teorema 1 O metodo das diferencas finitas para a equacao da difusao e instavel se s > 1/2 e estavelse s ≤ 1/2.

Iremos falar mais a respeito da estabilidade do metodo na subsecao 2.4.2

Exemplo 3 Tomamos ainda o mesmo problema ja estudado nos exemplos 1 e 2 e calculamos suasolucao aproximada ate t = 1/10 como nos exemplos anteriores. Fazemos agora ∆x = 1

50e tomamos

∆t = 15000

, o maior valor possıvel de modo a ainda ter estabilidade. Em vez de mostrar uma tabela,que seria grande demais e sem interesse (501 linhas e 51 colunas), optamos por ilustrar a solucao naforma de graficos obtidos a partir da tabela que nao mostramos. Na figura 2.2 sao mostrados algunsdos graficos de u em funcao de x para varios valores de t. Os calculos foram obviamente feitos porcomputador.

E razoavel que se quisermos diminuir o erro de discretizacao, podemos faze-lo diminuindo ∆x,ou seja, aumentando o numero de colunas numa tabela como 2.1. A condicao de estabilidade nosobriga a tambem diminuir ∆t, ou seja, tambem o numero de linhas na tabela deve ser aumentado.

16

x1 x2x

t1

t2

t

l1

l2

l3R

Figura 2.3: O retangulo R e os lados l1, l2 e l3 relativos ao teorema 2.

Note porem que que se o numero de colunas e aumentado por um fator α, o numero de linhas devecrescer por α2, ou seja, muito mais. Como o erro devido a discretizacao e O((∆x)2) = O(∆t), pode-se ver que para diminuir o erro por um fator β, o tempo computacional cresce como β3/2.

2.4.2 Estabilidade e princıpio do maximo e do mınimo

Citamos sem demonstracao o resultado seguinte sobre a equacao da difusao, conhecido como Princıpiodo maximo e do mınimo. A demonstracao pode ser encontrada em muitos dos livros padrao sobreEDP’s, por exemplo [3] ou [6]. A demonstracao e simples, usando somente a teoria basica de maximose mınimos de funcoes de varias variaveis. O princıpio do maximo e do mınimo e muito util na teoria daequacao da difusao, sendo usado por exemplo para demonstrar unicidade de solucao para problemasde difusao.

Teorema 2 Seja u(x, t) solucao da equacao da difusao ut = k uxx no retangulo fechado

R = (x, t) ∈ R2; x1 ≤ x ≤ x2, t1 ≤ t ≤ t2 .

Considere os seguintes lados de R, ver figura 2.3:

l1 = (x1, t) ∈ R2; t1 ≤ t ≤ t2 ,

l2 = (x, t1) ∈ R2; x1 ≤ x ≤ x2e

l3 = (x2, t) ∈ R2; t1 ≤ t ≤ t2 .

Entao os pontos de maximo e de mınimo de u em R encontram-se em l1 ∪ l2 ∪ l3.

Alem da utilidade matematica, o teorema acima possui ainda uma explicacao fısica intuitiva.Considere por exempo um problema de difusao e concentre sua atencao na parte do tubo comx ∈ [x1, x2]. Caso a maior parte do soluto se encontre nesta regiao quando t = t1, tendera a sair delacom a passagem do tempo e portanto o ponto de maximo de u em R estara em l2. Caso em t = t1a maior parte do soluto se encontre na regiao x < x1, entao com a passagem do tempo parte deste

17

soluto chegara ate x1. Neste caso, o maximo de u em R sera encontrado em l1. Similarmente, se amaior parte do soluto em t = t1 estiver em x > x2, o maximo de u em R estara em l3. Raciocıniosemelhante vale para o mınimo.

Concluımos portanto que o princıpio do maximo e do mınimo reflete uma propriedade fısica dassolucoes da equacao de difusao, propriedade que desejavelmente deve ser preservada nas solucoesaproximadas obtidas atraves de qualquer algoritmo numerico. Embora nao tenhamos demonstradoo teorema 1, e interessante que apresentemos aqui uma outra justificativa para a condicao de esta-bilidade 0 < s ≤ 1/2:

Teorema 3 Sejam R, l1, l2 e l3 como no teorema 2. Se u e uma solucao da equacao da difusaodiscreta (2.3), entao seu maximo e seu mınimo em R estarao em l1 ∪ l2 ∪ l3 independentemente dascondicoes iniciais se, e somente se, 0 < s ≤ 1/2.

Prova Como no exemplo 2, e facil construir exemplos de condicoes iniciais que violem o enunciadose s > 1/2. Resta entao a provar que se 0 < s ≤ 1/2 entao o maximo e o mınimo da solucao de (2.3)em R encontram-se em l1 ∪ l2 ∪ l3.

Sejam I1∆x = x1 e I2∆x = x2, I1, I2 ∈ Z, os ındices que definem os lados verticais de R.Sejam ainda J1 e J2 os ındices que definem os outros dois lados, ou seja, J1∆t = t1, J2∆t = t2.Defina finalmente Mi,j = maxui,j, ui+1,j, ui−1,j. Se 0 < s ≤ 1/2, entao 1 − 2s ≥ 0, logo, parai ∈ I1 + 1, I1 + 2, . . . , I2 − 1 tem-se

ui,j+1 = (1− 2s) ui,j + s (ui+1,j + ui−1,j)

≤ (1− 2s) Mi,j + s (Mi,j + Mi,j) = Mi,j . (2.11)

Defina agoraMj = max

i∈I1,I1+1,...,I2ui,j .

O resultado (2.11) prova entao que Mj+1 ≤ Mj, a nao ser que Mj+1 ocorra em l1∪ l3. Daı, o maximode u, maxj∈J1,J1+1,...,J2 Mj, tem que ocorrer em l1 ∪ l2 ∪ l3, demonstrando o teorema para o caso domaximo. O caso do mınimo e analogo.

Problemas do capıtulo 2

2.1 Em vez da situacao unidimensional em que um soluto se difunde por uma solucao contida emum tubo, considere a situacao bidimensional analoga em que o soluto se encontra em uma regiaotridimensional de profundidade desprezıvel. Neste caso, sendo z a coordenada que descreve a profun-didade, a posicao e dada aproximadamente pelas coordenadas x e y com (x, y) ∈ Λ ⊂ R2.

Repita o procedimento de modelagem da secao 2.1 e prove que se u(x, y, t) e a concentracao emfuncao da posicao e do tempo, entao a equacao da difusao bidimensional e

∂u

∂t= ∆u ,

onde aqui

∆u ≡ ∂2u

∂x2+

∂2u

∂y2

denota o laplaciano de u.Obtenha tambem a equacao de difusao tridimensional.

18

2.2 Considere uma situacao de difusao em um tubo com extremos em x = 0 e x = l em que umabomba e usada para injetar M > 0 gramas de soluto por unidade de tempo para dentro da solucao.

(a) Mostre que se a bomba estiver em x = 0, entao a condicao de contorno correspondente e umacondicao de Neumann do tipo ux(0, t) = −B, onde B e uma constante positiva. Obtenha B.

(b) Qual a condicao de contorno correspondente se a bomba estiver em x = l? Atencao ao sinal!

2.3 Mostre que se pelo menos uma das condicoes homogeneas de Dirichlet em (2.10) e substituıdapor uma condicao nao-homogenea, entao uma combinacao linear de solucoes do sistema de EDO’snao e necessariamente solucao do sistema de EDO’s.

2.4 Considere o problema de Dirichlet homogeneo (2.7).(a) Por derivacao direta prove que

un(x, t) = e−k n2π2

l2t sen

nπx

l

sao solucoes da equacao da difusao, onde n = 1, 2, . . ..(b) Verifique ainda que os un acima obedecem as condicoes de Dirichlet homogeneas.(c) Prove que qualquer combinacao linear da forma

∑Nk=0 ck uk(x, t), onde N e um inteiro positivo

qualquer, e solucao da equacao da difusao e obedece as condicoes de Dirichlet homogenea.(d) Prove que o conjunto de funcoes u1, u2, . . . , uN e linearmente independente para qualquer

N .Sugestao: Fixe t e prove que as un(x, t) sao ortogonais com relacao ao produto interno definido

por

〈f, g〉 ≡∫ l

0

f(x)g(x)dx .

(e) Conclua que o conjunto de todas as solucoes da equacao da difusao e que tambem obedecamas condicoes de Dirichlet homogeneas e um espaco vetorial de dimensao infinita.

2.5 Resolva numericamente ate t = 3/35 o problema de Dirichlet homogeneo com k = l = 1 e

f(x) =

1, se 2

5≤ x ≤ 4

5

0, se 0 ≤ x < 25

ou 45

< x ≤ 1.

Use para ∆x os valores 1/4, 1/8 e 1/16 e tome ∆t = 16/35(∆x)2. Compare as solucoes para osvarios valores de ∆x umas com as outras.

2.6 Utilizando o metodo de separacao de variaveis, a solucao exata do problema anterior e

u(x, t) =∞∑

n=1

an e−n2π2t sen nπx ,

onde

an =2

nπ

(cos

2nπ

5− cos

4nπ

5

).

(a) Faca graficos da soma parcial da serie acima ate n = 10 para t = 0, 1/35, 2/35, 3/35.(b) Coloque em um mesmo sistema de eixos os graficos produzidos em (a) e as solucoes aproxi-

madas correspondentes obtidas no problema anterior. Compare-as.

19

2.7 Considere conjuntamente dois problemas de difusao com condicoes de Dirichlet homogeneas eo mesmo valor de k. O primeiro em um tubo de comprimento l1 = 1 e o segundo em um tubo comcomprimento l2 = 2. A condicao inicial para o primeiro e

f1(x) = g(x) =

1, se 7

15≤ x ≤ 8

15

0, se 0 ≤ x < 715

ou 815

< x ≤ 1

e a condicao inicial para o segundo e f2(x) = g(x2). Semelhantemente, sejam ∆x1 e ∆x2 os respectivos

passos de discretizacao espacial.(a) Mostre que se se tomar para ambos os problemas o mesmo valor de s e se se fizer ∆x2 = 2∆x1,

entao as tabelas analogas a tabela 2.1 para ambos os problemas sao identicas. Observe que isto naodepende da forma particular das condicoes iniciais.

(b) Produza 5 linhas da tabela a qual se alude na letra (a) com ∆x1 = 1/20 e s = 1/5. Observeo “espalhamento”da mancha de soluto inicialmente concentrada em torno dos pontos medios dasbarras.

(c) Defina a largura da mancha de soluto como sendo a distancia entre os dois pontos ondea concentracao vale a metade do valor maximo. Observe que inicialmente a mancha se expanderapidamente e, a medida que passa o tempo, passa a espalhar-se mais devagar.

(d) Se no problema do primeiro tubo a mancha levou um tempo t1 para chegar ate a largura d1,quanto tempo leva no segundo problema para chegar a uma largura 2 vezes maior?

(e) Conclua que na letra (c) a largura da mancha de soluto e proporcional a raiz quadrada dotempo, desde que o soluto esteja longe o suficiente das extremidades do tubo para que a quantidadeperdida para os reservatorios de volume infinito nas extremidades seja desprezıvel.

2.8 Considere o problema de Neumann

ut = k uxx, ∀x ∈ (0, l),∀t > 0u(x, 0) = f(x), ∀x ∈ [0, l]ux(0, t) = g(t), ∀t > 0ux(l, t) = h(t), ∀t > 0

. (2.12)

Do ponto de vista de metodos numericos, a diferenca entre problemas de Dirichlet e Neumann refere-se somente a qual procedimento tomar com relacao ao calculo da primeira e da ultima colunas nastabelas analogas a tabela 2.1. Uma possibilidade e a de usar uma aproximacao para frente paraaproximar

ux(0, t) ≈ u1,j − u0,j

∆x,

de onde se tem que a primeira coluna da tabela e calculada como u0,j = u1,j−∆x g(j∆t), j = 1, 2, . . ..(a) Utilize a mesma ideia e uma aproximacao para tras e obtenha uma formula para os elementos

uN,j da ultima coluna da tabela.(b) Utilize valores convenientes para ∆x e ∆t e resolva o problema de Neumann homogeneo para

a mesma condicao inicial do exemplo 1 do texto.(c) Como explicado na secao 2.4, e de se esperar que o erro na solucao numerica de um problema

de Dirichlet pelo metodo la explicado seja O(∆t) = O((∆x)2). Sera que se pode esperar o mesmoerro para problemas de Neumann resolvidos da maneira indicada neste problema?

2.9 Em vez de usar a aproximacao para frente (ou para tras) para as derivadas ux(0, t) e ux(l, t)em um problema de Neumann, pode-se usar a aproximacao centrada para estas derivadas. Isto exige

20

que sejam introduzidos “pontos fantasmas”u−1,j e uN+1,j a cada linha na tabela. O valor de u−1,j ecalculado a partir da aproximacao

u1,j − u−1,j

2∆x≈ ux(0, j∆t) = g(j∆t) ,

enquantou0,j+1 = (1− 2s)u0,j + s(u−1,j + u1,j) ,

se j = 0, 1, 2, . . ..(a) Obtenha formulas analogas para uN+1,j e uN,j.(b) Resolva o problema de Neumann indicado na letra (b) do problema 2.8 utilizando o esquema

aqui sugerido com os pontos fantasmas.(c) Dentre o esquema sugerido no problema 2.8 e o presente, qual maneira deve produzir resultados

com o menor erro?

21

Capıtulo 3

A equacao da onda

3.1 Modelagem

Considere uma corda elastica tracionada de comprimento l e com extremidades fixadas. Um exemplopara ter em mente e o de uma corda de violao ou violino, cujas oscilacoes queremos estudar paraentender o mecanismo de producao do som. Estas cordas tem extremidades fixadas a pontos no bracodo instrumento e em sua caixa de ressonancia e as tarrachas utilizadas para “afinar”o instrumentofornecem a tracao necessaria. Em nosso modelo, denotaremos por T o valor da forca de tracaona corda quando em sua posicao de equilıbrio, ou seja, quando nao esta produzindo som. Maisadiante especificaremos melhor a tracao. Denotaremos por µ a massa por unidade de comprimento,ou densidade linear, da corda. Iremos desprezar forcas de atrito e a forca da gravidade, esta ultimaquase sempre desprezıvel frente a tracao. Suporemos ainda que cada pequeno trecho da corda semova somente na direcao transversal a corda (hipotese de transversalidade) e que nao se afaste muitode sua posicao de equilıbrio (hipotese de pequenas oscilacoes).

Um modelo discreto para a corda e pensa-la como N − 1 partıculas, cada uma com massa igual aµ∆x, situadas nas posicoes x = i∆x, i = 1, 2, . . . , N − 1 ao longo da corda. A extremidade esquerdada corda esta situada em x = 0 e e associada ao ındice i = 0. A extremidade direita e situada emx = l = N∆x. A elasticidade da corda e representada por molas que acoplam cada partıcula as suasvizinhas, ver figura 3.1.

A equacao da onda sera obtida aplicando-se as partıculas “discretas”a 2a. lei de Newton daMecanica. Da hipotese de transversalidade, a posicao da i-esima massa, localizada em x = i∆x,i = 1, 2, . . . , N − 1 e especificada pela coordenada transversal ui, ver figura 3.2, onde ui = 0 significa

x

u

Figura 3.1: Modelo discreto para a corda.

22

ui-1 ui ui+1

Fi-1

Fi

Figura 3.2: Forcas agindo sobre a i-esima partıcula.

que a i-esima partıcula esta em sua posicao de equilıbrio. Na figura 3.2 mostramos ainda os vetoresforca agindo sobre a partıcula i.

Denotamos por Fi−1 a forca sobre a partıcula i exercida pela partıcula i−1 e por Fi a forca sobrea partıcula i exercida pela partıcula i + 1. Definimos θi−1 como o angulo de rotacao anti-horariaentre o sentido positivo do eixo horizontal e o vetor ligando a partıcula i− 1 a i e, semelhantemente,θi e o angulo no sentido anti-horario entre o sentido positivo do eixo horizontal e o vetor ligando apartıcula i a i + 1. Para que o movimento da i-esima partıcula seja transversal, e necessario que aresultante na direcao longitudinal das forcas agindo sobre esta partıcula seja nula, ou seja,

Fi−1 cos θi−1 = Fi cos θi . (3.1)

Dito em palavras, a tracao na corda e independente do ındice i de posicao ao longo da corda.A tracao pode em princıpio ainda depender do tempo t. Se em um determinado instante todas

as massas estao afastadas da posicao de equilıbrio, a corda deve estar mais tracionada que em uminstante em que as massas estao todas proximas as suas posicoes de equilıbrio. Nossa hipotese depequenas oscilacoes e de que as amplitudes dos movimentos transversais das partıculas sao pequenaso suficiente para que se possa desprezar a variacao temporal da tracao na corda. Podemos portantoidentificar ambos os lados de (3.1) com a tracao T constante com relacao a posicao e tempo.

Analisando agora a componente transversal da forca resultante na partıcula i e usando que T =Fi cos θi, i = 1, 2, . . . , N − 1, temos (atencao aos sinais!) que esta componente vale

−Fi−1 sen θi−1 + Fi sen θi = − T

cos θi−1

sen θi−1 +T

cos θi

sen θi

= T (tan θi−1 − tan θi) = T

(ui+1 − ui

∆x− ui − ui−1

∆x

)

= Tui+1 + ui−1 − 2ui

∆x.

Pela 2a. lei de Newton, esta componente deve ser igual a massa µ∆x multiplicada pela componentetransversal da aceleracao d2ui

dt2. Portanto, temos

d2ui

dt2=

c2

(∆x)2(ui+1(t) + ui−1(t)− 2ui(t)) , (3.2)

23

onde a constante c e, em termos dos parametros fısicos,

c =

√T

µ. (3.3)

Esta equacao deve ser comparada a equacao diferencial na primeira linha de (2.10), obtida quandoremovemos somente a discretizacao temporal na equacao da difusao discreta (2.5). Observando oaparecimento de uma discretizacao de ∂2u

∂x2 no membro direito de (3.2), a equacao de movimento dacorda no limite contınuo e portanto

∂2u

∂t2= c2 ∂2u

∂x2. (3.4)

A equacao (3.4) acima e conhecida como equacao da onda unidimensional. A equacao da onda(3.4) ou sua versao em mais dimensoes aparece em diversos fenomenos fısicos, sempre associada afenomenos ondulatorios, por exemplo, no caso presente, ondas transversais em uma corda elasticatensionada, ou ainda ondas sonoras ou ondas eletromagneticas. No problema 3.1 daremos umajustificativa para o nome dado a equacao.

E interessante notar uma diferenca entre a deducao da equacao da difusao (2.6) no capıtulo 2e a presente deducao de outra EDP importante. Enquanto a primeira versao (2.3) que obtivemospara a equacao da difusao e discreta tanto na posicao quanto no tempo, no caso presente a primeiraversao para a EDP final (3.4) foi (3.2), discreta somente na posicao. Tal diferenca se deve ao fato determos usado na deducao da equacao da corda vibrante a 2a. lei de Newton, que ja e contınua navariavel tempo. Mais adiante, iremos discretizar a derivada temporal em (3.2) e obter uma versaodiscreta em ambas as variaveis para a equacao da onda. Tal versao sera usada, como (2.3) o foi, paraproduzir aproximacoes numericas para a solucao de problemas com a equacao da onda.

3.2 Problemas tıpicos de equacao da onda

Voltemos ao modelo discretizado de corda. Para determinar a posicao ui(t) da partıcula i, temosque resolver uma EDO de 2a. ordem (3.2). Como esta equacao envolve as posicoes das partıculasvizinhas, ou seja, ui−1 e ui+1, trata-se portanto de um sistema de EDO’s de 2a. ordem. Condicoesiniciais para este sistema sao a especificacao da posicao e da velocidade inicial para cada partıculada corda discreta, ou seja, ui(0) = f(i∆x) e dui

dt(0) = g(i∆x), i = 1, 2, . . . , N − 1, onde f e g sao

funcoes dadas definidas em [0, l].Observe ainda que as equacoes (3.2) para os ındices i = 1 e i = N − 1 envolvem respectivamente

funcoes u0 e uN . Estas ultimas nao sao funcoes incognitas, mas as posicoes das extremidades da corda,ou seja, u0 = uN = 0 no caso da corda de violao. Em outros problemas poderia ser interessanteestabelecer condicoes u0(t) = a(t) e uN(t) = b(t) que especifiquem como sao movidos os extremos dacorda por um agente externo. Um problema tıpico de corda vibrante discretizada e portanto

d2ui

dt2= c2

(∆x)2(ui+1(t) + ui−1(t) − 2ui(t)), i = 1, 2, . . . , N − 1

ui(0) = f(i∆x), i = 1, 2, . . . , N − 1dui

dt(0) = g(i∆x), i = 1, 2, . . . , N − 1

u0(t) = a(t), t > 0uN(t) = b(t), t > 0

, (3.5)

que deve ser comparado com seu analogo (2.10). Devido ao fato de as EDO’s neste caso serem de2a. ordem, temos,com relacao a (2.10), uma condicao inicial a mais para cada EDO.

24

Removendo-se a discretizacao espacial, temos o problema

utt = c2 uxx, ∀x ∈ (0, l),∀t > 0u(x, 0) = f(x), ∀x ∈ [0, l]ut(x, 0) = g(x), ∀x ∈ [0, l]u(0, t) = a(t), ∀t > 0u(l, t) = b(t), ∀t > 0

, (3.6)

analogo a (2.7). Em (3.6), a primeira linha e a EDP para o movimento da corda, a segunda e aterceira linhas sao as condicoes iniciais respectivamente para posicao e velocidade iniciais de todosos pontos da corda e as duas ultimas linhas sao condicoes de contorno de Dirichlet para a funcaoincognita u(x, t).

Analogamente ao que foi feito na secao 2.3, podemos convencer o leitor de que a solucao deum problema como (3.6) envolve espacos vetoriais de dimensao infinita. Para isto, volte a versaodiscreta (3.5) do mesmo. No caso em que as condicoes de contorno de Dirichlet sao homogeneas, osistema de EDO’s e linear e homogeneo. Da teoria geral das EDO’s lineares, o conjunto de todasas solucoes do sistema obedecendo as condicoes de Dirichlet homogeneas e um espaco vetorial dedimensao 2(N−1). No caso de condicoes de Dirichlet nao-homogeneas, a solucao geral do sistema deEDO’s sera a soma de uma solucao particular deste com uma combinacao linear de 2(N−1) solucoeslinearmente independentes do sistema homogeneo. Ao removermos a discretizacao espacial, N →∞e e razoavel esperarmos o aparecimento de espacos de dimensao infinita e solucoes dadas com series,por causa de combinacoes lineares com infinitos termos. No problema 3.2, o leitor sera levado aprovar que de fato o conjunto das solucoes da EDP da onda satisfazendo condicoes homogeneas deDirichlet e um espaco vetorial de dimensao infinita.

O leitor deve ver que esta mesma linha de raciocınio mostra que tambem a EDP da onda e “comose fosse”um sistema de infinitas EDO’s.

Para encerrar esta secao, mencionamos de passagem que problemas de equacao da onda comcondicoes de Neumann tambem sao interessantes. Um exemplo deste tipo e o problema da vibracaode uma coluna de ar dentro de um tubo com ambas as extremidades abertas, como numa flauta.Neste exemplo, u(x, t) esta relacionado a pressao do ar no ponto x e no instante t e obedece a equacaoda onda com condicoes de Neumann homogeneas.

3.3 Solucao numerica de problemas de onda

Partindo do problema (3.5) e discretizando da maneira usual (1.7) tambem a derivada com relacaoa posicao, obtemos

ui,j+1 + ui,j−1 − 2ui,j

(∆t)2= c2 ui+1,j + ui−1,j − 2ui,j

(∆x)2,

analoga a (2.5). Podemos tambem resolve-la para ui,j+1, obtendo

ui,j+1 = (2− 2s) ui,j + s (ui+1,j + ui−1,j) − ui,j−1 , (3.7)

onde i = 1, 2, . . . , N − 1, j = 1, 2, . . . e agora, no contexto de equacao da onda,

s = c2 (∆t)2

(∆x)2. (3.8)

Como no caso da equacao da difusao, (3.7) pode ser usada como base de um esquema paraaproximar numericamente solucoes de problemas de Dirichlet para a equacao da onda. Ha algumasdiferencas interessantes com relacao ao caso da difusao.

25

Primeiramente, devido ao fato de a equacao da onda ser de segunda ordem no tempo, enquanto ade difusao e de primeira ordem nesta mesma variavel, (3.7) somente consegue fornecer a linha j + 1de uma tabela analoga a tabela 2.1 a partir das duas linhas anteriores, j e j − 1. Em compensacao,ha um conjunto de condicoes iniciais a mais.

Em segundo lugar, os erros nas discretizacoes das derivadas que comparecem na equacao deonda sao O((∆t)2) e O((∆x)2). Para que o erro na aproximacao discreta de uma das derivadas naoprejudique a outra, e razoavel que se tome ∆t = O(∆x) e portanto s = O(1). Como no caso daequacao da difusao, e de se esperar que s deva obedecer a alguma condicao para que o esquemanumerico seja estavel. Esta condicao, provada em [5], e

0 < s ≤ 1 . (3.9)

Mais adiante iremos fornecer, como no caso da equacao da difusao, uma “interpretacao fısica”dacondicao de estabilidade acima.

Pelo exposto no penultimo paragrafo, para levar adiante a intencao de resolver numericamenteproblemas ondulatorios com condicoes de Dirichlet (3.6), precisamos comecar uma tabela analogaa tabela (2.1) com duas linhas obtidas das duas condicoes iniciais u(x, 0) = f(x) e ut(x, 0) = g(x).Obviamente, a linha j = 0 e dada por

ui,0 = f(i∆x) , (3.10)

i = 1, 2, . . . , N−1. Uma maneira simples porem ruim para obter a linha j = 1 e usar para a derivadaut(x, 0) uma discretizacao para frente

ut(i∆x, 0) ≈ ui,1 − ui,0

∆t,

a qual usarıamos para obter

ui,1 ≈ ui,0 + ut(i∆x, 0) ∆t = ui,0 + g(i∆x) ∆t .

De fato, o erro em tal discretizacao e O(∆t), muito maior que os erros O((∆t)2) = O((∆x)2) queseriam cometidos na continuacao. Desta maneira, estarıamos, antes mesmo de comecar o problema,contaminando sua solucao com um erro grande demais, motivo pelo qual esta proposta simples seradescartada. A situacao aqui e parecida com a do problema 2.8, onde usamos discretizacoes parafrente e para tras para aproximar condicoes de Neumann em problemas de difusao.

A solucao para este impasse e semelhante a do problema 2.9. Em vez de usar a aproximacao parafrente, usemos para ut(x, 0) a aproximacao centrada

ut(i∆x, 0) ≈ ui,1 − ui,−1

2∆t, (3.11)

na qual o erro cometido e O((∆t)2), como desejado. So que agora temos o problema de definir quemsao os valores a colocar na “linha fantasma”com j = −1. Para isto, usamos ut(i∆x, 0) = g(i∆x) em(3.11) para obter

ui,1 − ui,−1 = 2∆t g(i∆x) . (3.12)

Isto ainda nao resolve o problema, pois tampouco conhecemos os valores ui,1. Uma segunda relacaoligando ui,1 e ui,−1 e obtida extrapolando (3.7) para j = 0. Mais precisamente, temos

ui,1 + ui,−1 = (2− 2s) ui,0 + s (ui+1,0 + ui−1,0)

= (2− 2s) f(i∆x) + s [f((i + 1)∆x) + f((i− 1)∆x)] . (3.13)

26

0.2 0.4 0.6 0.8 1x

0.2

0.4

0.6

0.8

1

f HxL

Figura 3.3: Grafico da condicao inicial f(x) do exemplo 1. Observe que o pulso esta mais proximoda extremidade x = 0 da corda do que da extremidade x = 1.

Podemos encarar (3.12) e (3.13) como um sistema de equacoes para obter ui,1 e ui,−1. Como naonos interessam os valores ui,−1, assinalamos aqui somente a solucao deste sistema para ui,1:

ui,1 = (1− s) f(i∆x) +s

2[f((i + 1)∆x) + f((i− 1)∆x)] + ∆t g(i∆x) , (3.14)

i = 1, 2, . . . , N − 1.Conhecendo agora as linhas j = 0 e j = 1 da tabela, dadas respectivamente por (3.10) e (3.14),

podemos complementa-las com as condicoes de contorno

u0,j = a(j∆t)

euN,j = b(j∆t) ,

ambas a serem usadas com j = 0, 1, 2, . . . e temos todos os elementos para calcular linha por linhauma tabela analoga a tabela 2.1.

Para melhor esclarecer o metodo numerico acima exposto, estudemos dois exemplos.

Exemplo 1 Suponha um problema de Dirichlet (3.6) com c = l = 1, condicoes homogeneas (a(t) =b(t) ≡ 0) e com condicoes iniciais g(t) ≡ 0 e

f(x) =

1 − 36(x− 5

12)2 , se 1

4≤ x ≤ 7

12

0 , se 0 ≤ x < 14

ou 712

< x ≤ 1.

O fato de termos condicoes homogeneas significa que estamos tratando uma corda fixa em ambasextremidades, como uma corda de violao. A condicao inicial g(x) ≡ 0 significa que a corda e esticadaate uma conformacao inicial u(x, 0) = f(x) e solta com velocidade 0 a partir desta conformacao,exatamente como no caso de um violonista que primeiro puxa a corda e depois a solta. Em particular,a forma do grafico da conformacao inicial f(x) que iremos utilizar e a de um pulso centrado emx = 5/12 e largura 1/3, ver figura 3.3. Embora a forma do pulso que estamos empregando nao sejaum bom modelo para a maneira como um violonista dedilha uma corda, e interessante investigar oque acontece na presente situacao.

27

i 0 1 2 3 4 5 6j0 0 0 0,7500 0,7500 0 0 01 0 0,2109 0,5391 0,5391 0,2109 0 02 0 0,4878 0,1436 0,1436 0,4878 0,1187 03 0 0,2966 -0,0583 -0,0583 0,3634 0,3782 04 0 -0,2610 -0,0605 -0,0230 0,0101 0,4167 0

Tabela 3.1: Tabela de valores numericos ui,j da solucao do exemplo 1.

Para ter uma primeira ideia do que acontece e entender o algoritmo numerico, comecemos tomando∆x = 1/6 e ∆t = 1/8. Temos entao s = 3/4, dentro da condicao de estabilidade. A linha j = 0 databela com os valores aproximados da solucao e obtida atraves da discretizacao da condicao inicialf(x), ou seja, (3.10), acrescida de dois zeros nas extremidades x = 0 e x = 1 por causa das condicoesde Dirichlet homogeneas. A linha j = 1 e obtida a partir de (3.14) que, no caso presente, pode serreeescrita como

ui,1 =1

4f(i∆x) +

3

8[f((i + 1)∆x) + f((i− 1)∆x)] ,

i = 1, 2, . . . , 5 e complementada pelas condicoes de contorno u0,1 = u6,1 = 0. Os demais valores databela sao dados por u0,j = u6,j = 0, j = 2, 3, . . . e pela aplicacao recursiva de (3.7) que, neste caso,pode ser reescrita como

ui,j+1 =1

2ui,j +

3

4(ui+1,j + ui−1,j) − ui,j−1 ,

onde i = 1, 2, . . . , N − 1, j = 1, 2, . . .. Indicamos na tabela 3.1 os valores para a solucao do problemaate t = 1/2.

Da analise da tabela, podemos ver que o pulso inicial inicialmente se divide em dois pedacos, umse movendo para a esquerda e outro para a direita. A “crista”de cada um destes pulsos originadosde pulso original se move a cada intervalo de tempo ∆t. Como o pulso inicial estava mais proximoda extremidade x = 0, chega primeiro a esta extremidade o pulso que viajou para a esquerda. Nalinha j = 4 encontramos um valor negativo apreciavel para u1,4, indicando que aproximadamente emt = 4∆t = 1/2 o pulso chegou ate a extremidade esquerda da corda e foi ali refletido, originando ovalor negativo. Ate este instante, o pulso nao havia ainda sido refletido na extremidade x = 1, poisainda nao houve tempo suficiente para que ali chegasse.

Os resultados parecem bastante em acordo com o que se espera fisicamente de um pulso que viajeem uma corda elastica fixa em ambas as extremidades. Observe que o esquema numerico pareceestar levando em conta corretamente a propagacao dos pulsos a velocidades bem determinadas e atemesmo o fenomeno de reflexao do pulso quando este alcanca uma das extremidades da corda. Iremosconfirmar estas previsoes corretas no proximo exemplo, quando estudarmos o mesmo problema,usando agora melhores resolucoes espacial ∆x e temporal ∆t.

Exemplo 2 Voltemos ao mesmo problema do exemplo 1, agora com ∆x = 1/24 e ∆t = 1/32.Observe que ao melhorarmos a resolucao espacial dividindo o valor de ∆x do exemplo anterior porum fator 4, a resolucao temporal pode ser aumentada pelo mesmo fator e ainda manter a condicao deestabilidade. Como visto no capıtulo anterior, no caso de um problema de difusao, seria necessariodividir ∆t por 16 para manter a condicao de estabilidade, o que torna problemas de difusao maistrabalhosos do ponto de vista numerico do que problemas analogos para a equacao da onda.

28

-1

1t=1.5

-1

1t=1.625

-1

1t=1.75

-1

1t=1.875

-1

1t=1.

-1

1t=1.125

-1

1t=1.25

-1

1t=1.375

-1

1t=0.5

-1

1t=0.625

-1

1t=0.75

-1

1t=0.875

-1

1t=0

-1

1t=0.125

-1

1t=0.25

-1

1t=0.375

Figura 3.4: Instantaneos de u em funcao de x para varios instantes do tempo relativos ao exemplo2. Usamos ∆x = 1/24 e ∆t = 1/32.

Apresentamos na figura 3.4 alguns dos resultados dos calculos feito pelo computador para asolucao do problema. Cada grafico nesta figura corresponde a um quadro de um “filme”do movimentoda corda, ou seja, um retrato da posicao de cada ponto da corda, onde os retratos foram tomados aintervalos de 4∆t = 0, 125. Observe nesta figura como o pulso inicial se divide em dois pulsos, umviajando para cada lado, cada um com velocidade constante. Observe ainda as reflexoes dos pulsos acada vez que chegam a uma das extremidades, com mudanca de orientacao a cada reflexao. Observeainda a interferencia construtiva entre os dois pulsos em t = 1. Tal figura confirma flagrantemente asconclusoes retiradas do exemplo 1. Sugerimos ao leitor que procure construir seu proprio programade computador para reproduzir tal figura e outras analogas. E ainda interessante apresentar graficoscomo os da figura 3.4 como uma figura com animacao. Isto e relativamente facil de se fazer usandosoftwares como Mathematica ou Maple.

3.4 A condicao de estabilidade para a solucao numerica de

problemas com a equacao da onda

Mostramos na subsecao 2.4.2 que as solucoes numericas aproximadas para a equacao da difusaoobedecem a um importante requerimento fısico, o princıpio do maximo e do mınimo, se e somentese a condicao de estabilidade 0 < s ≤ 1/2 for satisfeita. Isto nao e uma prova da condicao deestabilidade, mas joga alguma luz sobre a existencia de uma condicao deste tipo.

Apresentamos aqui resultados analogos para a equacao da onda. O leitor observe primeiramenteque o algoritmo numerico descrito na secao anterior propaga a onda de regioes onde u 6= 0 a regioesonde u = 0 a taxa maxima de uma unidade de posicao discreta ∆x a cada passagem do tempodiscreto ∆t. Diremos portanto que a velocidade maxima de propagacao da onda pelo algoritmonumerico e ∆x

∆t.

29

Por outro lado, como o leitor pode ver na equacao (3.15) no problema 3.1, as informacoes saofisicamente transportadas ao longo da corda com velocidade c. O parametro c em (3.4) recebeportanto a interpretacao de velocidade com que o meio, por exemplo a corda, transmite informacaoatraves de ondas. Uma das mais importantes passagens da Fısica foi quando Maxwell conjecturouque a luz fosse uma onda eletromagnetica, provando que os campos eletricos e magneticos obedecema versoes tridimensionais de (3.4). Em particular, o parametro c na equacao de onda obtida porMaxwell estava em bom acordo com as medidas da epoca para a velocidade da luz, fornecendoexcelente apoio experimental a sua conjectura. Hoje a conjectura de Maxwell esta amplamentedemonstrada dos pontos de vista experimental e teorico.

Voltando a estabilidade do algoritmo numerico, para que as solucoes obtidas atraves do algoritmoconsigam acompanhar a velocidade com que a onda se move, e necessario que a velocidade numerica∆x∆t

seja pelo menos igual a velocidade c. Daı, segue facilmente que

s = c2

(∆t

∆x

)2

≤ 1 .

Problemas do capıtulo 3

3.1 O metodo de D’Alembert:(a) Usando a regra da cadeia para funcoes de varias variaveis, prove que se φ e ψ sao funcoes

arbitrarias duas vezes diferenciaveis de uma unica variavel, entao

u(x, t) = φ(x− c t) + ψ(x + c t) (3.15)

e sempre solucao da equacao da onda (3.4).(b) Mais do que isto, vamos mostrar ao leitor agora como provar, usando um metodo devido a

D’Alembert, que qualquer solucao da equacao da onda deve ser necessariamente da forma (3.15) parafuncoes convenientes φ e ψ. Inicialmente, introduza novas variaveis

ξ = x− c t (3.16)

η = x + c t . (3.17)

Usando a regra da cadeia e supondo que as derivadas parciais de 2a. ordem de u sao contınuas, proveque em termos das novas variaveis tem-se

ut = −c (uξ − uη) ,

ux = uξ + uη ,

utt = c2 (uξξ + uηη − 2uξη) ,

uxx = uξξ + uηη + 2uξη .

(c) Substituindo as expressoes acima na equacao da onda, prove que nas novas variaveis estatoma a forma

uξη = 0 .

(d) Prove, integrando separadamente com relacao a cada uma das variaveis, que as unicas solucoesda equacao acima sao da forma φ(ξ) + ψ(η) para funcoes arbitrarias φ e ψ.

(e) Observe usando (3.3) que o parametro c na equacao da corda vibrante possui dimensao develocidade.

30

(f) Definindo uma onda como uma forma que viaja pelo espaco com uma velocidade definida semalterar seu perfil, mostre que a solucao geral (3.15) pode ser interpretada como uma soma de duasondas de perfis arbitrarios φ e ψ, uma viajando para a direita com velocidade c e a outra viajandopara a esquerda com a mesma velocidade.

3.2 Considere um problema de Dirichlet homogeneo, ou seja, (3.6) com a(t) = b(t) ≡ 0.(a) Prove que se us(x, t) e solucao do problema homogeneo com f(t) ≡ 0 e se uc(x, t) e solucao

do mesmo problema com g(t) ≡ 0, entao us(x, t) + uc(x, t) e solucao do problema homogeneo com fe g quaisquer.

(b) Prove que o conjunto das funcoes us(x, t) que satisfaz a equacao da onda, as condicoes ho-mogeneas de Dirichlet e a condicao inicial u(x, 0) = 0 e um espaco vetorial. Da mesma forma,prove que o conjunto das funcoes uc(x, t) que satisfaz a equacao da onda, as condicoes homogeneasde Dirichlet e a condicao inicial ut(x, 0) = 0 e um espaco vetorial.

(c) Calcule as derivadas e mostre diretamente que as funcoes

un,s(x, t) = sennπct

lsen

nπx

l,

n = 1, 2, . . . sao solucoes da equacao da onda e satisfazem as condicoes de Dirichlet homogeneas etambem a condicao inicial u(x, 0) = 0.

(d) Analogamente, mostre que as funcoes

un,c(x, t) = cosnπct

lsen

nπx

l,

n = 1, 2, . . . sao solucoes da equacao da onda e satisfazem as condicoes de Dirichlet homogeneas etambem a condicao inicial ut(x, 0) = 0.

(e) Mostre que tanto o conjunto das funcoes em (c) quanto o das funcoes em (d) sao linearmenteindependentes. Portanto, o conjunto das solucoes da equacao da onda satisfazendo as condicoes deDirichlet homogeneas e um espaco vetorial de dimensao “duplamente”infinita.

Sugestao: Aplicar a sugestao do problema 2.4 separadamente a cada conjunto de solucoes.

3.3 A corda dedilhada.Suponha que o ponto onde um violonista dedilha a corda esteja a aproximadamente 1/3 de seu

comprimento a partir da extremidade esquerda. Uma boa aproximacao para as condicoes iniciaisusadas por um violonista e entao, na notacao de (3.6),

f(x) =

3x , se 0 ≤ x ≤ 1

332(1− x) , se 1

3< x ≤ 1

eg(x) ≡ 0 .

Resolva numericamente o problema de Dirichlet correspondente supondo c = l = 1 e valores apro-priados para ∆x e ∆t. Nas unidades aqui usadas, qual e aproximadamente a frequencia do somproduzido pela corda?

3.4 A corda de piano.Em um piano o som e produzido por cordas fixas em ambas as extremidades, como no violao.

Uma grande diferenca e que, em vez de dedilhadas como no problema 3.3, as cordas de um pianosao percutidas por um martelo. Supondo que o martelo possua 1/10 do comprimento de uma corda e

31

seu centro atinja-a num ponto a distancia 2/5 de seucomprimento a partir da extremidade esquerda,entao, na notacao de (3.6), boas aproximacoes para as condicoes iniciais sao

f(x) ≡ 0

e

g(x) =

1 , se 3

20≤ x ≤ 1

4

0 , se 0 ≤ x < 320

ou 14

< x ≤ 1.

Resolva numericamente o problema de Dirichlet correspondente supondo c = l = 1 e valores apro-priados para ∆x e ∆t. Nas unidades aqui usadas, qual e aproximadamente a frequencia do somproduzido pela corda?

3.5 RessonanciaVamos considerar aqui uma corda com c = l = 1. Em vez de extremos fixados, consideremos uma

corda em que o extremo x = 1 e fixado, enquanto o extremo x = 0 e movido periodicamente por umagente externo. Considere que o movimento do extremo esquerdo e dado por uma condicao do tipou(0, t) = a(t).

(a) Supondo que a corda inicialmente esta em repouso em sua posicao de equilıbrio, especifiquetodas as condicoes iniciais e de contorno a serem aplicadas.

(b) Suponha que a(t) = sen (πt). Calcule o movimento da corda ate t = 10, utilizando ∆x e ∆tapropriados.

(c) Repita o mesmo, agora com a(t) = 2sen t. Compare a amplitude do movimento da corda aofinal do intervalo de tempo nos dois casos.

(d) Repita ainda para a(t) = sen (nπt), n = 2, 3, . . ..(e) Observe que movimentos externos com frequencias apropriadas podem causar movimentos na

corda de grandes amplitudes, enquanto outras frequencias nao produzem tal efeito. Procure relacionareste fenomeno com o da ressonancia de um sistema massa-mola sob acao de uma forca externaperiodica, ver por exemplo a secao 3.9 de [1]

3.6 Modos normais ou ondas estacionarias em uma corda com entremos fixados.Suponha por simplicidade que l = 1 e considere uma corda discretizada com passo ∆x = l

Ne

extremos fixados.(a) Mostre que se, na notacao de (3.6) tivermos f(x) = sen (nπx) e g(x) ≡ 0, entao teremos, se

j = 0, 1, 2, . . .,

ui,j = Aj sennπi

N,

onde Aj satisfaz a equacao de diferencas

Aj+1 = [2− 2s(1− cosnπ

N)] Aj − Aj−1

com condicoes iniciais A0 = 1 e A1 = 1− s(1− cos nπN

).(b) Escolha valores para s, n e N e calcule numericamente os primeiros valores para Aj. Faca

um grafico destes valores.O resultado da letra (a) mostra que com as condicoes iniciais dadas, em cada instante de tempo

discreto a forma da corda e a mesma forma inicial com uma “amplitude”Aj dependente do tempoj. Na letra (b) voce deve ter visto que a amplitude oscila no tempo se s obedece a condicao deestabilidade 0 < s ≤ 1. Observe ainda o aparecimento de “problemas”quando s > 1.

Os movimentos da corda descritos neste problema, em que a corda preserva sua forma todo otempo tambem existem no caso contınuo; sao as funcoes un,c que aparecem no problema 3.2. Descubraquais condicoes iniciais devem ser usadas para que aparecam os analogos discretos das funcoes un,s

no mesmo problema?

32

Capıtulo 4

A equacao de Laplace

4.1 Aplicacoes da equacao de Laplace

Iremos agora tratar uma terceira EDP, a qual aparece em varias aplicacoes. Selecionamos duas dasaplicacoes em que a equacao de Laplace e importante. Conforme exemplificado pelas duas aplicacoesa seguir, a equacao de Laplace modela fenomenos estaticos, independentes do tempo. Desta forma,nao e de se esperar que os problemas relacionados possuam condicoes iniciais, ao contrario do queaconteceu nos casos das equacoes da difusao e da onda. Os problemas modelados pela equacao deLaplace sao puramente problemas de valores de contorno e a diferenca entre os varios problemas estanos diferentes tipos de condicoes de contorno.

4.1.1 Potencial eletrostatico

Pela lei de Coulomb, o campo eletrico no ponto ~x ∈ R3 causado por uma carga q localizada numponto ~x0 e dado por, a menos de uma constante multiplcativa,

~E(~x) =q (~x− ~x0)

||~x− ~x0||3 . (4.1)

Um calculo simples mostra que o rotacional do campo ~E nos pontos ~x 6= ~x0 e nulo.Por definicao, o campo eletrico no ponto ~x causado por um conjunto estatico de cargas estatico

e a soma dos campos eletricos causados por cada uma das cargas. Se C e uma curva fechada em R3

que nao passa sobre nenhuma das cargas e S e alguma superfıcie que tenha C como bordo, o teoremade Stokes mostra que ∫

C

~E · d~r =

∫

S

rot ~E · d ~A = 0 ,

e portanto o campo ~E e conservativo. Daı existe uma funcao φ, chamada potencial eletrostatico talque

~E = − grad φ . (4.2)

Em geral, o objetivo dos problemas de eltrostatica e encontrar, dadas informacoes sobre as cargasque o produzem, calcular o campo eletrico ~E. O potencial eltrostatico e portanto uma funcao auxiliarescalar que, por gradiente (ou seja, derivacao) permite obter o vetor ~E. Como veremos adiante, umaoutra vantagem de se trabalhar com o potencial eletrostatico em vez do campo eletrico, e que opotencial em muitos casos importantes pode ser obtido como solucao de uma EDP bem conhecida,a equacao de Laplace. Vejamos como esta EDP aparece.

33

Calculando o divergente do campo ~E dado em (4.1), ve-se facilmente que este e nulo, a nao ser em~x0, onde este nao e definido. Se S e uma superfıcie fechada orientada com o vetor normal apontandopara o exterior e q a carga total contida no interior de S, pode-se aplicar o teorema de Gauss dadivergencia, tendo-se o cuidado de isolar as singularidades de ~E com pequenas bolas centradas nospontos onde ha cargas. Obtem-se daı que

∫

S

~E · d ~A = 4π q , (4.3)

resultado conhecido com lei de Gauss da eletrostatica. Usando-se novamente o teorema da di-vergencia, chega-se a chamada forma diferencial da lei de Gauss da eletrostatica

div ~E = 4π ρ(~x) , (4.4)

onde ρ(~x) e a densidade de carga no ponto ~x.Substituindo-se (4.2) em (4.4), tem-se finalmente que o potencial eltrostatico e solucao da equacao

de Poisson∆φ = −4πρ , (4.5)

onde

∆φ = div grad φ =∂2φ

∂x2+

∂2φ

∂y2+

∂2φ

∂z2

e chamado o laplaciano da funcao φ. 1

Se numa regiao de R3 sabemos que nao existem cargas (portanto ρ ≡ 0), a equacao de Poissonreduz-se a equacao de Laplace

∆φ = 0 . (4.6)

E claro que se se conhece as posicoes de todas as cargas envolvidas em um problema, o campoeletrico pode ser encontrado somando-se (ou integrando-se) (4.1) sobre todas as cargas. A equacaode Laplace e um instrumento util para o calculo do campo eletrico, atraves do calculo do potencial,em problemas onde nao se conhecem as posicoes das cargas, mas outras informacoes. Ha 3 problemastıpicos:

(a) Problema de Dirichlet

Seja Ω ⊂ R3 regiao fechada e limitada que nao contenha cargas em seu interior int Ω. Em muitoscasos, o potencial φ e conhecido nos pontos da fronteira ∂Ω de Ω, por exemplo φ(~x) = f(~x),~x ∈ ∂Ω, onde f e funcao conhecida. O problema e entao encontrar φ resolvendo-se a equacaode Laplace com condicao de Dirichlet na fronteira:

∆φ(~x) = 0, ~x ∈ intΩφ(~x) = f(~x), ~x ∈ ∂Ω

. (4.7)

(b) Problema de Neumann

Em alguns problemas conhece-se o valor de ~E em ∂Ω e quer-se conhece-lo em int Ω, ondeΩ ⊂ R3 e regiao que nao contem cargas. Na verdade, basta conhecer a componente de ~E

1Cabe aqui lembrar que neste texto estamos usando o sımbolo ∆ com dois significados completamente distintos.Ate aqui vınhamos usando quase sempre este sımbolo para denotar o passo de discretizacao de uma variavel. Agorapassamos a usa-lo tambem para denotar o operador laplaciano. Uma vez que, devido aos contextos bastante distintos,ha pouca possibilidade de confusao entre os dois significados, continuaremos a usar esta notacao ambıgua, poremtradicional.

34

normal a superfıcie ∂Ω, ou seja grad φ · ~n, onde ~n e um vetor unitario normal a ∂Ω. Osproblemas de Neumann sao portanto da forma

∆φ(~x) = 0, ~x ∈ intΩgrad φ(~x) · ~n = g(~x), ~x ∈ ∂Ω

, (4.8)

onde g(~x) e uma funcao conhecida sobre ∂Ω.

(c) Problemas mistos

Podem aparecer problemas em que se quer encontrar o potencial φ em intΩ conhecendo-se ovalor de φ em alguns pontos da fronteira (condicao de Dirichlet) e o valor de grad φ(~x) · ~n nosdemais pontos (condicao de Neumann).

4.1.2 Solucoes estacionarias da equacao de difusao

No problema 2.1 o leitor deve ter mostrado que a versao multidimensional da equacao da difusao(2.6) e

∂u

∂t= k ∆u .

Em muitos casos e natural esperar-se que esta equacao possua solucoes estacionarias, ou seja,independentes do tempo. Como o primeiro membro na equacao acima se anula neste caso, entao assolucoes estacionarias da equacao da difusao sao solucoes da equacao de Laplace.

Vimos o aparecimento de solucoes estacionarias da equacao da difusao unidimensional no capıtulo2 por exemplo quando as condicoes de contorno sao as homogeneas de Dirichlet. Estas condicoesrepresentam fisicamente o caso em que nas extremidades do tubo estao colocados reservatorios devolume infinito com solucao a concentracao nula. Neste caso a concentracao em todos os pontos dotubo tende a zero quando t →∞. A funcao nula e uma solucao estacionaria da equacao de Laplaceunidimensional u′′(x) = 0. A convergencia da concentracao nos pontos do tubo para a solucaoestacionaria esta ilustrada em um exemplo na figura 2.2.

Solucoes estacionarias da equacao da difusao unidimensional possuem pouco interesse, uma vezque ela e exatamente soluvel de forma banal. Por isto e que nos interessaremos por solucoes esta-cionarias da equacao da difusao em dimensao maior que 1. Por simplicidade, restringir-nos-emos adimensao 2.

Elencamos a seguir os tipos de problemas para os quais e procurar-se solucoes estacionarias daequacao da difusao.

(a) Problemas de Dirichlet

No caso bidimensional, um problema de Dirichlet para a equacao da difusao especifica a con-centracao da solucao em todos os pontos de uma curva fechada fronteira de uma regiao Ω ⊂ R2.Gostarıamos de encontrar a concentracao da solucao nos pontos de int Ω dada a concentracaoinicial nestes pontos.

De um ponto de vista fısico, podemos pensar que a cada ponto da fronteira esta ligado umreservatorio de volume infinito, cada reservatorio contendo solucao a uma determinada con-centracao. Caso a concentracao nos reservatorios seja dependente do tempo, e pouco razoavelque a concentracao da solucao em int Ω atinja um estado estacionario. Por outro lado, se aconcentracao em cada um dos reservatorios for constante, e de se esperar que haja fluxo desoluto entre os reservatorios e int Ω e uma solucao estacionaria seja atingida quando t → ∞.Isto e o que acontece por exemplo no caso unidimensional quando os reservatorios possuem

35

solucao a concentracao nula e a concentracao em todos os pontos do tubo tende a zero quandot →∞.

(b) Problemas de Neumann

O leitor pode adaptar os argumentos na subsecao 2.2.2 e provar que o fluxo de massa, umvetor nos casos de dimensao maior que um, e proporcional ao gradiente da concentracao. Nosproblemas de Neumann, em vez de especificar a concentracao da solucao nos pontos da fronteira,o que se especifica e a componente normal a fronteira do fluxo de soluto em cada ponto desta.

Em alguns casos, torna-se claro que problemas de Neumann para a equacao de Laplace podemnao ter solucao. Por exemplo, se a componente normal do fluxo de massa for de mesmo sinalem todos os pontos da fronteira, entao temos soluto entrando em todos os pontos ou saindoem todos os pontos. Em ambos os casos nao e de se esperar que a equacao da difusao possuasolucao estacionaria, uma vez que a concentracao deve ou crescer ou decrescer todo o tempo.

(c) Problemas mistos Nestes deseja-se encontrar, caso exista, uma solucao estacionaria para aequacao da difusao especificando-se a concentracao em alguns pontos da fronteira e a compo-nente normal a fronteira do fluxo de massa em outros pontos.

Ou seja, os problemas interessantes de se estudar para a equacao de Laplace sao exatamente osmesmos em ambas as aplicacoes.

4.2 A equacao de Laplace discretizada

Iremos estudar aqui a questao de existencia e unicidade de solucoes para os problemas acima, poremdo ponto de vista da equacao discretizada. Boa parte do que aqui sera feito pode ser utilizado paraprovar teoremas sobre a equacao de Laplace sem discretizacao. Uma boa abordagem deste topicopode ser vista em [4].

4.2.1 O laplaciano discreto

Por simplicidade, vamos considerar daqui em diante problemas bidimensionais, sendo x e y, comousual, as coordenadas cartesianas. Os resultados que virao a seguir sao mais simples se discretizarmosambas as coordenadas com o mesmo passo de discretizacao h = ∆x = ∆y. Usando a discretizacaousual (1.7) para a derivada segunda, temos

∆u(ih, jh) ≈ ∆hui,j =1

h2(ui+1,j + ui−1,j + ui,j+1 + ui,j−1 − 4ui,j) , (4.9)

equacao esta que define o laplaciano discreto ∆hu.Podemos agora provar dois resultados que valem tambem para a equacao nao discretizada. Em

contraste com o caso nao-discreto, a prova destes resultados para a equacao discreta e quase trivial.

Teorema 4 (da media) Se u e uma solucao da equacao de Laplace discreta ∆hu = 0, entao ui,j e amedia aritmetica dos valores de u nos pontos vizinhos mais proximos da rede de discretizacao.

Prova E so igualar a zero o lado direito de (4.9).

36

Teorema 5 (do maximo e do mınimo) Seja Ω um polıgono cujos lados sao paralelos aos eixoscartesianos e com vertices sobre os pontos da rede de discretizacao e u uma solucao da equacao deLaplace discretizada definida em Ω. Entao tanto o maximo quanto o mınimo de u em Ω ocorrem emsua fronteira. Alem do mais, o maximo ou o mınimo ocorrerao tambem em um ponto interior de Ωsomente se u for constante.

Prova Vamos provar o teorema para o caso do maximo. A prova para o mınimo e analoga.Suponha que o maximo de u ocorra em um ponto do interior de Ω. Como este valor maximo e,

pelo teorema 4, a media aritmetica dos valores vizinhos, nenhum destes valores vizinhos pode sermenor que o valor maximo. Temos assim um conjunto de 5 pontos vizinhos tais que em todos eleso valor de u e constante igual ao valor maximo: o ponto original e seus 4 vizinhos. Repetindo oargumento da media para cada um dos 4 vizinhos que nao esteja na fronteira, apos um numero finitode passos teremos mostrado que u e constante.

4.2.2 Problemas de Dirichlet

Quando discutimos problemas de Dirichlet no contexto de solucoes estacionarias da equacao da di-fusao, afirmamos que e fisicamente intuitivo que problemas de Dirichlet para a equacao de Laplacepossuam solucao e que esta seja unica. Podemos provar esta afirmativa no caso da equacao dis-cretizada. Em [4] prova-se ainda que o resultado vale tambem no caso da equacao nao-discretizadae que a solucao do problema e o limite das solucoes discretas quando h → 0.

Teorema 6 Seja Ω polıgono como no enunciado do teorema 5 e f funcao definida nos pontos darede de discretizacao sobre a fronteira ∂Ω de Ω. Entao o problema de Dirichlet discreto

∆hui,j = 0, (ih, jh) ∈ int Ωui,j = f(ih, jh), (ih, jh) ∈ int ∂Ω

(4.10)

possui solucao e esta solucao e unica.

Prova Note primeiramente que para cada par (i, j) a equacao ∆hui,j = 0 e uma equacao algebricalinear envolvendo 4 valores de u. Os valores de u em pontos interiores sao incognitas, enquantoaqueles valores de u nos pontos de ∂Ω sao especificados pela condica ode contorno de Dirichlet.Portanto a questao de resolver (4.10) e reduzida a uma questao de sistemas de equacoes algebricaslineares da forma AX = B, onde A e matriz n×n, X e B matrizes n×1 e n, o numero de incognitas,e simplesmente o numero de pontos interiores de Ω sobre a rede de discretizacao.

O leitor verificara facilmente atraves de exemplos que os valores na matriz B sao simplesmenteos valores da funcao f com os sinais trocados. Em particular, no caso de condicao de Dirichlethomogenea f ≡ 0, tem-se B = 0 e o sistema de equacoes e homogeneo, ou seja AX = 0. A matriz Ae independente da funcao f .

Sabemos que o problema (4.10) tera solucao unica se e somente det A 6= 0. Para provar quedet A 6= 0, basta entao provar que o problema de Dirichlet homogeneo somente possui a solucaotrivial u ≡ 0. Isto porem e consequencia facil do teorema 5.

De fato, seja u solucao do problema homogeneo. Como u se anula na fronteira e possui maximona fronteira, entao ui,j ≤ 0 para todo (ih, jh) ∈ int Ω. Mas como o mınimo de u tambem esta nafronteira, entao ui,j ≥ 0 para todo (ih, jh) ∈ int Ω. Logo u ≡ 0.

37

4.2.3 Problemas de Neumann

Seja Ω um polıgono como o descrito no enunciado do teorema 5. Seja ni,j o vetor unitario normal ao

lado de Ω no ponto (ih, jh) ∈ ∂Ω apontando para o exterior de Ω. E natural definir-se a derivadanormal de u no ponto (ih, jh) como

dnui,j =1

h(ui,j − u(i,j)−ni,j

) , (4.11)

onde (i, j)− ni,j denota o ponto da rede de discretizacao mais proximo do ponto (ih, jh) alcancadoquando se caminha no sentido contrario ao de ni,j. O analogo discreto de uma condicao de Neumannseria portanto algo da forma

dnui,j = f(ih, jh) ,

onde f e uma funcao definida sobre a fronteira de Ω.Uma primeira propriedade de solucoes de problemas de Neumann para a equacao de Laplace,

igualmente facil de se provar em problemas discretos ou nao, e que se c e uma constante qualquer eu uma solucao, entao v = u + c e solucao do mesmo problema. De fato, a soma de uma constantenao altera a nulidade do laplaciano, seja ele discreto ou contınuo. Tampouco v deixa de obedecer ascondicoes de Neumann, pois elas dependem das derivadas da funcao (ou diferencas entre valores dafuncao em pontos vizinhos), ambas inalteradas com a soma de uma constante. Portanto, quando umproblema de Neumann para a equacao de Laplace possui solucao, esta nunca e unica.

Considere agora uma regiao limitada Ω e suponha que o vetor normal na fronteira de Ω estejaapontando para o exterior. Consideracoes fısicas nos dao uma ideia de uma condicao necessaria paraa existencia de solucao para problemas de Neumann para a equacao de Laplace.

Pensando em termos de solucoes estacionarias de problemas de difusao, e razoavel esperar-se queuma tal solucao somente seja possıvel se o fluxo total de soluto atraves da fronteira for nulo. Casoeste fluxo total seja negativo (mais soluto entrando em Ω do que saindo), e de se esperar que asconcentracoes crescam indefinidamente e, caso positivo, que estas decrescam ate que nao haja maissoluto na regiao e a condicao de fluxo total positivo nao possa mais se manter por falta de soluto.Como o fluxo normal a fronteira em um ponto desta e proporcional a grad u · n, a condicao e que aintegral de fluxo de grad u pela fronteira seja nula.

Pensando em termos eletrostaticos, para que o potencial em Ω seja solucao da equacao de Laplacee necessario que a carga total contida em Ω seja nula. Neste caso, a lei de Gauss (4.3) obriga que ofluxo do campo eletrico pela fronteira de Ω seja nulo. Levando em conta que o campo eletrico estarelacionado ao potencial por (4.2), temos novamente a mesma condicao.

No caso discreto, podemos provar facilmente este resultado. E interessante notar que na provaabaixo a identidade (4.14) e uma especie de versao discreta do teorema de Green da Analise Vetorial.

Teorema 7 Seja Ω polıgono como no enunciado do teorema 5 e f funcao definida nos pontos darede de discretizacao sobre a fronteira ∂Ω de Ω. Entao o problema de Neumann discreto

∆hui,j = 0, (ih, jh) ∈ int Ωdnui,j = f(ih, jh), (ih, jh) ∈ int ∂Ω

(4.12)

possui solucao somente se ∑

(ih,jh)∈∂Ω

f(ih, jh) = 0 .

38

Prova Da definicao (4.9) do laplaciano discreto, tem-se que

∆hui,j = − 1

h2

∑v

(ui,j − u(i,j)−v) , (4.13)

onde os vetores v sobre os quais se realiza a soma sao os quatro vetores unitarios paralelos aos eixoscoordenados e a notacao u(i,j)−v foi explicada em (4.11).

Vamos somar a equacao acima sobre todos os pares (i, j) tais que (ih, jh) ∈ int Ω. Sempre quetanto o ponto (ih, jh) quanto seu vizinho (ih, jh)+hv estiverem em int Ω, as parcelas correspondentesa diferenca dos valores de u nestes pontos irao anular uma a outra. Logo, quando realizarmos a somasobre todos os (ih, jh) ∈ int Ω de ∆hui,j, as unicas parcelas que sobrarao ao final sao as relativas adiferencas entre os valores de u em pontos no interior e seus vizinhos na fronteira. Ou melhor, temos

∑

(ih,jh)∈intΩ

∆hui,j =1

h

∑

(ih,jh)∈∂Ω

dnui,j . (4.14)

Se u e solucao de (4.12), entao o lado esquerdo na equacao acima e nulo e o lado direito e

1

h

∑

(ih,jh)∈∂Ω

f(ih, jh) ,

provando o resultado.

Problemas do capıtulo 4

4.1 (a) Resolva numericamente o problema de Dirichlet para a equacao de Laplace no quadrado[0, 1]× [0, 1] onde a funcao incognita obedece a condicao de ser nula nos lados y = 0, x = 0 e x = 1 evale 1 no lado y = 1. Use h = 1/4. Observe que, apesar do valor relativamente grande para h, vocetem que resolver um sistema de 9 equacoes lineares com 9 incognitas!

(b) A solucao exata do problema acima encontrada pelo metodo de separacao de variaveis e

u(x, y) =4

π

∞∑n=1

1

(2n− 1) senh ((2n− 1)π)sen ((2n− 1)πx) senh ((2n− 1)πy) .

Compare a solucao aproximada obtida em (a) com a soma de um bom numero de termos da solucaoexata. Faca graficos comparativos da solucao exata para valores fixos de y e dos valores correspon-dentes da solucao aproximada.

(c) Repita as letras (a) e (b) para h = 1/10.(d) Voce percebe como metodos eficientes para a resolucao de sistemas de equacoes lineares sao

importantes?

4.2 Considere o sistema de equacoes lineares AX = B obtido quando se discretiza um problema deDirichlet para a equacao de Laplace. Mostre que em cada linha de A aparecem no maximo 5 entradasnao-nulas. Em outras palavras, repare que quando a dimensao da matriz A e grande, entao a imensamaioria de seus elementos e nula. Uma matriz deste tipo e chamada esparsa. Procure conhecer emlivros de Algebra Linear Numerica algoritmos especiais para a solucao de sistemas lineares em que amatriz dos coeficientes e esparsa.

39

Referencias Bibliograficas

[1] W. E. Boyce and R. C. DiPrima. Equacoes Diferenciais Elementares e Problemas de Valores deContorno. LTC Editora, Rio de Janeiro, 6a. edicao, 1998.

[2] R. V. Churchill. Series de Fourier e Problemas de Valores de Contorno. Guanabara Dois, Riode Janeiro, 2a. edicao, 1978.

[3] D. G. Figueiredo. Analise de Fourier e Equacoes Diferenciais Parciais. Instituto de MatematicaPura e Aplicada, Rio de Janeiro, 1977.

[4] I. G. Petrovsky. Lectures on Partial Differential Equations. Dover, New York, 1991.

[5] R. D. Richtmyer and K. W. Morton. Difference Methods for Initial-Value Problems. PrenticeHall, Englewood Cliffs, 2a. edicao, 1967.

[6] W. A. Strauss. Partial Differential Equations - An Introduction. John Wiley and Sons, NewYork, 1992.

40