Calculo Numerico
35 pág.

Calculo Numerico


DisciplinaCálculo Numérico11.996 materiais247.718 seguidores
Pré-visualização35 páginas
gra´fica e´ mostrada na figura 7.3; note que quando h cresce o

erro de truncamento aumenta e, por outro lado, o erro de arredondamento diminui.

Logo, tal fato deve ser considerado para a escolha de h; nesse caso, h\u2217 e´ o melhor
espac¸amento.

7.5.1 Regia\u2dco de estabilidade de alguns me´todos

Uma maneira de analisar a estabilidade de um me´todo nume´rico consiste em esta-

belecer regio\u2dces do plano complexo constru´\u131das com base na seguinte definic¸a\u2dco:

Fund. de Ca´lculo Nume´rico para Engenheiros 175

Figura 7.3: Representac¸a\u2dco gra´fica dos erros na soluc¸a\u2dco nume´rica de EDOs

\u2022 Um me´todo de aproximac¸a\u2dco e´ dito ser absolutamente esta´vel num ponto \u3bbh do
plano complexo se a sequ¨e\u2c6ncia {yk} gerada pelo me´todo aplicado a` EDO de re-
fere\u2c6ncia

dy

dx
= \u3bby, \u3bb = constante qualquer (7.17)

com passo \u2206x = h, for limitada, isto e´, yk \u2192 0 quando xk \u2192\u221e.

Resulta, para os me´todos de Euler e Runge-Kutta:

a) Euler

Resolvendo por Euler a equac¸a\u2dco de refere\u2c6ncia (7.17) obte´m-se

yk+1 = yk + h(\u3bbyk) = yk(1 + \u3bbh)

yk+2 = yk(1 + \u3bbh)
2

...
...

...

yk+n = yk(1 + \u3bbh)
n

176 Cap´\u131tulo 7 - Soluc¸a\u2dco Nume´rica de EDO\u2019s

A sequ¨e\u2c6ncia {yk} e´ limitada se |1 + \u3bbh| \u2264 1. Enta\u2dco, a regia\u2dco de estabilidade e´ o
conjunto dos complexos \u3bbh tais que |1+\u3bbh| \u2264 1. Assim \u22122 \u2264 \u3bbh \u2264 0. Na figura 7.4
temos a representac¸a\u2dco geome´trica da estabilidade do me´todo de Euler.

Figura 7.4: Regia\u2dco de estabilidade para o me´todo de Euler

b) Runge-Kutta

A regia\u2dco de estabilidade do me´todo de Runge-Kutta para a equac¸a\u2dco de refere\u2c6ncia

(7.17) e´ o conjunto dos nu´meros complexos \u3bbh, tais que\u2223\u2223\u2223\u22231 + \u3bbh+ (\u3bbh)22
\u2223\u2223\u2223\u2223 \u2264 1 \u2192 Runge\u2212Kutta ordem 2\u2223\u2223\u2223\u22231 + \u3bbh+ (\u3bbh)22 + (\u3bbh)36
\u2223\u2223\u2223\u2223 \u2264 1 \u2192 Runge\u2212Kutta ordem 3\u2223\u2223\u2223\u22231 + \u3bbh+ (\u3bbh)22 + (\u3bbh)36 + +(\u3bbh)424
\u2223\u2223\u2223\u2223 \u2264 1 \u2192 Runge\u2212Kutta ordem 4

e graficamente resulta

Observe que o cruzamento das linhas com o eixo dos y (complexo) fornecera´ o

nu´mero de Courant-Friedrich-Lewwi - CFL limite que pode ser utilizado.

Fund. de Ca´lculo Nume´rico para Engenheiros 177

Figura 7.5: Regia\u2dco de estabilidade para o me´todo de RK2, RK3 e RK4, res-
pectivamente

7.6 Aplicac¸o\u2dces

As aplicac¸o\u2dces que sa\u2dco apresentadas a seguir se referem a problemas da f´\u131sica, en-

genharia e biomedicina, dentre outras. Muitas destas apresentam soluc¸a\u2dco anal´\u131tica e

visam introduzir o aluno ou leitor a` soluc¸a\u2dco nume´rica envolvendo problems te´cnicos.

