Desempenho dos métodos de Runge-Kutta 4-4 e Rosenbrock em função da

posição do cλ crítico em sistemas rígidos

Guilherme G. Loch

1, Joyce S. Bevilacqua,

Instituto de Matemática e Estatística da Universidade de São Paulo (IME-USP)

Departamento de Matemática Aplicada

Rua do Matão, 1010 - Butantã

05508-090, São Paulo, SP

E-mail: guiloch@ime.usp.br, joyce@ime.usp.br

Orlando Rodrigues Jr

Instituto de Pesquisas Energéticas e Nucleares (IPEN-CNEN/SP)

Av. Professor Lineu Prestes, 2242 - Butantã

05508-000, São Paulo, SP

E-mail: rodrijr@ipen.br

Resumo: Modelos n-compartimentais lineares são representados por um sistema de equações

diferenciais ordinárias com coeficientes constantes, que representam as taxas de transferência entre

os compartimentos. Quando uma das taxas de transferência, o c crítico, apresenta ordem de

grandeza muito diferente das demais, o problema torna-se rígido e se faz necessário o uso de

métodos numéricos específicos para esta classe de problemas. Este trabalho apresenta uma

comparação dos desempenhos do método clássico Runge-Kutta 4-4 e do método de Rosenbrock para

problemas compartimentais rígidos, considerando como exemplo a modelagem por compartimentos

da série de decaimento natural do Urânio-235. Por meio da avaliação do desempenho dos métodos

em função da posição do c crítico na cadeia de compartimentos foi possível perceber que quanto

mais próximo do final da cadeia ele estiver, maior será o número de iterações realizadas pelo

Rosenbrock para alcançar a precisão requerida. Porém, conforme previsto pela teoria, mesmo nos

casos mais críticos o método de Rosenbrock apresenta desempenho superior, já que é um método

projetado especificamente para problemas rígidos.

1. Modelos Compartimentais Rígidos

A modelagem matemática é uma ferramenta fundamental para a análise e compreensão de

problemas complexos oriundos de diversas áreas do conhecimento e, em muitos casos, a

representação do problema real é um sistema rígido de equações diferenciais ordinárias (EDO’s). A

rigidez de um sistema aparece quando algumas das componentes de sua solução decaem mais

rapidamente do que outras [2]. Problemas que envolvem sistemas de reações químicas, circuitos

elétricos, vibrações e oscilações configuram alguns exemplos que podem apresentar rigidez [1]. Além

destes, podemos citar os problemas que envolvem séries de decaimento radioativo e são modelados

por uma cadeia fechada de compartimentos.

Este trabalho considera um modelo n-compartimental linear, onde cada compartimento ix ,

ni ,...,1 , está relacionado aos seus adjacentes por meio de taxas de transferência 0i , ni ,...,1 ,

representado por um sistema de equações diferenciais ordinárias com coeficientes constantes:

00

'

xx

txAtx,

(1)

1 Bolsista de Mestrado CNPq

713

