apostila metodos

EQUAÇÃO DIFERENCIAL ORDINÁRIAEQUAÇÃO DIFERENCIAL ORDINÁRIAPROBLEMA DE VALOR INICIALPROBLEMA DE VALOR INICIAL

cayyxfdx

dy== )( :inicial cond.,),(

dx

Geralmente a variável x representa o tempo e a equação diferencialrepresenta a lei da natureza que descreve a taxa de variação de uma grandeza

Exemplo: A taxa de crescimento de uma população é proporcional à população:

h bit tdNKNdN

habitantes de numero :, NKNdt

=

É comum uma situação ser descrito por um sistema de equações diferenciais

siyyyxfdx

dysi

i ,,2,1,),,,,( 21 LL ==

ydEm notação vetorial: cyyf

y== )(,),( ax

dx

d

Equações diferenciais de ordem superior podem ser escritas comoum sistema de equações diferenciais de primeira ordem q ç p

⎞⎛ 23 ddd=′′=′=⎟⎟

⎠

⎞⎜⎜⎝

⎛=

2

3212

2

3

3

)0(,)0(,)0( :iniciais cond.,,,, γγγ

yddy

yyydx

yd

dx

dyyxg

dx

yd

⎪⎧ ==

===

1

2321

)0(

,,

γηηη

ηηη

ddx

yd

dx

dyy

⎪

⎪⎪⎪

⎨ ==

==

2232

112

)0(,

)0(,

γηηη

γηη

dx

ddx

⎪⎪⎪

⎩== 33321

3 )0(,),,,( γηηηηηxg

dx

ddx

cyyfy

== )(,),( axdx

d

dx

Se a variável x representa o tempo, a equação diferencial forneceo vetor velocidade f de uma partícula em função do vetor posição y.

A solução da equação diferencial descreve o movimento de umapartícula neste campo de velocidade

1 )0()(

:Exemplo

Ad

+η

22122

11211

)0(,)(

)0(,)(

Adt

d

Adt

=−−=

=+−=

ηηηη

ηηηη

dt

MÉTODO DE EULER

Método de Passo a Passo Explícito: O valor da função no instante k+1é calculado somente em função do valor da função no instante k

bé h id é d dTambém conhecido como Método da Tangente

Solução numérica: A função y(x) não será obtida para todos os valores de x

Definir pontos xi onde a função y(x) será calculada: Malha

y

2y

Série de Taylor de y(x) em xk:

)()()()(2h

h ′′′1y

0y),(

)(2

)()()(

1

1

kkkk

kkkk

yxfhyy

xyh

xyhxyxy

+=⇒

+′′+′+=

+

+ L

x0x 1x 2x 3x 4x

Exemplo:

dy x)(1)0(; :Resolver == yydx

dySolução exata:

xexy =)(

Método de Euler: kkkkk yhyyyxfhyy +=⇒+= ++ 11 ),(

xk y(xk) yk Erro0.0 1.000 1.000 0.0000 2 1 221 1 200 0 021

h = 0.2

xk y(xk) yk Erro0.0 1.000 1.000 0.0000 1 1 105 1 100 0 005

h = 0.1

0.2 1.221 1.200 -0.0210.4 1.492 1.440 -0.0520.6 1.822 1.728 -0.094

0.1 1.105 1.100 -0.0050.2 1.221 1.210 -0.0110.3 1.350 1.331 -0.0190.4 1.492 1.464 -0.0280 5 1 649 1 610 0 0390.5 1.649 1.610 -0.0390.6 1.822 1.771 -0.051

A solução apresenta errosA solução apresenta erros.

Aparentemente o erro é proporcional ao passo h

Dois tipos de Erros: Erro de Truncamento (Discretização)Erro de Arredondamento

Erro de Truncamento (Discretização) )(),( kkk xyyxhE −=

Erro que surge devido à aproximação da curva y(x)Erro que surge devido à aproximação da curva y(x)por uma reta no intervalo h

Propagação do erro: Solução numérica tende a aproximar a p g ç ç psolução exata que passa por (xk ,yk)

y

Família de curvas

Erro pode ser decrescido diminuindo-se o passo h

x

Aumento do número de intervalos leva a umaumento do número de cálculos necessários e consequentementeum aumento do tempo computacional e dos erros de arredondamento

Erro

T t l

Truncamento

Total

Arredondamento

hn 1=

Existe um tamanho de passo ótimo para produzir o erro mínimo

Análise do Erro de Truncamento para o Método de Euler

)()( xyyxhEERRO GLOBAL )(),( kkk xyyxhE −=ERRO GLOBAL

ERRO LOCAL: Diferença entre o valor yk e o valor da solução exataque passa pelo ponto (xk-1, yk-1)

y

Erro local

Erro global

x

),( 11 −− kk yx

x

Deseja-se obter uma estimativa do erro e saber se o erro cresce ou não

43421)(

1 ),(),(

kk yxy

kkkk yxfhyyyxfdx

dy

′

+ +=⇒=

0 Se >∂∂

y

f),( kk yxy

Í

0 Se <∂∂

y

f

FAMÍLIA DE CURVAS DIGERVE

y


y


))(,( 11 xyx

))(,( 22 xyx

))(,( 11 xyx

))(,( 22 xyx

x

),( 00 yx ),( 11 yx

x

),( 00 yx ),( 11 yx

x

))(,(),( 1111 xyxyyxy ′<′

x

))(,(),( 1111 xyxyyxy ′>′

ERRO AUMENTA EXPONENCIALMENTE ERRO LOCAL COMPENSA ERRO DA DISCRETIZAÇÃO

SE O ERRO LOCAL EM TODOS OS PONTOS FOR MENOR DO QUE ε , O ERRO GLOBAL SATISFAZ A SEGUINTE CONDIÇÃO:

y

fL

e

exyy

xxxLh

Lnh

nnn

max onde ;1

1)(

0 ∂∂

=−−

≤−<<

ε

Lh

y

:pequenofor Se

exyy

Lnh

nn

1)(

−≤− ε

Lhxyy nn )( ≤ ε

Erro pode crescer exponencialmente se L > 0

Erro pode ficar em uma faixa aceitável se h for muito pequeno

Análise do Estabilidade

ESTÁVEL: Produz solução limitada )( ∞<ESTÁVEL: Produz solução limitada

INSTÁVEL: Produz solução que tente ao infinito

)( ∞<

O Método de Euler é condicionalmente estável

ydy α=:Exemplo ydx

α :Exemplo

Solução exata:xeyxy α

0)( =

Solução pelo método de Euler:n

nnnn hyyyhyy )1(011 αα +=⇒+= ++

Observe que por serie de Taylor: L+++≈)(

12h

he h αααObserve que, por serie de Taylor: L+++≈2

1 he α

A solução pelo método de Euler produz os dois primeiros termos da serie de Taylor

:)1(2 Para 0 =−= yαMetodo de Euler:

α = -2.0000h = h = h =

0.2000 0.4000 1.0000Xk Exata Yk Yk Yk

0.00 1.0000 1.0000 1.0000 1.00000.20 0.6703 0.60000.40 0.4493 0.3600 0.20000.60 0.3012 0.21600.80 0.2019 0.1296 0.0400

1.0000

1.5000

1.00 0.1353 0.0778 -1.00001.20 0.0907 0.0467 0.00801.40 0.0608 0.02801.60 0.0408 0.0168 0.00161.80 0.0273 0.01012.00 0.0183 0.0060 0.0003 1.0000

-0.5000

0.0000

0.5000

0.00 1.00 2.00 3.00 4.00 5.00

Exata

h=0.2

h=0.4

h=1.0

2.20 0.0123 0.00362.40 0.0082 0.0022 0.00012.60 0.0055 0.00132.80 0.0037 0.0008 0.00003.00 0.0025 0.0005 -1.00003.20 0.0017 0.0003 0.0000

-1.5000

-1.0000

3.40 0.0011 0.00023.60 0.0007 0.0001 0.00003.80 0.0005 0.00014.00 0.0003 0.0000 0.0000 1.00004.20 0.0002 0.00004.40 0.0002 0.0000 0.0000

0 Para <α Solução exata decai com x

Solução aproximada pode crescer com x se 1h1 >+α

4.60 0.0001 0.00004.80 0.0001 0.0000 0.00005.00 0.0000 0.0000 -1.0000

Solução aproximada pode crescer com x se 1h1 >+α

Para α<0, o Método de Euler só é estável paraα

2−<h

Microsoft Excel

Worksheet

1)0(, problema oResolver =−=′ yyy

Comparar a solução exata com a obtida pelo método de Euler comp ç ph = 0.1, h = 0.5, h = 1.5 e h = 2. Metodo de Euler:

α = -1.0000h=0.1 h=0.5 h=1.5

0.00 1.0000 1.0000 1.0000 1.00000.10 0.9048 0.90000.20 0.8187 0.81000 30 0 7408 0 72900.30 0.7408 0.72900.40 0.6703 0.65610.50 0.6065 0.5905 0.50000.60 0.5488 0.53140.70 0.4966 0.47830.80 0.4493 0.43050.90 0.4066 0.38741.00 0.3679 0.3487 0.25001 10 0 3329 0 3138

0.6000

0.8000

1.0000

1.2000

exata

h=0.1

h=0.5

h=1 5 1.10 0.3329 0.31381.20 0.3012 0.28241.30 0.2725 0.25421.40 0.2466 0.22881.50 0.2231 0.2059 0.1250 -0.50001.60 0.2019 0.18531.70 0.1827 0.16681.80 0.1653 0.1501

-0.2000

0.0000

0.2000

0.4000

0.00 1.00 2.00 3.00 4.00 5.00

h 1.5

1.90 0.1496 0.13512.00 0.1353 0.1216 0.06252.10 0.1225 0.10942.20 0.1108 0.09852.30 0.1003 0.08862.40 0.0907 0.07982.50 0.0821 0.0718 0.03132.60 0.0743 0.0646

-0.6000

-0.4000

2.70 0.0672 0.05812.80 0.0608 0.05232.90 0.0550 0.04713.00 0.0498 0.0424 0.0156 0.25003.10 0.0450 0.03823.20 0.0408 0.03433.30 0.0369 0.03093.40 0.0334 0.02783.50 0.0302 0.0250 0.00783.60 0.0273 0.02253.70 0.0247 0.02033.80 0.0224 0.01823.90 0.0202 0.01644.00 0.0183 0.0148 0.0039

Microsoft Excel

Worksheet

MÉTODO DE EULER DE ORDEM ELEVADA2h

∂+

∂==′=′′

+′′+′+=+

),()( Onde

)(2

)()()( 1

k

kkkk

dyffyxf

dy

dxy

xyh

xyhxyxy L

⎥⎦

⎤⎢⎣

⎡∂∂

+∂∂

++=

∂∂

+ ),(2

),()()(

),()(

2

1 kkkkkk

k

yxfy

f

x

fhyxfhxyxy

dxyxyf

dxy

dxy

⎦⎣ ∂∂2 yx

1)(:ExataSol2)0(:Resolver ++ xexyyxydy x 1)( :ExataSol. 2)0(, :Resolver ++==−= xexyyxydx

Usar método de Euler de segunda ordem com h = 0.1

MÉTODO DE RUNGE-KUTTA

Método de Passo a Passo Explícito: O valor da função no instante k+1é calculado somente em função do valor da função no instante k

O fi i l i li h ál l d f i k 1O coeficiente que multiplica o passo h no cálculo da função no instante k+1é calculado de forma que a expansão coincida com o desenvolvimentoem série de Taylor até os termos de ordem N (ordem do método)

Derivação para Runge-Kutta de segunda ordem:

)(c i)( yxyxFdy

=

:ordem) (2Taylor de Serie

),( c.i.,),(

2

a

00

dFFh

yxyxFdx

⎤⎡ ∂∂

=

Aprox

2)()(

2

1kk

kkk dx

dy

y

F

x

FhFhxyxy +⎥

⎦

⎤⎢⎣

⎡∂

∂+

∂∂

++=+ L ( )

),(

Aprox.

21211 kkkkkk FhyhxFhFhyy μμλλ ++++=+

As constantes são determinadas de tal forma que asexpressões anteriores coincidam

2121 ,,, μμλλ

Por série de Taylor:

L+∂

+∂

+=++F

FhF

hyxFFhyhxF kk

kkkkkk 2121 )()( μμμμ +

∂+

∂+++

yFh

xhyxFFhyhxF kkkkkk 2121 ),(),( μμμμ

A solução aproximada fica:

( ) L+⎥⎦

⎤⎢⎣

⎡∂

∂+

∂∂

+++=⇒ + y

FF

x

FhFhyy k

kk

kkk 212

2211 μμλλλ

ç p

⎦⎣ ∂∂ y

( )Comparando com

1111 e

2

1

2

1;

2

1;1 2121221221 ====→===+ μμλλμλμλλλ

⎧⎞⎛ )(Ffff

⎩⎨⎧

++==

⎟⎠⎞

⎜⎝⎛ ++=+ ),(

),(;

22 12

1211 hfyhxFf

yxFfffhyy

kk

kkkk

RUNGE KUTTA DE QUARTA ORDEM

43211

)(

6336

Ff

ffffhyy kk ⎟

⎠⎞

⎜⎝⎛ ++++=+

12

1

2,

2

),(

fh

yh

xFf

yxFf

kk

kk

⎟⎠⎞

⎜⎝⎛ ++=

=

23 2,

2

22

fh

yh

xFf kk ⎟⎠⎞

⎜⎝⎛ ++=

⎠⎝

( )34 , hfyhxFf kk ++=

Método de Runge-Kutta mais utilizado

Boa combinação entre precisão e simplicidade de programaçãooa co b ação e e p ec são e s p c dade de p og a ação

)5.0 ordem, quarta Kutta-(Runge 1)0(, problema oResolver ==−=′ hyyy

( )432101 226

5.0),(

ffffyy

yyxF

++++=

−=

( )

102

01

2

5.0

6

fyf

yf

⎟⎠⎞

⎜⎝⎛ +−=

−= Metodo de Runge-Kutta (4 ordem):

h = 0.5000f1 f2 f3 f4 y - approx

0.00 1.0000 1.0000

( )304

203

5.02

5.0

fyf

fyf

+−=

⎟⎠⎞

⎜⎝⎛ +−=

⎠⎝0.50 0.6065 -1.0000 -0.7500 -0.8125 -0.5938 0.60681.00 0.3679 -0.6068 -0.4551 -0.4930 -0.3603 0.36821.50 0.2231 -0.3682 -0.2761 -0.2991 -0.2186 0.22342.00 0.1353 -0.2234 -0.1675 -0.1815 -0.1326 0.13552.50 0.0821 -0.1355 -0.1017 -0.1101 -0.0805 0.08223.00 0.0498 -0.0822 -0.0617 -0.0668 -0.0488 0.0499

M

3.50 0.0302 -0.0499 -0.0374 -0.0405 -0.0296 0.03034.00 0.0183 -0.0303 -0.0227 -0.0246 -0.0180 0.01844.50 0.0111 -0.0184 -0.0138 -0.0149 -0.0109 0.01115.00 0.0067 -0.0111 -0.0084 -0.0091 -0.0066 0.00685.50 0.0041 -0.0068 -0.0051 -0.0055 -0.0040 0.00416.00 0.0025 -0.0041 -0.0031 -0.0033 -0.0024 0.00256.50 0.0015 -0.0025 -0.0019 -0.0020 -0.0015 0.00157.00 0.0009 -0.0015 -0.0011 -0.0012 -0.0009 0.00097.50 0.0006 -0.0009 -0.0007 -0.0007 -0.0005 0.00068.00 0.0003 -0.0006 -0.0004 -0.0005 -0.0003 0.00038.50 0.0002 -0.0003 -0.0003 -0.0003 -0.0002 0.00029.00 0.0001 -0.0002 -0.0002 -0.0002 -0.0001 0.00019.50 0.0001 -0.0001 -0.0001 -0.0001 -0.0001 0.0001

10.00 0.0000 -0.0001 -0.0001 -0.0001 0.0000 0.0000

5.1 ordem, quarta Kutta-Runge =h

Metodo de Runge-Kutta (4 ordem):g ( )

h = 1.5000f1 f2 f3 f4 RK EULER

0.00 1.0000 1.0000 1.00001.50 0.2231 -1.0000 -0.2500 -0.8125 0.2188 0.2734 -0.50003 00 0 0498 -0 2734 -0 0684 -0 2222 0 0598 0 0748 0 25003.00 0.0498 -0.2734 -0.0684 -0.2222 0.0598 0.0748 0.25004.50 0.0111 -0.0748 -0.0187 -0.0607 0.0164 0.0204 -0.12506.00 0.0025 -0.0204 -0.0051 -0.0166 0.0045 0.0056 0.06257.50 0.0006 -0.0056 -0.0014 -0.0045 0.0012 0.0015 -0.03139.00 0.0001 -0.0015 -0.0004 -0.0012 0.0003 0.0004 0.0156

10.50 0.0000 -0.0004 -0.0001 -0.0003 0.0001 0.0001 -0.0078

1.0000

1.2000

0.2000

0.4000

0.6000

0.8000

Exata

RK

Euler

-0.6000

-0.4000

-0.2000

0.0000

0.00 2.00 4.00 6.00 8.00 10.00 12.00

RUNGE KUTTA PARA SISTEMA DE EDOS

d

0

0

)(;),,(

)(;),,(

BxzzyxGdz

AxyzyxFdx

dy

==

==

( )

0

22

)(;),,(

ffffh

yy

BxzzyxGdx

++++ ( )22 ggggh

zz ++++=( )

( )1

43211

,,

226

hhh

zyxFf

ffffyy

kkk

kk

⎟⎞

⎜⎛

=

++++=+( )

( )1

43211

,,

226

hf

hhG

zyxGg

ggggzz

kkk

kk

⎟⎞

⎜⎛

=

++++=+

223

112

2,

2,

2

2,

2,

2

gh

zfh

yh

xFf

gh

zfh

yh

xFf

kkk

kkk

⎟⎠⎞

⎜⎝⎛ +++=

⎟⎠⎞

⎜⎝⎛ +++=

223

112

2,

2,

2

2,

2,

2

gh

zfh

yh

xGg

gh

zfh

yh

xGg

kkk

kkk

⎟⎠⎞

⎜⎝⎛ +++=

⎟⎠⎞

⎜⎝⎛ +++=

( )334

223

,,

2,

2,

2

hgzhfyhxFf

gfyf

kkk

kkk

+++=

⎟⎠

⎜⎝

( )334 ,,

222

hgzhfyhxGg kkk +++=⎠⎝

MÉTODOS IMPLÍCITOS

Nos métodos vistos até agora, o valor da função no instante k+1é calculado somente em função de informações no instante k.

C i i é d lí i á iComo visto anteriormente, os métodos explícitos somente são estáveiscomo o tamanho do passo muito pequeno.

MÉTODO DE EULER IMPLÍCITO

y2y −=Δ hx

Série de Taylor para trás ao redor de xk+1:

MÉTODO DE EULER IMPLÍCITO

1y

2y

)(2

)()()( 1

2

111 ++++ +′′+′−=−

=Δ

kkkk xyh

xyhxyhxy

hx

L43421

x0x 1x 2x

),(

),(

111

111

+++

+++

+=⇒−=⇒

kkkk

kkkk

x

yxfhyy

yxfhyyk

Lado direito da equação não é conhecido !!

1)0(, problema oResolver =−=′ yyy

),( yyxF −=)(),(

),(

11111

yy

yhyyyxfhyy

yyxF

k

kkkkkkk

=

−+=⇒+= +++++

Em casos simples a equação pode ser facilmente rearrumada)1(1 h

yk +=+ Em casos simples, a equação pode ser facilmente rearrumada

para determinar a função no passo k+1Se F(x,y) for uma função não linear, a função no passo k+1

será calculada através da solução de uma equação

Metodo de Euler Implicito

α = -1.0000 1.2000

não linear (pelo Método de Newton, por exemplo).

α = 1.0000h=1.5 (EXP) h=0.5 (IMP)

0.00 1.0000 1.0000 1.00000.50 0.60651.00 0.36791.50 0.2231 -0.5000 0.4000 0 2000

0.4000

0.6000

0.8000

1.0000

Exata

Euler Explicito1.50 0.2231 0.5000 0.40002.00 0.13532.50 0.08213.00 0.0498 0.2500 0.16003.50 0.03024.00 0.0183 -0.6000

-0.4000

-0.2000

0.0000

0.2000

0.00 1.00 2.00 3.00 4.00 5.00

Euler Implicito

4.50 0.0111 -0.1250 0.0640

O método de Euler Implícito é incondicionalmente estável

1)0(, problema oResolver ==′ yyy x

x

xkkkkkkk

x

h

yhyyyxfhyy

yyxFk+=⇒+=

=+

+++++111111 ),(

),(

kxkk yyhy k =− +

++111

O valor de yk+1 não pode ser explicitada em função de yk

Para cada passo, deve-se determinar a raiz da equação não linear

conhecidos ,,;0 cbacayy b =−−)( bf

,,;yy

Cálculo pelo Método de Newton, por exemplo

)(

)0(

0

0)(

j

k

b

jyy

cayyyf

==

=−−=

( ) ( ))()()()1(

)(

fi1

repetir ,)(Enquantojjjj

j

jjyfyfyy

yf

+=′−=

<+

ε

)(1

ofim_equantj

k yy =+

MÉTODO IMPLÍCITO DE SEGUNDA ORDEM(MÉTODO DO TRAPÉZIO) y

2y

1y

2yTeorema do Valor Médio:

( ) )()()( que tal, 1

1+

+−

=′∈∃ kkkk xx

xyxyyxx ξξ

x0x 1x 2x

( ) ( )11

)(1

1

+

+

′′′

′+=⇒−

kk

xk

yhyyxx

ξ

( ) ( )),(),(2

1)()(

2

1)(Tomar 111 +++ +=′+′≈′ kkkkkk yxfyxfxyxyy ξ

( ) ( )[ ]++= yxfyxfh

yy ( ) ( )[ ]111 ,,2 +++ ++= kkkkkk yxfyxfyy

Mé d d d d ( d i h2)Método de segunda ordem (erro decai com h2)

Método Estável

Mesmas dificuldades do Método de Primeira Ordem (Euler)

MÉTODOS PREDITOR-CORRETOR

Na resolução de uma equação diferencial deve-se decidir entre o uso deét d lí it i fá il ã tá l dum método explícito, mais fácil e não estável, e o uso de um

método implícito, mais difícil, porém estáveis.

Quando a equação diferencial é não linear deve-se usar técnicas iterativasQuando a equação diferencial é não linear, deve-se usar técnicas iterativaspara resolver a equação não linear resultante em cado passo.

Para o Método de Newton, o chute inicial deve ser bom para o processoo é odo de ew o , o c u e c deve se bo p o p ocessoconvergir em poucas iterações.

Uma opção é usar um método explícito para obter o chute inicial doprocesso iterativo (Preditor) e um método implícito para obtera solução (Corretor).

é d d di l é iMétodo de Heun: Preditor: Euler e Corretor: Trapézio*

1 ),(+ += kkkk yxfhyy

( ) ( )[ ]111*

1 ,,2

:Inicial Chute ++++ ++=→ kkkkkkk yxfyxfh

yyy

MÉTODOS DE MÚLTIPLOS PASSOS

Nos métodos vistos até agora o valor y depende somente deNos métodos vistos até agora, o valor yk+1 depende somente de informações no instante anterior k no instante presente k+1.

Métodos de Múltiplos Passos: Usar informações em vários pontos anterioresMétodos de Múltiplos Passos: Usar informações em vários pontos anteriorespara obter maior precisão.

Passo Único: )()( 11kk

kkkk yxfyy

hOyydy

=−

≈+−

= ++ ),()( kk yxfh

hOhdx

=≈+=

Dois Passos: )(6

)(2

)()()(32

1 xyh

xyh

xyhxyxy kkkkk +′′′+′′+′+=+ L

)()(2)()(

)(6

)(2

)()()(

62

3

32

1

hOxyhxyxy

xyh

xyh

xyhxyxy kkkkk

+′=−⇒

+′′′−′′+′−=− L

),(2

)(2

11311kk

kkkk yxfh

yyhO

h

yy

dx

dy=

−≈+

−= −+−+

)()(2)()( 11 hOxyhxyxy kkk +=⇒ −+

),(2

22

11 kkkk yxfhyy

hhdx

+= −+ MÉTODO LEAPFROG

Idéia geral dos Métodos de Múltiplos Passos: Usar informações em vários pontos anteriores a xk para descrever como a função se comportapontos anteriores a xk para descrever como a função se comportaentre xk e xk+1.

( )∫∫++

=′=11

)()()()(kk xx

dxxyxfdxxyxyxy ( )

∫

∫∫+

=−⇒

==−+

1

)(

)(,)()()( 1

k

kk

xxx

kk

dxxpyy

dxxyxfdxxyxyxy

∫=⇒ + )(1

kx

kk dxxpyy

p(x) é um polinômio interpolador de grau N que aproxima f(x,y) e passapelo conjundo de dados Nkkkifx ii −−= ,,1,,),( L

Para N = 0: método de passo únicop

)(1x

x

kkk hfdxffxpk

k

=→= ∫+

),(1 kkkk

x

yxfhyyk

+=⇒ + Método de Euler

Para N = 1: método de Dois Passos

p(x) é um polinômio linear que interpola ),( e ),( 11 kkkk fxfx −−p( ) p q p )()( 11 kkkk ff

( ))( 11

11

1

ffh

xxff

xx

xxf

xx

xxxp kk

kkk

kk

kk

kk

k −−

−=−−

+−

−= −

−

−−

−

( ) ( )22

)()( 1

2

1 1

1

hffhf

h

xxffhfdxxp kkk

xx

kkkk

x

x

k

k

k

k

−−=−

−+=⇒ −− −

−

∫

Mé d d Ad B hf h d 2 d

( ) ( )1111 322 −+−+ −+=⇒−−+= kkkkkkkkk ffh

yyffh

fhyy

Método de Adams-Bashforth de 2a ordem

Para N = 2: método de Três Passos — Método de Adams-Bashforth 3a ordem

p(x) é um polinômio quadrático que interpola ),( e ),(,),( 1122 kkkkkk fxfxfx −−−−

( )h ( )211 5162324 −−+ +−+= kkkkk fffh

yy

Métodos de Adams-Bashforth são Métodos Explícitos

Pode-se formar os polinômios interpoladores utilizando-se pontos para frente.Métodos Implícitos

Situação mais comum é formar um polinômio de ordem N com os pontos Nkkkk xxxx −−+ L,,, 11

Método de Adams-Moulton de ordem NMétodo de Adams Moulton de ordem N

Para N = 0: método de passo único

Regra do Trapézio( ) ( )[ ]111 ,,2 +++ ++= kkkkkk yxfyxfh

yy

Para N = 3: Método de Adams-Moulton de 4a ordem

( ) ( ) ( ) ( )[ ]5199 +++ ffffh ( ) ( ) ( ) ( )[ ]2211111 ,,5,19,924 −−−−+++ +−++= kkkkkkkkkk yxfyxfyxfyxfyy

SOLUÇÃO DE EQUAÇÕES RÍGIDAS (STIFF)

EQUAÇÕES INSTÁVEISEQUAÇÕES INSTÁVEIS

=−′−′′ 01110 yyy Solução exata: )1()( xexy −=(1)

⎩⎨⎧

−=′=

1)0(1)0(

c.c.yy

⎧

=−′−′′ 01110 yyy Solução exata: )2(1212

111)( 11xx eexy

εε +⎟⎠⎞

⎜⎝⎛ += −(2)

⎩⎨⎧

−=′+=1)0(

1)0(c.c.

yy ε

8

10

2

4

6

Microsoft Excel0

0 0.2 0.4 0.6 0.8 1 1.2

Microsoft Excel

Worksheet

O problema (1) é instável.

Pequenas alterações na condição inicial podem produzir grandesalterações na solução do problema para x grande.

Esses problemas são chamados de mal condicionados.

A solução numérica é extremamente difícil, pois erros de arredondamentoA solução numérica é extremamente difícil, pois erros de arredondamentoe da discretização podem causar o mesmo efeito que a pequenamudança na condição inicial do problema e a solução tenderá adivergir para o infinito.

Estes problemas requerem método numéricos estáveis com passos bemmenores do que o usual.

As instabilidades são mais pronunciadas em problemas não-lineares.

)2( :Exemplo −=′ yxyySolução exata: )1(2)( =xy

2y(0)c.c. =

Este problema é instável. A solução para condição inicial é:0)0( yy =

)2(2

)( 0yxy = )2(

)2()( 2

00xeyy

xy−+

=

3

2

2.5

3

1

1.5

2

Microsoft Excel

Worksheet

0 0.2 0.4 0.6 0.8 1 1.2

MÉTODOS ESTÁVEIS E INSTÁVEIS

⎧ ′ 12 S l ã )1(11

)( 2 +− x

⎩⎨⎧

=+−=′

1)0(12

:Problemay

yy Solução exata: )1(22

)( 2 += xexy

E t bl é tá l i l ã ã d it lt dEste problema é estável, pois a solução não muda muito alterando-se a c.c.

2

1

2

1)(1)0( Se 2 +⎟

⎠⎞

⎜⎝⎛ +=⇒+= − xexyy εε

22)()(

⎠⎝yy

Aplicando-se o método de Leapfrog (segunda ordem, dois passos, explícito):

( ) hyhyyhyy kkkkk 24122 111 ++−=+−+= −−+

+− h 111 2 5

10

∞→∞→⇒

+==

ky

eyy

k

h

quando22

;1 210

-5

0

5

0 2 4 6 8

-10

Método implícitos são estáveis e portanto devem ser usados para problemas stiff

Resolver o problema: [ ] 0)0(;)sin(100 =−=′ yyxy

010)(010)i ( 100x

Solução exata:0001.1

01.0)cos(01.0)sin()(

100xexxxy

−+−=

Solução por Runge-Kutta:

Metodo de Runge-Kutta de quarta ordem:

h = 0.0300 y(0) = 0.0000

X f1 f2 f3 f4 Yk Yexato

0.00 0.0000 00.03 0.0000 1.4999 -0.7500 5.2495 0.033747047 0.0204960.06 -0.3752 1.6865 -1.4060 6.8397 0.068874772 0.0500020 09 0 8911 1 9421 2 3077 9 0234 0 105880707 0 0799120.09 -0.8911 1.9421 -2.3077 9.0234 0.105880707 0.0799120.12 -1.6002 2.2930 -3.5468 12.0236 0.14545912 0.1097730.15 -2.5747 2.7752 -5.2496 16.1467 0.188574804 0.1395360.18 -3.9137 3.4383 -7.5896 21.8144 0.236564525 0.1691740.21 -5.7535 4.3504 -10.8055 29.6059 0.291276497 0.198660.24 -8.2817 5.6055 -15.2252 40.3183 0.355262162 0.2279660 27 -11 7560 7 3323 -21 3001 55 0471 0 43203987 0 2570680.27 -11.7560 7.3323 -21.3001 55.0471 0.43203987 0.2570680.30 -16.5308 9.7080 -29.6503 75.2989 0.526457436 0.2859380.33 -23.0937 12.9765 -41.1288 103.1450 0.645190639 0.3145510.36 -32.1148 17.4727 -56.9085 141.4339 0.797428656 0.3428810.39 -44.5154 23.6576 -78.6019 194.0818 0.995816526 0.3709020.42 -61.5628 32.1644 -108.4264 266.4737 1.257751018 0.398590 45 -84 9991 43 8645 -149 4308 366 0140 1 607162455 0 4259180.45 -84.9991 43.8645 -149.4308 366.0140 1.607162455 0.425918

M eto do de R unge-Kutta de quarta o

h = 0.03 y(0) = 0h 0.03 y(0) 0

X f1 f2 f3 f4 Yk Yexato

0 =$E$3 =(SIN(A7)-0.01*COS(A7)+0.01*

=A7+$B$3 =100*(SIN(A7)-F7) =100*(SIN(A7+$B$3/2)-(F7+$B$3/2*B8)) =100*(SIN(A7+$B$3/2)-(F7+$B$3/2*C8)) =100*(SIN(A7+$B$3)-(F7+$B$3*D8)) =F7+$B$3/6*(B8+2*C8+2*D8+E8) =(SIN(A8)-0.01*COS(A8)+0.01*

=A8+$B$3 =100*(SIN(A8)-F8) =100*(SIN(A8+$B$3/2)-(F8+$B$3/2*B9)) =100*(SIN(A8+$B$3/2)-(F8+$B$3/2*C9)) =100*(SIN(A8+$B$3)-(F8+$B$3*D9)) =F8+$B$3/6*(B9+2*C9+2*D9+E9) =(SIN(A9)-0.01*COS(A9)+0.01*

=A9+$B$3 =100*(SIN(A9)-F9) =100*(SIN(A9+$B$3/2)-(F9+$B$3/2*B10)) =100*(SIN(A9+$B$3/2)-(F9+$B$3/2*C10)) =100*(SIN(A9+$B$3)-(F9+$B$3*D10)) =F9+$B$3/6*(B10+2*C10+2*D10+E10) =(SIN(A10)-0.01*COS(A10)+0.0

=A10+$B$3 =100*(SIN(A10)-F10) =100*(SIN(A10+$B$3/2)-(F10+$B$3/2*B11=100*(SIN(A10+$B$3/2)-(F10+$B$3/2*C11=100*(SIN(A10+$B$3)-(F10+$B$3*D11)) =F10+$B$3/6*(B11+2*C11+2*D11+E11) =(SIN(A11)-0.01*COS(A11)+0.0$ $ ( ( ) ) ( ( $ $ ) ( $ $ ( ( $ $ ) ( $ $ ( ( $ $ ) ( $ $ )) $ $ ( ) ( ( ) ( )

=A11+$B$3 =100*(SIN(A11)-F11) =100*(SIN(A11+$B$3/2)-(F11+$B$3/2*B12=100*(SIN(A11+$B$3/2)-(F11+$B$3/2*C12=100*(SIN(A11+$B$3)-(F11+$B$3*D12)) =F11+$B$3/6*(B12+2*C12+2*D12+E12) =(SIN(A12)-0.01*COS(A12)+0.0






Microsoft Excel

Worksheet

O método só é estável para passos muito pequenos h < 0.030

Solução por Trapézio:ç p p

( ) ( )[ ]

( ) ( )[ ])i (100)i (100

,,2 111

h

yxfyxfh

yy kkkkkk ++= +++

( ) ( )[ ]

( ))sin()sin(10022

1001

)sin(100)sin(1002

1

111

xxhh

y

yxyxh

yy

kkk

kkkkkk

++⎥⎦⎤

⎢⎣⎡ −

−+−+=

+

+++

1002

1

221 h

yk

+

⎥⎦⎢⎣=⇒ +

M eto do do T rapezio

h = 0.5

X Yk

0 =$E$3

=A7+$B$3 =(B7*(1-100*$B$3/2)+$B$3/2*100*(SIN(A7)+SIN(A8)))/(1+$B$3/2*100)

=A8+$B$3 =(B8*(1-100*$B$3/2)+$B$3/2*100*(SIN(A8)+SIN(A9)))/(1+$B$3/2*100)

=A9+$B$3 =(B9*(1-100*$B$3/2)+$B$3/2*100*(SIN(A9)+SIN(A10)))/(1+$B$3/2*100)

=A10+$B$3 =(B10*(1-100*$B$3/2)+$B$3/2*100*(SIN(A10)+SIN(A11)))/(1+$B$3/2*100)

=A11+$B$3 =(B11*(1-100*$B$3/2)+$B$3/2*100*(SIN(A11)+SIN(A12)))/(1+$B$3/2*100)

=A12+$B$3 =(B12*(1-100*$B$3/2)+$B$3/2*100*(SIN(A12)+SIN(A13)))/(1+$B$3/2*100)

=A13+$B$3 =(B13*(1-100*$B$3/2)+$B$3/2*100*(SIN(A13)+SIN(A14)))/(1+$B$3/2*100)

=A14+$B$3 =(B14*(1 100*$B$3/2)+$B$3/2*100*(SIN(A14)+SIN(A15)))/(1+$B$3/2*100)

Metodo do Trapezio

h = 0.5000 y(0) = 0.0000

=A14+$B$3 =(B14*(1-100*$B$3/2)+$B$3/2*100*(SIN(A14)+SIN(A15)))/(1+$B$3/2*100)

=A15+$B$3 =(B15*(1-100*$B$3/2)+$B$3/2*100*(SIN(A15)+SIN(A16)))/(1+$B$3/2*100)

=A16+$B$3 =(B16*(1-100*$B$3/2)+$B$3/2*100*(SIN(A16)+SIN(A17)))/(1+$B$3/2*100)

h 0.5000 y(0) 0.0000

X Yk Yexato

0.00 0.0000 00.50 0.460986095 0.4706026531.00 0.844567185 0.8359843631.50 0.988636033 0.9966879462.00 0.920867137 0.9133675582.50 0.59974723 0.6064229383.00 0.157533472 0.1510048333 50 -0 347014762 -0 3413845223.50 -0.347014762 -0.3413845224.00 -0.744664953 -0.750191044.50 -0.980244479 -0.9753246275.00 -0.95713432 -0.96166473

Microsoft Excel

Worksheet

MEC 1701 – Métodos Numéricos para Engenharia Mecânica MEC 2951 – Métodos Numéricos e Computacionais

Prof. Márcio Carvalho [email protected] Laboratório: Método de Euler Explícito Escreva uma rotina MatLab para solução de um problema de valor inicial usando o método de Euler explícito. Utilize a rotina desenvolvida para resolver o problema:

1)0( ==

=

tu

Kudt

du

Obtenha a solução para os seguintes valores de K: K = 4; K = -4; K = 0.5; K = -0.5 e diferentes passos de tempo. Comente o ocorrido.

4/17/11 9:44 PM C:\Marcio\Cursos-PUC\ENG1714\LABORATORIOS\lab7-e...\eulerexp.m 1 of 2

1 clc % Apaga tudo o escrito na janela de conandos 2 clear all % Apaga as variáveis usadas anteriornente 3 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 4 % PROGRAMA QUE RESOLVE UM PROBLEMA DE VALOR INICIAL PELO METODO DE EULER EXPLICITO 5 6 7 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 8 % Dados de entrada 9 10 disp( ' PROGRAMA QUE RESOLVE UM PROBLEMA DE VALOR INICIAL');11 disp( ' -------- --- ------- -- -------- -- ----- ------- ');12 disp( ' ');13 disp( ' du ');14 disp( ' A EQUAÇÃO DIFERENCIAL A SER RESOLVIDA É : ---- = c1*u ');15 disp( ' dt ');16 disp( ' ');17 18 c1=input('Entre o valor de " c1 " : '); 19 20 disp( ' ');21 disp( ' CONDIÇÃO INICIAL u(0)= c2');22 disp( ' ');23 disp( ' ');24 25 c2=input('Entre o valor de " c2 " : '); 26 27 disp( ' ');28 29 tm=input('Entre o tempo maximo 0 ---> : '); 30 31 clc32 disp( ' ');33 disp( ' ');34 disp( ' SOLUÇÃO ');35 disp( ' ------- ');36 disp( ' ');37 38 h=input('Entre o passo : '); 39 40 disp( ' ');41 42 t(1)=0;43 ue(1)=c2;44 45 k=1;46 47 while t(k)<tm 48 49 t(k+1)=t(k)+h;50 51 % METODO DE EULER

4/17/11 9:44 PM C:\Marcio\Cursos-PUC\ENG1714\LABORATORIOS\lab7-e...\eulerexp.m 2 of 2

52 53 ue(k+1)=ue(k)*(1+c1*h);54 55 k=k+1;56 end57 58 % SOLUÇÃO EXATA59 60 61 s=0;62 for ie = 1:k63 uex(ie) = c2*exp(c1*t(ie));64 s = s + abs( ue(ie) - uex(ie) )65 end66 E = s;67 68 disp('Erro ')69 disp(E)70 71 % GRAFICO 72 figure;73 plot(t,uex,t,ue); 74 legend('Exato','Euler');75 title(' Comparação do metodo de Euler e sol.exata ' );76 xlabel('t');77 ylabel('u');78 79 80 81 82 83 84 85


Prof. Márcio Carvalho [email protected] Laboratório 5 : Método de Runge-Kutta

Escreva uma rotina MatLab para solução de um problema de valor inicial usando o método de Runge-

Kutta de quarta ordem. Utilize a rotina desenvolvida para resolver o problema:

1)0( ==

−=

tu

udt

du t

Determine o valor de para diferentes valores do passo de tempo. Determine o valor do passo

de tal forma que a solução independa do passo escolhido. Faça o mesmo com o método de Euler

explícito, desenvolvido no último laboratório.

)( ∞→tu

4/17/11 9:46 PM C:\Marcio\Cursos-PUC\ENG1714\LABORATORIOS\lab...\rungekutta4.m 1 of 2

1 clc 2 clear all 3 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 4 % PROGRAMA QUE RESOLVE UM PROBLEMA DE VALOR INICIAL 5 % PELO METODO DE RUNGE-KUTTA de 4 ORDEM 6 7 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 8 % Dados de entrada 9 disp( ' PROGRAMA QUE RESOLVE UM PROBLEMA DE VALOR INICIAL');10 disp( ' -------- --- ------- -- -------- -- ----- ------- ');11 disp( ' ');12 13 disp( ' ');14 disp( ' ');15 U0 = input('Entre o valor initial u(t=0) ---> : '); 16 disp( ' ');17 tm = input('Entre o tempo maximo ---> : ');18 disp( ' ');19 h=input('Entre o passo : '); 20 disp( ' ');21 22 disp( ' ');23 t(1)=0;24 ue(1)=U0; % ue = solucao pelo metodo de Euler explicito25 ur(1)=U0; % ur = solucao pelo metodo de Runge Kutta26 k=1;27 28 while t(k)<tm 29 30 t(k+1)=t(k)+h;31 32 % METODO DE EULER33 34 ue(k+1)=ue(k)+h*F(t(k),ue(k));35 36 % METODO DE RUNGE-KUTTA37 38 f1=F(t(k),ur(k));39 f2=F(t(k)+h/2,ur(k)+h/2*f1);40 f3=F(t(k)+h/2,ur(k)+h/2*f2);41 f4=F(t(k)+h,ur(k)+h*f3);42 ur(k+1)= ur(k)+h*(f1/6+f2/3+f3/3+f4/6);43 k=k+1;44 45 end46 47 disp(' Valor em t = tm ');48 disp(' ');49 disp(' passo = ')50 h51 disp(' Euler explicito :');52 ue(k-1)53 disp(' Runge-Kutta :');54 ur(k-1)55

4/17/11 9:46 PM C:\Marcio\Cursos-PUC\ENG1714\LABORATORIOS\lab...\rungekutta4.m 2 of 2

56 % GRAFICO 57 58 plot(t,ue,t,ur); 59 legend('Euler','Runge-Kutta');60 title(' Comparação dos metodos de Euler e Runge-Kutta ' );61 xlabel('t');62 ylabel('u');63 64 65 66 67

4/17/11 9:46 PM C:\Marcio\Cursos-PUC\ENG1714\LABORATORIOS\lab8-edo-rk\F.m 1 of 1

1 function y=F(t,x)2 y=-x^t;


Prof. Márcio Carvalho [email protected] Laboratório 9 : Método Implícito Escreva uma rotina MatLab para solução de um problema de valor inicial usando o método de Euler implícito. Utilize a rotina desenvolvida para resolver o problema:

1)0( ==

−=

tu

udt

du t

Verifique a solução obtida para diferentes passos de tempo. Comente o ocorrido.

4/17/11 9:47 PM C:\Marcio\Cursos-PUC\ENG1714\LABORATORIOS\lab9-...\implicito.m 1 of 2

1 clc % Apaga tudo o escrito na janela de conandos 2 clear all % Apaga as variáveis usadas anteriornente 3 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 4 % PROGRAMA QUE RESOLVE UM PROBLEMA DE VALOR INICIAL POR METODO DE EULER IMPLICITO 5 6 7 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 8 % Dados de entrada 9 disp( ' PROGRAMA QUE RESOLVE UM PROBLEMA DE VALOR INICIAL');10 disp( ' -------- --- ------- -- -------- -- ----- ------- ');11 disp( ' ');12 disp( ' ');13 disp( ' du ');14 disp( ' A EQUAÇÃO DIFERENCIAL A SER RESOLVIDA É : ---- = -u^t ');15 disp( ' dt ');16 disp( ' ');17 disp( ' ');18 19 ci=input('Entre CONDIÇÃO INICIAL u(0) : '); 20 21 disp( ' ');22 disp( ' ');23 24 tm=input('Entre o tempo maximo 0 ---> : '); 25 Dt=input('Entre o passo : '); 26 27 clc28 29 disp( ' ');30 disp( ' SOLUÇÃO ');31 disp( ' ------- ');32 disp( ' ');33 34 35 t(1)=0;36 u(1)=ci;37 k=1;38 39 falha = 0;40 41 while (t(k)<tm) & (falha == 0)42 43 t(k+1)=t(k)+Dt;44 45 newton;46 47 u(k+1) = xraiz;48 49 k = k+1;50 51 end

4/17/11 9:47 PM C:\Marcio\Cursos-PUC\ENG1714\LABORATORIOS\lab9-...\implicito.m 2 of 2

52 53 u(k-1)54 55 % GRAFICO 56 57 plot(t,u); 58 xlabel('t');59 ylabel('u');60 61 62 63 64 65 66 67

4/17/11 9:47 PM C:\Marcio\Cursos-PUC\ENG1714\LABORATORIOS\lab9-edo...\newton.m 1 of 1

1 % SOLUÇÃO PELO MÉTODO DE NEWTON 2 3 maxiter = 10; 4 falha = 0; 5 6 a = t(k+1); 7 b = u(k); 8 9 x0 = u(k);10 11 erro = abs(x0-Dt*f1(a,x0)-b);12 13 i = 1;14 15 while (erro > 1e-4) & (i<maxiter) 16 17 Dx = -(x0-Dt*f1(a,x0)-b)/(1-Dt*df1dx(a,x0));18 x0 = x0 + Dx;19 20 i = i+1;21 22 erro = abs(x0-Dt*f1(a,x0)-b);23 24 end25 26 if (i == maxiter)27 28 disp(' Metodo de Newton nao convergiu');29 falha = 1;30 31 else32 33 xraiz = x0;34 35 end36

4/17/11 9:47 PM C:\Marcio\Cursos-PUC\ENG1714\LABORATORIOS\lab9-edo-eul...\f1.m 1 of 1

1 function y=f1(t,x)2 y=-x^t;3

4/17/11 9:48 PM C:\Marcio\Cursos-PUC\ENG1714\LABORATORIOS\lab9-edo-...\df1dx.m 1 of 1

1 function z=df1dx(t,x)2 z=-t*x^(t-1);

PROBLEMA DE VALOR DE CONTORNO

2

ByAyyyxfdx

yd==′= )1(;)0( :cont. cond.,),,(

2

2

ÉMÉTODO DE TIRO (SHOOTING METHOD)

Baseado nos Problemas de Valor Inicial.

Estima-se a derivada em x = 0:

Utilizando-se os métodos de valor inicial determina-se a solução em x = 1

)0(y′

Utilizando se os métodos de valor inicial, determina se a solução em x 1.

Verifica-se se o valor obtido em x = 1 é igual a B.

BRepetir até convergir.

A

B

Problema Linear

U P i í i d S i ãUsar o Princípio da Superposição.

)( e )( 21 xyxy : Funções que satisfazem a eq. Diferencial e a cond de contorno em x = 0cond. de contorno em x = 0.

)0()0()0()0(

21

21

yyAyy

′≠′==

)()()( 2211 xycxycxy += satisfaz a eq. Diferencial

ycycy += )0()0()0(

Aycc

ycycy

=⇒=++=

)0(1 Se

)0()0()0(

21

2211

)1()1()1( B+

)(

)1()1()1( 2211

xy

Bycycy =+=

é solução do problema se:

)1(1 yBcc⎧ =+12

21

21

2211

21 1)1()1(

)1(

)1()1(

1cce

yy

yBc

Bycyc

cc−=

−−

=⇒⎩⎨⎧

=+

=+

Problema Não-Linear

O Princípio da Superposição não pode ser usado.O Princípio da Superposição não pode ser usado.

Cada valor de corresponde um valor de .γ=′ )0(y )1(y

Utilizar métodos iterativos para determinar o valor de

)()1( γgy =

que fornece a condição de contorno correta em x = 1, isto édeterminar a raiz da equação: Bgf −= )()( γγ

Por exemplo: Método da Secante: )()(Ch t BfPor exemplo: Método da Secante:

( )

)()( Chute

)()( Chute

111

000

γγγγγ

γγγ

Bgf

Bgf

−

−=⇒•

−=⇒•

( )

M

)()()(

01

01112 γγ

γγγγγgg

Bg−−

−−=•

O Problema de Valor Inicial pode ser mal-condicionado mesmo queo Problema de Valor de Contorno seja bem-condicionado.

MÉTODO DE DIFERENÇAS FINITAS

A equação diferencial é transformada em um conjunto de equações algébricas

A função desconhecida é calculada apenas em N pontos – NÓS

Derivadas são aproximadas por diferença dos valores nodais

A equação aproximada é escrita em cada ponto nodal, gerando N equações

Exemplo:Lx

dx

dyyxf

dx

yd<<⎟

⎠⎞

⎜⎝⎛= 0;,,

2

2

LYLy

Yy

==

)(

)0( 0

Obter y(x).

i=1 i=2 i=Ni i+1i 1

Problema Linear

i i+1i-1

1x 1−ixix 1+ix

Nxhx =Δ

Aproximação das derivadas por uso de série de Taylor truncada:

22 hh)(

2)(

2

2

1

2

1 Byh

yhyyAyh

yhyy iiiiiiii KK +′′+′−=+′′+′+= −+

Diferentes aproximações para primeira derivada:Diferentes aproximações para primeira derivada:

DIFERENÇA PARA FRENTEh

yyyA ii

i

−=′⇒ +1)(

DIFERENÇA PARA TRÁSh

yyyB ii

i1)( −−

=′⇒

DIFERENÇA CENTRAL h

yyyBA ii

i 2)()( 11 −+ −

=′⇒−

Aproximação para segunda derivada (Deve-se eliminar o termo de 1a ordem):

DIFERENÇA CENTRALDIFERENÇA CENTRAL

211 2

)()(h

yyyyBA iii

i−+ +−

=′′⇒+h

Interpretação geométrica

dd11

2

2

hh

yyh

yy

h

dxdy

dxdy

dx

dy

dx

d

dx

ydiiii

ed

−+ −−

−

=−

≈⎟⎠⎞

⎜⎝⎛=

211

2

2 2

h

yyy

dx

yd

hhdxdxdx

iii −+ +−≈⇒

⎠⎝

i i+1i-1

x x 1ixe d1−ixix 1+ix

hx =Δ

d

EXEMPLO: Condução de calor em uma barra com convecção natural

AT BT∞Thc , Equação diferencial que descreve

a variação da temperatura na barra.

Lx

( )c TTkA

Ph

dx

Td=−− ∞ 0

2

2

B

A

TLxT

TxT

====

)(

)0(

i=1 i=2 i=Ni i+1i-1

1x 1−ixix 1+ix

Nx

Incógnitas do problema: Ni TTTT ,,,,, 21 KK

No ponto i: ( )∞ =−− TTkA

Ph

dx

Tdi

c

i

2

2

0

( )∞−+

⎟⎞

⎜⎛

⎟⎞

⎜⎛

⎟⎞

⎜⎛

=−−+−

⇒

PhPh

TTkA

Ph

h

TTTi

ciii2

11

121

02

∞+− −=⎟⎠⎞

⎜⎝⎛+⎟

⎠⎞

⎜⎝⎛ −−+⎟

⎠⎞

⎜⎝⎛⇒ T

kA

PhT

hT

kA

Ph

hT

hc

iic

i 12212

121

Para se obter as N equações necessárias para determinar as N incógnitas,a aproximação deve ser satisfeita em todos os nós:

(c.c.) 1 1 ATTi =⇒=

121

1212 322212

cc

PhPh

TkA

PhT

hT

kA

Ph

hT

hi

⎞⎛⎞⎛⎞⎛

−=⎟⎠⎞

⎜⎝⎛+⎟

⎠⎞

⎜⎝⎛ −−+⎟

⎠⎞

⎜⎝⎛⇒= ∞

1213 423222

cc TkA

PhT

hT

kA

Ph

hT

hi −=⎟

⎠⎞

⎜⎝⎛+⎟

⎠⎞

⎜⎝⎛ −−+⎟

⎠⎞

⎜⎝⎛⇒= ∞

M

(c.c.) BN TTNi =⇒=

O sistema pode ser escrito em forma matricial como:

⎥⎥⎥⎤

⎢⎢⎢⎡

−

⎥⎥⎥⎤

⎢⎢⎢⎡

⎥⎥⎥⎤

⎢⎢⎢⎡

⎟⎠⎞⎜

⎝⎛ +− ∞

c

A

c

Ph

TkAPhT

T

T

T

hkAPh

hhK

K

2

1

222 0012100001

⎥⎥⎥⎥⎥

⎢⎢⎢⎢⎢

−

−=

⎥⎥⎥⎥⎥

⎢⎢⎢⎢⎢

⎥⎥⎥⎥⎥

⎢⎢⎢⎢⎢

⎟⎠⎞⎜

⎝⎛ +−

⎠⎝

∞

∞

c

c

c

TkAPh

TkAPh

T

T

hkAPh

hh

hhh

MMM

K4

3

222 0121

⎥⎥⎥

⎦⎢⎢⎢

⎣⎥⎥⎦⎢

⎢⎣⎥⎦

⎢⎣

BN TT

MK 10000

Exemplo: Desenvolver o problema se os nós não são uniformemente distribuidos.

i=1 i=2 i=Ni i+1i-1

1−ixΔixΔ

1x 1−ixix 1+ix

Nx

EXEMPLO: Condição de contorno na derivada da função

=

TxTdx

Td

)0(

02

2

AT∞Thc ,

( ) ( )( )∞−===

==

TLxTK

hLx

dx

dT

TxT

e

A)0(L

x

Kdx

No ponto i:

2

0

11

2

2

+−

=

+ iii

i

TTT

dx

Td

0121

02

12212

211

=⎟⎠⎞

⎜⎝⎛+⎟

⎠⎞

⎜⎝⎛−+⎟

⎠⎞

⎜⎝⎛⇒

=+

⇒

+−

−+

iii

iii

Th

Th

Th

h

TTT

⎠⎝⎠⎝⎠⎝ hhh

Para se obter as N equações necessárias para determinar as N incógnitas,i ã d ti f it t d óa aproximação deve ser satisfeita em todos os nós:

=⇒= TTi A (c.c.) 1 1

⎞⎛⎞⎛⎞⎛

=⎟⎠⎞

⎜⎝⎛+⎟

⎠⎞

⎜⎝⎛−+⎟

⎠⎞

⎜⎝⎛⇒= T

hT

hT

hi

121

0121

2 322212

=⎟⎠⎞

⎜⎝⎛+⎟

⎠⎞

⎜⎝⎛−+⎟

⎠⎞

⎜⎝⎛⇒= T

hT

hT

hi 0

1213 423222

M

( ) ∞−∞− −=⎟

⎠⎞

⎜⎝⎛ −+⎟

⎠⎞

⎜⎝⎛−⇒−=

−⇒= T

K

h

K

h

hT

hTT

K

h

h

TTNi ee

NNeNN 11

11

O sistema pode ser escrito em forma matricial como:

( )( ) ⎥

⎥⎥⎤

⎢⎢⎢⎡

⎥⎥⎥⎤

⎢⎢⎢⎡

⎥⎥⎥⎤

⎢⎢⎢⎡

−

T

T

T

T

hhh

A

K

K

0

0000121000001

2

1

222

( )⎥⎥⎥⎥⎥

⎢⎢⎢⎢⎢

=

⎥⎥⎥⎥⎥

⎢⎢⎢⎢⎢

⎥⎥⎥⎥⎥

⎢⎢⎢⎢⎢

−T

T

h

hhh

hhh

MMM

K0

0

11

001214

3222

⎥⎥⎦⎢

⎢⎣−⎥

⎥⎦⎢

⎢⎣⎥⎥⎦⎢

⎢⎣ −−

∞TKh

TKh

hh eN

eK 110000


EXEMPLO: Problema de Convecção e Difusão

AV BVEquação diferencial que descreve

a velocidade em cada ponto.

Lx

p

dx

udK

dx

duu =− 0

2

2

B

A

VLxu

Vxudxdx

====

)(

)0(

B)(

i=1 i=2 i=Ni i+1i-1

1x 1−ixix 1+ix

Nx

Incógnitas do problema: Ni uuuu ,,,,, 21 KK

Aproximação por diferença central: 11

2h

uu

dx

du ii

i

−+ −≈

211

2

2 2

h

uuu

dx

ud iii

i

−+ +−≈

No ponto i:

02

2

=−i dx

udK

dx

duu

02

2 21111 =⎥⎦⎤

⎢⎣⎡ +−

−⎥⎦⎤

⎢⎣⎡ −

⇒ −+−+ iiiiii

ii

h

uuuK

h

uuu

dxdx

02

212

11212

=⎟⎠⎞

⎜⎝⎛−+⎥

⎦

⎤⎢⎣

⎡⎟⎠⎞

⎜⎝⎛ −

++⎟⎠⎞

⎜⎝⎛−⇒ +

−+− i

Ai

i

Bi

iii

Ai

uh

Ku

h

uu

h

Ku

h

K

321444 3444 21321AiBiAi

Para se obter as N equações necessárias para determinar as N incógnitas,a aproximação deve ser satisfeita em todos os nós:

02

2

(c.c.) 1

32213

212

1 A

uK

uuuK

uK

i

Vui

=⎟⎠⎞

⎜⎝⎛−+⎥

⎤⎢⎡

⎟⎠⎞

⎜⎝⎛ −

++⎟⎠⎞

⎜⎝⎛−⇒=

=⇒=

02

23

02

42324

222

322212

uh

Ku

h

uu

h

Ku

h

Ki

uh

uhh

uh

i

=⎟⎠⎞

⎜⎝⎛−+⎥

⎦

⎤⎢⎣

⎡⎟⎠⎞

⎜⎝⎛ −

++⎟⎠⎞

⎜⎝⎛−⇒=

⎟⎠

⎜⎝⎥

⎦⎢⎣

⎟⎠

⎜⎝

⎟⎠

⎜⎝

⇒

(c.c.)

2

BN VuNi

hhhh

=⇒=

⎠⎝⎦⎣ ⎠⎝⎠⎝M

⎥⎥⎥⎤

⎢⎢⎢⎡

⎥⎥⎥⎤

⎢⎢⎢⎡

⎥⎥⎤

⎢⎢⎡ AV

u

u

ABA K

K

0

000

000012

1

222

⎥⎥⎥⎥

⎢⎢⎢⎢

=

⎥⎥⎥⎥

⎢⎢⎢⎢

⎥⎥⎥⎥

⎢⎢⎢⎢

u

uABA

MMM

K0

00

4

3333

222

Função de

⎥⎥⎥

⎦⎢⎢⎢

⎣⎥⎥⎥

⎦⎢⎢⎢

⎣⎥⎥⎦⎢

⎢⎣

BN Tu

MMK 10000

u1 e u3

Os coeficientes da matriz dependem da solução do problemaOs coeficientes da matriz dependem da solução do problema.

Sistema de equações Não-Linear

Solução de Sistema de Equações Não-Linear

Método de Picard:Método de Picard:

1. Chute inicial; [ ])0()0(3

)0(2

)0(1

)0( ,,,, Nuuuuc K=

2. Calcular coeficientes da matriz usando o valor atual das incógnitas;

( ))(kcAA =

3. Resolver o sistema de equações e determinar o novo valor dasincógnitas;

( )cAA

g

( )( ) fcAc kk 1)()1( −+ =

• Comparar solução atual com anterior;

• Se não covergiu, voltar para 2.

Convergência Ruim

Método de Newton:

G li ã d Mé d d N 1 ã ã liGeneralização do Método de Newton para 1 equação não-linear

( ) ( ) ( ) ( )( ) ( ) ( )

xfxxfxxfxxf +′′+′+=+ ΔΔΔ 2 K

( ) ( ) ( )

( )xf

xfxxfxxf ′+≈=+ ΔΔ 0

( )( )xf

xfx

′−=Δ

:inicialChute xPROCEDIMENTO ITERATIVO

)(

)0(

)(do ,)( While

0 :inicialChute

>=

i

xfxf

ix

ε

)()1(

)(

)(

)(

)(

+ Δ+=′

−=Δ

ii

i

i

xxxxf

xfx

)1( :Raiz1

+

+=ixii

⎪⎪⎪

⎨

⎧

===

0)(

0),,,,(

0),,,,(

3212

3211

N

N

xxxxf

xxxxf

xxxxf

K

KSistema a ser resolvido:

⎪⎪⎪

⎩

⎨

=

=

0),,,,(

0),,,,(

321

3213

NN

N

xxxxf

xxxxf

K

M

K

Expansão por série de Taylor até termos de primeira ordem de cada equação:

NN xxxxxxf ΔΔΔ ≈=+++ K22111 0),,,(

NN

N

f

xx

fx

x

fx

x

fxxxf

ΔΔΔ

ΔΔΔ∂∂

++∂∂

+∂∂

+ KK 12

2

11

1

1211

0)(

),,,(

NN

NN

xx

fx

x

fx

x

fxxxf

xxxxxxf

ΔΔΔ

ΔΔΔ

∂∂

++∂∂

+∂∂

+

≈=+++

KK

K

22

21

2212

22112

),,,(

0),,,(

NNN

N

xxxxxxf

xxx

ΔΔΔ ≈=+++

∂∂∂

K

M

2211

21

0),,,(

NN

NNNNN

NNN

xx

fx

x

fx

x

fxxxf

f

ΔΔΔ∂∂

++∂∂

+∂∂

+ KK 22

11

21

2211

),,,(

),,,(

NN

N

fff

xx

fx

x

fx

x

fxxxf ΔΔΔ

∂∂∂∂∂

++∂∂

+∂∂

=− KK 12

2

11

1

1211 ),,,(

NN

N xx

fx

x

fx

x

fxxxf ΔΔΔ

∂∂

++∂∂

+∂∂

=−

M

KK 22

2

21

1

2212 ),,,(

NN

NNNNN x

x

fx

x

fx

x

fxxxf ΔΔΔ

∂∂

++∂∂

+∂∂

=− KK 22

11

21 ),,,(

f

fxfffx

f

x

f

x

f

N 11

222

1

2

1

1

1

⎥⎤

⎢⎡

⎥⎤

⎢⎡⎥⎥⎥⎤

⎢⎢⎢⎡

∂∂∂∂∂

∂∂

∂∂

ΔΔ

KSistema emForma matricial

f

f

x

x

fff

x

f

x

f

x

f

NN

N

222

2

2

1

2

⎥⎥⎥⎥

⎦⎢⎢⎢⎢

⎣

−=

⎥⎥⎥⎥

⎦⎢⎢⎢⎢

⎣⎥⎥⎥⎥⎥

⎢⎢⎢⎢⎢

∂∂∂

∂∂

∂∂

∂∂

Δ

ΔMM

M

K

{fx

x

f

x

f

x

ff

N

x

N

J

N

NNN

21

⎦⎣⎦⎣⎥⎥⎦⎢

⎢⎣ ∂

∂∂∂

∂∂ Δ

Δ321

4444 34444 21

K

if∂

fJx 1−=Δ

Matrix Jacobianaj

iij x

fJ

∂∂

=

PROCEDIMENTO ITERATIVO

)0(

0

:inicial Chute

i

c

( )1

)( do , While

0

>

=icf

i

ε

)()1(

1

+

−

+=

−=ii xcc

fJx

Δ

Δ Solução de um sistema linear

)1( :Raiz

1+

+=ic

ii

Convergência Quadrática

Voltando ao Problema Não-Linear:

dx

udK

dx

duu =− 0

2

2

B

A

VLxu

Vxu

====

)(

)0(

Sistema de equações algébricas não-linear resultante:

( )

( ) 02

0,,,,,,

1111

11111

=⎥⎤

⎢⎡ +−

−⎥⎤

⎢⎡ −

=

=−=

−+−+

+−

iiiii

ANiii

uuuK

uuuuuuuuf

Vuuuuuuf KK

( )

( ) 0

,,3,2

02

,,,,,,2111

=

=⎥⎦⎢⎣−⎥⎦⎢⎣

=+− iNiiii

Vf

Nih

Kh

uuuuuuf

K

KK

( ) 0,,,,,, 111 =−=+− BNNiiiN Vuuuuuuf KK

Solução pelo Método de Newton…

Cálculo da matriz Jacobiana:

2

0;;0;1 1

2

1

1

1

∂∂∂∂∂

=∂∂

=∂∂

=∂∂

N

ffKfKfKf

u

f

u

f

u

fK

;0;;2

;;0

0;;0;2

;2

2;

2

333243333

2

4

222

3

2132

2

222

1

2

=∂

+−=∂−

+=∂

−−=∂

=∂

=∂∂

=∂∂

+−=∂∂−

+=∂∂

−−=∂∂

N

fuKfuuKfuKff

u

f

u

f

h

u

h

K

u

f

h

uu

h

K

u

f

h

u

h

K

u

fK

;0;2

;2

;2

;05

24

23

221

∂∂∂

=∂

+=∂

+=∂

=∂

=∂

NNN fff

uhhuhhuhhuu

M

K

1;0;;011

=∂∂

=∂∂

=∂∂

− N

N

N

NN

u

f

u

f

u

fK


EXERCÍCIO: Escoamento desenvolvido de um fluido não Newtoniano emum tubo circular.

dd

Fx

02)()(

022

⇒=∑

ddpdp

rdxrdxxprxp

τ

πτππ

22

02)()( 22

⎟⎞

⎜⎛

⎟⎞

⎜⎛

⇒=−+−

x

r )(xp )( dxxp +

τ

drdx

pr

dx

p τ 22 =⎟⎠⎞

⎜⎝⎛−⇒=⎟

⎠⎞

⎜⎝⎛−

dx

d ⎞⎛⎞⎛ duddpLei de Newton de viscosidade

dr

duητ = ⎟⎠⎞

⎜⎝⎛=⎟

⎠⎞

⎜⎝⎛−

dr

du

dr

d

dx

dp η2

Condições de contorno:

00

0

=⇒=

=⇒=du

r

uRr

00 ⇒dr

r

Para fluido Newtoniano, a viscosidade é constante: μη =

2dddd ⎞⎛⎞⎛2

2

22dr

ud

dr

du

dr

d

dx

dp μμ =⎟⎠⎞

⎜⎝⎛=⎟

⎠⎞

⎜⎝⎛−

⎞⎛⎞⎛2Problema linear. Solução analítica:

⎟⎞

⎜⎛

+⎟⎠⎞

⎜⎝⎛−=⇒⎟

⎠⎞

⎜⎝⎛−=

2

12

2

22

1

dpr

Cdx

dpr

dr

du

dx

dp

dr

ud

μμ

⎥⎤

⎢⎡

⎟⎞

⎜⎛

⎟⎞

⎜⎛=

++⎟⎠⎞

⎜⎝⎛−=

22

21

1)(

4)(

rdpRru

CrCdx

dprru

μ

Perfil parabólico⎥⎥⎦⎢

⎢⎣

⎟⎠

⎜⎝

−⎟⎠

⎜⎝−= 1

4)(

Rdxru

μPerfil parabólico

⎟⎠⎞

⎜⎝⎛−=

dx

dpRuMAX μ4

2

Para fluido não Newtoniano, a viscosidade varia com a taxa de deformação.

( )2

12

0 1

−

∞∞⎥⎥⎦

⎤

⎢⎢⎣

⎡⎟⎠⎞

⎜⎝⎛+−+=

n

dr

duληηηη⎥⎦⎢⎣

( ) ⎟⎞

⎜⎛

⎟⎞

⎜⎛ dpdud ( ) ⇒⎟

⎠⎞

⎜⎝⎛−=⎟

⎠⎞

⎜⎝⎛

dx

dp

dr

du

dr

d γη &2

⎪⎫

⎪⎧

( ) ⎟⎠⎞

⎜⎝⎛−=

⎪⎪

⎪⎪

⎬

⎪⎪

⎪⎪

⎨⎥⎥⎥

⎦

⎤

⎢⎢⎢

⎣

⎡

⎥⎥⎦

⎤

⎢⎢⎣

⎡⎟⎠⎞

⎜⎝⎛+−+

−

∞∞ dx

dp

dr

du

dr

du

dr

dn

2

11

21

2

0 ληηη Problema Não-linear

⎪⎭

⎪⎩

⎥⎦⎢⎣ 444444 3444444 21η

0=⇒= uRr

00 =⇒=dr

dur

i=1 i=2 i=Ni i+1i-1

x x1x 1−ix

ix 1+ixNx

de

Ponto i:

( ) FFdudud ed

n

ληηη2

12

1

−

−≈⎪

⎪

⎬

⎫

⎪⎪

⎨

⎧

⎥⎥⎤

⎢⎢⎡

⎥⎤

⎢⎡

⎟⎞

⎜⎛+−+ ( )

rdrdrdr

F

Δληηη 0 1∞∞ ≈

⎪⎪

⎭

⎬

⎪⎪

⎩

⎨⎥⎥⎦⎢

⎢⎣ ⎥

⎥⎦⎢

⎢⎣

⎟⎠

⎜⎝

++

444444 3444444 21

( ) uuuuF ii

n

iid ληηη 1

21

2

10 1 +

−

+∞∞

−⎥⎥⎤

⎢⎢⎡

⎥⎥⎤

⎢⎢⎡

⎟⎠⎞

⎜⎝⎛ −

+−+= ( )

( ) uuuuF

rr

ii

n

ii

d

λ

ΔΔηηη

1

21

2

1

0

1

−

∞∞

−⎥⎤

⎢⎡

⎥⎤

⎢⎡

⎟⎞

⎜⎛ −

⎥⎥⎦⎢

⎢⎣ ⎥

⎥⎦⎢

⎢⎣

⎟⎠

⎜⎝

( )r

uu

r

uuF iiii

e ΔΔληηη 11

0 1 −−∞∞

⎥⎥⎥

⎦⎢⎢⎢

⎣ ⎥⎥⎦⎢

⎢⎣

⎟⎠⎞

⎜⎝⎛+−+=

Sistema de equações não-linear

n

uuuu

r

uuR

⎥⎤

⎢⎡ ⎤⎡ ⎞⎛

−=

−2

12

121 Δ

( )

n

iiiii r

uu

r

uuR

⎤⎡ ⎤⎡

−−

⎥⎥⎥

⎦⎢⎢⎢

⎣ ⎥⎥⎦

⎤

⎢⎢⎣

⎡⎟⎠⎞

⎜⎝⎛ −

+−+=

−

++∞∞ 1

1

211

0 ΔΔληηη

( ) ii

n

ii Nidx

dp

r

uu

r

uu−=⎟

⎠⎞

⎜⎝⎛−−

−

⎥⎥⎥

⎦

⎤

⎢⎢⎢

⎣

⎡

⎥⎥⎦

⎤

⎢⎢⎣

⎡⎟⎠⎞

⎜⎝⎛ −

+−+ −−∞∞ 1,,3,2,

2

11

21

21

2

10 K

ΔΔληηη

NN uR =

Matriz Jacobiana do método de Newton.

≥==−= kJr

Jr

J kΔΔ ,12,11,1 )3(0;1

;1

( ) +⎟⎠⎞

⎜⎝⎛−⎥⎥⎤

⎢⎢⎡

⎥⎥⎦

⎤

⎢⎢⎣

⎡⎟⎠⎞

⎜⎝⎛ −

+−+=

−

+∞∞ rr

uuJ

n

iiii ΔΔ

ληηη2

21

2

10,

11

( ) +⎟⎞

⎜⎛⎟⎞

⎜⎛ −

⎥⎤

⎢⎡

⎟⎞

⎜⎛ −

+−−

⎠⎝⎥⎥⎦⎢

⎢⎣ ⎥⎦⎢⎣ ⎠⎝

+

−−

++ uuuunuu

rr

ii

n

iiii λλλ

ΔΔ

1

12

12

11 211( )

( ) ⎟⎞

⎜⎛⎥⎤

⎢⎡

⎥⎤

⎢⎡

⎟⎞

⎜⎛ −

+⎟⎠⎞

⎜⎝⎛−⎟⎠⎞

⎜⎝⎛

⎥⎥⎦⎢

⎢⎣

⎟⎠⎞

⎜⎝⎛+−

−

++∞

+

uu

rrrr

n

iiiiii

ΔΔλ

Δληη

Δ

21

2

1

1102

1

1

212

( )

⎤⎡

+⎟⎠⎞

⎜⎝⎛

⎥⎥⎥

⎦⎢⎢⎢

⎣ ⎥⎥⎦

⎤

⎢⎢⎣

⎡⎟⎠⎞

⎜⎝⎛+−+

−−

−∞∞ rr

uu

n

ii

ΔΔληηη

11

21

0

11

( ) ⎟⎠⎞

⎜⎝⎛⎟⎠⎞

⎜⎝⎛ −

⎥⎥⎦

⎤

⎢⎢⎣

⎡⎟⎠⎞

⎜⎝⎛ −

+−

−− −−

∞−

rr

uu

r

uun

r

uu iiiiii

Δλ

Δλ

Δληη

Δ1

122

102

1 212

1

( ) 11

2

21

2

101

uuJ

n

iiii +⎟

⎠⎞

⎜⎝⎛⎥⎥⎤

⎢⎢⎡

⎥⎤

⎢⎡

⎟⎠⎞

⎜⎝⎛ −

+−+=

−

+∞∞+ ληηη ( )

( ) 1

1

12

12

201,

uuuunuu

rrJ

n

ii

⎞⎛⎞⎛ −⎤⎡ ⎞⎛ −−−

+⎟⎠

⎜⎝⎥⎥⎦⎢

⎢⎣ ⎥

⎥⎦⎢

⎢⎣

⎟⎠

⎜⎝

++

−−

∞∞+

λ

ΔΔληηη

( )

1

212

1

21

2

1102

1

rr

uu

r

uun

r

uu

n

iiiiii

⎞⎛⎥⎤

⎢⎡ ⎤⎡ ⎞⎛

⎟⎠⎞

⎜⎝⎛⎟⎠⎞

⎜⎝⎛ −

⎥⎥⎦

⎤

⎢⎢⎣

⎡⎟⎠⎞

⎜⎝⎛ −

+−

−−

−

++∞

+

Δλ

Δλ

Δληη

Δ

( ) 11

1

21

01, rr

uuJ

n

iiii +⎟

⎠⎞

⎜⎝⎛−⎥⎥⎥

⎦⎢⎢⎢

⎣ ⎥⎥⎦

⎤

⎢⎢⎣

⎡⎟⎠⎞

⎜⎝⎛ −

+−+=

−

−∞∞− ΔΔ

ληηη

( ) 212

1 1

12

12

102

1

rr

uu

r

uun

r

uu ii

n

iiii ⎟⎠⎞

⎜⎝⎛−⎟⎠⎞

⎜⎝⎛ −

⎥⎥⎦

⎤

⎢⎢⎣

⎡⎟⎠⎞

⎜⎝⎛ −

+−

−− −

−−

−∞

−

Δλ

Δλ

Δληη

Δ

)(0;1

)1,,1(0,

NkJJ

iiikJ ki

≠==

+−≠=

)(0;1 ,, NkJJ kNNN ≠==


Prof. Márcio Carvalho [email protected] Laboratório: Método das Diferenças Finitas para EDO linear Escreva uma rotina MatLab para solução de um problema de valor de contorno pelo método de diferenças finitas (diferença central). Utilize a rotina desenvolvida para resolver o problema:

2

2

2

2

1;4

20;25;50

/200

/10

)(,)0(

0)(

mAmP

CTCTCT

CmWk

CmWh

TLTTT

TTkA

hP

dx

Td

ooB

oA

o

o

BA

==

===

=

=

==

=−−

∞

∞

4/17/11 9:48 PM C:\Marcio\Cursos-PUC\ENG1714\LABORATORIOS\lab10-ed...\cond1d.m 1 of 2

1 clc % Apaga tudo o escrito na janela de conandos 2 clear all % Apaga as variáveis usadas anteriornente 3 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 4 % PROBLEMA DE TRANSFERENCIA DE CALOR UNI-DIMENSIONAL 5 6 % 28 / 10 / 2004 7 8 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 9 % ENTRADA DE DADOS10 11 disp(' PROGRAMA DE TRANSFERENCIA DE CALOR UNI-DIMENSIONAL ');12 disp( ' ');13 disp( ' ');14 disp(' DADOS DE ENTRADA ');15 disp( ' ----- -- ------- ');16 disp( ' ');17 TA = input(' Temperatura na parede da esquerda: ');18 TB = input(' Temperatura na parede da direita: ');19 Tinf = input(' Temperatura do ar: ');20 disp( ' ');21 h = input(' Coeficiente de Troca de Calor Convectivo: ');22 k = input(' Condutividade Termica da Barra: ');23 Area = input(' Area Transversal da Barra: ');24 P = input(' Perimetro Transversal da Barra: ');25 disp(' ');26 L = input(' Comprimento da Barra: ');27 N = input(' Número de nós N= : ');28 disp( ' ');29 disp( ' ');30 31 Dx=L/(N-1);32 33 % Geração da malha para o grafico34 35 for i=1:N36 x(i,1) = (i-1)*Dx;37 end38 % Nos no interior do dominio39 40 C = h*P/(k*Area);41 42 for i = 2:N-143 A(i,i) = -2/Dx^2-C;44 A(i,i-1) = 1/Dx^2;45 A(i,i+1) = 1/Dx^2;46 b(i,1) = -C*Tinf;47 end48 49 % Condicoes de Contorno50 51 A(1,1) = 1;52 b(1,1) = TA;53 A(N,N) = 1;54 b(N,1) = TB;55

4/17/11 9:48 PM C:\Marcio\Cursos-PUC\ENG1714\LABORATORIOS\lab10-ed...\cond1d.m 2 of 2

56 % Solução do sistema de equações57 58 T=inv(A)*b;59 60 61 % Grafico do perfil de temperatura62 63 plot(x,T,'b:o');64 xlabel('Posicao X')65 ylabel('Temperatura [oC]')66

MEC 1701 – Métodos Numéricos para Engenharia Mecânica

Prof. Márcio Carvalho [email protected] Laboratório: Método das Diferenças Finitas para EDO não-linear. Um modelo uni-dimensional de um escoamento leva a seguinte equação diferencial:

2

2

dx

udK

dx

duu = , com as condições de contorno:

01

/100

=→===→=

ux

smVux.

Determine a velocidade em cada ponto do domínio (entre 0 e 1) para 01.01;10 === kekk .

4/17/11 9:49 PM C:\Marcio\Cursos-PUC\ENG1714\LABORATORIOS\la...\edonaolinear.m 1 of 2

1 clc % Apaga tudo o escrito na janela de conandos 2 clear all % Apaga as variáveis usadas anteriornente 3 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 4 % ESCOAMENTO DE UM LIQUIDO NAO NEWTONIANO EM UM TUBO 5 6 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 7 % ENTRADA DE DADOS 8 9 L = 1;10 Va = 10;11 Vb = 0;12 13 disp(' EDO NAO LINEAR ');14 disp( ' ');15 disp( ' ');16 disp(' DADOS DE ENTRADA ');17 disp( ' ----- -- ------- ');18 disp( ' ');19 disp( ' ');20 K = input(' Coeficiente (K): ');21 22 disp( ' ');23 24 N = input(' Número de nós N : ');25 disp( ' ');26 disp( ' ');27 28 Dx=L/(N-1);29 30 % Chute Inicial para o Perfil de Velocidade:31 % Perfil Parabolico de um Liquid Newtoniano32 33 for i=1:N34 x(i)=(i-1)*Dx;35 v(i,1) = 0;36 % v(i,1)=Va-(Va-Vb)/L * x(i);37 end38 39 erro = 1;40 ERRO_MAX = 10^-5;41 ITER_MAX = 20;42 43 % Calculo do vetor residuo inicial44 45 [res] = calc_residuo(N, Dx, K, Va, Vb, v);46 norm_res = norm(res);47 48 iter = 1;49 50 % Metodo de Newton51 52 while (norm_res > ERRO_MAX) & (iter < ITER_MAX)53 [J] = calc_jac(N, Dx, K, Va, Vb, v);54 55 DV = pinv(J)*(-res');

4/17/11 9:49 PM C:\Marcio\Cursos-PUC\ENG1714\LABORATORIOS\la...\edonaolinear.m 2 of 2

56 57 v = v + DV;58 59 [res] = calc_residuo(N, Dx, K, Va, Vb, v); 60 61 iter = iter+162 norm_res = norm(res)63 end64 65 if (norm_res > ERRO_MAX) 66 disp('NAO CONVERGIU ');67 else68 plot(x,v);69 end70 71 72 73 74 75

4/17/11 9:49 PM C:\Marcio\Cursos-PUC\ENG1714\LABORATORIOS\la...\calc_residuo.m 1 of 1

1 function [res] = calcresiduo(N, Dx, K, Va, Vb, v) 2 3 % funcao que calcula a matriz jacobiana 4 5 res(1) = v(1)-Va; 6 7 res(N) = v(N)-Vb; 8 9 for i=2:N-110 res(i) = v(i)*(v(i+1)-v(i-1))/(2*Dx) - K*(v(i+1)-2*v(i)+v(i-1))/Dx^2;11 end

4/17/11 9:49 PM C:\Marcio\Cursos-PUC\ENG1714\LABORATORIOS\lab11-...\calc_jac.m 1 of 1

1 function [J] = calcjac(N, Dx, K, Va, Vb, v) 2 3 % funcao que calcula a matriz jacobiana 4 5 J(1,1) = 1; 6 7 J(N,N) = 1; 8 9 for i=2:N-110 J(i,i-1) = -v(i)/(2*Dx) - K/Dx^2;11 12 J(i,i) = (v(i+1)-v(i-1))/(2*Dx) + 2*K/Dx^2;13 14 J(i,i+1) = v(i)/(2*Dx) - K/Dx^2;15 end

EQUAÇÃO DIFERENCIAL PARCIALEQUAÇÃO DIFERENCIAL PARCIAL

PROBLEMA UNIDIMENSIONAL TRANSIENTE

TdK

dTc

2

ρ

TT

TxtTdx

Kdt

c

==

=

)0(

),0( 0

2ρ

Condição inicial0T

t

BT

AT

B

A

TLxtT

TxtT

====

),(

)0,(Condições de contorno

?)( =xtT

0T

x

?),( =xtT

( ) tjt Δ1−= jjiT ,

0t

tt Δ=

1j

2=jjiT ,

Ponto i Tempo j0=t 1=j

1 2 1−i i 1−N N1+i

Discretização no espaço:2

112

2 2

x

TTT

x

T iii

iΔ

−+ +−≈

∂∂

Equação diferencial deve ser satisfeita em todos os pontos i :

iiii Nix

TTT

c

k

dt

dT−=

+−= −+

211 1,,2;

2K

Δρ

BN

A

TT

TTxcdt

==1

Δρ

Uma vez discretizada as derivadas em relação a x, obtém-se umsistema de equações diferenciais ordinárias em t (prob. de valor inicial):

⎧ = ATT1

⎪⎪⎪⎪⎪⎧

+−=

A

TTTkdTx

TTT

c

k

dt

dT2

1232

1

2

2

Δρ Incógnitas do problema:

⎪⎪⎪

⎪⎪⎨

+−=

x

TTT

c

k

dt

dT

M

22343 2

Δρ )(,,)(,)( 21 tTtTtT NK

⎪⎪⎪

⎩

= BN TT

Método Explícito jiT ,

Ponto i Tempo jPonto i Tempo j

1,,2;2

2

,1,,1,11, −=+−

=−

≈ −++ NiTTTkTTdT jijijijjii K,,;

2xctdt ΔρΔ

Lado direito da EDO avaliada no instante anterior

[ ] 1,,2;2 ,1,,12,1, −=+−+= −++ NiTTTc

k

x

tTT jijijijiji K

ρΔΔ

Quando um método explícito é usado, as temperaturas em todos os pontos ino intante j+1 são calculadas diretamente em função das temperaturas nospontos i no instante j, conhecidas.

O método de Euler explícito é instável se2

2 2

1

2

1x

c

kt

c

k

x

t Δρ

ΔρΔ

Δ>⇒>

22 ccx ρρΔ

O passo de tempo tem que ser muito pequeno e função da discretização em x

TddT 2

Exemplo usando Excel:

TxtTdx

TdK

dt

dTc

==

=

),0( 0

2ρ

1;2;0

1.0;1;1 ===

TTT

xLc

k Δρ

B

A

TLxtT

TxtT

====

),(

)0,( 1;2;00 === BA TTT

[ ]ktΔ [ ] 1,,2;2 ,1,,12,1, −=+−+= −++ NiTTTc

k

x

tTT jijijijiji K

ρΔΔ

Planilha do Mi ft E lMicrosoft Excel

01.0;1.0;1 === txc

k ΔΔρ

6

8

10

t = 0

t=0.03

t=0 05

2

11

2>=

c

k

x

t

ρΔΔ

0

2

4

6

T

t=0.05

-6

-4

-20 0.2 0.4 0.6 0.8 1 1.2

x

005.0;1.0;1 === txc

k ΔΔρ

2.5

t = 0

2

12

=c

k

x

tc

ρΔΔρ

1.5

2

T

t 0

t=0.15

t=0.025

0.5

1

0

0 0.2 0.4 0.6 0.8 1 1.2

x

2.5t = 0

t=0.15

t=0 025Solução em regime permanente

2t=0.025

t=0.1

t=0.15

1

1.5

T

t=0.4

0.5

0

0 0.2 0.4 0.6 0.8 1 1.2

x

Método Implícito jiT ,

Ponto i Tempo jPonto i Tempo j

1,,2;2

2

1,11,1,1,1, −=+−

=−

≈ +−++++ NiTTTkTTdT jijijijijii K,,;

2xctdt ΔρΔ

Lado direito da EDO avaliada no instante atual

Conhecida as temperaturas no instante j, deseja-se determinar as temperaturas no instante j+1.

jijijiji

Aj

TTk

Tk

Tk

TT

=⎟⎟⎠

⎞⎜⎜⎝

⎛−+⎟⎟

⎠

⎞⎜⎜⎝

⎛++⎟⎟

⎠

⎞⎜⎜⎝

⎛−

=

+−+++

+

,1,121,21,12

1,1

;11211

ΔΔΔΔΔ

BjN

jijijiji

TT

Ni

txcxctxc

=−=

⎟⎠

⎜⎝

⎟⎠

⎜⎝

⎟⎠

⎜⎝

+

++++

1

,1,121,21,12

1,,2 K

ΔΔρΔρΔΔρ

BjN +1,

Sistema de equações linear.

Para cada instante de tempo, deve-se resolver um sistema de equações linear.

fxA =

Função das temperaturas no instante anteriorMatriz dos coeficientes

Vetor com as temperaturas em todos os nós no instante atual.

PROBLEMA BIDIMENSIONAL EM REGIME PERMANENTE

Parede isolada

)0(

02

2

2

2

=∂∂

+∂∂

TTy

T

x

TParede isolada

)0,(

),(

),0(

======

TyxT

TyLxT

TyxT

C

B

A

AT BT?),( =yxT

0),(

)0,(

==∂∂

Lyxy

T

TyxT C

TCT

jiT ,

Ny

j1+j

1−jCoord x Coord y

i 1+i1−i Nx1=i1=j

),( ji

)1,( +ji

),1( ji − )1( ji +2

,1,,1

2

2 2

x

TTT

x

T jijiji

jiΔ

−+ +−≈

∂∂

),( ji),( j ),1( ji +

)1,( −ji2

1,,1,

,

2

2

,

2

y

TTT

y

T jijiji

ji

ji

Δ−+ +−

≈∂∂

Equação algébrica resultante no ponto (i,,j) :1,2

1,,2

−=−=

Nyj

Nxi

K

K

011

21111

,221,21,2,12,12=⎟⎟

⎠

⎞⎜⎜⎝

⎛+−⎟⎟

⎠

⎞⎜⎜⎝

⎛+⎟⎟

⎠

⎞⎜⎜⎝

⎛+⎟

⎠⎞

⎜⎝⎛+⎟

⎠⎞

⎜⎝⎛

−+−+ jijijijiji Tyx

Ty

Ty

Tx

Tx ΔΔΔΔΔΔ

Condições de contorno:

1; == NyjTT

1,,2;

,,1;

,,1;

1

,

,1

−==

==

==

NxiTT

NyjTT

NyjTT

Ci

BjNx

Aj

K

K

K

1,,2;0

,,;

1,,

1,

−==− − NxiTT NyiNyi

Ci

K

As equações algébricas devem ser escritas em forma matricial

fxA =

Termo independenteMatriz dos coeficientes

⎥⎥⎥⎤

⎢⎢⎢⎡

==

TT

TT

2,12

1,11

Vetor com as temperaturas em todos os nós.

⎥⎥⎥⎥

⎦⎢⎢⎢⎢

⎣ =

==

TT

TTx 3,13

M

⎥⎦⎢⎣ =× NyNxNyNx TT ,

As incógnitas do problema devem ser numeradas de forma sequencialpara escrevermos o sistema de equações em forma matrical.

510

25

Exemplo de uma numeração sequencial:

TT

TT

212

1,11

==

4

10

924

NyNy TT ,1

2,12

=M

2

38

22

23 Ny

TT

TT 1,21

=

=+

M13

),( ji

1

2

6

7

21

22

NyNxNyNx

jijNyi

TT

TT

,

,)1(

=

=

×

+×−

M

6

Seguindo esta regra, a numeração sequencial do nó (i,j) é dada por:

i 1 l d ó d d ó (i j) d l i N ói-1 colunas de nós a esquerda do nó (i,j) ; cada coluna possui Ny nós:

{jNyiji +×−⇒43421

)1(),(

Número de nós nas colunas anteriores Número de nós na coluna i

Equação relativa ao nó # 8:

011

21111 ⎟

⎞⎜⎛

⎟⎞

⎜⎛

⎟⎞

⎜⎛

⎟⎞

⎜⎛

⎟⎞

⎜⎛

TTTTT

3,28 TT ⇔

02 3,2222,224,223,123,32

⎞⎛⎞⎛⎞⎛

=⎟⎟⎠

⎞⎜⎜⎝

⎛+−⎟⎟

⎠

⎞⎜⎜⎝

⎛+⎟⎟

⎠

⎞⎜⎜⎝

⎛+⎟

⎠⎞

⎜⎝⎛+⎟

⎠⎞

⎜⎝⎛

Tyx

Ty

Ty

Tx

Tx ΔΔΔΔΔΔ

011

21111

822729232132=⎟⎟

⎠

⎞⎜⎜⎝

⎛+−⎟⎟

⎠

⎞⎜⎜⎝

⎛+⎟⎟

⎠

⎞⎜⎜⎝

⎛+⎟

⎠⎞

⎜⎝⎛+⎟

⎠⎞

⎜⎝⎛

Tyx

Ty

Ty

Tx

Tx ΔΔΔΔΔΔ

251413121110987654321 K

=A

0000000000 KABCBALinha 8

Matrix é pentadiagonal

0)0(

20),(

10),0(

======

TT

TyLxT

TyxT

B

A

0),(

0)0,(

==∂∂

===

Lyxy

T

TyxT C

Gráfico de iso-linhas contourf

Gráfico 3D - superfície surf

20

251.6

1.8

2

5

10

15

20

0.8

1

1.2

1.4

1 521.5

2-5

0

0

0.2

0.4

0.6

00.5

11.5

0

0.5

1 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 20


Prof. Márcio Carvalho [email protected]

Laboratório: Equação Diferencial Parcial – Problema Transiente 1. Deseja-se analisar a torção de um eixo AB com um volante na extremidade B, submetido a

um torque M*. Parâmetros físicos e geométricos: eixo: • comprimento L = 1 m • momento de inércia polar J

4

2 barrarJ

π= (rbarra=5 cm)

• coeficiente de cisalhamento G G= E/[2(1+ν)] ; E = 200 106 Pa, ν=0,29

• massa específica ρ= 7800 kg/m3 volante: • momento de inércia polar massa – ID

ID= ρ Jvoltante ; 4

2 volanterJvolante

π=

• raio rvolante = 12 cm

M*

rL

Equilíbrio: 0 < x < L ⇒ 2

2

tJ

xJG

x ∂

∂=⎟⎟

⎠

⎞⎜⎜⎝

⎛∂∂

∂∂ θρθ

x = 0 ⇒ θ = 0

x = L ⇒ 2

2*

tI

xJGM D

∂

∂=

∂∂

−θθ

Determine numericamente θ (x,t), do instante inicial até t = 1 segundo,. para um momento M* = 10 Nm. Utilize o método implícito que desejar. Trace um gráfico com a variação de θ com o tempo para 4 posições ao longo do eixo (x=L/4; x=L/2; x=3 L/4 e x = L). Trace também, a distribuição de θ (x,t) para t=1 segundo. Verifique se sua solução independe da malha e passo de tempo utilizado (isto é, resolva para umas 2 malhas e 2 passos de tempo diferentes e compare as soluções).

mailto:[email protected]

4/17/11 9:51 PM C:\Marcio\Cursos-PUC\ENG1714\LABORATORIOS\l...\edptransiente.m 1 of 3

1 clc % Apaga tudo o escrito na janela de conandos 2 clear all % Apaga as variáveis usadas anteriornente 3 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 4 % PROGRAMA DE DIFERENÇAS FINITAS 5 6 % REGISTRO DE REVISÕES 7 % Fecha Programador Revisado por Descrição da mudança 8 % ----- ----------- ------------ -------------------- 9 % 30/11/01 O. Coronado M. Carvalho Codigo Original 10 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 11 % ENTRADA DE DADOS 12 13 disp(' PROGRAMA DE REGIME TRANSIENTE '); 14 disp( ' -------- -- ------ ---------- '); 15 disp( ' '); 16 disp( ' '); 17 disp(' DADOS DE ENTRADA '); 18 disp( ' ----- -- ------- '); 19 disp( ' '); 20 disp( ' '); 21 disp( ' '); 22 N=input('Entre o numero de nós N= : '); 23 disp( ' '); 24 IT=input('Entre o numero de intervalos de tempo IT= : '); 25 disp( ' '); 26 % Dados do problema 27 L=1; % Comprimento da barra 28 rb=0.05; % Raio da barra 29 rv=0.12; % Raio do volante 30 Rho=7800; % Massa especifica 31 E=200e6; 32 v=0.29; 33 % Calculo das constantes 34 J=pi*rb^4/2; % Momento de inercia polar 35 G=E/(2*(1+v)); % Coeficiente de cisalhamento 36 Id=Rho*pi*rv^4/2; % Momento de inercia polar massa 37 % Condicoes de contorno 38 disp( ' '); 39 disp( ' '); 40 disp( ' CONDIÇÕES DE CONTORNO'); 41 disp( ' --------- -- -------'); 42 disp( ' '); 43 disp( ' '); 44 t=1; % Tempo maximo 45 T0=input('Entre o angulo teta(x=0,t)= : '); 46 disp( ' '); 47 M=input('Entre o torque : '); 48 disp( ' ');


49 50 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 51 % SOLUÇÃO PELO METODO IMPLICITO 52 53 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 54 Dt=t/IT; %Passo do tempo 55 Dx=L/(N-1); %Passo do comprimento 56 57 x(1) = 0; 58 for i=2:N 59 x(i) = x(i-1)+Dx; 60 end 61 62 %Condições iniciais para t=0 63 64 for i=1:N 65 teta(i,1)=0; 66 teta(i,2)=0; 67 end 68 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 69 70 for j=3:IT+2 71 72 % Condicoes de contorno 73 74 A(1,1)=1; 75 f(1,1)=0; 76 A(N,N-1)=-1/Dx; 77 A(N,N)=1/Dx+Id/(G*J*Dt^2); 78 f(N,1)=M/(G*J)-Id/(G*J*Dt^2)*(-2*teta(N,j-1)+teta(N,j-2)); 79 80 for i=2:N-1 81 A(i,i+1)=1; 82 A(i,i)=-(2+Rho*Dx^2/(G*Dt^2)); 83 A(i,i-1)=1; 84 f(i,1)=Rho*Dx^2/(G*Dt^2)*(-2*teta(i,j-1)+teta(i,j-2)); 85 end 86 tetavetor=inv(A)*f; 87 88 for i=1:N 89 teta(i,j)=tetavetor(i); 90 end 91 92 end 93 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 94 plot(x,teta(:,IT+2)); 95 title(' Regime transiente'); 96 xlabel(' x ' ); 97 ylabel(' Teta ' ); 98 99 t(1)=0;100 for i=2:IT+2101 t(i)=Dt*(i-1);102 end103


104 figure105 a=round((N-1)/4);106 plot(t,teta(a+1,:),t,teta(2*a+1,:),t,teta(3*a+1,:),t,teta(N,:))107 title(' Regime transiente');108 figure(2);109 xlabel(' t ' );110 ylabel(' Teta ' );111 legend('L/4','L/2','3L/4','L');112 113 114 115 116 117 118 119 120 121 122 123

MEC 1701 Prof. Márcio Carvalho [email protected] As equações diferenciais que descrevem o escoamento tridimensional em regime permanente de um líquido Newtoniano em um duto retangular são apresentadas abaixo:

⎟⎟⎠

⎞⎜⎜⎝

⎛

∂∂

+∂∂

+∂∂

+∂∂

−=⎥⎦

⎤⎢⎣

⎡∂∂

+∂∂

+∂∂

⎟⎟⎠

⎞⎜⎜⎝

⎛

∂∂

+∂∂

+∂∂

+∂∂

−=⎥⎦

⎤⎢⎣

⎡∂∂

+∂∂

+∂∂

⎟⎟⎠

⎞⎜⎜⎝

⎛

∂∂

+∂∂

+∂∂

+∂∂

−=⎥⎦

⎤⎢⎣

⎡∂∂

+∂∂

+∂∂

2

2

2

2

2

2

2

2

2

2

2

2

2

2

2

2

2

2

Re

Re

Re

z

w

y

w

x

w

z

p

z

ww

y

wv

x

wu

z

v

y

v

x

v

y

p

z

vw

y

vv

x

vu

z

u

y

u

x

u

x

p

z

uw

y

uv

x

uu

Após o desenvolvimento, a pressão só varia na direção do escoamento (direção z) e as componentes da velocidade na direção x e y são nulas. Desta forma, a equação que descreve o movimento pode ser simplificada para

⎟⎟⎠

⎞⎜⎜⎝

⎛

∂∂

+∂∂

+∆∆

−=2

2

2

2

0y

w

x

w

z

p

Para um gradiente de pressão adimensional de 10−=∆∆

z

p, determine:

a) O perfil de velocidade w(x,y). Apresente o resultado na forma de um gráfico de w(x,y). b) A velocidade máxima VMax e a velocidade média VMED.

c) O fator de atrito X número de Reynolds; isto é:µ

ρ hDVf ×

O fator de atrito do escoamento pode ser definido como: 2

2

1V

DL

P

fh

ρ

∆

= , onde P

ADh

4= é o diâmetro

hidráulico e V a velocidade média do escoamento.

OBS: O produto µ

ρ hDVf × para o escoamento em um duto circular é igual a 64.

d) Utilize o programa desenvolvido para calcular o fator de atrito em função H / W. Em seus cálculos mantenha a área da seção reta constante e igual a 4=×WH . Determine a razão que leva ao menor valor do fator de atrito do escoamento.

x

zy

H=1

W=4

4/17/11 9:52 PM C:\Marcio\Cursos-PUC\ENG1714\LABORATORIOS\lab13-edp\labedp.m 1 of 3

1 clc % Apaga tudo o escrito na janela de conandos 2 clear all % Apaga as variáveis usadas anteriornente 3 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 4 % PROGRAMA DE DIFERENÇAS FINITAS 5 6 % REGISTRO DE REVISÕES 7 % Fecha Programador Revisado por Descrição da mudança 8 % ----- ----------- ------------ -------------------- 9 % 13/11/01 O. Coronado M. Carvalho Codigo Original 10 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 11 % ENTRA5DA DE DADOS 12 13 disp( ' PROGRAMA DE DIFERENÇAS FINITAS '); 14 disp( ' -------- -- ----------- ------- '); 15 disp( ' '); 16 disp( ' '); 17 disp(' DADOS DE ENTRADA '); 18 disp( ' ----- -- ------- '); 19 disp( ' '); 20 disp( ' '); 21 Wduto=input('Entre a largura do duto: '); 22 disp( ' '); 23 Hduto=input('Entre a altura do duto : '); 24 disp( ' '); 25 Nx=input('Entre o número de nós da direção x e y Nx=Ny= : '); 26 disp( ' '); 27 Ny=Nx; 28 disp( ' '); 29 Dx=Wduto/(Nx-1); 30 Dy=Hduto/(Ny-1); 31 dpdz=10 32 33 for ix=2:Nx-1 34 for jy=2:Ny-1 35 i=(ix-1)*Ny+jy; % Linha da matriz A 36 il=(ix+1-1)*Ny+jy; % Leste 37 iw=(ix-1-1)*Ny+jy; % Oeste 38 in=(ix-1)*Ny+jy+1; % Norte 39 is=(ix-1)*Ny+jy-1; % Sul 40 % Matriz dos coeficentes 41 A(i,i)=-2*(1/Dx^2+1/Dy^2); 42 A(i,il)=1/Dx^2; 43 A(i,iw)=1/Dx^2; 44 A(i,in)=1/Dy^2; 45 A(i,is)=1/Dy^2; 46 f(i,1)=-dpdz; 47 end 48 end


49 50 51 % CONDIÇÕES DE CONTORNO 52 53 for ix=1:Nx 54 jy = 1; 55 i=(ix-1)*Ny+jy; 56 A(i,i)=1; 57 f(i,1)=0; 58 59 jy = Ny; 60 i=(ix-1)*Ny+jy; 61 A(i,i)=1; 62 f(i,1)=0; 63 end 64 65 for jy=1:Ny 66 ix = 1; 67 i=(ix-1)*Ny+jy; 68 A(i,i)=1; 69 f(i,1)=0; 70 71 ix = Nx; 72 i=(ix-1)*Ny+jy; 73 A(i,i)=1; 74 f(i,1)=0; 75 end 76 77 78 % Solução do sistema de equações 79 80 WW=inv(A)*f; 81 82 % Geração da malha 83 x=0:Dx:Wduto; 84 y=0:Dy:Hduto; 85 [X,Y]=meshgrid(x,y); 86 87 % Armazenando os valores de WW na matriz W 88 89 for i=1:Nx 90 for j=1:Ny 91 W(i,j)=WW(j+(i-1)*Ny); 92 end 93 end 94 95 % Calculo da velocidade media 96 97 qtotal = 0; 98 99 for i = 1:Nx-1100 for j = 1:Ny-1101 102 w1 = W(i,j);103 w2 = W(i,j+1);


104 w3 = W(i+1,j);105 w4 = W(i+1,j+1);106 qelem = ((w1+w2+w3+w4)/4)*Dx*Dy;107 qtotal = qtotal + qelem;108 109 end110 end111 112 Wmed = qtotal/(Wduto*Hduto);113 114 115 % Geração do garfico 3D116 117 surf(X,Y,W);118 119 % Calculo do fator de atrito120 121 Dhid = 4*(Hduto*Wduto)/(2*(Hduto+Wduto));122 123 rho = 1.;124 mu = 1;125 126 fatrito = (-dpdz*Dhid)/(0.5*rho*Wmed^2);127 128 Rey = rho*Wmed*Dhid/mu;129 130 frey = fatrito*Rey;131 132 133 disp('W Medio =')134 disp(Wmed);135 136 137 disp('F*REYNOLDS =')138 disp(frey);139

Documents

apostila metodos