Cálculo Numérico
Prof: Reinaldo Haas
Equações Diferenciais Ordinárias
1 - Equações Diferenciais Ordinárias
Equações contendo derivadas são equações diferenciais.
Portanto, para compreender e investigar problemas envolvendo o
movimento de fluidos, o fluxo de corrente elétrica em circuitos, a
dissipação de calor em objetos sólidos, a propagação e detecção
de ondas sísmica, o aumento ou diminuição de populações, entre
muitos outros, é necessário saber alguma coisa sobre equações
diferenciais.
Vale lembrar que todo a parte do cálculo chamado de cálculo de
primitivas é nada mais nada menos que a determinação de
soluções de uma equação diferencial.
Como Resolver uma Equação Diferencial Ordinária (EDO)
Na solução de uma EDO dois caminhos podem ser seguidos. Isto é, o que tenta levar à solução exata do problema (método analítico) ou o que encontra uma solução
aproximada (método numérico). Do ponto de vista analítico, resolver uma EDO do tipo y’ = f ( x, y ) é encontrar
uma função y = F ( x ) que satisfaça a equação dada. Por exemplo, dada a equação
diferencial y’ = f ( x, y ) = 2 x + 3, sua solução é obtida por
Na verdade, temos uma família de soluções (para cada C R tem-se uma solução
particular). Na Figura 1 são mostradas algumas destas soluções. No caso para C = 0, C = 2
e C = 4.
y = ( 2x + 3) dx = x 2
+ 3x + C .
Representações de soluções particulares, para alguns valores de C, da
função
y= x 2 + 3 x + C.
Figura 1
C = 0
C = 2
C = 4
x
y
Para determinarmos uma solução específica é necessária a atribuição do valor de y em um dado x. Em outras palavras, deve ser dado um ponto ( x = a , y = s ) por onde a
solução particular deve obrigatoriamente passar.
O processo para encontrar esta solução específica y da equação y’ = f ( x, y ) com y ( a ) = s, onde a e s são dados numéricos, é chamado de problema de
condição inicial. Assim, podemos particularizar a solução do problema anterior atribuindo-lhe, por
exemplo, a seguinte condição:
Logo, a solução geral é dada por y = x 2 + 3 + C, e a particular será dada por
y ( 0 ) = 0 = 0 2
+ 3 x 0 + C C = 0. Ou seja, y = x 2 + 3 x .
0)0(
32
y
xdx
dy
Classificação de Equações
Diferenciais
Equações Diferenciais Ordinárias (EDO) -- se a função
desconhecida depende de uma única variável independente. Neste
caso, aparecem apenas derivadas simples.
Equações Diferenciais Parciais (EDP) -- se a função
desconhecida depende de diversas variáveis independentes. Neste
caso, aparecem as derivadas parciais.
Sistema de equações diferenciais -- se existem duas ou mais
funções que devem ser determinadas, precisamos de um
sistema de equações.
Ordem -- a ordem de uma EDO é a ordem da mais alta derivada que
aparece na equação.
Exemplos:
35 x
dx
dy 12
2
3
3
4
4
ydt
dy
dt
yd
dt
yd
dt
yd
Geralmente a equação F(y, y’, y”, ..., y(n)) = 0 é uma equação diferencial
de ordem n.
4'"2''' tyyyey t
Uma EDO dada para a maior derivada, obtendo-se
),...,",',,( 1 nn yyyytfy
Equações Lineares e não -lineares -- A equação diferencial
0),...,",',( )( nyyytF
É dita linear se F é uma função linear das varáveis y, y’, y”, ...
Assim a equação diferencial ordinária linear geral de ordem n é
)1()()()()( )1(
1
)(
0 tgytaytayta n
nn
A equação diferencial que não é da forma (1) é uma equação
não-linear.
Exemplo:
4'"2''' tyyyey t
Algumas questões relevantes
• Uma equação diferencial sempre tem solução?
(existência)
• Quantas soluções tem uma equação diferencial dada
que ela tem pelo menos uma? Que condições
adicionais devem ser especificadas para se obter
apenas uma única solução? (unicidade)
• Dada uma EDO, podemos determinar, de fato, uma
solução? E, se for o caso, como?
Um computador pode ser uma ferramenta extremamente útil no
estudo de equações diferenciais. Algoritmos já estão sendo usados
há muito tempo para solucioná-las. Entre eles podemos citar: o
método de Euler e Runge-Kutta.
Existem excelentes pacotes numéricos gerais que solucionam uma
gama de problemas matemáticos com versões para PC, estações,
etc. Entre eles temos: o Maple, o Mathematica e o Matlab.
Exemplo: Considere a equação diferencial
dy/dt + 2y = 3. Encontre sua solução.
Solução:
Temos que dy/dt = -2y + 3 ou dy/dt/(y-3/2) = -2
ln |y - 3/2 | = -2t + c
Logo,
y = 3/2 + c*e - 2t
Se g(t) = 0, então a equação é dita equação linear
homogênea e c=-3/2.
Definição 6.1.1: Uma equação que envolve derivadas até
ordem n, é chamada de equação diferencial ordinária (EDO)
de ordem n e pode ser escrita na forma:
(6.1.1)
Definição 6.1.2: A solução da equação (6.1.1) é qualquer função
y = F(x) que é definida em [a,b] e tem n derivadas neste
intervalo e que satisfaz (6.1.1).
y x f x y x y x y xn n( ) ( )( ) ( , ( ), ( ), , ( )) 1
bxa
Equações Diferenciais
A forma mais simples de uma EDO é
(6.1.2)
onde f é contínua para a < x < b.
A solução geral desta equação é:
(6.1.3)
com constante c determinada por
)(xfy
y x f x dx c( ) ( )
y x y( )0 0
De um modo geral temos
(6.1.4)
como por exemplo:
é conveniente reduzi-la a um sistema de EDO de primeira
ordem chamando :
n
n
nn
ayayay
bxa
yyyxfxy
)(,,)(,)(
),,,,()(
1
21
)1(
1)0(,,0)0(,1)0(
10
1)()()( 2
yyy
x
fxxyexyxxy x
y y x1 ( )
1)0(,,0)0(,1)0(
10
1)()()( 2
yyy
x
fxxyexyxxy x
y x y
y x y
y x xy x e y x x f
y y y
x
1 2
2 3
3 2 1
2
1 2 3
1
0 1 0 0 0 1
( )
( )
( ) ( ) ( )
( ) , ( ) , ( )
Para reduzir a uma EDO de primeira ordem
assumimos:
y y x1 ( )
10 x
)(2 xyy
)(3 xyy
isto é, se
3
2
1
~
y
y
y
y ~
y
y
y
y
1
2
3
~( )
( )
( )
( )
y
y
y
y
0
0
0
0
1
2
3
~
( , , , )
~( , )
y
y
y
f x y y y
F x y
2
3
1 2 3
de um modo geral (6.1.6)
dy
dtF x y x
y y
~ ~( , ~( ))
~( ) ~
0 0
17
Equações diferenciais de 1a ordem
Métodos numéricos são usados quando
não é possível obter uma solução geral,
ou a forma dela é tão complicada que
seu uso não é prático.
Uma equação diferencial de 1a ordem
tem a forma ,
e em geral podemos escrevê-la como:
),(
0),,(
yxfy
yyxF
Problema do valor inicial
- uma equação diferencial
- uma condição que deve
ser satisfeita pela solução
00 )(
),(
yxy
yxfy
18
Os métodos que estudaremos
partem da idéia de que o espaço da
variável independente (x) pode ser
discretizado, formando uma rede x0 x1= x0+h x2= x1+h.......
h é o passo .
)(x- 22
e),( yxyxf
O valor da função em cada
ponto da rede é calculado a
partir de expansões em série
de Taylor.
19
)(!2
)()()(
)(
0
2
000
00
xyh
xyhxyhxy
yxy
Método de Euler ou Euler-Cauchy
O valor de y para um passo h é dado
pela expansão de taylor:
Como em geral h é pequeno, suprimimos os
termos de ordem O(h2): h2, h3, .....
Resultando na aproximação
),()(
)()()()()()( 2
yxhfxy
xyhxyhOxyhxyhxy
20
O que resulta no processo iterativo
),(
),(
),(
1
1112
0001
nnnn yxhfyy
yxhfyy
yxhfyy
A omissão dos termos de ordem
superior a 2 causa erros de truncagem
(que podem ocorrer junto a erros de
arredondamento).
21
n xn yn 0,2(xn+yn) Exato erro
0 0,000 0,000 0,000 0,000 0,000
1 0,200 0,000 0,040 0,021 0,021
2 0,400 0,040 0,088 0,092 0,052
3 0,600 0,128 0,146 0,222 0,094
4 0,800 0,274 0,215 0,426 0,152
5 1,000 0,488 0,718 0,230
Exemplo:
passo h=0,2
O erro não é (em geral) conhecido. Podemos estimá-lo
utilizando um passo h´=2h
0)0(
y
yxy
)(2,01 nnnn yxyy
n xn yn 0,4(xn+yn) yn para h=0,2erro = yn(h=0,4) - yn(h=0,2)
0 0,000 0,000 0,000 0,000 0,000
1 0,400 0,000 0,160 0,040 0,040
2 0,800 0,160 0,384 0,274 0,114
3 1,200 0,544
22
Método de Euler melhorado (2a ordem)
Método chamado de preditor-corretor.
)],(),([2
1
),(
*
111
*
1
nnnnnn
nnnn
yxfyxfhyy
yxhfyy
)(2
1
),(
),(
21!
112
1
1
kkyy
kyxhfk
yxhfk
hxx
nn
nn
nn
nn
23
Exemplo:
o mesmo visto anteriormente
02,022,0)(2
2,0
2,02,02,0
2,0
2,0
211
2
1
nnnnn
nnnn
nn
yxykkyy
yxyxk
yxk
h
yxy
n xn yn 0,22(xn+yn) + 0,02 Exato erro
0 0,000 0,0000 0,0200 0,0000 0,0000
1 0,200 0,0200 0,0684 0,0214 0,0014
2 0,400 0,0884 0,1274 0,0918 0,0034
3 0,600 0,2158 0,1995 0,2221 0,0063
4 0,800 0,4153 0,2874 0,4255 0,0102
5 1,000 0,7027 0,3946 0,7183 0,0156
24
Método de Runge-Kutta (4a ordem)
),(
)2
1,
2
1(
)2
1,
2
1(
),(
)22(6
1
34
23
12
1
43211
1
kyhxhfk
kyhxhfk
kyhxhfk
yxhfk
kkkkyy
hxx
nn
nn
nn
nn
nn
nn
Se f(x,y) não depender
de y, o método reduz-se
à regra de integração
de Simpson
25
n xn yn 0,2214(xn+yn) + 0,0214 Exato erro
0 0,000 0,000000 0,021400 0,000000 0,000000
1 0,200 0,021400 0,070418 0,021403 0,000003
2 0,400 0,091818 0,130288 0,091825 0,000007
3 0,600 0,222106 0,203414 0,222119 0,000012
4 0,800 0,425521 0,292730 0,425541 0,000020
5 1,000 0,718251 0,401821 0,718282 0,000031
Comparação entre os métodos
xn y = ex - x - 1 Erro
Euler Euler melhorado Runge-Kutta
0,000 0,000000 0,000 0,0000 0,000000
0,200 0,021403 0,021 0,0014 0,000003
0,400 0,091825 0,052 0,0034 0,000007
0,600 0,222119 0,094 0,0063 0,000012
0,800 0,425541 0,152 0,0102 0,000020
1,000 0,718282 0,229 0,0156 0,000031
26
Qual o valor mais adequado para o passo h?
Se a função f varia muito com y, então h deve ser pequeno, para
evitar erros de truncagem. Em geral, adota-se a “proposta” de
que
12
232
05,001,0
kk
kkK
y
fhK
h h/2 se K 0,05
h 2h se 0,01 K
h não muda se
0,05 K 0,01
Estimativa de erro:
hy
hy
yy
2passooparaobtidoé~~onde
passooparaobtidoé~onde
~~~15
1
27
Métodos para eq. dif. de segunda ordem
P.V.I.
00
00
)(
)(
),,(
yxy
yxy
yyxfy
Novamente o problema é obter os
valores de yn e yn´ para a seqüência
x1 = x0 + h; x2 = x0 + 2h; ...
)(!3
)(!2
)()()(
)(!3
)(!2
)()()(
32
32
xyh
xyh
xyhxyhxy
xyh
xyh
xyhxyhxy
Começamos mais uma vez pelas expansões
em série de Taylor da função e de sua derivada:
28
O método mais simples consiste em desprezar
os termos em derivadas de ordem y´´´ ou
superiores
)()()(
)(2
)()()(2
xyhxyhxy
xyh
xyhxyhxy
001
0
2
001
0000
2
),,(
yhyy
yh
yhyy
yyxfy
1o passo: 2o passo:
112
1
2
112
1111
2
),,(
yhyy
yh
yhyy
yyxfy
Matlab
A biblioteca do Matlab de EDOs os seguintes métodos de valor
inicial:
ode23 método explícito de um passo, RK de ordem
baixa.
ode45 método explícito de um passo, RK de ordem
média. Este é geralmente, o primeiro
método a se tentar em um novo problema.
ode113 método de passo múltiplo, de Adams-
Bashforth-Moulton de ordens variadas.
ODE 23
Adequado para problemas que apresentam bruscas
variações na solução, para os quais é aceitável uma baixa
precisão, ou problemas em que f(t, y) não é suave, ou seja,
descontínua.
Equação da Mola
sujeita a uma força
)(tasin
function wd=spring(t,w);
a=2.0; m=2.0;
c=1.4;h=0.1;
wd = [w(2);-c/m*w(2)-h/m*w(1)+a/m*sin(t)];
)(122
21
tsinm
aw
m
hw
m
cw
ww
>> tspan=[0 100];
>> wo=[0;0];
>> [t,w] = ode23(@mola,tspan,wo);
>> plot(t,w(:,1));
>> grid
>> xlabel('tempo, segundos');
>> ylabel('Elongamento da mola, metros');
>> plot(w(:,1),w(:,2));
>> xlabel('Velocidade das oscilacoes da mola');
>> ylabel('Amplitude das oscilacoes da mola');
ODE 45
Adequado para problemas não-stiff que exijam precisão
moderada.
function yprime=vdpol(t,y)
mu=2;
yprime = [y(2);mu*(1-y(1)^2)*y(2)-y(1)];
Equação de Van der Pol
0)1( 2 xxxx
onde é um parâmetro positivo. Escolhendo
e a equação de Van der Pol se torna:
xy 1
xy 2
12
2
12
21
)1( yyyy
yy
>> tspan = [0 20];
>> yo = [2; 0];
>> [t,y] = ode45(@vdpol,tspan,yo);
>> plot(t,y(:,1),t,y(:,2),'--');
>> xlabel('tempo');
>> title('Solucao de Van der Pol');
37
Equações diferenciais parciais
0:parabólicaCalor)(
0:ahiperbólicOnda)(
0:elípticaLaplace)(0
222
2
2
22
2
2
22
bacTct
T
bacx
ct
bac
Uma equação é dita quasilinear se for
linear nas derivadas mais altas:
;;;onde
),,,,(2
2
2
2
x
uu
yx
uu
x
uu
uuuyxFcubuau
xxyxx
yxyyxyxx
38
Equações de diferenças para Eq. de
Laplace e Poisson
),(
0
2
2
2
2
222
yxfu
y
u
x
uuuuu yyxx
Laplace
Poisson
Vamos ver o caso mais simples em duas dimensões (x e y):
(x-h,y) (x,y) (x+h,y)
h h k
k
(x,y-k)
(x,y+k)
39
),(),(2
1),(
),(),(2
1),(
)(emtermosodesprezand,EE
),(6
1),(
2
1),(),(),(:E
),(6
1),(
2
1),(),(),(:E
3
21
2
2
2
1
kyxukyxuk
yxu
yhxuyhxuh
yxu
hO
yxuyxuhyxhuyxuyhxu
yxuyxuhyxhuyxuyhxu
y
x
xxxxxx
xxxxxx
(x-h,y) (x,y) (x+h,y)
h h k
k
(x,y-k)
(x,y+k)
40
)],(),(
),(),([4
1),(
),(),(2),(1
),(
),(),(2),(1
),(
),(),(2),(),(
2
2
2
kyhxukyhxu
kyhxukyhxuhk
yxu
kyxuyxukyxuk
yxu
yhxuyxuyhxuh
yxu
yxuhyxuyhxuyhxu
xy
yy
xx
xx
Para as derivadas segundas, desprezando os termos O(h4),
temos
Juntando as aproximações das derivadas primeiras e segundas,
fazendo h=k, obtemos a equação de diferenças correspondente à
equação de Poisson:
),(),(4),(),(),(),( 2 yxfhyxuhyxuyhxuhyxuyhxu
41
Para f(x,y) = 0 temos a equação de Laplace.
h é chamado de o comprimento da malha (mesh size).
Equações elípticas - em geral - devem levar em
conta problemas de contorno (condições
previamente definidas numa dada fronteira -
espacial, por exemplo). Casos mais comuns:
•Dirichlet: se u é definido na fronteira C
•Neumann: se un=u/n (derivada na
direção normal) é definida na fronteira.
Para resolver o problema, é necessário
criar uma malha.: nós da rede ou da malha (Pij)
Fronteira C
42
Exemplo
Uma placa de 12 cm de lado tem suas bordas mantidas às
temperaturas mostradas na figura. Quais os valores das
temperaturas no interior da placa? Será escolhido um
comprimento h = 4 cm.
12 x
y 12
u=0
u=100
u=100
u=100 R
u=0
u=100
u=100
P02
P10 P20
P01 P11 P21
P12
43
A equação de transferência de calor é
ut = c2(uxx+uyy)
Para o regime estacionário ut = 0, a
equação se reduz à de Laplace
uxx+uyy = 0
Para cada ponto da malha, temos a seguinte equação:
ui+1,j + ui-1,j + ui,j+1 + ui,j-1 -4 ui,j = 0
P11: - 4u11 + u21 + u01 + u12 + u10 = 0
- 4u11 + u21 + 100 + u12 + 100 = 0
- 4u11 + u21 + u12 = - 200
ui+1,j ui-1,j
ui,j+1
ui,j-1
ui,j
44
- 4u11 + u21 + u12 = -200
u11 - 4u21 + u22 = -200
u11 - 4u12 + u22 = -100
u21 +u12 - 4u22 = -100
Dando como resultados
u11 = u21 = 87,5 (88,1)
u12 = u22 = 62,5 (61,9)