Integração associativa em métodos Runge-Kutta

Universidade do Estado do Rio de JaneiroCentro de Tecnologia e Ciências

Instituto de Física Armando Dias Tavares

Alex Fraga Rocha

Integração associativa em métodos Runge-Kutta

Rio de Janeiro2016

Alex Fraga Rocha


Dissertação apresentada, como requisito par-cial para obtenção do título de Mestre, aoPrograma de Pós-Graduação em Física, daUniversidade do Estado do Rio de Janeiro.

Orientador: Prof. Dr. Luís Antônio Campinho Pereira da MotaCoorientador: Prof. Dr. Luiz Guilherme Silva Duarte

Rio de Janeiro2016

CATALOGAÇÃO NA FONTE UERJ/ REDE SIRIUS / BIBLIOTECA CTC/D

Autorizo, apenas para fins acadêmicos e científicos, a reprodução total ou parcial desta dissertação, desde que citada a fonte.

R672 Rocha, Alex Fraga. Integração associativa em métodos de Runge-Kutta / Alex Fraga Rocha. - 2016. 77 f.: il. Orientador: Luís Antônio Campinho Pereira da Mota. Coorientador : Luiz Guilherme Silva Duarte. Dissertação (mestrado) - Universidade do Estado do Rio de Janeiro, Instituto de Física Armando Dias Tavares.

1. Integração numérica -Teses. 2. Runge-Kutta, Fórmula de - Teses. 3. Equações diferencias - Soluções numéricas - Teses. I. Mota, Luís Antônio Campinho Pereira da. II. Duarte, Luiz Guilherme Silva. III. Universidade do Estado do Rio de Janeiro. Instituto de Física Armando Dias Tavares. IV. Título. CDU 517.91

Alex Fraga Rocha


Dissertação apresentada, como requisito par-cial para obtenção do título de Mestre, aoPrograma de Pós-Graduação em Física, daUniversidade do Estado do Rio de Janeiro.

Aprovada em 25 de fevereiro de 2016.Banca Examinadora:

Prof. Dr. Luís Antônio Campinho Pereira da Mota (Orientador)Instituto de Física Armando Dias Tavares – UERJ

Prof. Dr. Luiz Guilherme Silva Duarte (Coorientador)Instituto de Física Armando Dias Tavares – UERJ

Prof. Dr. Jair KoillerFundação Getúlio Vargas

Prof. Dr. Luís Fernando de OliveiraInstituto de Física Armando Dias Tavares – UERJ

Prof. Dr. Vitor Emanuel Rodino LemesInstituto de Física Armando Dias Tavares – UERJ

Prof. Dr. Jayr Avellar Costa FilhoFundação de Apoio à Escola Técnica

Rio de Janeiro2016

DEDICATÓRIA

À Ana Paula e João Victor

AGRADECIMENTOS

Primeiramente, a Deus que é a minha fonte de vida e o motivo da minha existência.Com carinho, à minha esposa Ana Paula e ao meu filho João Victor.Ao meu Orientador, Luís Antônio Campinho Pereira da Mota, por acreditar em

mim.Ao meu Coorientador, Luiz Guilherme Silva Duarte, pelo incentivo e motivação.Ao meu amigo Professor Luís Fernando de Oliveira, que foi quem me acolheu,

quando mais precisei, ainda na graduação, para concluir meu curso de física.Ao professor, Jayr Avellar Costa Filho, por ter me indicado o livro para o pontapé

inicial deste trabalho.A todos que contribuiram, de alguma forma, para essa conquista.

RESUMO

ROCHA, A. F. Integração associativa em métodos Runge-Kutta. 2016. 77 f. Dissertação(Mestrado em Física) – Instituto de Física Armando Dias Tavares, Universidade doEstado do Rio de Janeiro, Rio de Janeiro, 2016.

O presente texto apresenta um novo método de integração numérica que traz umamelhoria de desempenho. Esta melhoria será mais bem definida ao longo do texto. Basica-mente, o método trabalha com dois integradores de Runge-Kutta e consegue aperfeiçoar aprecisão obtida e/ou o tempo de processamento necessário. O presente trabalho apresentatambém uma implementação em Maple destas ideias para que possamos ver as ideias emação e quantificar esta melhoria.

Palavras-chave: Integração Numérica. Sistemas de Equações Diferenciais Ordinárias.Método Numérico de Runge-Kutta.

ABSTRACT

ROCHA, A. F. Associative integration Runge-Kutta methods. 2016. 77 f. Dissertação(Mestrado em Física) – Instituto de Física Armando Dias Tavares, Universidade doEstado do Rio de Janeiro, Rio de Janeiro, 2016.

This paper presents a new numerical integration method that brings an improve-ment to the performance in solving the differential equations system. This improvementwill be better defined as the reader progress on the text. Basicaly, the method workswith two Runge-kutta integration methods and improves either the precision obtainedor the computational time spent. This present work also brings the implementation ofthese ideas in Maple in order for the reader to see them in action and to quantify theimprovement above mentioned.

Keywords: Numerical Integration. Ordinary Differential Equations Systems. NumericRunge-Kutta Method.

LISTA DE ILUSTRAÇÕES

Figura 1 - Gráfico de y=y(x). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17Figura 2 - Gráfico de y=f(x). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17Figura 3 - Método de Euler . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22Figura 4 - Método de Runge-Kutta . . . . . . . . . . . . . . . . . . . . . . . . . . 24Tabela 1 - Desempenho do algoritmo AI no sistema de Lorenz. . . . . . . . . . . . 53Tabela 2 - Desempenho do algoritmo AI para sistema Hénon-Heiles. . . . . . . . . 57

LISTA DE ABREVIATURAS E SIGLAS

IFADT Instituto de Física Armando Dias TavaresPPG Programa de Pós-graduaçãoUERJ Universidade do Estado do Rio de Janeiro

SUMÁRIO

1 INTRODUÇÃO . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111.1 Este Trabalho . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 142 SISTEMAS DE EQUAÇÕES DIFERENCIAIS E A INTEGRA-

ÇÃO NUMÉRICA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 152.1 Introdução . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 152.2 Equações Diferenciais . . . . . . . . . . . . . . . . . . . . . . . . . . . 162.3 A importância das equações diferenciais . . . . . . . . . . . . . . . . 182.4 Classificação das Equações Diferenciais . . . . . . . . . . . . . . . . 182.5 Sistema de Equações Diferenciais . . . . . . . . . . . . . . . . . . . . 182.6 Ordem de uma Equação Diferencial . . . . . . . . . . . . . . . . . . . 192.7 Equações Diferenciais Lineares e Não Lineares . . . . . . . . . . . . 192.8 Métodos Numéricos de Integração . . . . . . . . . . . . . . . . . . . 202.9 Método de Euler ou Método da Reta Tangente . . . . . . . . . . . 202.10 Método de Runge-Kutta de quarta ordem em quatro estágios . . 212.11 Método de Passos Múltiplos . . . . . . . . . . . . . . . . . . . . . . . 232.11.1 Método de Adams-Bashforth . . . . . . . . . . . . . . . . . . . . . . . . . 232.11.2 Método de Adam-Moulton . . . . . . . . . . . . . . . . . . . . . . . . . . . 253 CONSTRUÇÃO DE UM ALGORITMO DE INTEGRAÇÃO . . 273.1 O Método de Taylor . . . . . . . . . . . . . . . . . . . . . . . . . . . . 273.1.1 A Base Matemática . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 273.1.2 Aplicando o Método de Taylor ao Sistema de Lorenz . . . . . . . . . . . . 293.2 Um Algoritmo de Integração Associativa (AI) . . . . . . . . . . . . 303.2.1 O Resultado Principal . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 313.2.2 A Base Estrutural de um AI . . . . . . . . . . . . . . . . . . . . . . . . . 333.2.3 O Algoritmo Propriamente Dito . . . . . . . . . . . . . . . . . . . . . . . 374 IMPLEMENTAÇÃO DE UM AI . . . . . . . . . . . . . . . . . . . . 404.1 Uma Implementação em Maple: O pacote AssISt . . . . . . . . . . . 414.1.1 Um resumo do pacote . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 414.1.2 Os comandos do Pacote . . . . . . . . . . . . . . . . . . . . . . . . . . . . 424.1.2.1 Comando: PTaylor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 424.1.2.2 Comando: Prk4maple . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 434.1.2.3 Comando: Prk78maple . . . . . . . . . . . . . . . . . . . . . . . . . . . . 434.1.2.4 Comando: Pcost . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 434.1.2.5 Comando: Dif . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 444.1.2.6 Comando: Pai . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 454.1.2.7 Comando: Pan . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46

4.1.3 Exemplificando o Uso do Pacote . . . . . . . . . . . . . . . . . . . . . . . 474.2 Uma Análise de Desempenho . . . . . . . . . . . . . . . . . . . . . . . 504.2.1 Sistema de Lorenz . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 504.2.1.1 Lorenz – Análise 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 504.2.1.2 Lorenz – Análise 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 514.2.1.3 Lorenz – Análise 3: Uma Análise mais Completa . . . . . . . . . . . . . . 524.2.2 Sistema de Hénon-Heiles . . . . . . . . . . . . . . . . . . . . . . . . . . . . 544.2.2.1 Henon-Heiles – Análise 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . 554.2.2.2 Henon-Heiles – Análise 2: Uma Análise mais Completa . . . . . . . . . . . 564.2.3 Algumas Observações . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58

CONCLUSÃO . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60REFERÊNCIAS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62APÊNDICE A – Pacote AssISt . . . . . . . . . . . . . . . . . . . . . . 63

11

1 INTRODUÇÃO

Desde os tempos de Newton e Liebnitz, as equações diferenciais tem ocupadolugar de destaque no entendimento e descrição de sistemas Físicos e, na verdade, emtodo o leque das ciências naturais, Economia, etc. Em suma, têm desempenhado papelpreponderante na descrição científica como regra geral.

Nosso grupo de pesquisa tem se dedicado a esta área de conhecimento, tanto emsua vertente algébrica quanto numérica. O trabalho aqui apresentado se encaixa nestetodo tratando de sistemas de equações diferencias ordinárias. Mais especificamente, nodesenvolvimento de uma nova abordagem à integração numérica deste tipo de sistema.

Em grande parte, nosso trabalho baseia-se em resultados apresentados para o casoparticular de análise de previsibilidade em séries temporais (BARBOSA et al., 2006)e que foi implementada em (CARLI; DUARTE; MOTA, 2014). Aqui nós adaptamos oconjunto de ideias e formulamos um algoritmo que se presta, a princípio muito bem, para aintegração de sistemas de equações diferenciais ordinárias. Brevemente, estas adaptaçõestem que ser feitas devido a peculiaridades inerentes ao caso da análise de séries temporais.Isto ficará claro ao tratarmos do método mais detalhadamente no corpo da Dissertação .

É claro que o tema da integração numérica de sistemas de equações diferenciasé por demais vasto para ser coberto em todos os seus aspectos aqui nesta Dissertação.Deixaremos tal pretensão de fora de nosso texto. Apenas cobriremos alguns casos paraexemplificar características e propriedades de certos métodos. Isto será feito até pararessaltar que, em certos termos, nosso método tem muito o que melhorar: é uma primeiraaproximação do uso de uma ideia que consideramos poderosa. Um pouco mais sobre istonos próximos parágrafos, ainda na Introdução.

Basicamente, em linhas incrivelmente gerais, a questão da eficiência (tempo com-putacional gasto versus precisão) na integração numérica pode ser dividida (muitos argu-mentam) no equilíbrio entre os erros de truncamento e arredondamento:

• erro de truncamento - essencialmente a ideia é que cada método de inte-gração equivale em precisão a uma ordem da série de Taylor correspon-dente, ou seja, se usamos até o terceiro termo da série de taylor (infinita),teremos erros na resposta da ordem do quarto termo, etc. Daí a origemdo nome.• erro de arredondamento - a questão aqui é a limitação oriunda do fato

que os computadores oferecem um número finito de casas decimais pararepresentar os números (de uma certa maneira, voltando à ideia grega dainexistência dos números reais, apenas números racionais são possíveis nocomputador). Isto também limita a precisão possível.

12

