Capítulo 10. Equações Diferenciais Parcias -...

214

Grupo: Marcelo Ricardo P. Pinto; Ricardo Rocha; Marcelo A. Carrasqueira; Felipe Aguiar; Danilo C. Pereira; Juliana de Miranda; Mariana Garcia.

Capítulo 10. Equações Diferenciais Parcias 10.1. Equações Diferenciais Parcias Elípticas

10.1.1. Introdução

As equações diferenciais se classificam em ordinárias e parciais.

No curso de análise numérica já foram tratados métodos numéricos para a

resolução das ordinárias.

Este trabalho tratará de métodos, numéricos, para solucionar equações

diferenciais parciais de dois tipos: elípticas e parabólicas.

Uma equação diferencial parcial (E.D.P.) é uma equação que contém uma

função incógnita de duas ou mais variáveis e suas derivadas parciais em relação a

essas variáveis.

A ordem de uma E.D.P. é a ordem da mais alta derivada presente.

Uma solução de uma E.D.P. é qualquer função que satisfaça à equação

identicamente.

Os tipos de E.D.P. que trataremos aqui são de segunda ordem, assim

resolvemos destacar também sua forma geral:

Au

xx y B

ux y

x y Cu

yx y D

ux

x y Euy

x y Fu x y G x y∂∂

∂∂ ∂

∂∂

∂∂

∂∂

2

2

2 2

2( , ) ( , ) ( , ) ( , ) ( , ) ( , ) ( , )+ + + + + =

Definindo ∆=B2-4AC temos:

215

• Se ∆<0 a equação será chamada de elíptica;

• Se ∆=0 a equação será chamada de parabólica;

• Se ∆>0 a equação será chamada de hiperbólica.

Um tipo especial de E.D.P. elíptica é a Equação de Poisson:

∇ =2u x y( , ) ∂∂

∂∂

2

2

2

2

ux

x yu

yx y f x y( , ) ( , ) ( , )+ = ,

e um tipo especial de E.D.P. parabólica é a Equação do Calor:

∂∂

α∂∂

ut

x tu

xx t( , ) ( , )= 2

2

2 .

Essas duas equações serão os nossos objetos de estudo.

10.1.2. Motivação

Em muitas aplicações de Física, Biologia, Química, Engenharia e outros,

mais especificamente na formulações de modelos matemáticos, as equações

diferencias são de grande importância, resolvendo a maioria dos problemas.

Ao contrário com que acontece com as equações diferenciais ordinárias,

não existem métodos de resoluções para as equações diferenciais parciais (com

exceção dos métodos numéricos), cada equação é um caso particular.

Os métodos numéricos são eficientes, de fácil entendimento e sendo de

iteração são facilmente programáveis.

10.1.3. Modelo

Primeiramente trataremos com as equações parciais diferenciais elípticas, e

como citado acima vamos considerar a Equação de Poisson.

216

10.1.4. Metodologia

Vamos considerar a equação de Poisson:

∇ =2u x y( , ) ∂∂

∂∂

2

2

2

2

ux

x yu

yx y f x y( , ) ( , ) ( , )+ = (1)

para (x,y) ∈ R e

u(x,y) = g(x,y) para (x,y) ∈ S

onde R = { (x,y) / a < x < b, c < y < d } e onde S denota o bordo de R.

Assumiremos que f e g são contínuas em seus domínios e que a solução a

encontrar será única.

O método usado é uma adaptação do método da diferença finita para

problemas de valores de borda. O primeiro passo é encontrar inteiros n e m de

modo que h = (b-a)/n e k=(d-c)/m. Particionaremos o intervalo [a,b] em n partes

iguais de largura h e o intervalo [c,d] em m partes iguais de largura k (ver figura 1)

fornecendo uma grade no retângulo R e desenhando linhas verticais e horizontais

através de pontos com coordenadas (xi,yj), onde:

xi = a + ih para i = 0,1,...,n

yj = c + jk para j = 0,1,...,m

217

As linhas x = xi e y = yj são chamadas linhas da grade e suas interseções