Para problemas mais complexos a ide´ia ba´sica aqui apresentada continua va´lida.

7.6.1 Sistema massa-mola

Considere um corpo oscilando no sentido vertical, preso a` extremidade de uma

mola. Admitindo que nenhuma forc¸a externa atue sobre o corpo para manter o

movimento, depois do mesmo iniciado, que a outra extremidade do mola esteja fixa

e que a massa do mesmo e a resiste\u2c6ncia de ar sejam desprezadas, o movimento do

corpo sera´ do tipo movimento harmo\u2c6nico simples:

y = A cos\u3c9t+B sin\u3c9t

178 Cap´\u131tulo 7 - Soluc¸a\u2dco Nume´rica de EDO\u2019s

onde y e´ o deslocamento do corpo no tempo t, a partir da sua posic¸a\u2dco inicial, conforme

figura 7.6. Para esse movimento:

Figura 7.6: Movimento harmo\u2c6nico simples

1. a amplitude ou deslocamento ma´ximo a partir da posic¸a\u2dco de equil´\u131brio e´
\u221a
A2 +B2,

porque quando dydt = 0, tg\u3c9t =
A
B e x =

\u221a
A2 +B2;

2. o per´\u131odo ou tempo necessa´rio para uma oscilac¸a\u2dco completa e´ 2pi\u3c9 s, porque quando

t varia de 2pi\u3c9 s os valores de x e
dx
dt permanecem invaria´veis, enquanto que qualquer

variac¸a\u2dco de t menor do que aquele valor acarreta variac¸a\u2dco em x ou dxdt ;

3. a frequ¨e\u2c6ncia ou nu´mero de oscilac¸o\u2dces por segundo e´ \u3c92pi ciclos/s;

4. a equac¸a\u2dco diferencial do movimento harmo\u2c6nico simples e´

m
d2x

dt2
= \u2212k(x), (7.18)

onde k e´ a constante ela´stica da mola e h(t) a forc¸a externa aplicada.

Fund. de Ca´lculo Nume´rico para Engenheiros 179

A equac¸a\u2dco diferencial ordina´ria (7.18) pode ser aproximada em diferenc¸as finitas

centrais conforme

m
yt+\u2206t \u2212 2yt + yt\u2212\u2206t

\u2206t2
+ kyt = h(t)

= sen\u3c9t

o que resulta em, no caso expl´\u131cito

yt+\u2206t = 2yt \u2212 yt\u2212\u2206t +\u2206t2
[
\u2212 k
m
yt +

1

m
sen\u3c9t

]
Note que para t = 0 teremos yt = yt\u22121 = 0.

Exemplo 7.8: Uma certa mola, cuja constante e´ k = 48kg/m, e´ mantida na verti-

cal, estando sua extremidade superior presa a um suporte. Um corpo, pesando 5kg,

e´ amarrado a` extremidade inferior da mola. Quando em repouso, o corpo e´ puxado

2cm para baixo e em seguida solto. Desprezando a resiste\u2c6ncia do ar e considerando

g = 10m/s2 (gravidade local), aproxime a soluc¸a\u2dco utilizando o me´todo de diferenc¸as

finitas centrais.

Soluc¸~ao: Considere a origem no centro de gravidade do corpo, com sistema em

repouso, e denomine de y o seu deslocamento no instante t; admita que y e´ positivo