Um outro aspecto importante nos integradores numéricos, uma classificação quese aplica a estes, é a divisão entre métodos com passo fixo, único, e aqueles com passovariável. Retomando a questão da equivalência entre um determinado integrador e umasérie de Taylor truncada em uma determinada ordem, se usamos uma equivalência atéa ordem n, com passo h, a série de Taylor equivalente vai até a ordem hn, e o erro, aimprecisão do cálculo, tem teoricamente erro da ordem hn+1. Vejam como o “tamanho”dopasso h é importante: Se tivermos um método numérico que equivalha a uma série deTaylor truncada em valores “altos” de n, o passo não precisa ser muito pequeno (postoque hn+1 será pequeno tendo em visto n ser grande, analogamente, se n for “pequeno”,h tem que ser menor para garantir a mesma precisão (tendo em visto que o erro serádado por uma potência menor). Isto é a essência do impacto que o tamanho do passo (h)tem na integração. E onde fica a questão do passo variável? Bem, grosso modo, certasfunções (a solução do problema que queremos resolver) comportam-se diferentementeem diferentes regiões de integração. Logo, os métodos que usam passo variável tentamavaliar a liberdade que têm de encontrar a mesma precisão em toda a região de integração,mudando o passo e, em caso do aumento dele, ganhando em tempo de processamento,ainda mantendo a precisão desejada. Lembrem-se que o termo desprezado (termo deordem n + 1 da série de Taylor correspondente) tem ainda a contribuição da derivadada função (a variação destes valores incorpora esta diferença entre as diversas regiões deintegração em termos de variabilidade para um mesmo passo h).

O acima exposto fala um pouco da questão do erro de truncamento e quanto aoerro de arredondamento (CELLIER; KOFMAN, 2006)? Bem, um aspecto já foi abordado.Se usarmos um passo “pequeno”, no sentido que foi mencionado acima, teremos queproceder a muitas etapas de integração para cobrir o intervalo desejado. Logo, muitasoperações, significando que os erros de arredondamento ocorrerão mais vezes, tornandomais “importante” este tipo de erro, ou seja, mais importante controlá-lo.

Um outro lado da questão do erro de arredondamento é relacionado à questão daprecisão simples (ou dupla) que esteja disponível na máquina de trabalho. Vejamos:

Supondo que estejamos trabalhando com uma máquina de precisão simples (e com32 bits) e desejemos utilizar um passo h = 0.001: Os computadores rodando com estaprecisão apresentam números com mantissa de 24 dígitos, deixando 8 para o expoente.Logo, o erro de arredondamento pode ser avaliado como sendo:

δarredondamento = 2−24 ∼ 0.0000000596 ∼ 10−6

Vamos analisar o seguinte: Se trabalharmos com um método com equivalência auma série de Taylor até ordem digamos, 4, e com o passo exemplificado acima (0.001), aexpansão terá termos com a seguinte ordem:

• primeiro termo → ordem 1

13

• segundo termo → ordem 0.001• terceiro termo → ordem 0.000001• quarto termo → ordem 0.000000001• quinto termo → ordem 0.000000000001• etc ...

Logo, se estivermos trabalhando com estes parâmetros, com esta máquina e preci-são, não faz muito sentido utilizarmos um método de integração de ordem 4 (por exemplo,rk4)...o erro de arredondamento seria um fator muito mais limitante à precisão do que ofato de usarmos um método de ordem menor. Apesar desta análise simples, muitas apli-cações comerciais usam este tipo de combinação, desnecessária, sem alertar ao usuáriosobre o problema.

Um outro aspecto que queremos mencionar sobre a questão da precisão (duplaou simples) é que, o erro de arredondamento, para uma precisão da mesma ordem dodiscutido acima, fica “resolvido” para uma gama de problemas. Não é mais uma grandequestão para este tipo de problema. Se utilizarmos, por exemplo, um integrador de ordem4 numa máquina rodando com precisão dupla (e com 32 bits). Nesta situação, em geral,a mantissa tem 52 bits e 12 bits são utilizados para o expoente. Se quisermos manter omesmo passo do ilustrado acima, teremos:

δarredondamento = 2−52 ∼ 2.22× 10−16 ∼ 10−13

Logo, se estivermos em problemas que se prestem a análise numérica, com passoda ordem do até aqui exemplificado, se utilizarmos processadores com 32 bits e precisãodupla, o problema de arredondamento está “controlado” (vide análise da progressão doserros com a ordem da expansão acima). É claro que o uso de precisão dupla é mais custosocomputacionalmente, inclusive existem técnicas para utilizar apenas alguns números comprecisão dupla para minimizar o custo computacional. Uma análise mais profunda podeser feita, mas acredito que o ponto principal está explicado.

Todos estes aspectos são relevantes na escolha de qual método (e com que parâ-metros, se de passo variável ou não, etc) numérico deve-se utilizar para certo tipo deproblema, que tipo de sistema de equações diferenciais ordinárias queremos integrar. Naverdade, no que toca esta Dissertação, todas estas análises são importantes para a criaçãode novos métodos e técnicas para integração numérica. Como esperamos deixar claro aolongo desta Dissertação, nosso método proposto é uma tentativa inicial de produção deuma nova abordagem. Logo, não refinamos o processo ao ponto de acolher a ideia de umpasso variável, por exemplo, este estudo ficará como um desenvolvimento posterior possí-vel. Mas todos estes aspectos, acima mencionados, são importantes para o entendimentodo porque nosso método funciona.

No capítulo seguinte, falaremos brevemente sobre alguns métodos de integração

14

numérica (bem conhecidos) onde estas questões são tratadas de maneiras variadas.Abaixo, detalhamos a estrutura deste trabalho.

1.1 Este Trabalho

Basicamente, esta dissertação divide-se da seguinte forma:

• No capítulo 2, apresentamos uma visão bem geral de integração numéricde sistemas de equações diferenciais ordinárias com o intuito de exem-plificar que o tópico tem sido exaustivamente explorado, indicando quehá real alpicação para qualquer boa ideia nesta área. Começamos comdefinições realmente básicas, o que são equações diferenciais, ordem deuma equação diferencial, sistemas delas, etc. Além disto, falamos breve-mente sobre a questão dos passo múltiplos 2.11. Damos então exemplosde algoritmos e métodos consagrados pelo seu uso e bem sucedidos emtermos de eficiência.

1. método de Newton

2. método de Euler 2.9

3. método de Runge-Kutta 2.10

4. método de Adams - Bashforth 2.11.1

5. método de Adam - Moulton 2.11.2

Basicamente, este capítulo é apenas uma rápida introdução, com algumasreferências básicas do assunto para o realmente iniciante. Se o leitor seconsiderar suficientemente “expert”, sugiro que pule este capítulo.• No capítulo 3 apresentaremos (em detalhe) uma extensão ao método de

nossas ideias apresentadas em (BARBOSA et al., 2006) para a criação deum integrador numérico de propósito geral.• No capítulo 4 introduzimos uma implementação computacional destas

ideias. Incluiremos também uma análise de desempenho de nossos al-goritmos e métodos, baseados em algumas aplicações para testarmos odesempenho do protótipo de integrador numérico.• Por fim, apresentaremos as conclusões e algumas linhas que estamos per-

seguindo para o futuro.

15

2 SISTEMAS DE EQUAÇÕES DIFERENCIAIS E A INTEGRAÇÃONUMÉRICA

2.1 Introdução

Antífon, que viveu na mesma época de Sócrates, foi um dos primeiros a manifestaro Cálculo Integral. O mesmo defendia que após sucessivas duplicações do número de ladosde um polígono regular inscrito num circulo, a diferença entre a área do circulo e a dospolígonos seria, ao final, nula. Quase analogamente, pensava Eudóxio, que tinha comobase a proposição: “Se de uma grandeza subtrai-se uma parte não menor que sua metade,do restante outra parte não menor que sua metade, e assim por diante, numa determinadaetapa do processo chega-se a uma grandeza menor que qualquer outra da mesma espéciefixada a priori”. Olhando bem, logo nos vem à cabeça a noção de soma de partes numasequência infinita, conhecida, na época, como método de exaustão, o qual Arquimedesdominava.

Embora se tivesse o conhecimento do método de exaustão, os grandes problemascientíficos do século xvii dependiam de um instrumento matemático mais rápido e maisamplo que o método de exaustão para a sua solução. Dentre os grandes problemas aserem resolvidos destacam-se quatro: O primeiro baseava-se em determinar a velocidadee a aceleração de um móvel; O segundo se referia à determinação de tangentes a curvas;O terceiro dizia a respeito do cálculo de máximos e mínimos e o quarto, por sua vez,referia-se à determinação de comprimentos, áreas, volumes e centros de gravidade, tarefasessas cujo método de exaustão não apresentava ser uma ferramenta adequada.

Diversos cientistas, muitos deles matemáticos, se depararam com esses problemas.Dentre eles, dois se destacaram, cada qual a sua maneira e forma: Newton e Leibniz.

Isaac Newton (1643 – 1727) nascido na Aldeia de Woolsthorpe, Inglaterra, graduou-se em Ciências no Trinity College em Cambridge e no período em que a peste bubônicaestava afligindo Londres o mesmo voltou a sua terra natal, onde ficou durante dois anos,período no qual desenvolveu técnicas científicas da Teoria da Gravitação Universal e doMétodo dos Fluxos, que mais tarde fora chamado de Cálculo Diferencial. Uma das pri-meiras noções de Cálculo expressas por Newton foi quando usou a expansão generalizadade (x + a)p e mostrou que a área sob a curva z = axp (p ∈ Q) é y = pax(p−1) (que hojesabemos ser a derivada de z). E, de forma “inversa”, mostrou que a área sob a curvay = pax(p−1) é z = axp.

Gottfried Wilhelm Leibniz (1646 – 1716) nascido em Leipzig, na Alemanha, e que seconfigurava num autodidata aprendendo Latim e Grego sozinho quando ainda era menino,graduou-se em Direito em Leipzig e em 1667 obteve o grau de Doutor em Filosofia naUniversidade de Altdorf com a tese Ars Combinatória (A arte das Combinações). Leibniz,

16

após viajar para Londres em 1673 e conhecer a obra de Isaac Barrow, ex-professor de IsaacNewton, teve o conhecimento da versão do Cálculo de Newton. Quando Leibniz esteveem Londres pela segunda vez, isso já em 1676, o mesmo já tinha desenvolvido as maisimportantes características e notações do seu Cálculo.

Para Newton, a taxa de variação representava a espinha dorsal do Cálculo, Já paraLeibniz, era a Diferencial. Como simbologia de Diferenciais de x e y, por exemplo, Leibnizusou dx e dy, respectivamente, assim como determinou as seguintes regras:

i) da = 0 se a for constanteii) d(u+ v) = du+dvii) d(uv) = udv+vdu

O Cálculo Integral de Leibniz tem uma relação com seu Cálculo diferencial atravésde somatórios de áreas infinitesimais, onde cada elemento de área está sob a curva y = y(x)de forma que seu valor é dado por ydx, como mostra a figura 1.

Para representar um somatório de todas estas áreas Leibniz inventou o símbolo∫.

Com isso, a área total sob a curva y = y(x) é dada por:∫ydx. De acordo com o gráfico1,

observa-se que a diferença entre as áreas S(OHG)-S(OEF) = ydx é a diferencial da área,e matematicamente fica: d

∫ydx = ydx.

É importante notar que a integral∫ydx recebeu o nome de integral indefinida pois

não possui limites de integração. Já a integral definida pode ser determinada da seguinteforma: Seja f(x) uma função contínua e integrável no intervalo [a,b] e F(x) sua primitiva.O número I =

∫ ba f(x)dx é a Integral definida de f(x) em [a,b], e seu valor é expresso

através da primitiva de f(x) como: I=F(b)-F(a). Algebricamente, esse valor representa aárea sob a curva do gráfico da figura 2 entre a e b.

2.2 Equações Diferenciais

Equações diferenciais são equações compostas de derivadas. Como exemplo deequação diferencial, temos a Eq.(1) que representa o movimento de um corpo em quedana atmosfera terrestre,

dv

dt= g − γt (1)

17

Figura 1 - Gráfico de y=y(x).

Legenda: Elemento de área ydx.Fonte: O autor, 2016.

Figura 2 - Gráfico de y=f(x).

Legenda: Área sob o gráfico no intervalo [a,b].Fonte: O autor, 2016.

18

2.3 A importância das equações diferenciais

Um problema físico, como o crescimento do volume de uma gota de chuva em quedalivre na atmosfera, ou a taxa de variação do fluxo sanguíneo na veia de um corpo humano,bem como o crescimento populacional de uma colônia de bactérias, são fenômenos quetêm como modelo matemático equações diferenciais.(BOYCE; DIPRIMA, 2015).

A solução dessas equações permite ao homem prever o comportamento desses fenô-menos, de forma que a sua ação sobre eles pode ser bem precisa. Muitos são os problemascientíficos que necessitam de um modelo matemático que os permita conhecê-los, mas aferramentas de cálculo nem sempre são possíveis. Entretanto, o crescimento de progra-mas com a finalidade de solucioná-los é um passo importante para o desenvolvimento dahumanidade, tais programas são chamados de métodos numéricos de solução de equaçõesdiferenciais (PRESS et al., 1992).

2.4 Classificação das Equações Diferenciais

As Equações Diferenciais são classificadas em equações diferenciais ordinárias,quando dependem de uma única variável independente, e em equações diferenciaisparciais, quando dependem de várias variáveis independentes. A Eq.(2), que mostracomo a carga elétrica varia num circuito simples formado por um capacitor, uma bateriae um resistor, é um exemplo de equação diferencial ordinária,

R dQdt

+ QC

=V. (2)

Já a Eq.(3), que representa a equação de calor, é um exemplo de equação diferencialparcial,

α2 ∂2u(x, t)∂x2 = ∂u(x, t)

∂t(3)

2.5 Sistema de Equações Diferenciais

Um sistema de equações diferenciais é formado quando precisamos determinarduas ou mais funções. As Eq.(4) e Eq.(5), que são conhecidas como equações de Lotka-Volterra(HAIRER; LUBICH; WANNER, 2006), ou predador-presa, são um exemplo desistema de equações diferenciais,

dx

dt= ax− αxy (4)

19

dy

dt= −cy + γxy (5)

2.6 Ordem de uma Equação Diferencial

A ordem de uma equação diferencial é a ordem da derivada de maior ordem queaparece na equação. A Eq.(6), que representa a lei de resfriamento de Newton, é umaequação diferencial ordinária de primeira ordem, enquanto a Eq.(7), que é a equaçãoque rege o movimento de um pêndulo simples para pequenos ângulos, é uma equaçãodiferencial ordinária de segunda ordem,

du

dt= −K(u− T ) (6)

d2θ

dt2+ g

ldθ = 0. (7)

De uma maneira geral, podemos expressar uma equação diferencial ordinária de ordem nde acordo com a Eq.(8),

F [t, u(t), u′(t), ..., u(n)(t)] = 0 (8)

Normalmente, a Eq.(8) é expressa pela Eq.(9), onde u(t) é substituído por y, e u′(t),u”(t), ..., u(n)(t) por y′, y”, ..., y(n), respectivamente,

F (t, y, y′, ..., y(n)) = 0 (9)

2.7 Equações Diferenciais Lineares e Não Lineares

A equação diferencial ordinária F (t, y, y′, ..., y(n)) = 0 é dita linear se F é umafunção linear das variáveis y, y′, ..., y(n); uma definição semelhante pode ser aplicada àsequações diferenciais parciais. Desta forma, a equação diferencial ordinária linear geral édada pela Eq.(10),

a0(t)y(n) + a1(t)y(n−1) + ...+ an(t)y = g(t). (10)

20

Uma equação que não é da forma da Eq.(10) é uma equação não linear. Um exemplo deequação diferencial ordinária linear é a Eq.(11),

t2d2y

dt2+ t

dy

dt+ 2y = sin(t) (11)

e um exemplo de equação diferencial ordinária não linear é a Eq.(12),

dy

dt+ ty2 = 0 (12)

2.8 Métodos Numéricos de Integração

Atualmente, sabemos como calcular integrais definidas de uma função sem o usode uma expressão analítica para a sua primitiva. Entretanto, para isso, foi necessário odesenvolvimento de alguns métodos de integração, denominados métodos numéricos.

Os métodos numéricos de integração são utilizados pelo fato de nem todas as fun-ções terem uma primitiva bem definida. Além disso, podemos ter uma função primitivamuito complexa de forma que sua integração, analiticamente, torna-se inviável. Os méto-dos numéricos de integração mais conhecidos são: Método de Euller ou Método da RetaTangente, o Método de Runge-Kutta e os Métodos de Adams-Bashforth e Adam-Moulton.

2.9 Método de Euler ou Método da Reta Tangente

Para entendermos o método de Euler, vamos considerar a equação diferencial deprimeira ordem, Eq.(13),

dy

dt= f(t, y) (13)

com condições iniciais tais que y(t0) = y0. Vamos considerar que as funções f e fy sejamcontínuas em algum retângulo no plano ty contendo o ponto (t0, y0). Então, pelo teoremada unicidade (única solução), existe uma única solução y = φ(t) do problema dado, emalgum intervalo em torno de t0. Se a Eq.(13) for não linear, então pode ser complicada adeterminação do intervalo de existência da solução, além de poder haver uma complicadarelação com a função f .

Agora, vamos escrever a Eq.(13) no ponto t0 na forma,

dφ

dt(tn) = f [tn, φ(tn)]. (14)

21

Aproximando a derivada na Eq.(14) pelo quociente de diferenças, obtemos a Eq.(15),

φ(tn+1)− φ(tn)tn+1 + tn

∼= f [tn, φ(tn)]. (15)

Agora, se substituirmos φ(tn+1) e φ(tn) pelos seus valores aproximados yn+1 e yn, respec-tivamente, e resolvermos para yn+1, obteremos a fórmula de Euler,

yn+1 = yn + f(tn, yn)(tn+1 − tn) n = 0, 1, 2, ..., . (16)

Se o tamanho do passo tn+1 − tn tiver valor uniforme h para todo n e se denotarmosf(tn, yn) por fn, então a Eq.(16) fica,

yn+1 = yn + hfn, n = 0, 1, 2, ..., . (17)

O método de Euler se resume em calcular, repetidamente, a Eq.(16) ou a Eq.(17),usando o resultado de cada passo para realizar o próximo passo. Dessa forma, obtemosuma sequência de valores y0, y1, y2, ..., yn, ... que aproximam os valores da solução φ(t) nospontos t0, t1, t2, ..., tn, ..., .

Se pensarmos num programa para utilizar o método de Euler, o mesmo terá umaestrutura parecida com a representada no esquema da figura 3.

Generalizando o método de Euler para sistemas de equações diferencias, temos quesua fórmula fica xn+1 = xn + hfn, onde X é um vetor e f é uma função vetorial.

2.10 Método de Runge-Kutta de quarta ordem em quatro estágios