são os pontos malhas da grade. Para cada ponto malha do interior da grade

usaremos a série de Taylor na variável x sobre xi para gerar a fórmula da diferença

central :

∂∂

∂∂

ε2

21 1

2

2 4

4

2

12u

xx y

u x y u x y u x y

hh u

xyi j

i j i j i ji j( )

( ) ( ) ( )( , ),,

, , ,=− +

−+ − (2)

onde ε ∈ (xi-1,xi+1);

E a série de Taylor na variável y sobre yj para gerar a fórmula da diferença

central :

∂∂

∂∂

η2

21 1

2

2 4

4

2

12u

yx y

u x y u x y u x y

hk u

yxi j

i j i j i ji j( )

( ) ( ) ( )( , ),,

, , ,=− +

−+ − (3)

onde η ∈ (yj-1,yj+1).

Usando essas fórmulas na Eq. (1), podemos expressar a equação de Poisson

nos pontos (xi,yj) como:

u x y u x y u x y

h

u x y u x y u x y

ki j i j i j i j i j i j( ) ( ) ( ) ( , ) ( , ) ( , ), , ,+ − + −− +

+− +

=1 12

1 12

2 2

= + +f x yh u

xy

k ux

x yi j i j i j( , ) ( , ) ( , )2 4

4

2 4

412 12∂∂

ε∂∂

para cada i = 1,2,...,(n-1) e j = 1,2,...,(m-1).

218

E as condições de bordo como:

u(x0,yj) = g(x0,yj), para cada j = 0,1,...,m;

u(xn,yj) = g(xn,yj), para cada j = 0,1,...,m;

u(xi,y0) = g(xi,y0) , para cada i = 0,1,...,n-1;

u(xi,ym)= g(xi,ym) , para cada i = 0,1,...,n-1;

Dessa forma, a equação diferencial, com o resultado do método da

diferença central, com erro local de truncamento de ordem O(h2 + k2), pode ser

escrita como : 2 12

1 1

2

1 12h

kw w w

hk

w w h f x yi j i j i j i j i j i j

�

��

�

�� +

�

�

�� − + −

�

��

�

�� + = −+ − + −, , , , ,( ) ( ) ( , )

(4)

para cada i = 1,2,...,n-1 e j = 1,2,...,m-1 e

w0,j = g(x0,yj) para cada j = 0,1,...,m

wn,j = g(xn,yj) para cada j = 0,1,...,m (5)

wi,0 = g(xi,yj) para cada i = 1,2,...,n-1

wi,m = g(xi,ym) para cada i = 1,2,...,n-1

onde wi,j aproxima u(xi,yj).

A equação típica em (4) envolve aproximações para u(x,y) nos pontos (xi-

1,yj), (xi,yj), (xi+1,yj), (xi,yj-1) e (xi,yj+1).

Reproduzindo uma porção de grades onde esses pontos estão situados (ver

figura 1), mostraremos que cada equação envolve aproximações na estrela formada

na região sobre (xi,yj).

219

Se usarmos as informações sobre as condições de bordo sempre que o

sistema dado por (4) for apropriado, isto é, todos os pontos (xi,yj) são adjacentes

aos pontos malhas do bordo, temos um sistema (n-1)(m-1) de equações lineares em

(n-1)(m-1) incógnitas, onde as incógnitas são as aproximações wi,j para u(xi,yj) no

interior dos pontos malhas.

O sistema linear envolvendo essas incógnitas é expresso por matriz

calculando mais eficientemente se uma reelaboração dos pontos malhas interiores

for feita. Uma recomendação de reelaboração desses pontos é dado por:

Pl = (xi,yj) e wl = wi,j,

onde l = i + (m-1-j)(n-1), para cada i = 1,2,...,n-1 e j = 1,2,...,m-1. Esta tabela de pontos

malhas consecutivamente deixa certo e levanta o sistema mais baixo.

10.1.5. Algoritmo de Equação de Diferença-Finita de

Poisson

Para aproximar a solução da equação de Poisson

∂∂

∂∂

2

2

2

2