para baixo. Com o sistema em repouso, a forc¸a da mola e´ igual e oposta a` da

gravidade. No instante t a forc¸a e´ \u2212ky, correspondente ao deslocamento y do corpo.
Enta\u2dco, tem-se

P

g

d2y

dt2
= \u2212ky \u21d2 d

2y

dt2
+ 96y = 0 (7.19)

Integrando, resulta y = C1sen
\u221a
96t+C2cos

\u221a
96t. Derivando em relac¸a\u2dco a t teremos

a velocidade do sistema dada por

dy

dt
=

\u221a
96(C1cos

\u221a
96t\u2212 C2sen

\u221a
96t).

180 Cap´\u131tulo 7 - Soluc¸a\u2dco Nume´rica de EDO\u2019s

Para t = 0, y = 1/6 e y\u2019= 0 tem-se C1 = 1/6, C2 = 0 e y = 1/6cos
\u221a
96t.

Discretizando a equac¸a\u2dco (7.19) por diferenc¸as finitas centrais obte´m-se

yt+\u2206t \u2212 2yt + yt\u2212\u2206t
\u2206t2

+ 96yt = 0

ou seja,

yt+\u2206t = 2yt \u2212 yt\u2212\u2206t \u2212 96\u2206t2yt

y(0) = 1/6

y\u2032(0) =
1

6
cos

\u221a
96t.

Utilizando passo de tempo \u2206t = 10\u22123, tem-se a soluc¸a\u2dco fornecida na tabela (7.2)
obtida pelo me´todo de diferenc¸as finitas centrais comparada com a anal´\u131tica (exata).

Tabela 7.2: Soluc¸a\u2dco do sistema massa-mola via diferenc¸as finitas centrais
t(s) y yexata t(s) y yexata
0.000 .1663 .1667 0.006 .1649 .1664
0.001 .1661 .1667 0.007 .1646 .1663
0.002 .1659 .1666 0.008 .1643 .1662
0.003 .1657 .1666 0.009 .1640 .1660
0.004 .1654 .1665 0.010 .1636 .1659
0.005 .1652 .1664

7.6.2 Vigas horizontais

A situac¸a\u2dco consiste em se determinar a deflexa\u2dco (flexa\u2dco) de uma viga sujeita as

cargas; considere somente vigas homoge\u2c6neas, quanto ao material, e uniformes. Ad-

mita que a viga seja formada por fibras longitudinais. Considere ainda que na flexa\u2dco,

Fund. de Ca´lculo Nume´rico para Engenheiros 181

vista na Fig. (7.7), as fibras da metade superior sa\u2dco comprimidas e as da metade

inferior sa\u2dco tracionadas; as duas partes sendo separadas por uma superf´\u131cie neutra

cujas fibras na\u2dco sofrem trac¸a\u2dco nem compressa\u2dco.

Figura 7.7: Representac¸a\u2dco esquema´tica da flexa\u2dco de uma viga

A fibra, que inicialmente coincide com o eixo da viga, encontra-se na superf´\u131cie

neutra apo´s a deformac¸a\u2dco (curva ela´stica ou curva das deflexo\u2dces) [2]. Para determinar

a equac¸a\u2dco desta curva considere uma sec¸a\u2dco transversal da viga, a uma dista\u2c6ncia x

de uma extremidade, conforme Fig. (7.7). Seja AB sua intersecc¸a\u2dco com a superf´\u131cie

neutra e P o trac¸o da curva ela´stica nessa sec¸a\u2dco. A meca\u2c6nica demonstra que o

momento fletor (torque) M , em relac¸a\u2dco a AB, de todas as forc¸as que agem em

qualquer das partes em que a viga foi dividida pela sec¸a\u2dco feita e´:

a) independente da parte considerada,

b) e dado por
EI

R
=M, (7.20)

onde E e´ o mo´dulo de elasticidade da viga, I o momento de ine´rcia da sec¸a\u2dco transver-

sal em relac¸a\u2dco a AB, e R e´ o raio de curvatura da curva ela´stica no ponto P .

182 Cap´\u131tulo 7 - Soluc¸a\u2dco Nume´rica de EDO\u2019s

Considere a origem na extremidade esquerda da viga, o eixo dos x na horizontal

e o ponto P com as coordenadas (x, y). Como a inclinac¸a\u2dco dydx da curva ela´stica, em

qualquer ponto, e´ uma quantidade necessariamente pequena, pode-se escrever

R =

[
1 +

(
dy
dx

)2]3/2
d2y
dx2

\u2248 1
d2y
dx2

,

e a equac¸a\u2dco 7.20 reduz-se a`

EI