14

MÉTODOS NUMÉRICOS PARA INTEGRAÇÃO

DE EQUAÇÕES DIFERENCIAIS ORDINÁRIAS

A integraçào numérica das equações de movimento pode ser realizada deduas formas:

• criação de seu próprio algoritmo de integraçào, utilizando algum método

numérico desenvolvido em uma linguagem qualquer de programação; e

• utilização de algum pacote de simulação comercialmente disponível.

Neste capítulo abordam-se técnicas para integração numérica de equações

diferenciais ordinárias, que permitem a simulação temporal de sistemas físicos

em computadores. Os métodos aqui apresentados sào gerais, à medida em que

se aplicam a modelos compostos por qualquer nÚmero de equações de movimcnto,lineares ou não.

14.1 SOLUÇÃO DE EQUAÇÕES DIFERENCIAIS ORDINÁRIAS. POR MÉTODOS NUMÉRICOS

As classesfundamcntais de problemas envolvendo equações diferencias ordi­

nárias são mostradas abaixo (Tao, 1988):

• problemas de valor inicial: são sistemas espacialmente homogêneos cujas pro­

priedades são assumidas como uniformes, tratados como sistemas a parâmetros

concentrados e que variam no tempo; e

394 • Nlodelagem e Simulaçáo de Processos Industriais e de Sistemas Eletromecânicos

• problemas de valores de contorno: envolvem sistemas em estado estacionário

que con têm gradien tes in ternos (variações espaciais).

Para resolver problemas de valor inicial há necessidade elovalor das variáveis

no instante inicial, ao passo que para resolver problemas de valor de contornosão necessários os valores elasvariáveis nas fronteiras do sistema.

14.1.1 Tipos de iVlétodos NwnéricosjJaTClResolver Equações Diferenciais Ordiná'lias

Os seguin tes métodos são utilizados na solução de equações diferencias ordi­

nárias (Tao, 1988):

• métodos de passo simples; e

• métodos de passo mÚltiplo.

Ambos os métodos estimam a solução como sendo uma série de valores a

intervalos específicos que geram uma função que satisfaz a equação diferencial.

CARACTERÍSTICAS DOS ivIÉTODOS DE PASSO SIMPLES

São métodos que necessitam apenas de um valor para estimar a solução em

cada intervalo de tempo. Essa estimativa é então usada para encontrar o valor da

função no próximo intervalo de tempo. Exemplo: a fórmula para encontrar )11/

depende explici tamente apenas de )1"1/-1'

CARACTERÍSTICAS DOS ivIÉTODOS DE PASSO MÚLTIPLO

São métodos em que o cálculo de )11/ depende explicitamente de dois ou

mais valores anteriores. Por exemplo, no método de passo duplo, o cálculo de)ln

depende dos seguintes valores: )1/7 = f( )l1/-l,)ln-2)' Esses métodos são também

conhecidos como predi tor-corretor.

COivlPARAÇi\.O ENTRE lvlÉTODOS DE PASSO SrrVIPLES/MÚLTIPLO

• ambos são precisos, embora os métodos de passo mÚltiplo normalmente apre­

sentem melhores resultados quando a função é descontínua em certos pontos

ou suas derivadas variam muito rapidamente (sistemas stiff);

• a maioria elos problemas de valor inicial não têm mÚltiplas condições iniciais,

sendo necessário gerá-Ias por um método de passo simples para inicializar o

método de passo mÚltiplo; e

Nlétodos Numéricos jJam Integraçâo de Equações Diferenciais Ordinárias • 395

• os requisitos computacionais (capacidade de memória) dos métodos de passo

mÚltiplo normalmente excedem os de métodos de passo simples.

Neste livro são abordados apenas os métodos de passo simples e se con­

centra nos problemas de simulação dinâmica a parâmetros concentrados (pro­blemas de valor inicial).

14.1.2 Nlétodos Numéricos de Integraç:ão de Passo Si1'njJles jJara Resolver Problel1w.s

de Valor Inicial

Discute-se aqui a solução para sistemas de 1ª ordem (Carnahan et al., 1969).

Sistemas de ordem superior são resolvidos conforme explicado no item 14.1.8.

Forma geral de sistemas de 1;;ordem:

y = f(t, y)

Objetivo:

com y( O ) = Yo