ux

uy

f x y+ = ( , ) a x b c y d≤ ≤ ≤ ≤,

220

Sejam as condições de fronteira:

u(x,y)=g(x,y) se x=a ou x=b e c y d≤ ≤

e

u(x,y)=g(x,y) se y=c ou y=d e a x b≤ ≤

As entradas são os extremos a,b,c,d ; os inteiros m,n ; a tolerância TOL ; e o número

máximos de iterações Lfront.

As saídas aproximadas w i j, para u x yi i( , ) para cada i=1,...,n-1 e j=1,...,m-1 ou uma

mensagem que diz que o número de iterações foi excedido.

Passo 1) Faça h=(b-a)/n

k=(d-c)/m

Passo 2) Para i=1,...,n-1 faça x a ihi = + (os passos 2 e 3 constrõem

os pontos de malha)

Passo 3) Para j=1,...,m-1 faça y c jki = +

Passo 4) Para i=1,...,n-1

para j=1,...,m-1 faça w i j, = 0

Passo 5) Faça

r h kq rl

== =

=

2 2

2 11

/( )

Passo 6) Enquanto I ≤ Lfront faça o passo 7 ao 20.

(Esses passos executam as iterações de Gauss-Seidel)

Passo 7) Faça

221

( ) ( ) ( )( )z h f x y g a y rg x d rw w q

norm z w

w z

m m m m

m

m

= − + + + +

= −

=

− − − −

−

−

21 1 1 1 1 2 2 1

1 1

1 1

, , , /, ,

,

,

Passo 8) Para i = 2,...,n-2 Faça

( ) ( )( )z h f x y rg x d w w rw q

se

w z norm

i m i i m i m i m

i m

= − + + + +

− >

− − − + − −

−

21 1 1 1 1 2

1

, , /, , ,

,

então faça norm w zi m= −−, 1

faça w zi m, − =1

Passo 9) Faça

( ) ( ) ( )( )z h f x y g b y rg x d w rw q

se

w z norm

n m m n n m n m

n m

= − + + + +

− >

− − − − − − − −

− −

21 1 1 1 2 1 1 2

1 1

, , , /, ,

,

então faça

faça

norm w zn m= −− −1 1,

w zn m− − =1 1,

Passo 10) Para j = m-2,...,2 faça o passo 11 ao 13

Passo 11) Faça ( ) ( )( )z h f x y g a y rw rw w qj j j j j= − + + + ++ −2

1 1 1 1 1 2, , /, , ,

se w z normj1, − > então faça norm w zj= −1,

faça w zj1, =

222

Passo 12) Para i=2,...,n-2 faça

z h f x y w rw w rw qi j i j i j i j i j= − + + + +− + + −( ( , ) ) /, , , ,

21 1 1 1

se w zi j, − >norm então faça norm w zi j= −,

faça w zi j, =

Passo 13) Faça

z=( ( , ) ( , ) ) /, , ,− + + + +− − − + − −h f x y g b y w rw rw qn j j n j n j n j2

1 2 1 1 1 1

Passo 14) Faça z h f x y g a y rg x c rw w q= − + + + +( ( , ) ( , ) ( , ) ) /, ,2

1 1 1 1 1 2 2 1

se w z norm1 1, − > então faça norm w z= −1 1,

faça w z1 1, =

Passo 15) Para i=2,...,n-2 faça

z h f x y g x c w rw w qi i i i= − + + + +− +( ( , ) ( , ) ) /, , ,2

1 1 1 1 2 1 1

se w z normi ,1 − > então faça norm w zi= −,1

faça w zi ,1 =

Passo 16) Faça

z h f x y g b y rg x c w rw qn n n n= − + + + +− − − −( ( , ) ( , ) ( , ) ), ,2

1 1 1 1 2 1 1 2

se w z normn− − >1 1, então faça norm w zn= −−1 1,

faça w zn− =1 1,

Passo 17) Se norm � TOL então faça passos 18 e 19

Passo 18) Para i=1,...,n-1

para j=1,...,m-1 imprima ( , , ),x y wi j i j

