Introducao a`s 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 EDPs 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 princpio do maximo e do mnimo . . . . . . . . . . . . . . . . . 17

3 A equacao da onda 223.1 Modelagem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 223.2 Problemas tpicos 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, EDPs sao estudadas a nvel basico apos um primeiro curso sobre Equacoes DiferenciaisOrdinarias (EDOs), 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 EDOs lineares, por exemplo ao nvelde [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 diferentea 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 EDPs. Por que os metodos para resolver EDPs sao tao diferentes dos metodospara resolver EDOs? Por que aparecem as series de Fourier em EDPs e nao em EDOs? De ondevieram estas tais condicoes de contorno que nao apareciam antes em EDOs? Por que usualmenteas solucoes de EDPs aparecem como series? Por que nao existe para EDPs uma teoria geral deexistencia e unicidade tao poderosa quanto a correpondente para EDOs?

Nao consegui responder a estas perguntas naquela epoca. Embora as respostas sejam simplese estejam espalhadas pelas presentes notas, a maior parte dos livros sobre EDPs nao as fornece.Quando, como professor, comecei a ensinar EDPs 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 EDPs.Estes metodos, assim como o presente texto, partem de uma abordagem discreta a`s EDPs. Conformeacima explicado, meu objetivo principal nao e o de ensinar metodos numericos para EDPs e simensinar EDPs. Porem, ensinar EDPs da maneira presente leva-nos tao perto de aprender algo sobremetodos numericos que seria um grande desperdcio 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 EDPs 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 cobaiasde 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 OCarroll e Giovanni Gallavotti, que em estagios diversos de minha formacao me ofereceramexemplos da utilidade de discretizar sistemas contnuos. 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

Captulo 1

Discretizacao

1.1 EDPs 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=

2ux2

, 2ut2

= u + 2ux2

ou cos u ut

= 2u

xtsao exemplos de EDPs.

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 comoEDOs.

Antes de introduzir as primeiras EDPs 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 EDOs: a discretizacao de funcoes.