sendo [,[ 0 tt , as componentes do vetor x(t) representando cada compartimento, as condições

iniciais dadas por Tnxxxx 0

2

0

1

00 e a matriz A dos coeficientes, que representam as taxas de

transferência entre os compartimentos, contendo elementos somente na diagonal e nas posições

imediatamente abaixo dela. Havendo ramificações no modelo, a matriz A continua sendo triangular

inferior, porém surgem elementos não nulos fora das posições anteriormente citadas.

Um sistema linear com coeficientes constantes é rígido quando todos os autovalores i ,

ni ,...,1 , da matriz dos coeficientes A tiverem parte real 0)Re( i e a sua taxa de rigidez , dada

por:

,)Re(

)Re(

1

1

ini

ini

min

max

(2)

for muito grande [2].

Observe que para o modelo compartimental considerado, as taxas de transferência são

exatamente os autovalores da matriz dos coeficientes A. Dessa forma, quando uma das taxas de

transferência, a qual será chamada de c crítico, apresentar ordem de grandeza muito diferente das

demais, o problema é considerado rígido.

A comparação do desempenho dos métodos numéricos de Runge-Kutta 4-4 e Rosenbrock para

problemas rígidos aparece como escopo principal deste trabalho. Apresentamos um exemplo de

aplicação à área nuclear, considerando a série de decaimento natural do Urânio-235, que é modelada

por um sistema n-compartimenal, onde cada compartimento representa a concentração de um dado

nuclídeo e as taxas de transferência, que têm ordens de grandeza entre 10-12

e 1001

, são as constantes

de decaimento. Também realizamos uma análise das soluções numéricas em função da posição do c

crítico na cadeia de compartimentos.

2. Solução Analítica

Supondo que no instante inicial 00 t existe apenas uma quantidade 0

11 )0( xx e não são

consideradas recirculação, perdas e acréscimos de quantidades durante a dinâmica do sistema (1), é

possível escrever a solução do problema analiticamente da seguinte forma:

.1,1 ,1

10

1 )()(

n

n

i

n

ijj

ijj

t

i

n

nie

xtx

(3)

onde os 0i , ni ,...,1 , são exatamente as taxas de transferência entre os compartimentos e

correspondem aos autovalores da matriz A.

Para modelos onde existem ramificações, ainda é possível utilizar a solução analítica dada pela

equação (3) realizando o desdobramento e considerando cada ramo da cadeia separadamente como

um modelo compartimental linear.

Mesmo que a existência e a unicidade da solução exata sejam garantidas, quando pelo menos

duas taxas de transferência são iguais ou muito próximas, surgem indeterminações no cálculo da

solução analítica dada pela equação (3), de modo que alguns denominadores tornam-se nulos ou

aproximam-se de zero. Neste caso, se faz necessário o uso de estratégias para que estas

indeterminações sejam removidas, tornando-se imperativo o uso de técnicas semi-analíticas ou

numéricas.

714

3. Soluções Numéricas

Independente da existência de soluções analíticas, soluções numéricas permitem adaptações na

modelagem ou em ajustes de dados com maior facilidade. No entanto, sistemas rígidos usualmente

geram instabilidades numéricas, exigindo cautela na escolha do método numérico. Métodos implícitos

como os Rosenbrock e os BDF (Backward Differentiation Formulae) apresentam desempenho muito

superior do que os métodos explícitos [1]. Em geral, para manter a estabilidade dos algoritmos

explícitos é necessário um passo de integração muito pequeno, gerando, consequentemente,

simulações muito lentas. Porém, existem casos em que os métodos explícitos não apresentam

dificuldades na resolução desta classe de problemas.

A proposta deste estudo é comparar o desempenho do método clássico de Runge-Kutta 4-4 com o

do Rosenbrock, que é mais adequado para sistemas rígidos.

Os métodos numéricos foram implementados em Fortran 90 e os algoritmos executados em uma

máquina com processador Intel® Core™ i5 e Windows® 7 64-bit.

3.1 Método de Runge-Kutta 4-4 e a estimativa do passo de integração

O método de Runge-Kutta 4-4 é um método explícito de passo único, de 4ª ordem e 4-estágios,

obtido pelo truncamento da série de Taylor. O método foi implementado de forma que sua

estabilidade e convergência fossem garantidas para uma precisão pré-fixada e o tempo de integração

em aberto.

O tamanho do passo de integração h para que o erro absoluto cometido pelo método de Runge-

Kutta 4-4 no tempo final ft seja inferior a uma precisão 0 pré-fixada foi estimado através da

delimitação do seu erro de discretização global. Sabendo que não existe recirculação, perda e

acréscimo de quantidades durante a dinâmica do sistema, realizando algumas maximizações e

substituições adequadas ao problema e considerando valores para h dentro da região de estabilidade

do método, o tamanho do passo de integração deve satisfazer: 41

0

11

5

120

xAth

f

. (4)

3.2 Método de Rosenbrock

Os métodos de Rosenbrock, também conhecidos como métodos de Runge-Kutta linearmente

implícitos, surgiram como uma alternativa para resolver equações implícitas, que geralmente eram

resolvidas por processos iterativos [1]. O método Rosenbrock implementado para este trabalho é de 3ª

ordem e 4-estágios, utiliza uma mudança de variáveis para reduzir os gastos computacionais e uma

fórmula embutida de 2ª ordem para controle do passo de integração.

4. Simulações e Resultados

4.1 Decaimento Natural do Urânio-235

Considerando os sete primeiros elementos da série de decaimento natural do Urânio-235,

conforme apresentado na Figura 1, é possível visualizar um exemplo de problema rígido, onde a taxa

de rigidez é da ordem de 1013

.

Enquanto a constante de decaimento do 235

U é da ordem de 10-12

, a do 223

Fr é da ordem de 1001

.

Essa grande diferença entre as ordens de grandeza se traduz em uma matriz dos coeficientes A quase

singular, cujo valor do determinante é aproximadamente 2510140,9 .

715

A solução analítica e as aproximações numéricas para o problema proposto foram construídas

considerando que no instante inicial existem 100 unidades de 235

U e que o inventário dos demais

elementos da série é nulo.

Figura 1. Série de decaimento natural do Urânio-235 com taxas de decaimento em dia-1

[3].

A Tabela 1 apresenta a solução analítica para o problema proposto após terem sido percorridos

1006

dias a partir do instante inicial e os erros absolutos cometidos pelos métodos, dada uma precisão 0410 . O tamanho do passo de integração utilizado para o método de Runge-Kutta 4-4 foi

05105 h dia, estimado pela equação (4). Este mesmo valor foi utilizado como passo inicial para o

método de Rosenbrock e os passos mínimo e máximo utilizados, calculados por meio da estratégia de

controle do passo, foram, respectivamente, 05105 dia e aproximadamente 3.781 dias.

Elemento Solução Analítica Erro Absoluto do

Runge-Kutta 4-4

Erro Absoluto do

Rosenbrock

235U 9,99997300003644E+01 1,4E-05 9,9E-14

231Th 4,14196627908533E-10 5,9E-17 2,0E-26

231Pa 2,62322611868427E-04 6,2E-11 5,1E-17

227Ac 1,72402403970940E-07 4,1E-14 3,4E-20

227Th 4,05986467670791E-10 5,5E-12 5,5E-12

223Fr 1,40090679162464E-13 1,4E-13 1,4E-13

223Ra 3,53799074572466E-10 1,1E-10 1,1E-10

Tabela 1. Solução analítica para o problema e erros absolutos após 1006

dias.

Para que a solução aproximada atingisse a precisão pré-fixada, foram necessárias 10102 e 415

iterações para os métodos de Runge-Kutta 4-4 e Rosenbrock, respectivamente. O método de

Rosenbrock mostrou um desempenho superior ao método de Runge-Kutta 4-4 em função da

velocidade com a qual as soluções são calculadas, uma vez que apresentou soluções convergentes e

consistentes com um número de iterações cerca de 1008

vezes inferior.

716

4.2 Desempenho dos métodos numéricos com relação à posição do cλ crítico

O impacto da posição do c crítico foi estudado para uma cadeia linear composta por 7

compartimentos, na qual todas as taxas de transferência são iguais a 1, com exceção do c , que varia

de 10-06

a 10-12

. Para a posição do c crítico na cadeia foram consideradas três situações: na posição

inicial ( 1 ), na posição intermediária ( 4 ) e na posição final (7 ). As aproximações numéricas para o

problema proposto foram construídas para o tempo final 0610ft , considerando Tx 001000 .

Observe que a taxa de rigidez do sistema varia entre 1006

e 1012

conforme a variação do c , o que

comprova a rigidez do problema.

Em todas as situações consideradas, sabendo que 321

5 A e fixando uma precisão 0410 , o

passo de integração que garante a convergência e a estabilidade do método de Runge-Kutta 4-4,

estimado por meio da equação (4), é 0310h , o que corresponde a 1009

iterações do método.

A análise do desempenho do método de Rosenbrock de acordo com a posição Pi do c foi

realizada usando 0310h como passo de integração inicial. Os tamanhos mínimo (hmin) e máximo

(hmax) para os passos de integração utilizados e calculados pela estratégia de controle do passo, bem

como o número de iterações (n) necessárias para que o método encontrasse soluções convergentes e

consistentes para o problema são apresentados na Tabela 2.

Pi 0610i 0910i 1210i

hmin hmax n hmin hmax n hmin hmax n

1 5,7E-04 3,6E+01 35.084 1,0E-03 3,6E+04 770 1,0E-03 5,0E+04 98

4 4,8E-06 3,6E+01 360.062 4,8E-06 3,6E+04 332.401 4,8E-06 5,0E+04 332.393

7 4,8E-06 3,6E+01 447.734 4,8E-06 3,6E+04 420.074 4,8E-06 5,0E+04 420.066

Tabela 2. Desempenho do Rosenbrock em função da variação da posição Pi e do valor do c crítico.

A Figura 2 mostra a variação do tamanho do passo de integração e o número de iterações

realizadas pelo método de Rosenbrock de acordo com a posição do c crítico no caso em que ele vale

10-06

.

É possível perceber a partir das situações consideradas que o método de Rosenbrock é bastante

sensível a mudanças na posição do c crítico. Quando o elemento crítico está na posição P1, o

número de iterações do Rosenbrock é significativamente inferior ao número de iterações nos casos em

que ele está nas posições P4 e P7.

Quanto mais próximo do final da cadeia estiver posicionado o c crítico, maior será o número de

iterações necessárias para o método de Rosenbrock obter soluções convergentes e consistentes, já que

para manter a estabilidade do método, o tamanho do passo de integração é reduzido. Porém, mesmo

na pior das situações, em que o c crítico vale 10

-06 e está na posição P7, o método de Rosenbrock

continua em vantagem com relação ao método de Runge-Kutta 4-4 no que diz respeito ao número de

passos realizados, já que para gerarem soluções com mesma precisão, o primeiro requer 447.734

iterações enquanto o segundo necessita de 1009

iterações.

Novamente, o método Rosenbrock teve um desempenho superior ao Runge-Kutta 4-4, atingindo

a mesma precisão com um número significativamente inferior de passos, devido ao fato de ser um

método especificamente projetado para equações rígidas.

717

Figura 2. Tamanho e número de passos do Rosenbrock em função da posição Pi do c crítico.

5. Conclusão

Foram comparados os desempenhos dos métodos de Runge-Kutta 4-4 com passo fixo e

Rosenbrock de 3ª ordem com passo adaptativo para a resolução de problemas rígidos.

A modelagem da série de decaimento natural do Urânio-235 foi considerada como exemplo de

aplicação para comparação dos métodos. O método de Rosenbrock mostrou desempenho superior,

apresentando soluções convergentes e consistentes com um número de iterações cerca de 1008

vezes

menor.

O desempenho do método de Rosenbrock fica condicionado à posição, na cadeia, do

compartimento que tem taxa de transferência com ordem de grandeza muito diferente das demais, de

modo que quanto mais próximo do final da cadeia ele estiver, maior será o número de iterações

realizadas pelo Rosenbrock para alcançar a precisão requerida. Porém, mesmo nos casos mais críticos,

o Rosenbrock apresenta soluções com mesma precisão do que o método de Runge-Kutta 4-4

realizando um número de iterações cerda de 1004

vezes inferior.

Para trabalhos futuros, a proposta é analisar os casos em que ocorrem retiradas e/ou acréscimos

de quantidades nos compartimentos durante a dinâmica do sistema, o que permite considerar

aplicações na produção de radiofármacos.

Referências

[1] E. Hairer, G. Wanner, “Solving Ordinary Differential Equations II (Stiff and Differential-

Algebraic Problems)”, Springer, Berlin, 1996.

[2] J. D. Lambert, “Numerical methods for ordinary differential systems: the initial value problem”,

John Wiley & Sons, New York, 1991.

[3] G. F. Thomas, D. H. Barber, Stiffness in radioactive decay chains, Ann. Nucl. Energy, vol 21 (5),

pp.309–320, (1994).

718