223

Passo 19) Fim. (programa terminado com sucesso)

Passo 20) Faça i=i+1

Passo 21) Imprima (“O número máximo de iterações foi atingido”)

Fim.

10.1.6. Fluxograma

Início

DadosIniciais

Cálculo de h e k

i = 1

Cálculo de xi

i = n-1 i = i+1No

j = 1

Yes

Cálculo de yi

j = m-1 j = j+1No

i = 1

Yes

j = 1

wi,j = 0

j = m-1 j = j+1No

i = i+1

224

i = n-1 No

Cálculo de r,q e i=1

Yes

Cálculo de z,norm e"w"

norm é menorou igual TOL i = i+1No

Imprima osresultados

Yes

FIM

i menor queLFront

Número de iteraçõesfoi excedido

No

Yes

10.1.7. Exemplo

Considere a equação de Poisson

( ) ( )∂∂

∂∂

2

2

2

2

ux

x yu

yx y xey, , ,+ = 0<x<2, 0<y<1,

com condições de contorno

( )( )

u y

u x x

0 0

0

, ,

, ,

=

=

( )( )

u y e

u x ex

y2 2

1

, ,

, ,

=

=

0 10 2

≤ ≤≤ ≤

yx

,.

usaremos o algoritmo para aproximar a solução exata ( )u x y xey, = com n = 6 e m = 5.

O critério de parada para o método Gauss-Seidel no passo 17 requer

( ) ( )w wi jl

i jl

, , ,− ≤− −1 1010

225

para cada i = 1,...,5, e j = 1,...,4; então a solução para a equação diferencial foi obtida de

modo preciso, e o procedimento de parada para l = 61. O resultado com os valores

corretos são apresentados na tabela abaixo:

j xi yi wi,j(61) u(xi,yi) |u(xi,yj) - wi,j

(61)|

1 0,33330 0,20000 0,40726 0,40713 1,3 x 10-4

2 0,33330 0,40000 0,49748 0,49727 2,08 x 10-4

3 0,33330 0,60000 0,60760 0,60737 2,23 x 10-4

4 0,33330 0,80000 0,74201 0,74185 1,6 x 10-4

1 0,66670 0,20000 0,81452 0,81427 2,55 x 10-4

2 0,66670 0,40000 0,99496 0,99455 4,08 x 10-4

3 0,66670 0,60000 1,21520 1,21470 4,37 x 10-4

4 0,66670 0,80000 1,48400 1,48370 3,15 x 10-4

1 1,00000 0,20000 1,22180 1,22140 3,64 x 10-4

2 1,00000 0,40000 1,49240 1,49180 5,8 x 10-4

3 1,00000 0,60000 1,82270 1,82210 6,24 x 10-4

4 1,00000 0,80000 2,22600 2,22550 4,51 x 10-4

1 1,33330 0,20000 1,62900 1,62850 4,27 x 10-4

2 1,33330 0,40000 1,98980 1,98910 6,79 x 10-4

3 1,33330 0,60000 2,43020 2,42950 7,35 x 10-4

4 1,33330 0,80000 2,96790 2,96740 5,4 x 10-4

1 1,66670 0,20000 2,03600 2,03570 3,71 x 10-4

2 1,66670 0,40000 2,48700 2,48640 5,84 x 10-4

3 1,66670 0,60000 3,03750 3,03690 6,41 x 10-4

4 1,66670 0,80000 3,70970 3,70920 4,89 x 10-4

10.1.8. Implementação

#include<stdio.h>

#include<conio.h>

#include<stdlib.h>

#include<complex.h>

float exemplo(float a1,int a2);

void main(void)