Considere uma funcao real f definida em um intervalo [a, b]. Se f e contnua, entao para dar umadescricao bastante boa do comportamento de f nao e necessario especificar o valor de f em todosos pontos do seu domnio; basta especificar este valor em um conjunto X = {x0, x1, x2, . . . , xN} depontos do domnio, onde N e um inteiro positivo suficientemente grande e os elementos razoavelmentebem espalhadospor 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 conjuntocontnuo de informacoes (todos os valores de f) por um conjunto finito (o conjunto dos valores{f(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 domnio 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 = ba

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

xi = a + ix, 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, pi]. Toda a informacao sobre osvalores de f no intervalo considerado esta contida em seu grafico. Podemos discretizar a funcao futilizando um passo x = pi

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

pi 5

2 pi 5

3 pi 5

4 pi 5

pi

0.2

0.4

0.6

0.8

1

Figura 1.1: Uma discretizacao do seno.

pi 5

2 pi 5

3 pi 5

4 pi 5

pi

-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 = pi

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| < Mx , (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(xi1)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) + 12!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 Os, 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 = xi1, x = xi xi1, 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) 12

[f(xi+1) f(xi)

x+

f(xi) f(xi1)x

]=

f(xi+1) f(xi1)2x

. (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(xi1) = f(xi) f (xi)x + 12f (xi) (x)2 1

6f (3)(2) (x)

3 .

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

f (xi) =f(xi+1) f(xi1)

2x 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(xi1)

x

]=

f(xi+1) + f(xi1) 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 captulo 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 possvel que coincidisseexatamente com o seno nos pontos da discretizacao utilizada e tambem tal que p(pi

2) = 1

2. Encontre

a formula para p(x).

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

4, pi5, pi6, pi7, pi8. 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 concludo 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 = xi1, c = xi e n = 3.Some os resultados.

7

Captulo 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 espalhamentodo soluto e chamado difusao.

Microscopicamente, a difusao e devida ao movimento termico desordenado das moleculas dosoluto. Gostaramos 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 construdasa partir da explicacao microscopica qualitativa acima.

Descrevemos inicialmente a situacao fsica.

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 impossvel 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 contnuas, 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(ix, jt) e fi em vez de f(ix), 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 j

em vez de posicao ix e instante jt.A discretizacao do tempo significa que desprezaremos a passagem contnua do tempo e olharemos

para 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 = 2t, 3t, . . ., 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 stios de tamanho x iniciandoem x = 0, x = x, x = 2x, . . . , x = (N 1)x = l x e especificaremos somente o ndicei = 0, 1, . . . , N 1 do stio em que se localiza cada molecula.

A cada passagem do tempo discreto, havera uma redistribuicao das moleculas do soluto entreos stios. Dada uma regra de passagem do soluto de um stio 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 stio para um dos stios imediatamentevizinhos a` direita ou a` esquerda. Portanto, uma molecula que no instante j esteja no stio i, noproximo instante discreto, j + 1, ou permanecera onde estava, ou entao passara para um dentreos stios i 1, ou i + 1. Como a passagem de soluto de um stio para o outro e consequencia domovimento termico desordenado, entao a passagem do stio i para o stio vizinho a` direita i + 1 etao provavel quanto a passagem para o stio vizinho a` esquerda i 1. Alem do mais, a massa desoluto que passa do stio i para i 1 e proporcional a` massa de soluto no proprio stio i e portantoproporcional a` concentracao em i. A massa de soluto que passa de i para i1 sera ainda proporcionalao intervalo temporal t, a` area A da secao transversal do tubo e inversamente proporcional a x.

Em resumo, se mii1,j representa a massa de soluto que no instante j passa do stio i parai+1 (e tambem para i 1), entao a nossa hipotese fundamental sobre a qual esta baseada a equacaoda difusao e

mii1,j =k At

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 stios vizinhos, assim como acrescentara massa que veio dos stios vizinhos i 1:

mi,j+1 = mij (mii+1,j +mii1,j) + mi+1i,j +mi1i,j .Usando agora (2.1), temos

mi,j+1 = mij k Atx

2ui,j +k At

x(ui+1,j + ui1,j) .

Para fazermos aparecer concentracoes no lugar de massas na ultima equacao, dividimo-la pelo volumeAx de cada stio, obtendo

ui,j+1 = ui,j +kt

(x)2(ui+1,j + ui1,j 2ui,j) (2.2)

= (1 2s)ui,j + s (ui+1,j + ui1,j) , (2.3)onde o paramtero de estabilidade s > 0 e definido por

s = kt

(x)2(2.4)

9

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

ui,j+1 ui,jt

= kui+1,j + ui1,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 fsica, uma vez que o tempo e o espaco sao contnuos. A continuidade do espacoe do tempo sao recuperadas ao fazermos o limite do contnuo x 0 e t 0. Neste limite, deacordo com os resultados da secao anterior sobre a discretizacao de derivadas, temos

ui,j+1 ui,jt

ut

eui+1,j + ui1,j 2ui,j

(x)2

2u

x2.

No limite contnuo, (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 partculas (moleculas). No incio 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 raciocnio contnuo 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 EDOs, quase sempre problemas de EDPs 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 fsico 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 realsticos 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 fsica 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

l0f(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

(ix, jt) =mii+1,j mi+1i,j

t

= kA ui+1,j ui,jx

limite contnuo kA ux

(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 realsticos 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 substitudas 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 sada 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:

duidt

= k(x)2

(ui+1(t) + ui1(t) 2ui(t)), i = 1, 2, . . . , N 1ui(0) = f(ix), i = 1, 2, . . . , N 1u0(t) = 0, t > 0uN(t) = 0, t > 0

, (2.10)

onde ui(t) = u(ix, 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 EDOs de 1a. ordem. Cada uma das N 1 EDOsescritas na primeira linha de (2.10) envolve alem da concentracao ui que comparece no lado direito

12

com sua derivada, as concentracoes nos stios vizinhos ui1 e ui+1. As N 1 equacoes na 2a. linhasao as condicoes iniciais para cada uma das EDOs 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 EDOs. 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 EDOs. Observeporem que o numero de EDOs 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 EDOs.

Alem do mais, as EDOs em (2.10) sao todas lineares e, por causa da homogeneidade das condicoesde Dirichlet, homogeneas. Combinacoes lineares de solucoes deste sistema de EDOs 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(ix),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 EDPs. Paramos aqui com a compreensao de que estes problemas possuemanalogia com sistemas de infinitas EDOs 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 EDOs deve ter ouvido que nem todas as EDOs 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 EDPs possivelmente se encontram emsituacao pior que a das EDOs. Para resolver problemas que nao possam ser resolvidos de maneiraexata existem tambem metodos numericos para EDPs. 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 EDPs dometodo de Euler para EDOs.

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 fsica 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 a`s 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 = 14iremos escolher, com base em nossa experiencia, t = 1

40. O

parametro 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(ix), 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 = 2t = 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 fsica 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 fsico. 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. Usamosx = 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 possvel 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 princpio do maximo e do mnimo

Citamos sem demonstracao o resultado seguinte sobre a equacao da difusao, conhecido como Princpiodo maximo e do mnimo. A demonstracao pode ser encontrada em muitos dos livros padrao sobreEDPs, por exemplo [3] ou [6]. A demonstracao e simples, usando somente a teoria basica de maximose mnimos de funcoes de varias variaveis. O princpio do maximo e do mnimo 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 x2}e