A fórmula de Runge-Kutta é composta de uma média ponderada de valores def(t, y) em pontos diferentes no intervalo tn ≤ t ≤ tn+1 e dada pela Eq.(18),

yn+1 = yn + h(kn1 + 2kn2 + 2kn3 + kn4

6 ), (18)

em que,

kn1 = f(tn, yn),

kn2 = f(tn + 12h, yn + 1

2hkn1),

kn3 = f(tn + 12h, yn + 1

2hkn2),

kn4 = f(tn + h, yn + hkn3).

22

Figura 3 - Método de Euler

Legenda: Programa para cálculo do Método de EulerFonte: BOYCE; DIPRIMA, 2015, p.377.

23

A soma (kn1+2kn2+2kn3+kn46 ) pode ser considerada como um coeficiente angular médio.

Observa-se que kn1 é o coeficiente angular no extremo esquerdo do intervalo, kn2 é ocoeficiente angular no ponto médio usando a fórmula de Euler para ir de tn a tn + h

2 ,kn3 é a segunda aproximação do coeficiente angular no ponto médio e kn4 é o coeficienteangular em tn + h usando a fórmula de Euler e o coeficiente angular kn3, para ir de tn atn + h.

Da mesma forma que o método de Euler, o método de Runge-Kutta pode ser ge-neralizado para sistemas de equações diferenciais, e sua fórmula fica,

xn+1 = xn + h(kn1+2kn2+2kn3+kn46 ), em que:

kn1 = f(tn, xn),

kn2 = f(tn + 12h, xn + 1

2hkn1),

kn3 = f(tn + 12h, xn + 1

2hkn2),

kn4 = f(tn + h, xn + hkn3)

Analogamente ao método de Euler, podemos atribuir um programa simples parao método de Runge-Kutta de acordo com o quadro da figura 4.

2.11 Método de Passos Múltiplos

Até agora os métodos vistos foram métodos de partida ou de passo único, onde ovalor de y = φ(t) em qualquer ponto da partição depende, apenas, dos dados no pontoanterior da partição. O método de passos múltiplos (ROMA; NÓS, 2012) é um métodoque utiliza informação em mais do que o último ponto da partição.

2.11.1 Método de Adams-Bashforth

Sabendo que a Eq.(19) é,

φ(tn+1)− φ(tn) =∫ tn+1

tnφ′(t)dt (19)

24

Figura 4 - Método de Runge-Kutta

Legenda: Programa para o cálculo do Método deRunge-kutta

Fonte: BOYCE; DIPRIMA, 2015, p.391.

25

em que φ(t) é a solução do problema de valor inicial y′ = f(t, y) com y(t0) = y0, o métodode Adams tem como objetivo aproximar φ′(t) por um polinômio Pk(t) de grau k e usar opolinômio para calcular a integral na Eq.(19). Os coeficientes de Pk(t) são determinadoscom a utilização dos k + 1 dados calculados anteriormente. Para ilustrar, suponha quedesejamos usar um polinômio de grau P1(t) = At+B. para isso, necessitamos somente, dedois pontos de dados, (tn, yn) e (tn−1, yn−1). Para P1 ser uma aproximação de φ′, é precisoque P1(tn) = f(tn, yn) e P1(tn−1) = f(tn−1, yn−1). É importante saber que denotamosf(tj, yj) por fj para j inteiro. Temos, portanto, que A e B tem que satisfazer as equações,

Atn +B = fn

Atn−1 +B = fn−1

Resolvendo para A e B, obtemos,

A = fn−fn−1h

B = fn−1fn+fnfn−1h

.

Substituindo φ′(t) por P1(t) e calculando a integral na Eq.(19), vemos que

φ(tn+1)− φ(tn) = A2 (t2n+1 − t2n).

Finalmente, substituindo φ(tn+1) e φ(tn) por yn+1 e yn, respectivamente, e fazendo al-gumas simplificações algébricas e considerando um passo constante h, obtemos a Eq.(20)

yn+1 = yn + 32hfn −

12hfn−1, (20)

que é a chamada fórmula de Adams-Bashforth de segunda ordem.

2.11.2 Método de Adam-Moulton

Outro método da fórmula de Adams é o método de Adam-Moulton (GARCIA,2012), que é obtido através de uma variação na fórmula de Adams-Bashforth. Usandonovamente um polinômio de primeiro grau, Q1 = αt + β, determinamos os coeficientesusando os pontos (tn, yn) e (tn+1, yn+1). Dessa forma, α e β têm que satisfazer as equações,

26

fn = αtn + β

fn+1 = αtn+1 + β

que dão como resultado,

α = fn+1−fn

h

β = fntn+1−fn+1tnh

.

Substituindo φ′(t) na Eq.(19) por Q1 e simplificando, obtemos a Eq.(21),

yn+1 = yn + 12hfn + 1

2hf(tn+1, yn+1) (21)

que é a fórmula de Adam-Moulton de segunda ordem.

27

3 CONSTRUÇÃO DE UM ALGORITMO DE INTEGRAÇÃO

Neste capítulo vamos apresentar um algoritmo que desenvolvemos para aumentara precisão de um integrador numérico usando um integrador mais preciso em um númeroreduzido de passos.

• Em primeiro lugar vamos apresentar o método de integração numérica ba-seado no desenvolvimento das soluções de um sistema em séries de Taylormostrando, com mais detalhes, a estrutura matemática envolvida no pro-cesso. Em seguida, apresentaremos um exemplo prático da aplicação dométodo.• Em seguida vamos mostrar um resultado importante envolvendo a ideia

contida no método de séries de Taylor e, a partir desse resultado, construirum procedimento para melhorar a precisão de um integrador numérico.A ideia é que possamos fazer isso usando um integrador mais preciso emum número limitado de pontos e, portanto, com uso reduzido de tempo/ memória.

3.1 O Método de Taylor

Nesta seção vamos mostrar o funcionamento do método de Taylor.

3.1.1 A Base Matemática

Considere um sistema dinâmico autônomo em n variáveis dado por

xi = dxidt

= fi(x1, x2, · · · , xn) , (i = 1, · · · , n). (22)

Considere que xi = φi(t) é uma curva-solução (em forma paramétrica) passando peloponto ~x0 em t = 0. Desenvolvendo as funções φi(t) em série de Taylor em torno do pontot = 0, vamos obter:

xi = φi(t) = φi(0) + dφidt

(0) t+ d2φidt2

(0) t2

2! + · · · , (23)

Como o sistema é definido por (22), temos que (sobre a curva-solução xi = φi(t))

xi = φi(t) = fi(~x(t)) ⇒ φi(0) = fi(~x(0)) = fi(~x0), (24)

28

implicando que, em segunda ordem, ficamos com:

d

dt

(dφidt

)= dfidt

=∑j

∂fi∂xj

dxjdt

=∑j

∂fi∂xj

fj =∑j

fj∂

∂xj[fi]. (25)

Uma vez que podemos escrever fi como dxidt

, temos (omitindo o ∑ para índices repetidos)

d

dt

(dφidt

)= fj

∂

∂xj[fi] = fj

∂

∂xj

[∂xi∂xk

dxkdt

]= fj

∂

∂xj

[fk

∂

∂xk[xi]

]. (26)

Assim, definindo

X ≡ fi∂

∂xi, (27)

podemos escrever

du

dtu= Xu. (28)

Finalmente, temos:

xi(t) = xi(0) + tX [xi](0) + t2

2!X2[xi](0) + · · · =

=(

1 + tX + t2

2!X2 + · · ·

)[xi](0) =

∞∑k=0

tk

k!Xk[xi](0). (29)

Baseados em (29) podemos entender o método de séries de Taylor para integrarum sistema dinâmico. Olhando com atenção para o resultado final

xi(t) =∞∑k=0

tk

k!Xk[xi](0), (30)

é fácil perceber que, se quisermos calcular xi(δt), onde δt � 1, podemos aproximar oresultado truncando a série infinita em uma dada ordem N, isto é,

xi(δt) ≈N∑k=0

δtk

k! X k[xi](0) = Fi(~x0; δt). (31)

Quanto maior o N, maior (teoricamente) a aproximação1. Dessa maneira, a partir de um

1 Em se tratando de uma integração numérica de um sistema dinâmico específico, a melhora na aproxi-mação dependeria do número de dígitos empregados pelo computador para fazer as contas. A melhoriacontínua só aconteceria se o número de dígitos empregados fosse arbitrário – veja o capítulo 1 sobreerros de arredondamento.

29

dado ponto inicial ~x0 podemos calcular um ponto subsequente da trajetória (que podemoschamar de ~x1). A partir do ponto ~x1 podemos (usando o mesmo processo) calcular umponto subsequente ~x2, onde ~x2 = ~F (~x1; δt). E assim sucessivamente.

Na seção seguinte vamos apresentar um exemplo prático de integração numéricafazendo uso do método de séries de Taylor.

3.1.2 Aplicando o Método de Taylor ao Sistema de Lorenz

Nesta seção vamos exemplificar o uso do método de Taylor usando um sistemadinâmico caótico bastante famoso, importante e (apesar disso) um dos mais simples queexistem, o conhecido sistema de Lorenz (LORENZ, 1963).

x = σ(y − x),

y = −y − xz +Rx, (32)

z = xy − bz,

onde σ, R e b são parâmetros constantes do problema (com comportamento caótico paraR > 24, 74). Por que mais simples? Percebam que esse sistema possui o número mínimode equações diferenciais autônomas (o tempo não aparece explicitamente) para que hajacaos: três (em duas dimensões não há caos, uma vez que as trajetórias não podem secruzar). Além disso, como o caos é um fenômeno que só se manifesta em sistemas não-lineares, a ‘menor parcela’ de não-linearidade que poderíamos pensar em acrescentar aum sistema linear para torná-lo caótico, seria um termo quadrático. Podemos notar entãoque o sistema de Lorenz possui apenas dois termos não-lineares quadráticos. Considereuma condição inicial x(0) = x0, y(0) = y0, z(0) = z0. Podemos expandir a soluçãocorrespondente a essa condição inicial em série de Taylor, obtendo:

x = φ(t) = φ(0) + dφ

dt(0) t+ d2φ

dt2(0) t

2

2! + · · · ,

y = γ(t) = γ(0) + dγ

dt(0) t+ d2γ

dt2(0) t

2

2! + · · · , (33)

z = η(t) = η(0) + dη

dt(0) t+ d2η

dt2(0) t

2

2! + · · · .

Como o sistema é definido pelas equações

dx

dt= f(x, y, z),

dy

dt= g(x, y, z), (34)

dz

dt= h(x, y, z),

30

temos que x = φ = f, y = γ = g, z = η = h, implicando que, em segunda ordem, ficamoscom:

d

dt

(dφ

dt

)= df

dt= ∂f

∂x

dx

dt+ ∂f

∂y

dy

dt+ ∂f

∂z

dz

dt= ∂f

∂xf + ∂f

∂yg + ∂f

∂zh,

d

dt

(dγ

dt

)= dg

dt= ∂g

∂x

dx

dt+ ∂g

∂y

dy

dt+ ∂g

∂z

dz

dt= ∂g

∂xf + ∂g

∂yg + ∂g

∂zh, (35)

d

dt

(dη

dt

)= dh

dt= ∂h

∂x

dx

dt+ ∂h

∂y

dy

dt+ ∂h

∂z

dz

dt= ∂h

∂xf + ∂h

∂yg + ∂h

∂zh,

que podemos escrever:

φ = fx f + fy g + fz h,

γ = gx f + gy g + gz h, (36)

η = hx f + hy g + hz h,

e assim por diante até a ordem desejada. Usando o resultado em (33) vamos obteruma formulação que, para um δt � 1 determinado, gera um mapeamento que, apli-cado a um dado ponto (xn, yn, zn) de uma curva-solução obtém outro ponto subsequente(xn+1, yn+1, zn+1):

xn+1 = F (xn, yn, zn; δt),

yn+1 = G(xn, yn, zn; δt), (37)

zn+1 = H(xn, yn, zn; δt),

Em outras palavras, o mapeamento (37) ‘anda’ sobre a curva-solução do sistema de Lorenz.É importante perceber que o mapeamento não é exato uma vez que a série foi truncadaem alguma ordem. Esse fato está intimamente conectado à ideia que nos permitiu geraro algoritmo de integração associativa (descrito na próxima seção).

3.2 Um Algoritmo de Integração Associativa (AI)

Nesta seção vamos mostrar como construir um algoritmo de integração usando doisintegradores com precisão diferenciada. Esta seção está organizada da seguinte maneira:

• Para clarear a ideia principal sobre a qual se baseia a construção donosso algoritmo, vamos apresentar alguns resultados envolvendo mapea-mentos que representam soluções de sistemas de equações diferenciais –veja eq.(30).• A seguir, vamos construir um procedimento para melhorar a precisão de

31

um integrador. A ideia é que possamos fazer isso usando um integradormais preciso em um número limitado de pontos e, portanto, com usoreduzido de tempo / memória.

3.2.1 O Resultado Principal

Considere o sistema dinâmico definido por

xi = fi(~x), (i = 1, · · · , n), (38)

cuja solução é

xi(t) = Fi(~x0, t), (i = 1, · · · , n), (39)

onde fi(~x) ≡ ∂Fi(~x,t)∂t

|t=0 e xi ≡ dxi

dt. Vimos que a solução (39) do sistema dinâmico (38)

pode ser obtida a partir do operador X ≡ ∑ni=1 fi

∂∂xi

:

xi(t) = Fi(~x0, t) = xi(0) + tX [xi](0) + t2

2!X2[xi](0) + · · · =

∞∑k=0

tk

k!Xk[xi](0). (40)

Dessa maneira, começando de um ponto genérico P0 (de coordenadas ~x(P0)) e escolhendoum intervalo de tempo δt, podemos gerar um mapeamento M que leva um ponto de umadada curva-solução do sistema para outro ponto (da mesma curva-solução) que corres-ponde a um incremento δt no tempo, isto é, ao ponto que corresponde à curva-solução dosistema após a passagem de um tempo δt a partir da posição P0:

xi(P+δP ) = Fi(~x(P ), δt) =∞∑k=0

δtk

k! X k[xi](P ). (41)

Como vimos, o processo de resolver numericamente o sistema pode ser resumido emescolher um pequeno intervalo de tempo (δt � 1) e truncar o mapeamento M (41) emalguma ordem N, obtendo o mapeamento M dado por:

xi(P+δP ) = F i(~x(P ), δt) =N∑k=0

δtk

k! X k[xi](P ), (42)

onde xi(P+δP ) se aproxima de xi(P+δP ) quando δt → 0. Vamos agora a um resultadointeressante:

32

Definição 3.1 Considere as funções εi e δrεi definidas como

εi(~x)(P ) ≡ xi(P+δP ) − xi(P+δP ) =∞∑

k=N+1

δtk

k! X k[xi](P ) (43)

δεi(~x)(P ) ≡ εi(~x)(P+δP ) − εi(~x)(P ) (44)

δrεi(~x)(P ) ≡ δr−1εi(~x)(P+δP ) − δr−1εi(~x)(P ) (r = 2, . . .) (45)

Teorema 3.1 Considere o sistema dinâmico xi = fi(~x) cuja solução é xi(t) = Fi(~x0, t) eo mapeamento M dado por xi(P+δP ) = Fi(~x(P ), δt) = ∑∞

k=0δtk

k! X k[xi](P ), onde o operadorX é definido por X ≡ ∑n

i=1 fi∂∂xi

. Considere ainda as funções εi e δrεi como definidasacima. Se as funções fi e suas derivadas existem em um ponto P (de coordenadas ~x(P ))e em uma vizinhança não nula de P então temos que:

limδt→0

δr+1εi(~x)(P )

δrεi(~x)(P )= 0, (46)

onde r é um inteiro positivo.

Prova 3.1 Da definição (3.1) podemos escrever

δεi(~x)(P ) = εi(~x)(P+δP ) − εi(~x)(P ) =∞∑

k=N+1

δtk

k! X k[xi](P+δP ) −∞∑

k=N+1

δtk

k! X k[xi](P ). (47)

Como xi(P+δP ) = xi(P ) + δxi, temos (definindo Φki(~x)(P ) = X k[xi](P )):

δεi(~x)(P ) =∞∑

k=N+1

δtk

k!(X k[xi](P+δP ) − X k[xi](P )

)=

=∞∑

k=N+1

δtk

k!(Φk

i(~x) |(P+δP ) −Φki(~x) |(P )

)=

=∞∑

k=N+1

δtk

k!

n∑j=1

∂ Φki(~x)

∂xj|(P ) δxj + O(δxj2)

, (48)

δ2εi(~x)(P ) =∞∑

k=N+1

δtk

k!

n∑j=1

(∂ Φk

i(~x)∂xj

|(P+δP ) −∂ Φk

i(~x)∂xj

|(P )

)δxj + O(δxj2)

=

=∞∑

k=N+1

δtk

k!

n∑j=1

∂2 Φki(~x)

∂xl∂xj|(P ) δxl δxj + O(δx3)

, (49)

33

· · ·

δrεi(~x)(P ) =∞∑

k=N+1

δtk

k!

n∑j1,··· ,jr=1

∂r Φki(~x)

∂xj1 · · · ∂xjr|(P ) δxjl · · · δxjr + O(δxr+1)

. (50)

Logo,

limδt→0

δr+1εiδrεi

=

∑∞k=N+1

δtk

k!

(∑nj1,··· ,jr+1=1

∂r+1 Φki(~x)

∂xj1 ···∂xjr+1|(P ) δxjl · · · δxjr+1 + O(δxr+2)

)∑∞k=N+1

δtk

k!

(∑nj1,··· ,jr=1

∂r Φki(~x)

∂xj1 ···∂xjr|(P ) δxjl · · · δxjr + O(δxr+1)

) . (51)

Como estamos em pontos regulares do sistema, então limδt→0δxδt

= K, onde K éuma n−upla de números reais finitos e, portanto, δt → 0 implica que δx → K δt o que,por sua vez, implica que o limite acima é zero.

Na seção seguinte vamos mostrar, baseados nesse resultado, uma maneira de con-jugar dois integradores de diferentes precisões para realizar o trabalho com precisão com-parável à do integrador mais preciso com menos tempo de uso de CPU.

3.2.2 A Base Estrutural de um AI

Vamos considerar um sistema dinâmico de baixa dimensionalidade, isto é, umsistema ~x = ~f(~x) onde a dimensão do vetor ~x é um inteiro n não muito grande2. Vimosque podemos desenvolver as soluções do sistema em série de Taylor e, para um dadointervalo de tempo δt, construir um mapeamento M dado por

xi(P+δP ) =∞∑k=0

δtk

k! X k[xi](P ) = Fi(~x(P ); δt),

que, a partir de um dado ponto P0 (de coordenadas ~x0), calcula os pontos subsequentes dacurva-solução (que passa por P0) separados por um intervalo de tempo δt. Teoricamente,o mapeamento M faria isso com uma precisão infinita. Uma vez que isso não é possívelno mundo real (pois o número de termos na série é infinito), o que fazemos para calcularos pontos na prática é truncar a série em alguma ordem N finita. Considere que fazemosisso e escolhemos um δt � 1 (fazendo com que as ordens mais altas sejam desprezíveis).Obtemos assim um mapeamentoM , com uma precisão correspondente a um certo númeroD de dígitos. Considere que escolhemos um ponto ‘inicial’ P0 e que podemos ter acessoa um certo número de pontos da curva-solução com uma precisão arbitrária PrA �

2 Por não muito grande queremos dizer, na prática, não maior que 12.

34

D. Considere, por fim, as funções δrεi (a função εi corresponde a r = 0) nos pontossubsequentes ao ponto P0, gerados pelo mapeamento M . Podemos observar:

Observação 3.1 Como δt� 1 e o mapeamento M corresponde ao truncamento do ma-peamento M em uma ordem N, existe um número inteiro positivo máximo L tal que

δLε

δL−1ε< 1. (52)

Observação 3.2 O número L é tanto maior quanto maior é o número N (a ordem detruncamento) e quanto menor é o número δt (o intervalo de tempo entre dois pontosconsecutivos da sequência de pontos gerados pelo mapeamento M).

Observação 3.3 Como consequência do fato que δt � 1 e do teorema 3.1, os valoresdos módulos das funções δrεi são � 1.

Vamos examinar agora as funções δrεi no pontos P0 e nos pontos subsequentesgerados pelo mapeamento M . Em P0, temos:

δ0εi(~x)(P0) = εi(~x)(P0) = xi(P0+δP0) − xi(P0+δP 0) =∞∑

k=N+1

δtk

k! X k[xi](P0). (53)

Chamaremos (para simplificar a notação) P1 ≡ P0 + δP 0 e P 1 ≡ P0 + δP 0, e, de maneirasimilar, Pj+1 ≡ Pj + δP j e P j+1 ≡ Pj + δP j. Considere agora os p pontos após P0:

δ0εi(~x)(P1) = εi(~x)(P1) = xi(P1+δP1) − xi(P1+δP 1) =∞∑

k=N+1

δtk

k! X k[xi](P1). (54)

· · ·

δ0εi(~x)(Pp−1) = εi(~x)(Pp−1) = xi(Pp−1+δPp−1) − xi(Pp−1+δP p−1) =∞∑

k=N+1

δtk

k! X k[xi](Pp−1). (55)

Para as funções δεi, temos:

δεi(~x)(P0) = εi(~x)(P0+δP0) − εi(~x)(P0), (56)

δεi(~x)(P1) = εi(~x)(P1+δP1) − εi(~x)(P1), (57)

· · ·

δεi(~x)(Pp−1) = εi(~x)(Pp−1+δPp−1) − εi(~x)(Pp−1), (58)

35

Podemos exprimir, para as funções δrεi,

δrεi(~x)(Pj−1) = δr−1εi(~x)(Pj−1+δPj−1) − δr−1εi(~x)(Pj−1). (59)

Agora, estamos prontos para apresentar a ideia central que permite a construçãodo algoritmo: Considere que temos p pontos P1, P2, . . . , Pp, gerados a partir de P0 pelomapeamentoM (com precisão arbitrária que, para efeitos práticos, vamos supor infinita).Podemos observar que:

Observação 3.4 Como δt � 1, o ponto P1 é muito próximo3 do ponto P0, ou seja adistância de P0 a P1 é muito menor que 1 (dd(P0, P1) � 1). Isso implica que |δxi| � 1.Além disso, como

xi(P+δP ) =∞∑k=0

δtk

k! X k[xi](P ) = xi(P ) +∞∑k=1

δtk

k! Φki(~x)(P ),

⇒ xi(P+δP ) − xi(P ) = δxi = Φ1i(~x)(P )δt+

∞∑k=2

δtk

k! Φki(~x)(P )

⇒ O(δxi) ≈ O(δt). (60)

Observação 3.5 De acordo com a definição das funções δrεi, temos que:

• De (53) temos εi(~x)(P ) = ∑∞k=N+1

δtk

k! Φki(~x)(P ) ≈ O(δtN+1).

• De (47,56) temos que

δεi(~x)(P ) = εi(~x)(P+δP )−εi(~x)(P ) =∞∑

k=N+1

δtk

k!

n∑j=1

∂ Φki(~x)

∂xj|(P ) δxj + O(δxj2)

,⇒ δεi(~x)(P ) ≈ O(δtN+2).• Finalmente, de (50) temos que

δrεi(~x)(P ) =∞∑

k=N+1

δtk

k!

n∑j1,··· ,jr=1

∂r Φki(~x)

∂xj1 · · · ∂xjr|(P ) δxjl · · · δxjr + O(δxr+1)

,⇒ δrεi(~x)(P ) ≈ O(δtN+r+1).

Suponha que, a partir de um dado ponto P0, usamos o mapeamentoM para calcu-larmos p pontos P1, . . . , Pp e, a partir desses usamos o mapeamento M para calcularmos

3 Estamos considerando a métrica usual, isto é, a distância de P0 a P1 como dd(P0, P1) ≡√∑

(δxi(P0))2.

36

p + 1 pontos P 1, . . . , P p+1. Suponha agora que queremos calcular o ponto Pp+1, aproxi-madamente, sem usar o mapeamento M , ou seja, queremos calcular xi(Pp+1). Sabemosde (55) que εi(~x)(Pp) = xi(Pp+1)− xi(P p+1), o que implica que xi(Pp+1) = xi(P p+1) + εi(~x)(Pp).Portanto, se soubéssemos εi(~x)(Pp), poderíamos calcular xi(Pp+1). Contudo, de (58), sabe-mos que δεi(~x)(Pp−1) = εi(~x)(Pp)− εi(~x)(Pp−1) e, assim, εi(~x)(Pp) = εi(~x)(Pp−1) + δεi(~x)(Pp−1).Isso nos leva a que

xi(Pp+1) = xi(P p+1) + εi(~x)(Pp−1) + δεi(~x)(Pp−1). (61)

Notem que os dois primeiros termos do lado direito de (61) são conhecidos e, além disso,sua soma (xi(P p+1) + εi(~x)(Pp−1)) representa uma aproximação melhor que xi(P p+1) para ovalor de xi(Pp+1). Para entender porque, basta dar uma olhada na análise que fizemos nasobservações 3.4 e 3.5:

εi(~x)(P ) ≈ O(δtN+1), δεi(~x)(P ) ≈ O(δtN+2). (62)

Assim, o erro na aproximação por xi(P p+1) é da ordem O(δtN+1), ao passo que o errona aproximação por xi(P p+1) + εi(~x)(Pp−1) é da ordem O(δtN+2). Esse raciocínio podeser estendido: não conhecemos o valor de δεi(~x)(Pp−1), mas sabemos que δ2εi(~x)(Pp−2) =δεi(~x)(Pp−1) − δεi(~x)(Pp−2) e, portanto, δεi(~x)(Pp−1) = δεi(~x)(Pp−2) + δ2εi(~x)(Pp−2), o queimplica em

xi(Pp+1) = xi(P p+1) + εi(~x)(Pp−1) + δεi(~x)(Pp−2) + δ2εi(~x)(Pp−2). (63)

Analogamente à situação anterior, não conhecemos δ2εi(~x)(Pp−2) mas, sendo δ2εi(~x)(P ) ≈O(δtN+3), a aproximação de xi(Pp+1) por

xi(Pp+1) ≈ xi(P p+1) + εi(~x)(Pp−1) + δεi(~x)(Pp−2), (64)

é melhor que a aproximação

xi(Pp+1) ≈ xi(P p+1) + εi(~x)(Pp−1). (65)

Podemos continuar com esse raciocínio, o que nos leva ao seguinte resultado:

xi(Pp+1) = xi(P p+1)+εi(~x)(Pp−1)+δεi(~x)(Pp−2)+δ2εi(~x)(Pp−3)+· · ·+δrεi(~x)(Pp−r−1)+· · · . (66)

Na seção seguinte vamos mostrar como usar o resultado (66) para construir umalgoritmo de integração.

37

3.2.3 O Algoritmo Propriamente Dito

O mapeamento M definido na seção anterior não pode ser usado na prática poisenvolve infinitas operações. Podemos, no entanto, usar praticamente o mapeamento Mque é definido por xi(P+δP ) = ∑N

k=0δtk

k! X k[xi](P ), uma vez que N é um inteiro positivo.Nesta seção vamos construir um algoritmo de integração associando dois mapeamentosde precisões diferenciadas.

Definição 3.2 Considere os mapeamentos M+ e M− definidos por:

M+ ≡ xi(P+δP ) =N+∑k=0

δtk

k! X k[xi](P ), (67)

M− ≡ xi(P+δP ) =N−∑k=0

δtk

k! X k[xi](P ), (68)

onde N+ e N− são dois inteiros positivos tais que N+ > N−.

A ideia é utilizar o mapeamento mais preciso (M+) de maneira análoga ao ma-peamento M e o mapeamento menos preciso (M−) de maneira análoga ao mapeamentoM . Notem que, nessa situação ‘real’, ambos os conjuntos de n−uplas que representam ospontos P+ e P− possuem um número limitado de dígitos e, dessa maneira, a aproxima-ção análoga à aproximação representada por (66), também só valeria até algum r finito.Adicionando esse resultado ao que apresentamos na seção anterior, podemos construir umalgoritmo que conjuga o uso dos mapeamentos M+ e M−. Para isso, vamos definir asfunções ∆kεi (análogas às funções (45)) como:

Definição 3.3 Considere um ponto P0 e, usando o mapeamento M+ definido acima,calculamos um número p de pontos que imediatamente sucedem o ponto P0. Vamosdenominá-los P+

1 , P+2 , . . . , P

+p (onde P+

1 ≡M+[P0] e P+j+1 ≡M+[P+

j ]). Usando o mape-amento M− definido acima, podemos, usando os pontos P0, P

+1 , P

+2 , . . . , P

+p , determinar

p+ 1 pontos definidos por:

P−1 ≡M−[P0], P−2 ≡M−[P+1 ], P−3 ≡M−[P+

2 ], . . . , P−p+1 ≡M−[P+p ]. (69)

A partir dos pontos P0, P+ e P− podemos definir funções ∆rεi (análogas às funções (45))

38

como:

∆0εi(~x)(P+j ) ≡ xi(P+

j ) − xi(P−j ), (j = 1, . . . , p),

∆1εi(~x)(P+j ) ≡ ∆0εi(~x)(P+

j ) −∆0εi(~x)(P+j−1), (j = 2, . . . , p),

... ...

∆rεi(~x)(P+j ) ≡ ∆r−1εi(~x)(P+

j ) −∆r−1εi(~x)(P+j−1), (j = r + 1, . . . , p). (70)

(71)

A partir das funções ∆rεi definidas acima e dos pontos P0, P+ e P− podemos calcularuma aproximação para o ponto P+

p+1 sem usar o mapeamento M+. Começaremos comuma aproximação ‘de ordem zero’:

xi(P+p+1) ≈ xi(P−

p+1), (72)

onde a ordem da aproximação é, de acordo com o que mostramos na observação 3.5,O(δt(N+−N−+1)). Mas sabemos que xi(P+

p+1) − xi(P−p+1) = ∆0εi(~x)(P+

p+1) e ∆1εi(~x)(P+p+1) =

∆0εi(~x)(P+p+1) −∆0εi(~x)(P+

p ), portanto,

xi(P+p+1) = xi(P−

p+1) + ∆0εi(~x)(P+p ) + ∆1εi(~x)(P+

p+1), (73)

levando à aproximação ‘de ordem um’:


p+1) + ∆0εi(~x)(P+p ), (74)

onde a ordem da aproximação é O(δt(N+−N−+2)). Mas sabemos que ∆2εi(~x)(P+p+1) =

∆1εi(~x)(P+p+1) −∆1εi(~x)(P+

p ) e, assim,

xi(P+p+1) = xi(P−

p+1) + ∆0εi(~x)(P+p ) + ∆1εi(~x)(P+

p ) + ∆2εi(~x)(P+p+1), (75)

levando à aproximação ‘de ordem dois’:


p+1) + ∆0εi(~x)(P+p ) + ∆1εi(~x)(P+

p ), (76)