{

226

float m,n,r,q,l,a,b,c,d,i,j,w[4][4],LFront,norm,TOL;

float h,k,x[200],y[200],aux,aux1,aux2,aux3,aux4,aux5,z;

m=5;

n=6;

printf("entre com o valor de a: ");

scanf("%f",&a);

printf("entre com o valor de b: ");

scanf("%f",&b);

printf("entre com o valor de c: ");

scanf("%f",&c);

printf("entre com o valor de d: ");

scanf("%f",&d);

printf("entre com o valor de tolerƒncia: ");

scanf("%f",&TOL);

printf("entre com o m ximo de itera‡äes : ");

scanf("%f",&LFront);

h=(b-a)/n;

k=(d-c)/m;

for (i=1;i<=(n-1);i++) //passo 2

x[i]=a+(i*h);

for (j=1;j<=(m-1);j++) //passo 3

y[i]=c+(j*k);

227

for (i=1;i<=(n-1);i++) //passo 4

for (j=1;j<=(m-1);j++)

w[i][j]=0;

r=h*h/k*k; //passo 5

q=2*(1+r);

l=1;

while (l<=LFront) //passo 6

{

aux=exemplo(x[1],y[m-1]); // passo 7

aux1=exemplo(a,y[m-1]);

aux2=exemplo(x[1],d);

aux3=exemplo(x[1],y[m-2]);

aux4=exemplo(x[2],y[m-1]);

aux5=exemplo(x[1],y[m-1]);

z=(-(h*h)*aux + aux1 + r*aux2 + r*aux3 + aux4)/q;

norm=fabs(z-aux5);

aux5=z;

for (i=2;i<=n-2;i++) //passo 8

{

aux=exemplo(x[i],y[m-1]);

aux1=exemplo(x[i],d);

aux2=exemplo(x[i-1],y[m-1]);

aux3=exemplo(x[i+1],y[m-1]);

228

aux4=exemplo(x[i],y[m-2]);

aux5=exemplo(x[i],y[m-1]);

z=(-(h*h)*aux+ r*aux1 + aux2 + aux3 +r*aux4)/q;

if( abs(aux5-z) > norm)

norm=fabs(aux5-z);

aux5=z;

}

// fim de 8

{

aux=exemplo(x[n-1],y[m-1]);

aux1=exemplo(b,y[m-1]);

aux2=exemplo(x[n-1],d);

aux3=exemplo(x[n-2],y[m-1]);

aux4=exemplo(x[n-1],y[m-2]);

aux5=exemplo(x[n-1],y[m-1]); // passo 9

z=(-(h*h)*aux + aux1 +r*aux2 + aux3 + r*aux4)/q;

if( fabs(aux5-z) > norm)

norm=fabs(aux5-z);

aux5=z;

}

{

for(j=m-2;j<=2;j--) //passo 10

{

229

aux=exemplo(x[1],y[j]); // passo 11

aux1=exemplo(a,y[j]);

aux2=exemplo(x[1],y[j+1]);

aux3=exemplo(x[1],y[j-1]);

aux4=exemplo(x[2],y[j]);

aux5=exemplo(x[1],y[j]);

z=(-(h*h)*aux + aux1 +r*aux2 + r*aux3 + aux4)/q;

if( fabs(aux5-z) > norm)

norm=fabs(aux5-z);

aux5=z;

}

for(i=2;i<=n-2;i++)

{

aux=exemplo(x[i],y[j]); // passo 12

aux1=exemplo(x[i-1],y[j]);

aux2=exemplo(x[i],y[j+1]);

aux3=exemplo(x[i+1],y[j]);

aux4=exemplo(x[i],y[j-1]);

aux5=exemplo(x[i],y[j]);

z=(-(h*h)*aux + aux1 +r*aux2 + aux +r*aux4)/q;

if( fabs(aux5-z) > norm)

norm=fabs(aux5-z);

aux5=z;

230

}

//fim do 12

{

aux=exemplo(x[n-1],y[j]); // passo 13

aux1=exemplo(b,y[j]);

aux2=exemplo(x[n-2],y[j]);

aux3=exemplo(x[n-1],y[j+1]);

aux4=exemplo(x[n-1],y[j-1]);

aux5=exemplo(x[n-1],y[j]);

z=(-(h*h)*aux + aux1 + aux2 +r*aux3 +r*aux4)/q;

if( fabs(aux5-z) > norm)

norm=fabs(aux5-z);

aux5=z;

}

}

//fim do 10

{

aux=exemplo(x[1],y[1]); // passo 14

aux1=exemplo(a,y[1]);

aux2=exemplo(x[1],c);

aux3=exemplo(x[1],y[2]);

aux4=exemplo(x[2],y[1]);

aux5=exemplo(x[1],y[1]);

231

z=(-(h*h)*aux + aux1 + r*aux2 +r*aux3 + aux4)/q;

if( fabs(aux5-z) > norm)

norm=fabs(aux5-z);

aux5=z;

}

for(i=2;i<=n-2;i++)

{

aux=exemplo(x[i],y[1]); // passo 15

aux1=exemplo(x[i],c);

aux2=exemplo(x[i-1],y[1]);

aux3=exemplo(x[i],y[2]);

aux4=exemplo(x[i+1],y[1]);

aux5=exemplo(x[i],y[1]);

z=(-(h*h)*aux + aux1 +r*aux2 + aux3 +r*aux4)/q;

if( fabs(aux5-z) > norm)

norm=fabs(aux5-z);

aux5=z;

}

//fim do 15

{

aux=exemplo(x[n-1],y[1]); // passo 16

aux1=exemplo(b,y[1]);

aux2=exemplo(x[n-1],c);

232

aux3=exemplo(x[n-2],y[1]);

aux4=exemplo(x[n-1],y[2]);

aux5=exemplo(x[n-1],y[1]);

z=(-(h*h)*aux + aux1 +r*aux2 + aux3 +r*aux4)/q;

if( fabs(aux5-z) > norm)

norm=fabs(aux5-z);

aux5=z;

}

if(norm<=TOL) //passo 17

{

for(i=1;i<=(n-1);i++)

for(j=1;j<=(m-1);j++) //passos 18 e 19

{

printf("%f","os valores de x[i] y[j] e

aux5[i][j]:",x[i],y[j],aux5);

getch();

}

printf("Procedimento completado com sucesso!");

}

//fim do if norm

l=l+1; //passo 20

}

printf("Procedimento completado sem sucesso!");

233

}//fim do main

float exemplo(float a1,int a2) { float prod; /* definicao da funcao */

prod=a1*(exp(a2));

return(prod);

}

10.2. Equações Diferenciais Parciais Parabólicas

10.2.1. Modelo

Cuidaremos, agora das equações parciais diferenciais parabólicas., considerando a

Equação de propagação do Calor.

10.2.2. Metodologia

Consideremos a Equação de propagação do Calor:

∂∂

α∂∂

ut

x tu

xx t( , ) ( , )= 2

2

2 (1)

0<x<1 , t>0

com condições iniciais

u(t,0) = 0

u(1,t) = 0, t>0 e u(x,0) = f(x), 0 ≤ x ≤ 1

234

Usaremos a aproximação da solução do problema de diferença finita e o método

similar usado em equações diferenciais elípticas.

Primeiramente, selecionaremos duas malhas constantes, h e k, tal que m = 1/h seja

um inteiro. Os pontos dessa grade serão da forma (xi,tj), onde xi = ih para i= 0,1,...,m e tj

= jk, para j = 0,1,... .

Obteremos o método de diferença usando a série de Taylor em t para formar o

quociente de diferença.

∂∂

∂∂

µut

x tu x t k u x t

kk

tu xi j

i j i ji j( , )

( , ) ( , )( , )=

+ −−

2

2

2 (2)

para algum µj ∈ (tj,tj+1) e a série de Taylor em x para formar o quociente de diferença:

∂∂

∂∂

ε2

2 2

2 4

4

2

12u

xx t

u x h t u x t u x h t

hh u

xti j

i j i j i ji j( )

( , ) ( ) ( , )( , ),,

, , ,=+ − + −

− (3)

onde ε ∈ (xi-1,xi+1).

A equação diferencial parcial (1) implica que para todo ponto interior da grade

(xi,tj), com i = 1,2,...,m-1 e j = 1,2,..., teremos:

∂∂

α∂∂

ut

x tu

xx ti j i j( , ) ( , )− =2

2

2 0

então o método da diferença usando quocientes de diferença (2) e (3) é dado por:

w w

k

w w w

hi j i j i j i j i j, , , , ,+ + −−

−− +

=1 2 1 12

20α (4)

onde wij aproxima u(xi,tj).

O erro local truncado da equação diferencial é:

τ∂∂

µ α∂∂

εi j i j i j

k ut

xh u

xt, ( , ) ( , )= −

2 12

2

22

2 4

4 (5)

Resolvendo a eq. (--0 para wi,j+1) dados, teremos:

235

wk

hw

kh

w wi j i j i j i j, , , ,( )+ + −= −�

��

�

�� + +1

2

22

2 1 112α

α (6)

para cada i= 1,2,...,(m-1) e j = 1,2,... .

Assim a condição inicial u(x,0) = f(x), para 0≤x≤1, implica que wi,0=f(xi), para i

=

0,1,...,m, esse valor pode ser usado na eq. (6) para achar o valor de wi,j para i= 1,2,...,(m-

1). A condição inicial u(0,t) = 0 e u(1,t)= 0 implica que w0,1 = wm,1 = 0, então todas as

entradas da forma wi,1 pode ser determinada.

Se o procedimento é reaplicado uma vez em todas as aproximações wi,j1 os

valores de wi,2, ...,wi,m-1 podem ser obtidos de maneira similar.

Uma explicitação natural do método da diferença na eq. (6) implica que a matriz

(m-1) por (m-1) associada com seus sistemas pode ser escrita na forma tridiagonal :

( )( )

( )

A =

−−

−

�

�

�

�������

1 2 0 01 2

00

0 0 1 2

λ λλ λ

λλ λ

� �

� � � �

� � � � �

� � � � �

� � � � �

� �

onde λ = α2(k/h2).

Se tivermos:

w(0) = (f(x1), ... ,f(xm-1))t

e

w(j) = (wi,j, ... , wm-1,j)t para cada j = 1,2,...

então uma aproximação para a solução é dada por:

w(j) = Aw(j-1), para cada j = 1,2,... .

236

Este método é conhecido como método Avançado de Diferença. Se a solução da

equação diferencial parcial tem quatro derivadas parciais contínuas em x e duas em t,

então a eq. (5) implica que o método é de ordem O(k + h2).

10.2.3. Algoritmo

Para aproximar solução de equações parciais parabólicas

( ) ( )∂∂

α∂∂

ut

x tu

xx t, , ,− =2

2

20 0 1< <x , 0 < <t T ,

sujeito as condições de contorno

( ) ( )u t u t0 1 0, , ,= = 0 < <t T ,

e condições iniciais

( ) ( )u x f x, ,0 = 0 1≤ ≤x :

Entrada ponto final l; tempo máximo T; constante α; inteiros m≥3,N≥1.

Saída aproximações w i j, para ( )u x ti j, para i = 1,...,m-1 e j = 1,...,N.

Passo 1) Faça h l m= ;

