11.1 Método de diferenças finitas
Nesta seção, discutimos os fundamentos do método de diferenças finitas (MDF) para problemas de valores de contorno (PVC). Este método consiste na reformulação do problema contínuo em um problema discreto usando fórmulas de diferenças finitas tomadas sobre uma malha apropriada.
Para introduzir os conceitos principais, consideramos o seguinte problema de valor de contorno (PVC) −uxx=f(x,u),a<x<b,(11.1) u(a)=ua,(11.2)u(b)=ub ,(11.3)
onde ua e ub são dados. Por ter fixados os valores da variável u nos contornos, este é chamado de PVC com condições de Dirichlet1 .
A resolução de um tal problema pelo método de diferenças finitas consiste em quatro etapas fundamentais: 1. construção da malha, 2. construção do problema discreto, 3. resolução do problema discreto e 4. visualização e interpretação dos resultados.
Figura 11.1: Malha uniforme de N pontos em um intervalo [a,b ].
1. Construção da malha. A malha consiste em uma representação discreta do domínio [a,b]. Como veremos, sua construção tem impacto direto nas próximas etapas do método. Aqui, vamos construir a malha mais simples possível, aquela que consiste de N pontos igualmente espaçados, isto é, a chamada malha uniforme.
Para tanto, seja N∈ℕ dado e, então, tomamos o seguinte conjunto discreto PN ={x1,x2,…,xN} (a malha), onde
| xi=a+(i−1)h,i=1,2,…,N, | (11.4) |
com
o qual é chamado de tamanho (ou passo) da malha (veja a Figura 11.1).
2. Construção do problema discreto. A segunda etapa consiste na discretização das equações, no nosso caso, das equações (11.1)-(11.3).
Vamos começar pela Equação (11.1). Em um ponto da malha xi, i=2,3,…,N−1, temos
| −uxx(xi)=f(xi,u(xi )). | (11.6) |
Usando a fórmula de diferenças finitas central de ordem 2 para a segunda derivada, temos
| −u(xi−h)−2u(xi)+u( xi+h)h2+O(h2)=f(xi,u(x i)). | (11.7) |
Rearranjando os termos, obtemos
| −u(xi−h)−2u(xi)+u(xi +h)h2=f(xi,u(xi))+O(h2) . | (11.8) |
Agora, denotando por ui a aproximação numérica de u(xi), a equação acima nos fornece
| 1h2ui−1−2h2ui+ 1h2ui+1=−f(xi,ui), | (11.9) |
para i=2,3,…,N−1. Observamos que trata-se de um sistema de N incógnitas, a saber ui, e de N−2 equações, isto é, um sistema subdeterminado.
Para obtermos um sistema determinado, aplicamos as condições de contorno. Da condição de contorno dada na Equação (11.2), temos
| u(a)=ua⇒u1=ua. | (11.10) |
Analogamente, da condição de contorno dada na Equação (11.2), temos
| u(b)=ub⇒uN=ub. | (11.11) |
Por fim, as equações (11.11), (11.9) e (11.10) determinam o problema discreto associado u1=ua,(11.12)1h2 ui−1−2h2ui+1h2ui+1= −f(xi,ui),i=2,…,N−1,(11.13)uN =ub.(11.14)
Este é um sistema de equações de N incógnitas e N equações.
3. Resolução do sistema discreto. Esta etapa consiste em resolver o sistema discreto construído na etapa anterior.
Para o PVC (11.1)-(11.3), construímos o problema discreto (11.12)-(11.14). Este é um problema de N equações e N incógnitas. Observamos que se f(x,u) é uma função linear, o sistema será linear e podemos resolver o sistema usando de técnicas numéricas para sistema lineares. Agora, se f(x,u) é uma função não linear, podemos usar, por exemplo, do método de Newton para sistemas.
4. Visualização e interpretação dos resultados. A solução do problema discreto consiste dos valores ui, isto é, de aproximações dos valores de u nos pontos da malha. Para visualizarmos a solução podemos, por exemplo, construir o gráfico do conjunto de pontos {(xi,ui)}. Ainda, para obtermos aproximações da solução em outros pontos que não fazem parte da malha, podemos usar de técnicas de interpolação e/ou ajuste.
Exemplo 11.1.1. Use o método de diferenças finitas para resolver o seguinte problema de valor de contorno com condições de Dirichlet homogêneas: −u xx=100(x−1)2,0<x<1,(11.15)u(0) =0,(11.16)u(1)=0.(11.17)
Use a fórmula de diferenças finitas central de ordem 2 para discretizar a derivada em uma malha uniforme de 11 pontos. Calcule, também, a solução analítica deste problema, faça um esboço das soluções numérica e analítica e compute o erro absoluto médio definido por
| E:=1N∑i=1Nu(xi )−ui, | (11.18) |
onde xi é o i-ésimo ponto da malha, i=1,2,…,N e N é o número de pontos na mesma. Por fim, repita seus cálculos para uma malha com 101 pontos. O que ocorre com o erro absoluto médio?
Solução. Vamos seguir as etapas conforme acima.
1. Construção da malha. Tomando N=11, definimos os pontos da malha no domínio [0, 1] por:
| xi=(i−1)h,i=1,2,…,N, | (11.19) |
com h=1∕(N−1).
No Scilab, podemos construir a malha da seguinte forma:
a = 0
b = 1
N = 11
h = (b-a)/(N-1)
x = linspace(a,b,N)
2. Construção do problema discreto. Usando a fórmula de diferenças finitas central de ordem 2 para aproximar a derivada na Equação (11.15), obtemos o seguinte sistema de equações:
| −ui−1−2ui+ui+1h2=100 (xi−1)2,i=2,…,N−1. | (11.20) |
Completamos este sistema com as condições de contorno dadas nas equações (11.16) e (11.17), donde
Ou seja, obtemos o seguinte problema discreto: u1=0,(11.22)−1h2 ui+1−2ui+ui+1=100(xi−1)2 ,i=2,…,N−1,(11.23)uN=0.(11.24)
Observamos que este é um sistema linear N×N, o qual pode ser escrito na forma matricial Au=b, cujos matriz de coeficientes é
| A=10000 ⋯01−2100⋯001−210⋯0⋮⋮⋮⋮⋮⋮⋮00000⋯1, | (11.25) |
o vetor das incógnitas e o vetor dos termos constantes são
| u=u1 u2u3⋮uNeb= 0−100h2(x2−1)2−100h2(x3−1)2⋮0. | (11.26) |
No Scilab, podemos construir o problema discreto a seguinte forma:
A = zeros(N,N)
b = zeros(N)
A(1,1) = 1
b(1) = 0
for i = 2:N-1
A(i,i-1) = 1
A(i,i) = -2
A(i,i+1) = 1
b(i) = -100 * h^2 * (x(i)-1)^2
end
A(N,N) = 1
b(N) = 0
3. Resolução do problema discreto. Neste caso, o problema discreto consiste no sistema linear Au=b e, portanto, a solução é
No Scilab, podemos computar a solução do sistema Au=b com:
u = A\b
Figura 11.2: Esboço dos gráficos das soluções analítica (linha) e numérica (pontos) do PVC dado no Exemplo 11.1.1.
4. Visualização e interpretação dos resultados. Tendo resolvido o problema discreto Au=b, obtemos os valores da solução numérico de u nos pontos da malha, isto é, obtivemos o conjunto de pontos {(xi ,ui)}i=1N. Neste exemplo, queremos comparar a solução numérica com a solução analítica.
A solução analítica pode ser obtida por integração. Temos:
| −uxx=100(x−1)2⇒−ux +c1=100(x−1)33⇒−u+c2x+c1= 100(x−1)412, | (11.28) |
ou seja, u(x)=−(x−1)412+c2x+c1 . As constantes são determinadas pelas condições de contorno dadas pelas equações (11.16) e (11.17), isto é:
| u(0)=0⇒c1=10012,u( 1)=0⇒c2=−10012. | (11.29) |
Portanto, a solução analítica é:
| u(x)=−100(x−1)412−100x12 +10012 | (11.30) |
A Figura 11.2 mostra o esboço dos gráficos das soluções analítica (11.30) e a da solução numérica (11.27).
No Scilab, podemos fazer o esboço das soluções analítica e numérica da seguinte forma:
//def. sol. analitica
deff(’y = ue(x)’,’y = -100.0*(x-1).^4/12 - 100*x/12 + 100.0/12’)
//grafico
xx = linspace(0,1)
yy = ue(xx)
plot(x,u,’ro’,xx,yy,’b-’)
Tabela 11.1: Erro absoluto médio das soluções numéricas com N=11 e N=101 do PVC dado no Exemplo 11.1.1.
| N | h | E |
| 11 | 0,1 | 1,3×10−2 |
| 101 | 0,01 | 1,4×10−4 |
Por fim, computamos o erro absoluto médio das soluções numéricas com N=11 e N=101. A Tabela 11.1 mostra os resultados obtidos. Observamos, que ao diminuirmos 10 vezes o tamanho do passo h, o erro absoluto médio diminui aproximadamente 100 vezes. Este resultado é esperado, pois o problema discreto (11.22)-(11.24) aproxima o problema contínuo (11.15)-(11.17) com erro de truncamento de ordem h2. Verifique!
No Scilab, podemos computar o erro absoluto médio da seguinte forma:
E = sum(abs(ue(x)’ - u))/N
♢
Exercícios resolvidos
ER 11.1.1. Use o método de diferenças finitas para resolver o seguinte problema de valor de contorno: −uxx+u=e−x,0<x<1,(11.31) u(0,5)=1,(11.32)u(1,5)=2.(11.33)
Para tanto, use a fórmula de diferenças finitas central de ordem 2 para discretizar a derivada em uma malha uniforme com passo h=0,1. Faça, então, um esboço do gráfico da solução computada.
Solução. O passo h é uma malha uniforme com N pontos no domínio [0,5,1,5] satisfaz:
| h=(b−a)N−1⇒N=(b−a)h+1. | (11.34) |
Ou seja, a malha deve conter N=11 pontos igualmente espaçados. Denotamos os pontos na malha por xi, onde xi= 0,5+(i−1)h.
| x | u | x | u |
| 0.50 | 1.000000 | 1.00 | 1.643900 |
| 0.60 | 1.143722 | 1.10 | 1.745332 |
| 0.70 | 1.280661 | 1.20 | 1.834176 |
| 0.80 | 1.410269 | 1.30 | 1.908160 |
| 0.90 | 1.531724 | 1.40 | 1.964534 |
| 1.00 | 1.643900 | 1.50 | 2.000000 |
Tabela 11.2: Solução numérica do Exercício 11.1.1.
Agora, a equação diferencial dada no i-ésimo ponto da malha é:
| −uxx(xi)+u(xi)=ex i,i=2,3,…,N−1. | (11.35) |
Denotando ui≈u(xi) e usando a fórmula de diferenças finitas central de ordem dois para a derivada uxx , obtemos:
| −ui−1−2ui+ui+1h2 +ui=exi, | (11.36) |
para i=2,3,…,N−1. Rearranjando os termos e aplicando as condições de contorno, temos o problema discretizado como segue:
| u1=1−ui−1+(2+h2) ui−ui+1=h2exi,i=2,…,N−1, uN=2. | (11.37) |
Figura 11.3: Esboço do gráfico da solução numérica do Exercício 11.1.1.
O problema discreto obtido é um sistema linear N×N. Resolvendo este sistema, obtemos a solução discreta apresentada na Tabela 11.2. A Figura 11.3 mostra um esboço do gráfico da solução computada.
No Scilab, podemos computar a solução numérica e graficá-la com o seguinte código:
//malha
a = 0.5
b = 1.5
N = 11
h = (b-a)/(N-1)
x = linspace(a,b,N)’
//sistema
A = zeros(N,N)
b = zeros(N,1)
A(1,1) = 1
b(1) = 1
for i = 2:N-1
A(i,i-1) = -1
A(i,i) = 2 + h^2
A(i,i+1) = -1
b(i) = h^2 * exp(x(i))
end
A(N,N) = 1
b(N) = 2
//solucao
u = A\b
//grafico
plot(x,u,’b-o’)
♢
Exercícios
E 11.1.1.Considere o seguinte problema de valor de contorno para a equação de calor no estado estacionário:
| −uxx=32,0<x<1. u(0)=5u(1)=10 | (11.38) |
Defina uj=u(xj) onde xj=(j−1) h e j=1,…,5. Aproxime a derivada segunda por um esquema de segunda ordem e transforme a equação diferencial em um sistema de equações lineares. Escreva este sistema linear na forma matricial e resolva-o. Faça o mesmo com o dobro de subintervalos, isto é, com malha de 9 pontos.
Resposta.
| 10000− 12−1000−12−1000−12−100001u1u2u3u4 u5=522210 | (11.39) |
Solução: [5, 9.25, 11.5, 11.75, 10]
| 10000 0000−12−100000 00−12−1000000 0−12−10000000 −12−10000000−1 2−10000000−12 −10000000−12−1 000000001u1u2u3u4 u5u6u7u8u9 =50.50.50.50.5 0.50.50.510 | (11.40) |
Solução: [5,7.375,9.25,10.625,11.5,11.875,11.75, 1.125,10]
E 11.1.2.Considere o seguinte problema de valor de contorno para a equação de calor no estado estacionário:
| −uxx=200e−(x−1) 2,0<x<2.u(0)=120u(2)=100 | (11.41) |
Defina uj=u(xj) onde xj=(j−1) h e j=1,…,21. Aproxime a derivada segunda por um esquema de segunda ordem e transforme a equação diferencial em um sistema de equações lineares. Resolva o sistema linear obtido.
Resposta. 120. 133.56 146.22 157.83 168.22 177.21 184.65 190.38 194.28 196.26 196.26 194.26 190.28 184.38 176.65 167.21 156.22 143.83 130.22 115.56 100.
E 11.1.3.Considere o seguinte problema de valor de contorno para a equação de calor no estado estacionário:
| −uxx=200e−(x−1) 2,0<x<2.u′(0)=0u(2)=10 0 | (11.42) |
Defina uj=u(xj) onde xj=(j−1) h e j=1,…,21. Aproxime a derivada segunda por um esquema de segunda ordem, a derivada primeira na fronteira por um esquema de primeira ordem e transforme a equação diferencial em um sistema de equações lineares. Resolva o sistema linear obtido.
Resposta. 391.13 391.13 390.24 388.29 385.12 380.56 374.44 366.61 356.95 345.38 331.82 316.27 298.73 279.27 257.99 234.99 210.45 184.5 157.34 129.11 100.
E 11.1.4.Considere o seguinte problema de valor de contorno para a equação de calor no estado estacionário com um termo não linear de radiação:
| −uxx=100−u41000 0,0<x<2.u(0)=0u(2)=10 | (11.43) |
Defina uj=u(xj) onde xj=(j−1) h e j=1,…,21. Aproxime a derivada segunda por um esquema de segunda ordem e transforme a equação diferencial em um sistema de equações não lineares. Resolva o sistema obtido. Expresse a solução com dois algarismos depois do separador decimal. Dica: Veja problema 38 da lista 2, seção de sistemas não lineares.
Resposta. 0., 6.57, 12.14, 16.73, 20.4, 23.24, 25.38, 26.93 , 28, 28.7, 29.06, 29.15, 28.95, 28.46, 27.62 , 26.36, 24.59, 22.18, 19.02, 14.98, 10.
E 11.1.5.Considere o seguinte problema de valor de contorno para a equação de calor no estado estacionário com um termo não linear de radiação e um termo de convecção:
| −uxx+3ux=100−u 410000,0<x<2.u′(0)=0u(2 )=10 | (11.44) |
Defina uj=u(xj) onde xj=(j−1) h e j=1,…,21. Aproxime a derivada segunda por um esquema de segunda ordem, a derivada primeira na fronteira por um esquema de primeira ordem, a derivada primeira no interior por um esquema de segunda ordem e transforme a equação diferencial em um sistema de equações não lineares. Resolva o sistema obtido.
Resposta. u(0)=31.62, u(1)=31,50, u(1,9)=18,17.
E 11.1.6.Considere o seguinte problema de valor de contorno:
| −u″+2u′=e−x− u2100,1<x<4.u′(1)+u(1)=2 u′(4)=−1 | (11.45) |
Defina uj=u(xj) onde xj=1+(j−1 )h e j=1,…,101. Aproxime a derivada segunda por um esquema de segunda ordem, a derivada primeira na fronteira por um esquema de primeira ordem, a derivada primeira no interior por um esquema de segunda ordem e transforme a equação diferencial em um sistema de equações não lineares. Resolva o sistema obtido.
Resposta. u(1)=1,900362, u(2,5)=1.943 681, u(4)=1,456517.