l3 = {(x2, t) R2; t1 t t2} .Entao os pontos de maximo e de mnimo de u em R encontram-se em l1 l2 l3.

Alem da utilidade matematica, o teorema acima possui ainda uma explicacao fsica 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. Raciocniosemelhante vale para o mnimo.

Conclumos portanto que o princpio do maximo e do mnimo reflete uma propriedade fsica 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 mnimo 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 mnimo da solucao de (2.3)em R encontram-se em l1 l2 l3.

Sejam I1x = x1 e I2x = 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, J1t = t1, J2t = t2.Defina finalmente Mi,j = max{ui,j, ui+1,j, ui1,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 + ui1,j) (1 2s)Mi,j + s (Mi,j +Mi,j) = Mi,j . (2.11)

Defina agoraMj = max

i{I1,I1+1,...,I2}ui,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 mnimo e analogo.

Problemas do captulo 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 desprezvel. 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 substitudapor uma condicao nao-homogenea, entao uma combinacao linear de solucoes do sistema de EDOsnao e necessariamente solucao do sistema de EDOs.

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

un(x, t) = ek n2pi2

l2t sen

npix

l

sao solucoes da equacao da difusao, onde n = 1, 2, . . ..(b) Verifique ainda que os un acima obedecem a`s 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 a`s 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 l0

f(x)g(x)dx .

(e) Conclua que o conjunto de todas as solucoes da equacao da difusao e que tambem obedecama`s 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 en2pi2t sennpix ,

onde

an =2

npi

(cos

2npi

5 cos 4npi

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 = 2x1,

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 espalhamentoda 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 desprezvel.

2.8 Considere o problema de Neumannut = k uxx, x (0, l),t > 0

u(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,jx

,

de onde se tem que a primeira coluna da tabela e calculada como u0,j = u1,jx g(jt), 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 fantasmasu1,j e uN+1,j a cada linha na tabela. O valor de u1,j ecalculado a partir da aproximacao

u1,j u1,j2x

ux(0, jt) = g(jt) ,

enquantou0,j+1 = (1 2s)u0,j + s(u1,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

Captulo 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 afinaro instrumentofornecem a tracao necessaria. Em nosso modelo, denotaremos por T o valor da forca de tracaona corda quando em sua posicao de equilbrio, 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 desprezvel 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 equilbrio (hipotese de pequenas oscilacoes).

Um modelo discreto para a corda e pensa-la como N 1 partculas, cada uma com massa igual ax, situadas nas posicoes x = ix, 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 = Nx. A elasticidade da corda e representada por molas que acoplam cada partcula a`s suasvizinhas, ver figura 3.1.

A equacao da onda sera obtida aplicando-se a`s partculas discretasa 2a. lei de Newton daMecanica. Da hipotese de transversalidade, a posicao da i-esima massa, localizada em x = ix,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-1Fi

Figura 3.2: Forcas agindo sobre a i-esima partcula.

que a i-esima partcula esta em sua posicao de equilbrio. Na figura 3.2 mostramos ainda os vetoresforca agindo sobre a partcula i.

Denotamos por Fi1 a forca sobre a partcula i exercida pela partcula i1 e por Fi a forca sobrea partcula i exercida pela partcula i + 1. Definimos i1 como o angulo de rotacao anti-horariaentre o sentido positivo do eixo horizontal e o vetor ligando a partcula 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 apartcula i a` i + 1. Para que o movimento da i-esima partcula seja transversal, e necessario que aresultante na direcao longitudinal das forcas agindo sobre esta partcula seja nula, ou seja,

Fi1 cos i1 = 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 princpio ainda depender do tempo t. Se em um determinado instante todas