k T N= ;

λ α= 2 2k h .

Passo 2) Para i = 1,...,m-1 faça ( )w f ihi = . (valores iniciais)

Passo 3) Faça l1 1 2= + λ;

u l1 1= −λ .

Passo 4) Para i = 2,...,m-2 faça l u i1 11 2= + + −λ λ ;

u li = −λ 1 .

237

Passo 5) Faça l um m− −= + +1 21 2λ λ .

Passo 6) Para j = 1,...,N faça passos 7 - 11.

Passo 7) Faça t = jk;

z w l1 1 1= .

Passo 8) Para i = 2,...,m - 1 faça ( )z w z li i i i= + −λ 1 .

Passo 9) Faça w zm m− −=1 1 .

Passo 10) Para i = m - 2,...,1 faça w z u wi i i i= − +1.

Passo 11) SAÍDA(t);

Para i = 1,...,m - 1 faça x = ih; SAÍDA ( )x w i, .

Passo 12) FIM.

10.2.4. Fluxograma

Início

Dados iniciais

Cálculo de h,k e λ

i = 1

wi = f(ih)

Cálculo de l 1 e u1

i = m - 1

No

Yes

i = 2

Cálculo de l 1 e u1

i = m - 2 No

i = i + 1

i = i + 1

Yes