• estimar y( t); ou

• estimar t corresponden te a um dado valor de y.

o método básico aplicado neste tipo de algoritmo é descrito a seguir.

Iniciando nas condições iniciais (to), é feita uma estimativa de y(t) no

próximo valor de t (tI) correspondente a tI == to + h, onde h== passo de integração.

O novo valor é usado para estimar o próximo ponto e assim por diante.

A base matemática permeando este método consiste em expandir a função

y( t) em série de Taylor. Estima-se o valor da função a uma pequena distância deum valor conhecido, usando as derivadas dafunção, calculadas no valor conhecido

da mesma.

Assim, dada a função g(t) com derivadas contínuas na região em torno de

t == a, o valor de g( a + h), onde h é pequeno, é dado por:

Repetindo este cálculo seqüencialmente para valores de tu == a + '1'/ •• h, onde

n == 1, 2, 3 ..., resulta em uma série de estimativas de g( t). Se em cada cálculo um

nÚmero infinito de termos é incluído, a solução exata é obtida. Na prática, no

entanto, o nÚmero de termos na série deve ser ±lnito. Assim, algum erro é aceito.

A ordem de magnitude desse erro de truncamento é definida pela potência do

\

396 • iVIodeLagem e Simula[rlo de Processos Industriais e de Sistemas Eletromecânicos

Último termo incluído. Assim, se o Último termo é de ordem h5, o erro resultante

é no máximo de ordem h5 e a ordem de magnitude desse erro é abreviada porO(h5). A inclusão de mais termos reduz o erro de truncamento.

São descritos a seguir dois métodos numéricos de integração que se enquadram

nessa filosofia: Euler e Runge-Kutta (Carnahan et a!., 1969).

14.1.3 lYlétodode Elder

Consiste em utilizar a 1C! derivada Uá conhecida) e usar um passo pequeno

de integração h, truncando termos de ordem 2": 2:

)I(a+h) = )I(a)+h. jJ(a) = )I(a)+h'f[a,y(a)]

onde jl(t) = f(t,)I)

ou, equivalentemente:

Lil = h . f (tn ' )ln )

)1'11+1 = )ln + Ul

Avalia-se a equação com as condições iniciais conhecidas e se realiza a extra­

polação segundo a tangente à curva nesse ponto.

Mostra-se, na figura 1"1.1,a interpretação gráfica do método de Euler.

)'

tn+l = tn + h

)I = Yn quando t = t'll ; procura-se )In+1 quando t = ln+h.

)'11+1

Yn

I

---r------ I-l--

I II II I UII II II I------4---

l~h~1

)'11+1 (real)I, )'*n+1 (calculado)II1

IIII .

)'n p.!:------------l E=)'n+l-y"'n+1I II I

tn [n+ I h

Figura H.L Interpretação grMica elo métoelo de Euler.

_-\..- uéwdo de Euler:

Dados:

• tamanho do passo de integração h

• y(O) = )10

• nÚmero de iterações N·5<t)=f(t,)I)

Para 11, = Oa N - 1 Ülça:

• )1'11 + 1 = )ln + h . jUn,)ln)

·t.II+l=tn+h

• saída = )ln + 1

o erro de truncamento é O(h), sendo proporcional ao tamanho li do passo.

Como o nÚmero de cálculos é inversamente proporcional ao tamanho do pas­so, uma solução precisa requer um grande nÚmero de cômputos.

A relação entre amagnitude do erro e o tamanho do passo no método de Euler

é ilustrada pela figura 1"1.2(Franks, 1972), que mostra a melhora que ocorre quan­

do um passo único de tI a t2 é comparado com dois meio-passos: tI a tm e tm a '2'

)'1

ti

Figllra 1'1.2 1\iétoelo ele Ell/er com passo redllzielo;, metade.

Para o caso com dois meio-passos, começando com a derivada 5'1 em 11,

avança-se meio passo até tm, onde a derivada jlm é reavaliada. Verifica-se que o

resultado final )I';'2m está mais próximo ao valor correto )12 que o valor Y"'2 obtidocom um passo Único.

Para um método de integração de 1il ordem, os erros numéricos são dire­

tamente proporcionais ao tamanho do passo (E =]C li), de forma que a toleráncia

desejada, isto é, o máximo erro aceitável, determina o tamanho do passo a ser usado.