onde a ordem da aproximação é O(δt(N+−N−+3)). Podemos continuar esse raciocínio e, deforma geral, temos:

xi(P+p+1) = xi(P−

p+1) + ∆0εi(~x)(P+p ) + · · ·+ ∆r−1εi(~x)(P+

p ) + ∆rεi(~x)(P+p+1), (77)

levando à aproximação ‘de ordem r’:


p+1) + ∆0εi(~x)(P+p ) + · · ·+ ∆r−1εi(~x)(P+

p ), (78)

39

onde a ordem da aproximação é O(δt(N+−N−+r+1)) e r ≤ p.

Observação 3.6 Usando a aproximação (78) podemos (usando o mapeamento M− e asfunções ∆rεi) calcular o ponto xi(P+

p+1) com uma precisão O(δtr) melhor que a precisãousual fornecida pelo mapeamento M−.

Observação 3.7 Dependendo do número de dígitos usados para representar os pontosda curva-solução do sistema e do número p de pontos usados para calcularmos as funções∆rεi, podemos, em teoria, chegar até uma ordem de aproximação semelhante à obtida pelomapeamento M+.

Passos do Algoritmo:

1. Escolha dois mapeamentos de precisões diferenciadas M+ e M−.

2. Escolha uma condição inicial P0 e um número de pontos Npoints a serem calculados.

3. Escolha um número inteiro positivo p e, usando o mapeamento mais preciso M+,calcule p pontos a partir da condição inicial P0 (os pontos P+).

4. A partir desses p + 1 pontos, calcule, usando o mapeamento menos preciso M−,p+ 1 pontos a partir dos p+ 1 pontos {P0, P

+1 , . . . , P

+p } (os pontos P−).

5. Usando o conjunto de pontos {P0} ∪ {P+} ∪ {P−} calcule as funções ∆rεi, (r =0, . . . , p).

6. Calcule a aproximação ao ponto xi(P+p+1) usando a relação xi(P+

p+1) ≈ xi(P−p+1) +

∆0εi(~x)(P+p ) + · · ·+ ∆r−1εi(~x)(P+

p ).

7. Calcule o ponto xi(P−p+2) usando o mapeamentoM− na aproximação ao ponto xi(P+

p+1).

8. Calcule as funções ∆rεi, r = (0, . . . , p) para o ponto xi(P+p+1).

9. Calcule a aproximação ao ponto xi(P+p+2) usando a relação xi(P+

p+2) ≈ xi(P−p+2) +

∆0εi(~x)(P+p+1) + · · ·+ ∆r−1εi(~x)(P+

p+1).

10. Repita os passos 7, 8 e 9 para os pontos P+p+2, . . . , P

+Npoints.

No próximo capítulo vamos apresentar uma implementação computacional paraum estudo do desempenho desse tipo de algoritmo.

40

4 IMPLEMENTAÇÃO DE UM AI

Um Estudo de Desempenho

Neste capítulo vamos implementar um algoritmo AI em um ambiente computaci-onal para estudar o seu desempenho quando aplicado a casos concretos. Nesse primeiroestudo vamos observar como se comporta um algoritmo AI no caso em que os mapeamen-tos M+ e M− são integradores do tipo Runge-Kutta (veja o capítulo 2). O capítulo seráorganizado da seguinte maneira:

• Primeiro vamos discutir uma implementação computacional de um AI:Escolhemos como mapeamentos M+ e M− os integradores numéricos dotipo Runge-Kutta (de sétima-oitava e quarta ordens, respectivamente).Como linguagem de programação escolhemos a plataforma Maple de com-putação simbólica4.• Em seguida, vamos comparar a aplicação de quatro integradores a dois

problemas clássicos de sistemas dinâmicos de baixa dimensionalidade queapresentam caos: O famoso sistema de Lorenz (um sistema caótico dis-sipativo) e o (também bastante conhecido) sistema de Hénon-Heiles (umsistema caótico conservativo). Os quatro integradores (usados para com-paração do desempenho) são:

– Taylor: Integrador usando o método de séries de Taylor para obtermosuma base segura de resultados com alta precisão.

– RK78: Integrador usando o método de Runge-Kutta de sétima-oitavaordem.

– RK4: Integrador usando o método de Runge-Kutta de quarta ordem.– AI(RK78−RK4): Integrador usando o nosso método (que conjuga os in-tegradores RK4 e RK78).

4 O ponto é que, embora a integração numérica não seja o forte do Maple, para um primeiro estudode desempenho comparativo, essa plataforma (pelo fato do Maple ser uma linguagem interpretada)permite a implementação muito mais simples de ferramentas de análise do que em linguagens maisadequadas para a integração numérica como o C, Fortran etc (que são linguagens compiladas).

41

4.1 Uma Implementação em Maple: O pacote AssISt

Nesta seção vamos mostrar um implementação computacional em Maple de umalgoritmo de integração associativa (AI), na qual usamos os integradores numéricos Runge-Kutta de quarta e sétima-oitava ordens (RK4 e RK78) como os integradores associados.

4.1.1 Um resumo do pacote

Dados Técnicos:

Título do Pacote: AssIST – Associative Integrator Study.

Sistema Operacional no qual o Pacote foi testado: Windows 7, Windows 8.

Plataforma Computacional usada: Maple 17.

Memoria requerida para executar os testes: Em torno de 200 Megabytes.

Número de linhas do programa: 637.

Palavras-chave: Implementação em Maple, Integradores de Runge-Kutta, Integração usando o Método

de Séries de Taylor, Integração Associada, Análise de Precisão e Desempenho.

Natureza do problema matemático

Análise de possíveis melhorias no desempenho / precisão de integradores numéricos.

Algoritmo Proposto:

Um algoritmo de integração que conjuga dois integradores de diferentes níveis de precisão.

Tempo de execução típico:

Intermediário entre o tempo de execução do integrador numérico mais simples e o tempo do integrador

mais preciso.

Sumário dos comandos:

• PTaylor - Essa rotina gera uma lista de pontos obtidos por integraçãonumérica do sistema de equações diferenciais em questão (o sistema di-nâmico estudado) usando o método de séries de Taylor.• Prk4maple - Essa rotina gera uma lista de pontos obtidos por integração

numérica do sistema de equações diferenciais em questão usando o métodode Runge-Kutta de quarta ordem.• Prk78maple - Essa rotina gera uma lista de pontos obtidos por integração

numérica do sistema de equações diferenciais em questão usando o método

42

de Runge-Kutta de sétima-oitava ordem.• Pcost - Essa rotina é a responsável pela ‘costura’ entre os pontos gerados

pelos mapeamentos RK4 e RK78.• Dif - Essa rotina calcula as diferenças (em diversos níveis – veja o capí-

tulo anterior na definição das funções ∆rεi) entre os pontos gerados porPrk78maple e Pcost.• Pai - Essa rotina gera uma lista de pontos obtidos por integração nu-

mérica do sistema de equações diferenciais em questão usando o nossométodo de associação de integradores.• Pan - Essa rotina faz diversas análises baseadas nas listas de pontos ge-

radas pelos integradores.

4.1.2 Os comandos do Pacote

Vamos agora a uma descrição mais detalhada dos comandos do Pacote.

4.1.2.1 Comando: PTaylor

Funcionalidade: Esse comando produz uma lista de listas com as coordenadas dos pontoscalculados pelo método de Taylor.

Chamada: 5