Cálculo de l m-1

j = 1

Cálculo de t e z 1

i = 2

238

Cálculo de zi

i = m -1

No

i = i + 1

Yes

Cálculo de wm-1

i = m - 2

Cálculo de wi

i = 1 No i = i - 1

Yes

SAÍDA (t)

i = 1

Cálculo de x

SAÍDA (x,wi)

i = m -1 No i = i + 1

Yes

j = N

Yes

FIM

No j = j + 1

10.2.5. Exemplo

Para h = 0.1 e k = 0.01, usado para aproximar a equação do calor

( ) ( )∂∂

α∂∂

ut

x tu

xx t, , ,− =2

2

20 0 1< <x , 0 < t,

( ) ( )u t u t0 1 0, , ,= = 0 < t, e ( )u x x, sen ,0 = π 0 1≤ ≤x ,

Os resultados de w i ,50 até ( )u x i , .0 5 onde i = 0,1,...,10 estão listados na tabela

abaixo:

239

x i w i , 50 u(x i, 0.5) | w i , 50 - u(x i, 0.5) |

0 0 0 -

0,1 0,00289802 0,00222241 6,756 x 10-4

0,2 0,00551236 0,00422728 1,285 x 10-3

0,3 0,00758711 0,00581836 1,769 x 10-3

0,4 0,00891918 0,00683989 2,079 x 10-3

0,5 0,00937818 0,00719188 2,186 x 10-3

0,6 0,00891918 0,00683989 2,079 x 10-3

0,7 0,00758711 0,00581836 1,769 x 10-3

0,8 0,00551236 0,00422728 1,285 x 10-3

0,9 0,00289802 0,00222241 6,756 x 10-4

1 0 0 -

10.2.6. Implementação

#include<stdio.h>

#include<conio.h>

#include<stdlib.h>

#include<complex.h>

float f(float a1);

void main(void)

{

//declaração de variáveis

float l[4],lfinal,T;

float i,j,t,u[4],w[4]; float lam,h,k,x,z[4]; //constantes const A = 5 ; const m=5; const N=6; //entrada de dados printf("entre com o valor de l: ");

240

scanf("%f",&lfinal); printf("entre com o tempo m ximo T: "); scanf("%f",&T); // PASSO 1 h=lfinal/m; k=T/N; lam = A*A*(k/h*h); w[m]=0; //PASSO 2 for (i=1;i<=(m-1);i++) w[i]=f(i*h); // passo 3-11 resolve um sistema linear tridimensional usando // o algoritmo 6.7.) //PASSO 3 l[1] = 1 + lam; u[1] = -lam/(2*l[1]); //PASSO 4 for (i=2;i<=(m-2);i++) { l[i] = 1+lam+((lam*u[i-1])/2); u[i] = -lam/(2*l[i]); } //PASSO 5 l[m-1] = 1+lam+((lam*u[m-2])/2); //PASSO 6 for (j=1;i<=N;j++) //passos 7-11 { //PASSO 7 t = j*k; z[1] = (((1-lam)*w[1])+((lam/2)*w[2]))/l[1]; //PASSO 8 for (i=2;i<=(m-1);i++) z[i] = (((1-lam)*w[i])+((lam/2)*w[i+1]+w[i-1]+z[i-1]))/l[i]; //PASSO 9 w[m-1] = z[m-1]; //PASSO 10 for (i=(m-2);i<=1;i++)

241

w[i] = z[i]-u[i]*w[i]+1; //PASSO 11 // OUTPUT for(i=1;i<=(m-1);i++) { x = i*h; printf("%f","os valores de x w[i]:",x,w[i]); getch(); } } //PASSO 12 //STOP printf("Procedimento completado !"); }//fim do main float f(float a1) { float exponencial; /* definição da função */ exponencial = exp(a1); return(exponencial); }