as massas estao afastadas da posicao de equilbrio, a corda deve estar mais tracionada que em uminstante em que as massas estao todas proximas a`s suas posicoes de equilbrio. Nossa hipotese depequenas oscilacoes e de que as amplitudes dos movimentos transversais das partculas 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 partcula i e usando que T =Fi cos i, i = 1, 2, . . . , N 1, temos (atencao aos sinais!) que esta componente vale

Fi1 sen i1 + Fi sen i = Tcos i1

sen i1 +T

cos isen i

= T (tan i1 tan i) = T(ui+1 ui

x ui ui1

x

)= T

ui+1 + ui1 2uix

.

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

2uidt2

. Portanto, temos

d2uidt2

=c2

(x)2(ui+1(t) + ui1(t) 2ui(t)) , (3.2)

23

onde a constante c e, em termos dos parametros fsicos,

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

2ux2

no membro direito de (3.2), a equacao de movimento dacorda no limite contnuo 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 fsicos, 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 captulo 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 contnua 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 tpicos de equacao da onda

Voltemos ao modelo discretizado de corda. Para determinar a posicao ui(t) da partcula i, temosque resolver uma EDO de 2a. ordem (3.2). Como esta equacao envolve as posicoes das partculasvizinhas, ou seja, ui1 e ui+1, trata-se portanto de um sistema de EDOs de 2a. ordem. Condicoesiniciais para este sistema sao a especificacao da posicao e da velocidade inicial para cada partculada corda discreta, ou seja, ui(0) = f(ix) e

duidt(0) = g(ix), 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 tpico de corda vibrante discretizada e portanto

d2uidt2

= c2

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

ui(0) = f(ix), i = 1, 2, . . . , N 1duidt(0) = g(ix), 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 EDOs 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 problemautt = c

2 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 EDOs e linear e homogeneo. Da teoria geral das EDOs lineares, o conjunto de todasas solucoes do sistema obedecendo a`s condicoes de Dirichlet homogeneas e um espaco vetorial dedimensao 2(N1). No caso de condicoes de Dirichlet nao-homogeneas, a solucao geral do sistema deEDOs sera a soma de uma solucao particular deste com uma combinacao linear de 2(N1) 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 raciocnio mostra que tambem a EDP da onda e comose fosseum sistema de infinitas EDOs.

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,j1 2ui,j(t)2

= c2ui+1,j + ui1,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 + ui1,j) ui,j1 , (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 fsicadacondicao 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(ix) , (3.10)

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

ut(ix, 0) ui,1 ui,0t

,

a qual usaramos para obter

ui,1 ui,0 + ut(ix, 0)t = ui,0 + g(ix)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, estaramos, 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(ix, 0) ui,1 ui,12t

, (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 fantasmacom j = 1. Para isto, usamos ut(ix, 0) = g(ix) em(3.11) para obter

ui,1 ui,1 = 2t g(ix) . (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 + ui1,0)= (2 2s) f(ix) + 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(ix) + s2[f((i+ 1)x) + f((i 1)x)] + t g(ix) , (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(jt)

euN,j = b(jt) ,

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 tomandox = 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(ix) +

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 + ui1,j) ui,j1 ,

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 cristade 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 = 4t = 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 captulo 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 filmedo movimentoda corda, ou seja, um retrato da posicao de cada ponto da corda, onde os retratos foram tomados aintervalos de 4t = 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 fsico, o princpio do maximo e do mnimo, 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 Fsica 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 numericaxt

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

s = c2(t

x

)2 1 .

Problemas do captulo 3

3.1 O metodo de DAlembert:(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

DAlembert, 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 contnuas, 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) = sennpict

lsen

npix

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) = cosnpict

lsen

npix

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 duplamenteinfinita.

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) 0e

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 equilbrio, especifiquetodas as condicoes iniciais e de contorno a serem aplicadas.

(b) Suponha que a(t) = sen (pit). 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 (npit), 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 (npix) e g(x) 0, entao teremos, se

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

ui,j = Aj sennpii

N,

onde Aj satisfaz a equacao de diferencas

Aj+1 = [2 2s(1 cos npiN

)]Aj Aj1com condicoes iniciais A0 = 1 e A1 = 1 s(1 cos npiN ).

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

O resultado da letra (a) mostra que com as condicoes iniciais dadas, em cada instante de tempodiscreto a forma da corda e a mesma forma inicial com uma amplitudeAj 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 problemasquando s > 1.

Os movimentos da corda descritos neste problema, em que a corda preserva sua forma todo otempo tambem existem no caso contnuo; 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,sno mesmo problema?

32

Captulo 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 R3que nao passa sobre nenhuma das cargas e S e alguma superfcie 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 cargas

que 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 superfcie 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 = 4pi 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 = 4pi (~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 = 4pi , (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 problemastpicos:

(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 smbolo com dois significados completamente distintos.Ate aqui vnhamos usando quase sempre este smbolo 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 ambgua, poremtradicional.

34

normal a` superfcie , ou seja grad ~n, onde ~n e um vetor unitario normal a . Osproblemas de Neumann sao portanto da forma{

(~x) = 0, ~x intgrad(~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= ku .

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 captulo2 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 ilustra