[> PTaylor(syst_dvd,ic,passo,Npoints);

Parâmetros:

syst_dvd = sistema de equações diferenciais (em um formato específico);

ic = condição inicial;

passo = o passo de integração;

Npoints = o número de pontos que se deseja obter.

5 No que se segue o símbolo [> representa a entrada para um comando em uma sessão de Maple.

43

4.1.2.2 Comando: Prk4maple

Funcionalidade: Esse comando produz uma lista de listas com as coordenadas dos pontoscalculados pelo método RK4.

Chamada:

[> Prk4maple(syst_dvd,ic,passo,Npoints);

Parâmetros: Os mesmos usados pelo comando PTaylor.

4.1.2.3 Comando: Prk78maple

Funcionalidade: Esse comando produz uma lista de listas com as coordenadas dos pontoscalculados pelo método RK78.

Chamada:

[> Prk78maple(syst_dvd,ic,passo,Npoints);

Parâmetros: Os mesmos usados pelo comando PTaylor.

4.1.2.4 Comando: Pcost

Funcionalidade: Esse comando produz uma lista de listas com as coordenadas dos pontosgerados com um mapeamento RK4 a partir de pontos gerados com um mapeamentoRK78.

Chamada:

[> Pcost(syst_dvd,ListT,passo,Npoints_cost);

Parâmetros:


ListT = Uma lista de listas com as coordenadas dos pontos gerados pelomapeamento RK78;


Npoints_cost = o número de pontos que se deseja ‘costurar’.

44

Detalhamento:

Pcost é a rotina que constrói a base para a comparação entre os níveis de precisãoatingidos pelos mapeamentos RK78 e RK4. A partir de um conjunto PP de pontossequenciais PP = {PP1, PP2, . . . , PPs} gerados pelo mapeamento RK78 ele constrói umoutro conjunto PPcost de pontos calculado pela aplicação do mapeamento RK4 a cada umdos pontos do conjunto PP , ou seja,PPcost = {MRK4[PP1],MRK4[PP2], . . . ,MRK4[PPs]}.

4.1.2.5 Comando: Dif

Funcionalidade: Esse comando produz uma lista de listas com as coordenadas dos pontoscorrigidos pelas funções ∆εi.

Chamada:

[> Dif(syst_dvd,passo,Npoints,Npoints_dif,ListT,ListC);

Parâmetros:


passo = o passo de integração.

Npoints = o número de pontos que se deseja obter.

Npoints_dif = o número de níveis de diferenças.

ListT = Uma lista de listas com as coordenadas dos pontos gerados pelomapeamento RK78;

ListC = Uma lista de listas com as coordenadas dos pontos gerados pelocomando Pcost a partir de ListT.

Detalhamento:

A rotina Dif é o ‘coração’ do algoritmo de integração associativa: ele recebe umconjunto de pontos calculados pelo mapeamento RK78 (ListT) e um outro conjunto(ListC) que foi gerado pelo comando Pcost a partir do conjunto ListT. A partir dessesdois conjuntos ele gera, inicialmente, um conjunto ({∆0εi(~x)PP2 = xi(PP2) − xi(PP 2), . . .})com as diferenças das coordenadas dos pontos dos dois conjuntos. Em seguida, calculaas diferenças das diferenças ({∆1εi(~x)PP3 = ∆0εi(~x)(PP3)−∆0εi(~x)(PP2), . . .}), e assim pordiante até atingir o número de níveis de diferença (Npoints_dif). Vamos ver isso num

45

exemplo simples: considere que

ListT = [PP1, PP2, PP3, PP4] ,

ListC =[PP 2, PP 3, PP 4, PP 5

],

Npoints_dif = 3,

onde

PP 2 = MRK4[PP1], PP 3 = MRK4[PP2], PP 4 = MRK4[PP3], PP 5 = MRK4[PP4].

Dif calcula os ∆0ε:

DD[0,2] = PP2 − PP 2, DD[0,3] = PP3 − PP 3, DD[0,4] = PP4 − PP 4;

os ∆1ε:

DD[1,3] = DD[0,3] −DD[0,2], DD[1,4] = DD[0,4] −DD[0,3];

e os ∆2ε:

DD[2,4] = DD[1,4] −DD[1,3].

Após isso, Dif calcula a aproximação para o ponto PP 5:

PP 5∼ = PP 5 +DD[0,4] +DD[1,4] +DD[2,4].

A partir do ponto PP 5∼, calcula PP 6 = MRK4[PP5∼] e as diferençasDD[0,5], DD[1,5], DD[2,5]

e a aproximação para o ponto PP 6:

PP 6∼ = PP 6 +DD[0,5] +DD[1,5] +DD[2,5];

e assim por diante.

4.1.2.6 Comando: Pai

Funcionalidade: Esse comando gerencia os comandos Prk78maple, Prk4maple, Pcost eDif para fazer a integração associada (isto é, ele representa o algoritmo como um todo).

Chamada:

[> Pai(syst_dvd,ic,passo,Npoints,Npoints_dif);

46

Parâmetros:


ic = a condição inicial;


Npoints = o número de pontos que se deseja obter;

Npoints_dif = o número de níveis de diferenças.

Detalhamento:

A rotina Pai representa o algoritmo: a partir da condição inicial ic ele calcula oconjunto de pontos ListT usando a rotina Prk78maple e, a partir desse, o conjunto depontos ListC com a rotina Pcost (que usa a rotina Prk4maple). A partir desses doisconjuntos ele gera os pontos subsequentes usando a rotina Dif.

4.1.2.7 Comando: Pan

Funcionalidade: Esse comando auxilia na análise do algoritmo para diversos sistemas/ condições iniciais / número de pontos integrados / número de pontos ‘costurados’ /precisão / tempo de CPU.

Chamada:

[> Pan(syst_dvd,ic,passo,Npoints,Npoints_dif,pop);

Parâmetros:


ic = a condição inicial;


Npoints = o número de pontos após os quais desejamos fazer a análise;

Npoints_dif = o número de níveis de diferenças;

pop = determina as opções de saída da análise produzida: uma lista comdois operandos – o primeiro relacionado à coordenada que desejamos analisar e o segundoao aspecto que desejamos analizar.

Detalhamento:

A rotina Pan apresenta opções de saída definidas pelo parâmetro pop = [n1, n2]. Onúmero n1 se refere à coordenada que será analisada e n2 ao desempenho dos integradores.

47

Dessa maneira, se o sistema tem, por exemplo, 3 dimensões (x1, x2, x3), o valor n1 = 1 sereferirá à coordenada x1; se n1 6= 1, 2, 3, a análise (designada por n2) será sobre todas ascoordenadas. Em relação à n2 temos 3 opções:

1. n2 = 1 - Análise do tempo gasto por cada um dos integradores.

2. n2 = 2 - Análise do tempo gasto e do número de dígitos significativos por cada umdos integradores.

3. n2 = 3 - Análise do tempo gasto, do número de dígitos significativos de todas ascoordenadas e um ‘printing’ da coordenada designada por n1.

Variáveis Globais:

A rotina Pan faz uso de algumas variáveis globais:

• Nheuristic - Esta variável se destina a determinar o número de vezesque os integradores calcularão os pontos determinados pelo parâmetroNpoints – isto se deve ao fato que, para poucos pontos, o contador doMaple (comando time()) não apresenta precisão suficiente.• Nremember - Quando queremos fazer uma análise diferente (por exem-

plo, de uma outra coordenada) usando a mesma entrada básica, estavariável (Nremember = 1) evita que os cálculos sejam refeitos.• PointsTay, PointsRK78, PointsRK4, PointsAI ; tTAY, tRK78, tRK4 ,tAI ; DD - Estas variáveis guardam os valores dos pontos calculados pe-los integradores e os respectivos tempos gastos no cálculo. A variável DDguarda as diferenças entre os valores do integrador que usa o método deTaylor e os demais.

4.1.3 Exemplificando o Uso do Pacote

Para exemplificar o uso dos comandos do pacote AssISt vamos usar o sistemade Hénon-Heiles

x = Az, (79)

y = w, (80)

z = −(A+ 2 k

Ay

)x, (81)

w = −k x2

A+ y + σ y2. (82)

48

Apresentaremos os comandos na mesma ordem em que foram descritos na seção anteriore simularemos a maneira como eles aparecem em uma sessão de Maple. Considere osseguintes parâmetros:

[> sigma := 5.000000: k:=0.999900: A:= sqrt(2*k/sigma +1):[> passo:=0.001:[> syst_dvd := [dvd[1]=A*vd[3],dvd[2]=vd[4],dvd[3]=-(A+2*k*vd[2]/A)*vd[1],

dvd[4]=-k*vd[1]^2/A + vd[2] + sigma*vd[2]^2]:[> Npoints:=80: Npoints_dif := 13: Npoints_cost := Npoints_dif+1:[> Digits := 55:[> ic := [1.000093761322944,-0.0006000959740288238,

-0.000040043010069843160, -0.8968057300859730197438631]:

Podemos usar os comandos PTaylor, Prk78maple e Prk4maple para integrar o sistema apartir da condição inicial dada. Por exemplo:

[> pontosTAYLOR := PTaylor(sysdvd,ic,passo,Npoints):[> pontosRK78 := Prk78maple(sysdvd,ic,passo,Npoints):[> pontosRK4 := Prk4maple(sysdvd,ic,passo,Npoints):

O sinal “:” após a linha de comando evita a impressão (na tela) do resultado. Se quisermosver algum dos pontos (por exemplo, o ponto 80 integrado pelo método de Taylor) podemosusar o comando

[> pontosTAYLOR[80];

[0.079, 0.9958778706064519485208857839805919534750878504070506864,

−0.07414440924614320256849352884246604440456934326251522196,

−0.08846938823439947934809235258932828727811020561028739221,

−0.9655906415488712273170132317689523963571242521298090207]

Para utilizarmos o nosso algoritmo (através do comando Pai), precisamos acrescentar oparâmetro Npoints_dif, que determina o número de níveis de diferenças que devemoscalcular. A rotina Pai chama (internamente) as rotinas Pcost e Dif que fazem a ‘costura’e calculam as diferenças (respectivamente).

[> pontosAI := Pai(sysdvd,ic,passo,Npoints,Npoints_dif):

Para ver o último ponto podemos usar

[> pontosAI[80];

49

[0.079, 0.9958778706064519485208857839771050926504558431325745645,

−0.07414440924614320256849352885361459284454270447896652664,

−0.08846938823439947934809235258963359126537211094725698514,

−.9655906415488712273170132317445272670566540907309817945]

Para fazer uma análise do desempenho dos integradores podemos usar o comando Pan.Se quisermos comparar o tempo gasto pelos integradores e fazer uma análise da precisãoem uma das coordenadas (por exemplo a coordenada x – coordenada 1) usamos:

[> Nheuristic := 50:[> Nremember := 0:[> pop := [1,2]:[> Pan(sysdvd,ic,passo,Npoints,Npoints_dif,pop);

Time spent (in miliseconds) :

Time RK78 =, 151, T ime AI =, 72, T imeRK4 =, 24

—————————–

Number of accurate digits (Nad) :

Nad RK78 =, 29, Nad AI =, 29, Nad RK4 =, 14

Por fim, se quisermos uma análise mais completa podemos escrever

[> Nremember := 1:[> pop := [1,3]:[> Pan(sysdvd,ic,passo,Npoints,Npoints_dif,pop);



—————————–


Nad RK78 = [29, 30, 29, 28], Nad AI = [29, 28, 30, 28], Nad RK4 = [14, 15, 14, 14]

PointsTay[80][1] = 0.9958778706064519485208857839805919534750878504070506864

Pointsrk78[80][1] = 0.9958778706064519485208857839768777550194781749259410418

PointsAI[80][1] = 0.9958778706064519485208857839771050926504558431325745645

Pointsrk4[80][1] = 0.9958778706064539125304520815792525636620511943215558880

50

4.2 Uma Análise de Desempenho

Nesta seção vamos analisar o desempenho do algoritmo AI aplicado a dois sitemasdinâmicos:

• O sistema de Lorenz – um sistema caótico dissipativo.• O sistema de Hénon-Heiles – um sistema caótico conservativo (Hamilto-

niano).

4.2.1 Sistema de Lorenz

Já falamos brevemente sobre o sistema de Lorenz na seção 3.1.2. O sistema deLorenz foi introduzido por Edward Lorenz em 1963, que o derivou a partir das equaçõessimplificadas de rolos de convecção que ocorrem nas equações da atmosfera. É um mapacaótico que mostra como o estado de um sistema dinâmico evolui no tempo num padrãocomplexo, não-repetitivo e cuja forma é conhecida por se assemelhar a uma borboleta.

Trata-se de um sistema não-linear, tridimensional e determinístico que exibe com-portamento caótico e demonstra aquilo a que hoje se chama um atrator estranho. Asequações que governam o sistema de Lorenz são:

x = σ(y − x),

y = −y − xz +Rx, (83)

z = xy − bz.

em que σ se chama o número de Prandtl e R se chama o número de Rayleigh. Todos osparâmetros σ,R e b > 0, mas usualmente σ = 10, b = 8/3, enquanto R varia. O sistemaexibe comportamento caótico para R > 28 mas tem órbitas periódicas para outros valoresde R. Usaremos a configuração

[> sigma := 10.0: R:=28.0: b:= 8/3:[> passo:=0.001:[> syst_dvd := [dvd[1]=10*vd[2]-10*vd[1],dvd[2]=-vd[1]*vd[3]+

28*vd[1]-vd[2],dvd[3]=vd[1]*vd[2]-8/3*vd[3]]:

para fazer algumas análises:

4.2.1.1 Lorenz – Análise 1

[> Npoints:=37: Npoints_dif := 17: Npoints_cost := Npoints_dif+1:

51

[> Digits := 55:[> ic := [5,5,5]:[> Nheuristic := 50: Nremember := 0:

[> pop := [1,2]:[> Pan(sysdvd,ic,passo,Npoints,Npoints_dif,pop);



—————————–



Quanto às outras coordenadas, temos




—————————–


Nad RK78 = [21, 20, 22], Nad AI = [21, 20, 21], Nad RK4 = [9, 8, 9]

PointsTay[37][1] = 5.635697329741590793410533085232013836758578119077257351

Pointsrk78[37][1] = 5.635697329741590793410395908246223183156290746748275911

PointsAI[37][1] = 5.635697329741590793410064396074898331210980530401826086

Pointsrk4[37][1] = 5.635697330694193316461462394186049286259678541533919190

4.2.1.2 Lorenz – Análise 2

[> Npoints := 74: Npoints_dif := 27: Npoints_cost := Npoints_dif+1:[> Digits := 50:[> ic := [5,5,5]:[> Nheuristic := 50: Nremember := 0:

[> pop := [1,2]:[> Pan(sysdvd,ic,passo,Npoints,Npoints_dif,pop);

52



—————————–



Outras coordenadas:




—————————–



PointsTay[73][1] = 7.359435594101039328467573045347355735002079257041890945

Pointsrk78[73][1] = 7.359435594101039328468751508117837528422704235205877825

PointsAI[73][1] = 7.359435594101039328470234415302385251586701759528755262

Pointsrk4[73][1] = 7.359435594963909052588038640434925917161406151659933287

4.2.1.3 Lorenz – Análise 3: Uma Análise mais Completa

As análises 1 e 2 mostradas exemplificam o uso / desempenho do nosso algoritmo.Para um estudo mais detalhado, precisamos examinar as correlações entre os temposgastos / número de dígitos significativos e as configurações dos parâmetros Npoints_dife Digits. Para uma melhor apreciação dos resultados encontrados veja a tabela 1 com asinformações principais.

Podemos perceber que o algoritmo AI gasta um pouco menos de tempo que oRK78para conseguir o mesmo nível de precisão. Já no caso de uma precisão intermediária entreas precisões apresentadas pelo RK78 e RK4, o tempo gasto é também intermediário(entre os tempos gastos pelos integradores RK78 e RK4).

Podemos também mudar a condição inicial e observar se o padrão se mantém.Em outras palavras, podemos tentar responder à questão:

53

Tabela 1 - Desempenho do algoritmo AI no sistema de Lorenz.

Time spent (ms) Number of acurate digitsNpts N_dif Digits RK78 AI RK4 RK78 AI RK433 17 55 45 38 8 [21, 20, 22] [21, 20, 22] [9, 8, 9]33 12 55 45 30 7 [21, 20, 22] [18, 17, 18] [9, 8, 9]33 7 55 44 22 7 [21, 20, 22] [13, 13, 13] [9, 8, 9]67 29 55 93 81 16 [21, 20, 22] [21, 20, 22] [9, 8, 10]67 22 55 91 78 16 [21, 20, 22] [19, 19, 18] [9, 8, 10]67 15 55 90 56 15 [21, 20, 22] [14, 14, 14] [9, 8, 10]100 37 63 138 126 23 [20, 21, 20] [20, 20, 22] [9, 8, 9]100 28 60 134 101 22 [20, 21, 20] [16, 16, 16] [9, 8, 9]100 19 60 135 79 22 [20, 21, 20] [13, 13, 12] [9, 8, 9]

Legenda: Estudo comparativo no sistema de Lorenz. Npts e N_dif significam Npoints eNpoints_dif, respectivamente.

Fonte: O autor, 2016.

Questão 4.1 Se recomeçarmos a integração em um ponto posterior à condição inicial (emuma dada trajetória do sistema) mantendo os outros parâmetros, a precisão do integradorAI vai se manter estável, assim como acontece com os integradores RK78 e RK4?

Para o sistema de Lorenz, podemos encontrar um começo de resposta para a ques-tão 4.1. Na análise que se segue vamos omitir os dados de saída com os valores dascoordenadas (apenas para não colocar dados desnecessários nesta análise em particular):

[> Npoints := 37: Npoints_dif := 17: Npoints_cost := Npoints_dif+1:[> Digits := 55:[> ic := [5,5,5]:[> Nheuristic := 50: Nremember := 0:[> pop := [1,3]:[> Pan(sysdvd,ic,passo,Npoints,Npoints_dif,pop):



—————————–



Agora, vamos tomar como condições iniciais as coordenadas dos pontos Pointsrk78[37] ePointsAI[37] e refazer uma integração de 37 pontos a partir delas:

54

[> Npoints := 37: Npoints_dif := 17: Npoints_cost := Npoints_dif+1:[> Digits := 55:[> ic1 := [op([2..4],Pointsrk78[Npoints])]:[> Nheuristic := 50: Nremember := 0:[> pop := [1,3]:[> Pan(sysdvd,ic1,passo,Npoints,Npoints_dif,pop):



—————————–



[> Npoints := 37: Npoints_dif := 17: Npoints_cost := Npoints_dif+1:[> Digits := 55:[> ic2 := [op([2..4],PointsAI[Npoints])]:[> Nheuristic := 50: Nremember := 0:[> pop := [1,3]:[> Pan(sysdvd,ic2,passo,Npoints,Npoints_dif,pop):



—————————–



O que nos leva a supor que o integrador AI é estavel no nível de equivalência com ointegrador RK78. Resta fazer uma bateria de testes para verificar se este resultado seestende para outros níveis de precisão, para um conjunto denso de trajetórias, para quenúmero de re-integrações, para uma outra configuração, para outros sistemas dissipativosetc.

4.2.2 Sistema de Hénon-Heiles

Enquanto em Princeton, em 1962, Michel Hénon e Carl Heiles trabalharam no mo-vimento não-linear de uma estrela em torno de um centro galáctico onde o movimento é

55

restrito a um plano. Em 1964, eles publicaram um artigo (HÉNON; HEILES, 1964) inti-tulado “A aplicabilidade da terceira integral de movimento: Alguns experimentos numéri-cos”. Sua ideia original era encontrar uma terceira integral de movimento numa dinâmicagaláctica. Para isso, eles tomaram um potencial axissimétrico (não linear) bidimensionalsimplificado e descobriram que a terceira integral existe apenas para um número limitadode condições iniciais. Na perspectiva moderna as trajetórias cujas condições iniciais nãoapresentam a terceira integral de movimento são chamadas órbitas caóticas.

O potencial de Hénon-Heiles é expresso por

V (x, y) = 12(x2 + y2) + λ(x2y − y3

3 ). (84)

e, assim, a Hamiltoniana de Hénon-Heiles pode ser escrita como

H = 12m(p2

x + p2y) + k

2(x2 + y2) + λ (x2y − y3

3 ), (85)

que resultam (após uma mudança de escala) nas quatro equações que definem o sistemade Hénon-Heiles:

x = Az, (86)

y = w, (87)

z = −(A+ 2 k

Ay

)x, (88)

w = −k x2

A+ y + σ y2. (89)

Usaremos a configuração

[> sigma:=5.000000: k:=0.999900: A:= sqrt(2*k/sigma +1):[> passo:=0.001:[> syst_dvd := [dvd[1]=A*vd[3],dvd[2]=vd[4],dvd[3]=-(A+2*k*vd[2]/A)*vd[1],

dvd[4]=-k*vd[1]^2/A + vd[2] + sigma*vd[2]^2]:

para nossa análise.

4.2.2.1 Henon-Heiles – Análise 1

[> Npoints := 140: Npoints_dif := 17: Npoints_cost := Npoints_dif+1:[> Digits := 60:[> ic := [1.000093761322944,-0.0006000959740288238,

-0.000040043010069843160,-0.8968057300859730197438631]:[> Nheuristic := 50: Nremember := 0:

56

[> pop := [1,3]:[> Pan(sysdvd,ic,passo,Npoints,Npoints_dif,pop):



—————————–


Nad RK78 =, [29, 29, 29, 28], Nad AI =, [29, 29, 29, 28], Nad RK4 =, [14, 15, 14, 13]

Vemos que, no caso do sistema de Hénon-Heiles, o integrador AI tem um desempe-nho muito bom: o tempo gasto é menos da metade do tempo gasto pelo integrador RK78(para a mesma precisão). Vamos analisar, na próxima seção, o desempenho do nossointegrador (tal como no caso do sistema de Lorenz) para vários valores dos parâmetros deconfiguração.

4.2.2.2 Henon-Heiles – Análise 2: Uma Análise mais Completa

Tal como no caso do sistema de Lorenz, a análise 1 nos serviu para termos umaideia da relação precisão / tempo gasto no uso do nosso algoritmo para a integraçãodo sistema de Hénon-Heiles. Para um estudo mais detalhado, vamos examinar (comofizemos para o sistema de Lorenz) as correlações entre os tempos gastos / número dedígitos significativos e as configurações dos parâmetros Npoints_dif e Digits. Vejam atabela 2 com os detalhes da análise.

Podemos perceber que, para o sistema de Hénon-Heiles, o algoritmo AI gasta, emgeral, menos da metade do tempo gasto pelo RK78 (no mesmo nível de precisão).

Podemos examinar, como no caso do sistema de Lorenz, se o nível de precisão semantém ao mudarmos a condição inicial, ou seja, continuar respondendo à questão 4.1:

[> Npoints := 100: Npoints_dif := 15: Npoints_cost := Npoints_dif+1:[> Digits := 60:[> ic := [1.000093761322944,-0.0006000959740288238,

-0.000040043010069843160,-0.8968057300859730197438631]:[> Nheuristic := 50: Nremember := 0:[> pop := [1,3]:[> Pan(sysdvd,ic,passo,Npoints,Npoints_dif,pop):

57

Tabela 2 - Desempenho do algoritmo AI para sistema Hénon-Heiles.

Time spent (ms) Number of acurate digitspts dif dig RK78 AI RK4 RK78 AI RK4120 16 60 228 108 38 [29, 29, 29, 28] [29, 29, 29, 28] [14, 15, 14, 13]185 19 60 355 169 57 [28, 29, 29, 28] [28, 29, 29, 28] [14, 14, 14, 13]210 21 65 420 199 68 [28, 29, 29, 28] [28, 29, 29, 28] [14, 14, 14, 13]210 21 60 430 202 70 [28, 29, 29, 28] [25, 27, 27, 26] [14, 14, 14, 13]220 21 70 431 203 69 [28, 29, 28, 28] [28, 29, 28, 28] [14, 14, 14, 14]230 21 75 444 216 71 [28, 29, 28, 28] [28, 29, 28, 28] [14, 14, 14, 14]320 25 90 605 306 92 [28, 28, 28, 27] [28, 28, 29, 27] [13, 14, 14, 13]320 16 90 608 240 92 [28, 28, 28, 27] [23, 21, 22, 22] [13, 14, 14, 13]320 19 90 611 262 94 [28, 28, 28, 27] [25, 24, 24, 24] [13, 14, 14, 13]320 22 90 612 289 99 [28, 28, 28, 27] [28, 26, 26, 26] [13, 14, 14, 13]

Legenda: Estudo comparativo no sistema de Hénon-Heiles. pts e dif significam Npoints eNpoints_dif, respectivamente.

Fonte: O autor, 2016.



—————————–



Agora, vamos tomar como condições iniciais as coordenadas dos pontos Pointsrk78[100]e PointsAI[100] e refazer uma integração de 100 pontos a partir delas:


58



—————————–






—————————–



O que confirma a estabilidade do integrador AI no nível de precisão do integrador RK78.

4.2.3 Algumas Observações

Nesta seção vamos reunir observações sobre as características do nosso integrador.Embora a maioria delas já deva ter sido notada pelo leitor mais atento, nada como umarepetição do óbvio para o assentamento dos novos conceitos e para a geração de novasperguntas / ideias que levem a uma melhor adequação / aplicação desse tipo de algoritmo.

Observação 4.1

• A primeira observação é que todos os resultados numéricos obtidos queenvolvem o tempo de integração são, na realidade, ‘semi’- quantitativos,uma vez que o pacote AssISt (criado para estudar o nosso integrador)é implementado em Maple – Uma linguagem de computação simbólica.Para uma comparação quantitativa mais acurada dos tempos de integra-ção, precisaríamos implementar o nosso algoritmo em uma liguagem de‘baixo nível’ como o C.

59

• Outra observação interessante é que não esperávamos obter a mesma pre-cisão do melhor integrador (o RK78) para nenhuma configuração dosparâmetros, uma vez que a base inicial é o próprio RK78.• Também não esperávamos, para uma configuração que reproduz (por que?)a mesma precisão do RK78, que o tempo de integração do AI fosse me-nor.• O item anterior nos permite pensar em usar o AI como integrador nolugar do RK78, uma vez que o tempo de integração é menor (para amesma precisão).• Os resultados apresentados nas tabelas mostradas nas seções anterioresnos sugerem que podemos, com uma simples escolha de parâmetros, cons-truir integradores de precisão intermediária.

60

CONCLUSÃO

Como dissemos na Introdução deste texto, nossa intenção aqui foi desenvolver umconjunto de ideias que possibilitariam uma integração mais eficiente, pelo menos em umagama de problemas.

Em que sentido “mais eficiente”? Para melhor definir o termo, poderíamos levantaralguns considerações.

Suponha que tenhamos um método de integração numérica X e que este leve umtempo TX para integrar um determinado sistema, dentro de um intervalo predefinido, comuma precisão PX . E que um outro método Y leve um tempo TY para executar o mesmotrabalho, com uma precisão PY . Supondo também que TX > TY e que PX > PY .

De uma maneira simples, chamaremos de um método mais eficiente de integraçãonumérica M aquele que ou alcançar uma precisão PM maior (ou igual) do que PY empre-gando um tempo TM menor do que TY ou alcançar uma precisão PM maior (ou igual) doque PX empregando um tempo TM menor do que TX .

Nós usaremos também uma outra (menos óbvia) definição de “mais eficiente”.Supondo que nosso método M apresente uma precisão que satisfaça:

PX > PM > PY

e tempos de integração que obedeçam a

TX > TM > TY ,

diremos que o métodoM é uma maneira “mais eficiente” de lidar com o problema numéricoem questão.

Por que a inclusão do terceiro e último item acima?Como tentamos deixar claro ao longo da Dissertação, a integração numérica de sistemas deequações diferenciais lida com uma constante batalha entre precisão e tempo de integração(resumindo, uma das questões mais importantes). Logo, usando a terminologia acimaintroduzida, podemos não “precisar” da acurácia PX mas precisarmos de uma precisãomaior do que PY . Logo, nossa abordagem se prova uma ferramenta importante nestecenário: ⇒ Podemos “controlar” a procura da precisão obtida e teremos (porconstrução) a princípio PY ≤ PM ≤ PX . E, como vimos no corpo da Dissertação, esteresultado foi obtido com TY ≤ TM ≤ TX .

Na verdade, conseguimos, para os dois exemplos utilizados, obter a situação naqual temos: TM < TX com PM = PX .

O que seria, na verdade, uma melhoria em um dos dois primeiros critérios quemencionamos acima, no começo da Conclusão. Mas acredito que deixamos claro o que

61

queremos dizer, podemos “mexer” na relação precisão/tempo de processamento commuito mais facilidade do que antes e obter soluções (mesmo quanto não conseguirmosPM = PX) “melhoradas” no (terceiro) sentido que mencionamos acima.

Esperamos ter deixado claro que nosso método proposto é uma melhoria, no sentidoque explicamos, sobre a integração numérica utilizando o Método de Runge-Kutta. Poristo demos nome Integração associativa em métodos de Runge-Kutta. Tentamosdeixar claro as ideias por detrás desta melhoria e ter, convincentemente, demonstradoos efeitos (reais) que uma implementação destes algoritmos e conceitos tem em umaintegração de fato (no caso, de dois sistemas bem conhecidos: Lorenz e Hénon-Heiles).

Na verdade, acreditamos que o mesmo conjunto de ideias pode ser aplicado amuitos outros pares de integradores, não apenas aos dois integradores de Runge-Kuttautilizados aqui. Isto será estudado no futuro.

Outro passo óbvio a ser tomado é realizar uma implementação destes algoritmosem conjunção com plataformas numéricas (por exemplo C++) para uma efetiva utilizaçãona prática, já que a eficácia da abordagem ficou verificada na implementação em Maple.

Acima, mencionamos dois possíveis caminhos para uma extensão do presente tra-balho. Existem outros, um pouco mais técnicos.

1. Não criamos ainda um algoritmo de análise mais completa de um dado sistema deequações diferenciais ordinárias, como costumamos fazer. Esta análise serve paraestimarmos, com maior precisão, os parâmetros a serem usados para este específicosistema, em que partes deste específico sistema, etc.

2. Durante o processo de cálculo de diferenças, podemos melhorar o algoritmo.

Resumindo, achamos que nosso método se mostrou bem sucedido e gera muitaslinhas de futuro trabalho.

62

REFERÊNCIAS

BARBOSA, L. M. C. R. et al. Improving the global fitting method on nonlinear timeseries analysis. Physical Review E, [S.l.], v. 74, n. 2, p. 026702, Aug. 2006.

BOYCE, W. E.; DIPRIMA, R. C. Equaceções Diferenciais Elementares e Problemas deValores de Contorno. 10. ed. Rio de Janeiro: LTC, 2015.

CARLI, H.; DUARTE, L.G.S.; MOTA, L.A.C.P. A maple package for improved globalmapping forecast. Computer Physics Communications, [S.l.], v. 185, n. 3, p. 1115 – 1129,Mar. 2014.

CELLIER, F. E.; KOFMAN, E. Continuous System Simulation. USA: Springer, 2006.

GARCIA, T. M. C. Solução numérica de equações diferenciais ordinárias. 2012. 114 f.Monografia (Especialização em Ciências) — Programa de Pós-graduação em Matemática,Universidade Tecnológica Federal do Paraná, Campo Mourão, 2012.

HAIRER, E.; LUBICH, C.; WANNER, G. Geometric Numerical Integration: Structurepreserving algorithms for ordinary differential equations. 2. ed. Berlin: Springer, 2006.

HÉNON, M.; HEILES, C. The applicability of the third integral of motion: Somenumerical experiments. The Astrophysical Journal, [S.l.], v. 69, n. 1, p. 73–79, Feb. 1964.

LORENZ, E. N. Deterministic nonperiodic flow. Journal of the Atmospheric Sciences,[S.l.], v. 20, n. 2, p. 130–141, Jan. 1963.

PRESS, W. H. et al. Numerical Recipes in C. 2. ed. Cambridge, New York: CambridgeUniversity Press, 1992.

ROMA, A. M.; NÓS, R. M. Tratamento Numérico de Equações Diferenciais. São Paulo:USP, 2012. Disponível em: <https://linux.ime.usp.br/~daniloss/antes-2012/tnedo.pdf>.Acesso em: 15 nov. 2015.

https://linux.ime.usp.br/~daniloss/antes-2012/tnedo.pdf

63

APÊNDICE A – Pacote AssISt

Neste apêndice vamos apresentar a implementação propriamente dita, ou seja, alistagem dos programas.

A.1 Cabeçalho

############################################################################# These are routines to study our Associative Integrator#############################################################################print(______________________________________);print(‘AssISt Package Development - V.03.16‘);print(______________________________________);print(‘J Avellar, LACP da Mota & LGS Duarte‘);print(______________________________________);print(‘LO Pereira & AF Rocha‘);############################################################################ Index:#### PTaylor - A routine to integrate using the Taylor series.#### Prk4maple - A routine to integrate using the RK4 method.#### Prk78maple - A routine to integrate using the RK78 method.#### Pcost - The ’sew’ routine.#### Dif - The routine that makes the difference(s).#### Pai - The algorithm AI.#### Pan - A routine to analize the AI method.########################################################################

64

A.2 Rotinas:

A.2.1 PTaylor

PTaylor := proc(syst_dvd,ic,passo,Npoints)local i,nv,sys,xxx,t,mapa,dsys1,dsol1,listi,rhssys;

## Comment# - Entries:# syst_dvd - system of differential equations.## ic - initial conditions.## passo - the integration step.## Npoints - the number of points you want to performe.## - Features:# Integrates using Taylor method.## - Output:# A list with the points (list of lists).##

nv := nops(syst_dvd);rhssys := map(u->rhs(u),syst_dvd);mapa := [seq(vd[i]=xxx[i](t),i=1..nv)];for i to nv dosys[i] := diff(xxx[i](t),t)=subs(mapa,rhssys[i]);end do;

dsys1 := [seq(sys[i],i=1..nv),seq(xxx[i](0)=ic[i],i=1..nv)];dsol1 := dsolve(dsys1,numeric,method=taylorseries,

output=procedurelist):listi:=[seq(i*passo,i=0..Npoints)]:[seq(map(u->rhs(u),dsol1(i)),i in listi)]:

end proc:

65

A.2.2 Prk4maple

Prk4maple := proc(syst_dvd,ic,passo,Npoints)local i,j,k,x,y,yt,dy,dyt,dym,nv,h,hh,h6,Pointsrk4;

## Comment# - Entries:# syst_dvd - equations of the system.## ic - initial conditions.## passo - the integration step.## Npoints - the number of points you want to performe.## - Features:# Integrates using RK4 method.## - Output:# A list with the points (list of lists).##

nv := nops(syst_dvd);h := passo;hh := h*0.5;h6 := h/6.0;Pointsrk4 := [[0,op(ic)]];x := 0;for i to nv do y[i] := ic[i] end do;

for k to Npoints dofor i to nv do dy[i] := subs([seq(vd[j]=y[j],j=1..nv)],rhs(syst_dvd[i])) end do;for i to nv do yt[i] := y[i]+hh*dy[i] end do;for i to nv do dyt[i] := subs([seq(vd[j]=yt[j],j=1..nv)],rhs(syst_dvd[i])) end do;for i to nv do yt[i] := y[i]+hh*dyt[i] end do;for i to nv do dym[i] := subs([seq(vd[j]=yt[j],j=1..nv)],

66

rhs(syst_dvd[i])) end do;for i to nv do

yt[i] := y[i]+h*dym[i];dym[i] := dym[i] + dyt[i];

end do;for i to nv do dyt[i] := subs([seq(vd[j]=yt[j],j=1..nv)],rhs(syst_dvd[i])) end do;for i to nv do y[i] := y[i]+h6*(dy[i]+dyt[i]+2*dym[i]) end do;x := x + h;Pointsrk4 := [op(Pointsrk4),[x,seq(y[j],j=1..nv)]];

end do;

return(Pointsrk4);

end proc:

A.2.3 Prk78maple

Prk78maple := proc(syst_dvd, ic, passo, Npoints)

## Comment# - Entries:# syst_dvd - equations of the system.## ic - initial conditions.## passo - the integration step.## Npoints - the number of points you want to performe.## - Features:# Integrates using RK78 method.## - Output:# A list with the points (list of lists).##

local i,j,k,c1,c2,c3,c4,c5,c6,c7,c8,c9,c10,c11,c12,c13,b2_1,b3_1,b3_2,

67

b4_1,b4_2,b4_3,b5_1,b5_2,b5_3,b5_4,b6_1,b6_2,b6_3,b6_4,b6_5,b7_1,b7_2,b7_3,b7_4,b7_5,b7_6,b8_1,b8_2,b8_3,b8_4,b8_5,b8_6,b8_7,b9_1,b9_2,b9_3,b9_4,b9_5,b9_6,b9_7,b9_8,b10_1,b10_2,b10_3,b10_4,b10_5,b10_6,b10_7,b10_8,b10_9,b11_1,b11_2,b11_3,b11_4,b11_5,b11_6,b11_7,b11_8,b11_9,b11_10,b12_1,b12_2,b12_3,b12_4,b12_5,b12_6,b12_7,b12_8,b12_9,b12_10,b12_11,b13_1,b13_2,b13_3,b13_4,b13_5,b13_6,b13_7,b13_8,b13_9,b13_10,b13_11,b13_12, Pointsrk78,dydx,x,y,ytemp,nv,h,ak2,ak3,ak4,ak5,ak6,ak7,ak8,ak9,ak10,ak11,ak12,ak13;

c1 := 4.17474911415302462220859284685e-2;c2 := 0.0;c3 := 0.0;c4 := 0.0;c5 := 0.0;c6 := -5.54523286112393089615218946547e-2;c7 := 2.39312807201180097046747354249e-1;c8 := 7.0351066940344302305804641089e-1;c9 := -7.59759613814460929884487677085e-1;c10 := 6.60563030922286341461378594838e-1;c11 := 1.58187482510123335529614838601e-1;c12 := -2.38109538752862804471863555306e-1;c13 := 2.5e-1;

b2_1 := 5.55555555555555555555555555556e-2;b3_1 := 2.08333333333333333333333333333e-2;b3_2 := 6.25e-2;b4_1 := 3.125e-2;b4_2 := 0.0e+0;b4_3 := 9.375e-2;b5_1 := 3.125e-1;b5_2 := 0.0e+0;b5_3 := -1.171875e+0;b5_4 := 1.171875e+0;b6_1 := 3.75e-2;b6_2 := 0.0e+0;b6_3 := 0.0e+0;b6_4 := 1.875e-1;b6_5 := 1.5e-1;

68

b7_1 := 4.79101371111111111111111111111e-2;b7_2 := 0.0e+0;b7_3 := 0.0e+0;b7_4 := 1.12248712777777777777777777778e-1;b7_5 := -2.55056737777777777777777777778e-2;b7_6 := 1.28468238888888888888888888889e-2;b8_1 := 1.6917989787292281181431107136e-2;b8_2 := 0.0e+0;b8_3 := 0.0e+0;b8_4 := 3.87848278486043169526545744159e-1;b8_5 := 3.59773698515003278967008896348e-2;b8_6 := 1.96970214215666060156715256072e-1;b8_7 := -1.72713852340501838761392997002e-1;b9_1 := 6.90957533591923006485645489846e-2;b9_2 := 0.0e+0;b9_3 := 0.0e+0;b9_4 := -6.34247976728854151882807874972e-1;b9_5 := -1.61197575224604080366876923982e-1;b9_6 := 1.38650309458825255419866950133e-1;b9_7 := 9.4092861403575626972423968413e-1;b9_8 := 2.11636326481943981855372117132e-1;b10_1 := 1.83556996839045385489806023537e-1;b10_2 := 0.0e+0;b10_3 := 0.0e+0;b10_4 := -2.46876808431559245274431575997e+0;b10_5 := -2.91286887816300456388002572804e-1;b10_6 := -2.6473020233117375688439799466e-2;b10_7 := 2.84783876419280044916451825422e+0;b10_8 := 2.81387331469849792539403641827e-1;b10_9 := 1.23744899863314657627030212664e-1;b11_1 := -1.21542481739588805916051052503e+0;b11_2 := 0.0e+0;b11_3 := 0.0e+0;b11_4 := 1.66726086659457724322804132886e+1;b11_5 := 9.15741828416817960595718650451e-1;b11_6 := -6.05660580435747094755450554309e+0;b11_7 := -1.60035735941561781118417064101e+1;b11_8 := 1.4849303086297662557545391898e+1;b11_9 := -1.33715757352898493182930413962e+1;

69

b11_10 := 5.13418264817963793317325361166e+0;b12_1 := 2.58860916438264283815730932232e-1;b12_2 := 0.0e+0;b12_3 := 0.0e+0;b12_4 := -4.77448578548920511231011750971e+0;b12_5 := -4.3509301377703250944070041181e-1;b12_6 := -3.04948333207224150956051286631e+0;b12_7 := 5.57792003993609911742367663447e+0;b12_8 := 6.15583158986104009733868912669e+0;b12_9 := -5.06210458673693837007740643391e+0;b12_10 := 2.19392617318067906127491429047e+0;b12_11 := 1.34627998659334941535726237887e-1;b13_1 := 8.22427599626507477963168204773e-1;b13_2 := 0.0e+0;b13_3 := 0.0e+0;b13_4 := -1.16586732572776642839765530355e+1;b13_5 := -7.57622116690936195881116154088e-1;b13_6 := 7.13973588159581527978269282765e-1;b13_7 := 1.20757749868900567395661704486e+1;b13_8 := -2.12765911392040265639082085897e+0;b13_9 := 1.99016620704895541832807169835e+0;b13_10 := -2.34286471544040292660294691857e-1;b13_11 := 1.7589857770794226507310510589e-1;b13_12 := 0.0e+0;

nv:=nops(syst_dvd);h:=passo;Pointsrk78:=[[0,op(ic)]];x:=0;for i to nv do y[i]:= ic[i] end do;

for k to Npoints dofor i to nv do dydx[i]:=subs([seq(vd[j]=y[j], j=1..nv)],

rhs(syst_dvd[i])) end do;for i to nv do ytemp[i]:=y[i]+b2_1*h*dydx[i] end do;

for i to nv do ak2[i]:=subs([seq(vd[j]=ytemp[j], j=1..nv)],rhs(syst_dvd[i])) end do;

for i to nv do ytemp[i]:=y[i]+h*(b3_1*dydx[i]+b3_2*ak2[i]) end do;for i to nv do ak3[i]:=subs([seq(vd[j]=ytemp[j], j=1..nv)],

70

rhs(syst_dvd[i])) end do;for i to nv do ytemp[i]:=y[i]+h*(b4_1*dydx[i]+b4_2*ak2[i]+

b4_3*ak3[i]) end do;for i to nv do ak4[i]:=subs([seq(vd[j]=ytemp[j], j=1..nv)],

rhs(syst_dvd[i])) end do;for i to nv do ytemp[i]:=y[i]+h*(b5_1*dydx[i]+b5_2*ak2[i]+b5_3*ak3[i]+b5_4*ak4[i]) end do;

for i to nv do ak5[i]:=subs([seq(vd[j]=ytemp[j],j=1..nv)], rhs(syst_dvd[i])) end do;

for i to nv do ytemp[i]:=y[i]+h*(b6_1*dydx[i]+b6_2*ak2[i]+b6_3*ak3[i]+b6_4*ak4[i]+b6_5*ak5[i]) end do;


for i to nv do ytemp[i]:= y[i]+h*(b7_1*dydx[i]+b7_2*ak2[i]+b7_3*ak3[i]+b7_4*ak4[i]+b7_5*ak5[i]+b7_6*ak6[i]) end do;


for i to nv do ytemp[i]:= y[i]+h*(b8_1*dydx[i]+b8_2*ak2[i]+b8_3*ak3[i]+b8_4*ak4[i]+b8_5*ak5[i]+b8_6*ak6[i]+b8_7*ak7[i])end do;


for i to nv do ytemp[i]:= y[i]+h*(b9_1*dydx[i]+b9_2*ak2[i]+b9_3*ak3[i]+b9_4*ak4[i]+b9_5*ak5[i]+b9_6*ak6[i]+b9_7*ak7[i]+b9_8*ak8[i])end do;


for i to nv do ytemp[i]:= y[i]+h*(b10_1*dydx[i]+b10_2*ak2[i]+b10_3*ak3[i]+b10_4*ak4[i]+b10_5*ak5[i]+b10_6*ak6[i]+b10_7*ak7[i]+b10_8*ak8[i]+b10_9*ak9[i])end do;


for i to nv do ytemp[i]:= y[i]+h*(b11_1*dydx[i]+b11_2*ak2[i]+b11_3*ak3[i]+b11_4*ak4[i]+b11_5*ak5[i]+b11_6*ak6[i]+b11_7*ak7[i]+b11_8*ak8[i]+b11_9*ak9[i]+b11_10*ak10[i]) end do;


for i to nv do ytemp[i]:=y[i]+h*(b12_1*dydx[i]+b12_2*ak2[i]+b12_3*ak3[i]+b12_4*ak4[i]+b12_5*ak5[i]+b12_6*ak6[i]+b12_7*ak7[i]+b12_8*ak8[i]+b12_9*ak9[i]+b12_10*ak10[i]+b12_11*ak11[i]) end do;

71


for i to nv do ytemp[i]:= y[i]+h*(b13_1*dydx[i]+b13_2*ak2[i]+b13_3*ak3[i]+b13_4*ak4[i]+b13_5*ak5[i]+b13_6*ak6[i]+b13_7*ak7[i]+b13_8*ak8[i]+b13_9*ak9[i]+b13_10*ak10[i]+b13_11*ak11[i]+b13_12*ak12[i]) end do;for i to nv do ak13[i]:=subs([seq(vd[j]=ytemp[j],j=1..nv)], rhs(syst_dvd[i])) end do;

for i to nv do y[i]:=y[i]+h*(c1*dydx[i]+c2*ak2[i]+c3*ak3[i]+c4*ak4[i]+c5*ak5[i]+c6*ak6[i]+c7*ak7[i]+c8*ak8[i]+c9*ak9[i]+c10*ak10[i]+c11*ak11[i]+c12*ak12[i]+c13*ak13[i]) end do;

x:=x+h;Pointsrk78:=[op(Pointsrk78),[x, seq(y[j], j=1..nv)]];end do;

return (Pointsrk78);

end proc;

A.2.4 Pcost

Pcost := proc(syst_dvd,ListT,passo,Npoints_cost)local i,j,nv,ic,dsol,pont,listi,listT;

## Comment# - Entries:# syst_dvd - system of differential equations.## ListT - list with the points calculated with the better# integrator.## passo - the integration step.## Npoints_cost - the number of points you want to ’sew’.## - Features:# Performes the ’sew’ (for now using only the rk4 as backup# integrator).

72

## - Output:# A list with the points (list of lists).##

nv := nops(syst_dvd);listT := [op(1..Npoints_cost,ListT)];for j to Npoints_cost do

ic := [op([2 .. nv+1],op(j,listT))];dsol:=Prk4maple(syst_dvd,ic,passo,1);pont[j]:=[j*passo,op(2..nv+1,dsol[2])];

end do:

[ListT[1],seq(pont[j],j=1..Npoints_cost)];

end proc:

A.2.5 Dif

Dif := proc(syst_dvd,passo,Npoints,Npoints_dif,ListT,ListC)local i,j,k,n,p,m,h,ic,dd,d1,d2,result,listi,pont,dsol,LT,LC,nn,nv,vd;

## Comment# - Entries:# syst_dvd - system of differential equations.## passo - the integration step.## Npoints - the number of points you want to performe.## Npoints_dif - the number ’dif’ levels.## ListT - list Taylor (the preciser series).## ListC - list Sew (the sewed).## - Features:# Integrates using Our method.

73

## - Output:# A list with the points (list of lists).##

p := Npoints_dif;LT := [op(1..p+1,ListT)];LC := [op(1..p+2,ListC)];nv := nops(syst_dvd);vd := [1..(nv+1)];m := p+1;for i from 2 by 1 to m do dd[1,i]:=(op(vd,op(i,LT))-op(vd,op(i,LC))); end do:for k from 2 to p do assign([seq(dd[k,j]=dd[k-1,j]-dd[k-1,j-1],j=3..m)]); end do:result:=[op(vd,op(m+1,LC))+op(convert([seq([dd[i,m]],i=1..p)],‘+‘))];LT:=[op(LT),result];h:=m+1;ic := [op([2..nv+1],op(h,LT))];dsol:=Prk4maple(syst_dvd,ic,passo,1);pont[h]:=[h*passo,op([2..nv+1],dsol[2])];LC:=[op(LC),pont[h]];for i to p do d1[i]:=dd[i,m] end do;

for n from 2 to Npoints-p-1 dom := n+p;d2[1]:=(op(vd,op(m,LT))-op(vd,op(m,LC)));for k from 2 to p do d2[k]:=d2[k-1]-d1[k-1];d1[k-1]:=d2[k-1] end do:d1[p]:=d2[p];result:=[op(vd,op(m+1,LC))+op(convert([seq([d2[i]],i=1..p)],‘+‘))];LT:=[op(LT),result];h:=m+1;ic := [op([2..nv+1],op(h,LT))];dsol:=Prk4maple(syst_dvd,ic,passo,1);pont[h]:=[h*passo,op([2..nv+1],dsol[2])];LC:=[op(LC),pont[h]];

end do:LT;end proc:

74

A.2.6 Pai

Pai := proc(syst_dvd,ic,passo,Npoints,Npoints_dif)local Pointsrk78,Pointscost,LF;

## Comment# - Entries:# syst_dvd - system of differential equations.## passo - the integration step.## Npoints - the number of points you want to performe.## Npoints_dif - the number ’dif’ levels.## ListT - list Taylor (the preciser series).## ListC - list Sew (the sewed).## - Features:# Integrates using Our method.## - Output:# A list with the points (list of lists).##

Pointsrk78 := Prk78maple(syst_dvd,ic,passo,Npoints_dif+1);Pointscost := Pcost(syst_dvd,Pointsrk78,passo,Npoints_dif+1);LF := Dif(syst_dvd,passo,Npoints,Npoints_dif,Pointsrk78,Pointscost);

end proc:

A.2.7 Pan

Pan := proc(syst_dvd,ic,passo,Npoints,Npoints_dif,pop)local n1,n2,nv,conjv,pt,p78,p4,pai,t0,i,nh;global tTAY,tRK78,tRK4,tAI,Pointsrk78,Pointsrk4,PointsTay,PointsAI,DD,Nheuristic,Nremember;

75

## Comment# - Entries:# syst_dvd - system of differential equations.## ic - the initial condition.## passo - the integration step.## Npoints - the number of points you want to performe.## Npoints_dif - the number ’dif’ levels.## pop - print options.## - Features:# Performance comparison.## - Output: printings on the screen.##

nh := Nheuristic;if Nremember<>1 thent0 := time[real]();PointsTay := PTaylor(syst_dvd,ic,passo,Npoints);tTAY := trunc((time[real]()-t0)*1000);t0 := time[real]();for i to nh doPointsrk78 := Prk78maple(syst_dvd,ic,passo,Npoints):end do:tRK78 := trunc((time[real]()-t0)/nh*1000);t0 := time[real]();for i to nh doPointsAI := Pai(syst_dvd,ic,passo,Npoints,Npoints_dif):end do:tAI := trunc((time[real]()-t0)/nh*1000);t0 := time[real]();for i to nh doPointsrk4 := Prk4maple(syst_dvd,ic,passo,Npoints):

76

end do:tRK4 := trunc((time[real]()-t0)/nh*1000);end if;

n1 := pop[1];n2 := pop[2];nv := nops(syst_dvd);conjv := [seq(i,i=1..nv)];

pt := [op([2..nv+1],op(Npoints,PointsTay))];p78 := [op([2..nv+1],op(Npoints,Pointsrk78))];p4 := [op([2..nv+1],op(Npoints,Pointsrk4))];pai := [op([2..nv+1],op(Npoints,PointsAI))];

for i to nv doDD[i,intrk78] := abs(op(2,pt[i]-p78[i]))-length(op(1,pt[i]-p78[i]));DD[i,intrk4] := abs(op(2,pt[i]-p4[i]))-length(op(1,pt[i]-p4[i]));DD[i,intai] := abs(op(2,pt[i]-pai[i]))-length(op(1,pt[i]-pai[i]));

end do;

if not has([seq(i,i=1..nv)],n1) then print(‘pop[1] isthe coordinate number.‘);return end if;

if n2 = 1 thenprint(‘Time spent (in miliseconds):‘);print(‘Time RK78 = ‘,tRK78,‘ Time AI = ‘,tAI,‘Time RK4 = ‘,tRK4);elif n2=2 thenprint(‘Time spent (in miliseconds):‘);print(‘Time RK78 = ‘,tRK78,‘ Time AI = ‘,tAI,‘Time RK4 = ‘,tRK4);print(____________________);print(‘Number of accurate digits (Nad):‘);print(‘Nad RK78 = ‘,DD[n1,intrk78],‘ Nad AI = ‘,DD[n1,intai],‘ Nad RK4 = ‘,DD[n1,intrk4]);elif n2=3 then

77

print(‘Time spent (in miliseconds):‘);print(‘Time RK78 = ‘,tRK78,‘ Time AI = ‘,tAI,‘Time RK4 = ‘,tRK4);print(____________________);print(‘Number of accurate digits (Nad):‘);print(‘Nad RK78 = ‘,[seq(DD[i,intrk78],i=1..nv)],‘Nad AI = ‘,[seq(DD[i,intai],i=1..nv)],‘ Nad RK4 = ‘,[seq(DD[i,intrk4],i=1..nv)]);print(____________________);print(‘PointsTay[‘,Npoints,‘][‘,n1,‘] = ‘,pt[n1]);print(‘ Pointsrk78[‘,Npoints,‘][‘,n1,‘] = ‘,p78[n1]);print(‘ PointsAI[‘,Npoints,‘][‘,n1,‘] = ‘,pai[n1]);print(‘ Pointsrk4[‘,Npoints,‘][‘,n1,‘] = ‘,p4[n1]);else print(‘pop[2] must have one of the values: 1, 2, 3.‘);end if;

end proc:

Documents

Integração associativa em métodos Runge-Kutta