h

~l ~ --1- - )'::~~ll~

2

I T------ - - - - - ---------1- - - --- --I

I

)'

)'n

)'.11+] = )1. + ul +U'). '11 -

2

U2 é avaliado no ponto definido pelo método de I ª ordem (Euler)

y(O) = yo

14.1.4 iVlétoeloeleRunge-Kutta

Parte-se direto para a 2" ordem, porque o método de Runge-Kutta de I"

ordem equivale ao método de Euler. A equação é avaliada nos pontos extremosdo in tervalo h.

MÉTODO DE RUNGE-KUTTA DE 2" ORDEM

Este método imita os termos da série de Taylor sem, no entanto, derivar a

equação original.

398 • iVIodelagem e Simulaçâo de Processos Industriais e de Sistemas Eletromecánicos

figura 14.3 Interpretaç~1o gr<1fica elo método ele Runge-Kutla ele 2" orelem.

Visto que o erro de truncamento é O(h2), um método de Runge-Kutta de

mais alta ordem é necessário se um resultado mais preciso é desejado.

Métodos NuméTicos para Integração de lJquações Difirrenâais OTdinárias • 399

)'

:\IÉTODO DE RUNGE-KUTTA DE 4" ORDEM

)'

)'4

')'4,,

,,,

,,

li +h/2li

Figura HA Inlerpretação gráfica do método ele Rllnge-Klllta ele 4" ordem.

É o mais utilizado. Consiste em avaliar as derivadas no início, meio e fim do

intervalo de integração. O Últimopasso é fazer uma soma ponderada dessas derivadas.

Passos:

• a derivada YI é avaliada em tI e, usando o método de Euler, o valor da função

é calculado em (tI + h/2), resultando Y2;

• a derivada é avaliada em (tI + h/2) resultando Y2;

• iniciando em (YI,tI), a função em tI + h/2 é recalculada usando a derivada Y2

para gerar Y3;

• a derivada )'3 é avaliada em tI + h/2;

• iniciando em (YI,tI), a função é calculada em t2 = tI + h usando a derivada )13

para gerar Y4;

• a derivada Y4 é avaliada em t2; e

• usando as derivadas YI a Y4, o valor da função )I (t2) é calculado por:

o algoritmo de Runge-Kutta de 4ª ordem é:

Dados:

• passo de integração h

• y(O) = Yo

• nÚmero N de iterações

• y(t)=f(t,y)

400 • Nfodelagem e Simulaçâo de ProcessosIndustriais e de Sistemas Eletmmecânicos

Para n = O a N - 1, faça:

'Ul = h . f([lI' )111)

.( h '11'1 )u2 = h . j til + 2" ,)ln + 2

(h 'lI'))u3 = h .f til + "2' )ln + 2U4 = h . .r (til + h,)l1I + us)

tn + 1 = til + h

saída = )111 + 1

Um esquema alternativo para encerrar o algoritmo é especificar t e não N

Visto que o erro de truncamento é O(h4), o erro é muito reduzido quando com­

parado ao método de Euler. Mas são necessários quatro cálculos da função por

passo, resultando que uma redução no tamanho do passo aumen ta muito o nÚmerode cálculos.

Outra forma para reduzir o nÚmero de cálculos é ajustar o tamanho do

passo para manter o erro de truncamento especificado. Em regiões de baixas taxas

de variação de )I(t), o tamanho do passo pode crescer, reduzindo o nÚmero de

cômputos, o contrário ocorreildo para altas taxas de variação de )J( t). Um esquema

muito utilizado para verificação do erro de truncamento, atribuído a Fehlberg, émostrado abaixo.

MÉTODO DE RUNGE-KUTTA-FEI-ILBERG (MÉTODO RKF-45)

Usa uma fórmula para cálculo do erro de truncamento O(h5) baseado em dife­

rentes escolhas para os coeficientes do método de 4;]ordem (Tao, 1988). A diferença

entre os resultados de O(h4) e O(h5) dá uma estimativa do erro de truncamento.

Diversas estratégias são disponíveis para alterações no tamanho do passo

para reduzir o erro de truncamento. Por simplicidade, o algoritmo apresentado

a seguir reduz à metade ou duplica o tamanho do passo de integração.

Algoritmo RKF-45: