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)