Estimação Bayesiana de Parâmetros envolvidos em Modelos Determinísticos - PDF
Estimação Bayesiana de Parâmetros envolvidos em Modelos Determinísticos
Download "Estimação Bayesiana de Parâmetros envolvidos em Modelos Determinísticos"
Isabela Sales Cordeiro
1 Estimação Bayesiana de Parâmetros envolvidos em Modelos Determinísticos por Josiane da Silva Cordeiro Universidade Federal do Rio de Janeiro Instituto de Matemática Departamento de Métodos Estatísticos 2009
2 Estimação Bayesiana de Parâmetros envolvidos em Modelos Determinísticos Josiane da Silva Cordeiro Dissertação submetida ao Corpo Docente do Instituto de Matemática - Departamento de Métodos Estatísticos da Universidade Federal do Rio de Janeiro - UFRJ, como parte dos requisitos necessários à obtenção do grau de Mestre em Estatística. Aprovada por: Prof a. Alexandra M. Schmidt. PhD - IM - UFRJ - Orientadora. Prof. Cláudio J. Struchiner PhD - PROCC - FIOCRUZ - Co-Orientador. Prof. Hélio S. Migon PhD - IM - UFRJ. Prof. Dirceu Silveira Reis Jr. PhD - DEHID - FUNCEME. ii
3 FICHA CATALOGRÁFICA Cordeiro, Josiane da Silva. Estimação Bayesiana de Parâmetros envolvidos em Modelos Determinísticos / Josiane da Silva Cordeiro. Rio de Janeiro: UFRJ, IM, DME, Dissertação - Universidade Federal do Rio de Janeiro, IM, DME. 1. Introdução. 2. Modelos Determinísticos. 3. Estimação Bayesiana de Parâmetros envolvidos em Modelos Determinísticos. 4. Estimação Bayesiana de Parâmetros envolvidos em Equações a Diferenças Finitas: Uma Aplicação à Modelagem de Chuva-Vazão. 5. Conclusões e Trabalhos Futuros. (Mestrado-UFRJ/IM/DME) I. Schmidt, Alexandra II. Universidade Federal do Rio de Janeiro III. Título. iii
4 iv Ao Criador, meu Pai Celestial.
5 Agradecimentos A Deus, pelo seu grande amor. Em particular, este amor tem sido revelado através do seu cuidado, das oportunidades surgidas e a clareza necessária para enxergá-las e, principalmente, o fortalecimento para prosseguir. Agradeço a Deus também, por cada pessoa que faz parte da minha vida. À minha família pelo amor. Aos meu pais, Josias e Maria, pelo cuidado e carinho demonstrados de forma muito peculiar, segundo o jeito de ser de cada um. Ao meu irmão Welington, pela paciência e por me dar a certeza de que fiz a escolha certa quanto ao meu curso de mestrado. E ao meu irmão Welton, também pela paciência e pelos vários momentos descontraídos. Ao Felipe, pelo amor e apoio. E também, por ter contribuído para que meu ano de dissertação fosse um ano mais feliz. Às minhas avós (in memorian), pelo amor que, enquanto possível, foi dado a mim. Aos meus tios, tias, primos e primas que de diversas maneiras me amaram e incentivaram ao longo da minha vida. Aos meus amigos, aos atualmente presentes e aos ausentes fisicamente, mas que de maneiras distintas me apoiaram muito. Aos meus companheiros e amigos do DME, que compartilharam comigo momentos de angústia, de satisfação e de alegria, nestes dois anos de intensa convivência. E, é claro, obrigada pelos grandes momentos no laboratório casamenteiro, o LPGE. Em especial, a minha turma Denise, Mariana, Nícia, Patrícia, Vera e Alexandre, e a Valmária, João, Targino e Vinícius. À minha orientadora Alexandra Schmidt, por todo incentivo desde quando ainda não v
6 era aluna da UFRJ. Agradeço pela dedicação como orientadora e, pelo compartilhamento de conhecimento e experiências científicas. Ao meu atual co-orientador Cláudio Struchiner, pela orientação acadêmica dada, desde o ensino médio, com liberaridade e, também pela transmissão de um pouco, que foi muito para mim, do seu conhecimento científico. Aos meus professores do DME-IM/UFRJ e do DEMAT/UFRuralRJ que, em diferentes etapas, foram essenciais para meu crescimento através do conhecimento passado e da constante motivação. Ao Hélio Migon e ao Dirceu Reis pela cordialidade e colaboração científica. A todos os vinculados a UFRJ, pela colaboração para promover um ambiente de estudo de qualidade. A todos os profissionais do PROCC/FIOCRUZ, pelo ambiente agradável e pela disponibilidade em diversos momentos. À FUNCEME pela concessão dos dados utilizados nesta dissertação. À CAPES e à FAPERJ, pela bolsa especial, que possibilitaram o prosseguimento dos meus estudos com mais entusiasmo. vi
7 Resumo Neste trabalho, temos interesse em investigar procedimentos de inferência bayesiana de parâmetros originalmente envolvidos em modelos determinísticos que descrevem problemas reais. A fim de acoplar incerteza ao estudo determinístico, assumimos que a solução determinística é a média de uma realização da distribuição normal. Desta forma, especificamos a função de verossimilhança e, sob o paradigma bayesiano, atribuímos uma distribuição a priori para os parâmetros, assim obtendo a distribuição a posteriori. Devido à complexidade destes modelos, surge a necessidade de métodos computacionais para obtermos amostras desta distribuição. Neste contexto, escolhemos algoritmos baseados em dois métodos específicos, a saber, o Monte Carlo via cadeias de Markov e o Monte Carlo sequencial. Particularmente, objetivamos comparar estes métodos para a estimação dos parâmetros relacionados a um modelo determinístico que descreve interações entre espécies de presas e de predadores. Os resultados atingidos via todos os algoritmos utilizados mostraram que tais procedimentos são razoáveis neste contexto de estimação, embora a descrição da incerteza, associada à distribuição a posteriori dos parâmetros, reportada por ambos os métodos tenham sido distintas. Outro interesse específico neste trabalho, assumindo um modelo determinístico da vazão de uma bacia hidrográfica, é ajustar o modelo estatístico proposto a conjuntos de dados de três bacias do estado do Ceará, cedidos pela FUNCEME. Neste caso, obtivemos um ajuste razoável do modelo aos diferentes conjuntos de dados atingindo distribuições a posteriori informativas para os parâmetros de interesse. Palavras Chaves: Modelos Determinísticos, Estimação Bayesiana, Monte Carlo via cadeias de Markov, Monte Carlo sequencial. vii
8 Abstract In this work, we are interested in investigating Bayesian inference procedures of parameters related to deterministic models which describe real problems. In order to account for uncertainty estimation into the deterministic study, we assume the deterministic solution as the mean of a normal distribution. In this way, we specify the likelihood function and, under the Bayesian paradigm, we assign a prior distribution to the parameters, hence obtaining a posterior distribution. Due to the complexity of these models, it emerges the need of computational methods to obtain samples of the resultant posterior distribution. In this context, we choose algorithms based on two specific methods, Markov chain Monte Carlo and sequential Monte Carlo algorithms. Particularly, we aim to compare these methods to estimate parameters involved in a deterministic model that describes interactions between prey species and predator species. The achieved results through all the algorithms showed that these procedures are reasonable in this context of estimation, although the description of uncertainty, associated with the posterior distribution of the parameters, reported by both methods have been distinct. Another specific aim of this work is to model the runoff of a hydrographic basin through a deterministic model. The proposed statistical model is fitted to datasets from three basins in the State of Ceará, made available by FUNCEME. In this case, we obtained a reasonable fit of the model to the different datasets reaching informative posterior distributions to the parameters of interest. Key Words: Deterministic Models, Bayesian Estimation, Markov Chain Monte Carlo, Sequential Monte Carlo. viii
9 Sumário 1 Introdução 1 2 Modelos Determinísticos Equações a Diferenças Finitas Equações Diferenciais Ordinárias Métodos Numéricos de Aproximação de Solução de EDO Método de Euler Método de Runge-Kutta Métodos de Passos Múltiplos Método da Colocação Calibração e Estimação de Parâmetros no contexto de Modelos Determinísticos Estimação Bayesiana de Parâmetros envolvidos em Modelos Determinísticos MCMC MCS CBA-MCS CBA-CRP Estudo Simulado: Uma Aplicação em um Modelo Clássico de EDO Função de Verossimilhança Geração dos Dados Artificiais Formulação Bayesiana ix
10 3.3.4 Implementação Computacional Resultados Estimação Bayesiana de Parâmetros envolvidos em Equações a Diferenças Finitas: Uma Aplicação à Modelagem de Chuva-Vazão Modelo SMAP Função de Verossimilhança Formulação Bayesiana Estudo Simulado Geração dos Dados Artificiais Resultados através do MCMC Dados Reais Bacia Bacia Bacia Conclusões e Trabalhos Futuros 68 x
11 Lista de Tabelas 3.1 Parâmetros do CBA-MCS ajustados SMAP/Bacia Estimativas pontual (mediana a posteriori) e intervalar (intervalo de credibilidade de 90%) para cada parâmetro SMAP/Bacia Estimativas pontual (mediana a posteriori) e intervalar (intervalo de credibilidade de 90%) para cada parâmetro SMAP/Bacia Estimativas pontual (mediana a posteriori) e intervalar (intervalo de credibilidade de 90%) para cada parâmetro xi
12 Lista de Figuras 3.1 Traços das cadeias dos parâmetros α e β (colunas), para as diferentes observações com N=16, N=51 e N=151 (linhas), e supondo a priori 3 via MCMC. A linha horizontal preta é o valor verdadeiro do parâmetro Mediana (ponto cheio) e intervalos de credibilidade de 90% da distribuição a posteriori de α e β supondo a priori 1 (P1), priori 2 (P2) e priori 3 (P3), tanto via MCMC (em azul) quanto via CBA-MCS (em vermelho), com N=16, N=51 e N=151 (linhas) I.C. de 90% (linhas tracejadas) e a mediana (linhas cheias) da distribuição preditiva das populações de presas e predadores, obtidos a partir dos valores das cadeias dos parâmetros via MCMC (em azul) e via CBA-MCS (em vermelho), assumindo a priori 3, e para as diferentes amostras com N=16, N=51 e N=151 (linhas). Os círculos representam os dados observados Gráficos de Dispersão entre α e β, no caso da priori 3, via MCMC (primeira coluna) e CBA-MCS (segunda coluna), e para os diferentes conjuntos de dados gerados N=16, N=51 e N=151 (linhas 1, 2 e 3, respectivamente) Mediana (ponto cheio) e intervalos de credibilidade de 90% da distribuição a posteriori de α e de β via CBA-MCS (em azul) e CBA-MCS com B = 10 (em vermelho), a partir das amostras de tamanhos N=16 (N16), N=51 (N51) e N=151 (N151) xii
13 3.6 I.C. de 90% (linha tracejada) e a mediana (linha cheia) da distribuição preditiva das populações de presas e de predadores, obtidos a partir dos valores das cadeias dos parâmetros via CBA-MCS, com B = 10, assumindo a priori 3, e para as diferentes amostras com N=16, N=51 e N=151 (linhas). Os círculos representam os dados observados Mediana (ponto cheio) e intervalos de credibilidade de 90% da distribuição a posteriori de α e de β via CBA-MCS (em azul) e CBA-CRP (em vermelho), a partir das amostras de tamanhos N=16 (N16), N=51 (N51) e N=151 (N151) I.C. de 90% (linha tracejada) e a mediana (linha cheia) da distribuição preditiva das populações de presas e de predadores, obtidos a partir dos valores das cadeias dos parâmetros via CBA-CRP assumindo a priori 3. Os círculos representam os dados observados com N=16, N=51 e N=151 (linhas) Mediana (ponto cheio) e intervalos de credibilidade de 90% da distribuição a posteriori de α e de β via CBA-MCS (em azul) e CBA-CRP com o núcleo de transição para trás sub-ótimo (em vermelho), a partir das amostras de tamanhos N=16 (N16), N=51 (N51) e N=151 (N151) I.C. de 90% (linha tracejada) e a mediana (linha cheia) da distribuição preditiva das populações de presas e de predadores, obtidos a partir dos valores das cadeias dos parâmetros via CBA-CRP, com o núcleo de transição para trás sub-ótimo, assumindo a priori 3. Os círculos representam os dados observados com N=16, N=51 e N=151 (linhas) Diagrama do modelo SMAP SMAP/Dados Artificiais - Traços das cadeias dos parâmetros (via formulação 1). A linha horizontal preta é o valor verdadeiro do parâmetro SMAP/Dados Artificiais - Traços das cadeias dos parâmetros (via formulação 2). A linha horizontal preta é o valor verdadeiro do parâmetro.. 48 xiii
14 4.4 SMAP/Dados Artificiais - Histogramas dos parâmetros (via formulação 1). A linha vertical vermelha é o valor verdadeiro do parâmetro SMAP/Dados Artificiais - Histogramas dos parâmetros (via formulação 2). A linha vertical vermelha é o valor verdadeiro do parâmetro SMAP/Dados Artificiais - I.C. de 90% (linha tracejada) e a mediana (linha cheia em azul) da distribuição preditiva da vazão obtidos a partir dos valores das cadeias dos parâmetros via MCMC e assumindo a formulação 1. Os círculos representam os dados observados SMAP/Dados Artificiais - I.C. de 90% (linha tracejada) e a mediana (linha cheia em azul) da distribuição preditiva da vazão obtidos a partir dos valores das cadeias dos parâmetros via MCMC e assumindo a formulação 2. Os círculos representam os dados observados SMAP/Bacia Traços das cadeias dos parâmetros SMAP/Bacia Histogramas dos parâmetros SMAP/Bacia Mediana (linha tracejada preta) da distribuição preditiva da vazão obtidos a partir dos valores das cadeias dos parâmetros via MCMC e os dados observados da bacia (linha cheia azul) SMAP/Bacia Quantis de 25% e de 75% da distribuição preditiva da vazão a partir das cadeias dos parâmetros, e os dados observados da bacia (linha cheia azul) SMAP/Bacia Histogramas dos parâmetros SMAP/Bacia Mediana (linha tracejada preta) da distribuição preditiva da vazão obtidos a partir dos valores das cadeias dos parâmetros via MCMC e os dados observados da bacia (linha cheia azul) SMAP/Bacia Quantis de 25% e de 75% da distribuição preditiva da vazão a partir das cadeias dos parâmetros, e os dados observados da bacia (linha cheia azul) SMAP/Bacia Histogramas dos parâmetros xiv
15 4.16 SMAP/Bacia Mediana (linha tracejada preta) da distribuição preditiva da vazão obtidos a partir dos valores das cadeias dos parâmetros via MCMC e os dados observados da bacia (linha cheia azul) SMAP/Bacia Quantis de 25% e de 75% da distribuição preditiva da vazão a partir das cadeias dos parâmetros, e os dados observados da bacia (linha cheia azul) xv
16 Capítulo 1 Introdução Modelos matemáticos determinísticos são um conjunto de equações ou inequações matemáticas, organizadas de forma que, conhecidas algumas condições deste sistema, é possível obter sua solução em um dado momento. Estes modelos são utilizados em diversas áreas de pesquisas científicas, tais como biologia, química e física. Equações do tipo diferenciais (ordinárias, parciais ou com retardo), integro-diferenciais, e a diferenças finitas, têm sido amplamente utilizadas nestes contextos. Em modelagem matemática, podemos também estar interessados na calibração dos parâmetros que definem o modelo determinístico adotado. A calibração dos parâmetros de um modelo determinístico consiste em buscar um conjunto de valores para estes parâmetros até que se alcance uma melhor representação do problema de estudo. Desta forma, ajustamos o modelo determinístico aos dados em relação aos parâmetros que o definem. O crescente uso da modelagem determinística e o aumento da complexidade dos modelos, para descrever problemas de diversas áreas, vêm acentuando a necessidade de métodos de calibração mais robustos neste contexto. Particularmente, no caso de parâmetros que definem equações diferenciais ordinárias não-lineares, a calibração pode ser feita linearizando-se estas estruturas. Neste contexto, o método dos mínimos quadrados não-lineares têm sido amplamente aplicado. Por outro lado, diversos procedimentos estatísticos mais elaborados, clássicos e bayesianos, para estimar parâmetros envolvidos em modelos determinísticos têm sido propos- 1
17 tos na literatura. Veja, por exemplo, Poole e Raftery (2000); Cancré et al. (2000); Huang et al. (2006); Ramsay et al. (2007); Campbell (2007); Toni et al. (2008). Neste trabalho, temos interesse em revisar e investigar métodos de estimação bayesiana em modelos originalmente expressos de forma determinística. Neste contexto, adicionamos incerteza ao estudo determinístico. Desta forma, assumimos que a solução determinística entra como, por exemplo, a média de uma realização de uma distribuição de probabilidade bem definida. Além disto, seguindo o paradigma de Bayes, obtemos como resultado da análise uma distribuição de probabilidade a posteriori para os parâmetros de interesse. Dessa forma a incerteza inerente ao procedimento de estimação e a informação observada é naturalmente descrita. Em diversos casos, devido a complexidade dos modelos, a distribuição a posteriori resultante não tem forma analítica fechada. Portanto, surge a necessidade de métodos eficazes para a geração de amostras desta distribuição. Em particular, selecionamos métodos de amostragem baseados em Monte Carlo via cadeias de Markov e em Monte Carlo sequencial, sob a perspectiva da estatística bayesiana. Nosso objetivo específico é comparar a eficiência destes dois métodos em diversos cenários. A seguir, descrevemos a organização deste texto. No capítulo 2, fazemos uma revisão de literatura sobre modelos determinísticos. Em particular, apresentamos conceitos e resultados matemáticos relacionados as duas classes de modelos determinísticos estudados aqui, equações a diferenças finitas e equações diferenciais ordinárias. Além disto, discutimos alguns dos métodos mais conhecidos para aproximar soluções de equações diferenciais ordinárias. Finalmente, na seção 2.4, fornecemos uma visão geral sobre os métodos de calibração de parâmetros pertencentes a modelos determinísticos, propostos na literatura. No capítulo 3 concentra-se o foco deste trabalho, que são os métodos computacionais bayesianos para estimar parâmetros envolvidos em modelos determinísticos. Primeiro, discutimos brevemente o que tem sido proposto na literatura. Em seguida, nas seções 3.1 e 3.2, apresentamos os métodos bayesianos de amostragem que temos interesse em investigar, inclusive descrevemos os algoritmos para implementá-los. Mais especificamente, elegemos os métodos de Monte Carlo via cadeias de Markov e o de Monte Carlo 2
18 sequencial, adaptados ao contexto do nosso estudo. Na seção 3.3, utilizamos os métodos discutidos num estudo simulado. Escolhemos um modelo especificado por um sistema de duas equações diferenciais ordinárias não-lineares, que visa descrever dinâmicas populacionais entre uma espécie de presas e seus predadores. Já no capítulo 4, aplicamos a metodologia estudada, através de um dos métodos de estimação de interesse, a um modelo hidrológico descrito por um sistema de equações a diferenças finitas no tempo, proposto por Lopes et al. (1981). Este modelo tem sido amplamente utilizado por especialistas em hidrologia, e é conhecido como Soil Moisture Accounting Procedure (SMAP). Ajustamos o modelo estatístico proposto a um conjunto de dados artificialmente gerados. Além disto, três conjuntos de dados coletados no estado do Ceará, cedidos pelo grupo de pesquisa da Fundação Cearense de Meteorologia e Recursos Hídricos (FUNCEME), também foram ajustados. Estes conjuntos de dados diferem principalmente em relação à quantidade de observações que compõem suas séries temporais. Finalmente, no capítulo 5, concluímos nosso estudo apresentando uma discussão sobre os métodos e resultados obtidos. Além disto, indicamos possíveis extensões deste trabalho, apontando tópicos de nosso interesse de prosseguimento nesta linha de pesquisa. 3
19 Capítulo 2 Modelos Determinísticos Modelos matemáticos determinísticos são um conjunto de equações ou inequações matemáticas, organizadas de forma que, conhecidas algumas condições deste sistema, é possível obter sua solução em um momento desejado. Portanto, dadas as mesmas condições do sistema, a solução é sempre igual. Mais especificamente, aqui definimos um modelo determinístico M como uma função que, para cada t em [0, T ], relaciona um conjunto de parâmetros θ Θ R m a um conjunto de variáveis de saída M(θ, t) R p, dada as condições de entrada deste sistema. Desta forma, modelos determinísticos não levam em conta nenhuma aleatoriedade. Isto é, dado um conjunto de entradas do modelo, sua saída é unicamente determinada. Em particular, são modelos determinísticos sistemas de equações a diferenças finitas (discreto) e sistemas de equações diferenciais ordinárias (EDO s) (contínuo), os quais temos interesse neste trabalho. Exemplos de condições de entrada, que não são parâmetros do modelo, são as condições iniciais do sistema, M(θ, 0), no caso de EDO s. Estas condições são necessárias para obtermos uma solução única. 2.1 Equações a Diferenças Finitas Bassanezi (2002) sugere que o uso de equações a diferenças finitas em modelagem é apropriado quando, por exemplo, o crescimento populacional, entre gerações sucessivas, ocorre em etapas discretas e sem sobreposição de gerações da espécie analisada. Várias 4
20 aplicações desta classe de equações são apresentadas em Bassanezi (2002), assim como um estudo introdutório desta teoria é visto em Lima (2006). Segue abaixo a definição matemática destas equações. Definição Uma equação a diferenças de 1 a ordem é uma equação da forma x t+1 = f(x t ), (2.1) onde f é uma função determinada e t {0, 1,..., N,...}. A solução da equação a diferenças finitas em (2.1) é uma sequência, cujos termos satisfazem a relação definida por f, da forma (x 0, x 1,..., x N,...). Além disto, segundo Lima (2006), supondo conhecido o ponto inicial x 0, ou primeiro termo da sequência, então a equação em (2.1) possui uma única solução. 2.2 Equações Diferenciais Ordinárias A vasta aplicabilidade dos conceitos do cálculo diferencial para resolver problemas práticos tem movido vários pesquisadores, desde o final do século XVII quando I. Newton e G. W. Leibnitz, motivados por problemas físicos e geométricos, o desenvolveram. A seguir, apresentamos alguns conceitos indispensáveis para uma melhor compreensão desta classe de equações, inclusive sua definição. Definição Seja f uma função definida num aberto Ω que associa cada par (t, x) Ω (R R n ) a um ponto pertencente ao R n com dx dt = f(t, x), (2.2) então, dizemos que (2.2) é uma equação diferencial ordinária em R n definida por f. Uma solução da EDO em (2.2) é uma função derivável φ : I R n com (t, φ(t)) Ω, para todo t I R, I intervalo. Além disto, supondo conhecido um ponto da solução, digamos x(t 0 ) = x 0, temos um problema de valor inicial (PVI), também chamado de problema de Cauchy. 5
21 A fim de garantir a existência de soluções de EDO s definidas no R n, basta que a função f seja contínua. Entretanto, para garantir unicidade de solução de um PVI é preciso assumir outras hipóteses além da continuidade da função f. O teorema abaixo trata deste problema, e sua demonstração pode ser apreciada em Doering e Lopes (2007). Teorema 2.1 (Teorema de Existência e Unicidade de Solução) Seja f : Ω R n uma função contínua definida num aberto Ω R n+1. Suponhamos que a derivada parcial espacial, f x, seja também contínua neste aberto. Então, para cada (t 0, x 0 ) Ω existem um intervalo aberto I contendo t 0 e uma única função diferenciável φ : I R n com (t, φ(t)) Ω, para todo t I, que é solução do problema de valor inicial dx dt = f(t, x), x(t 0 ) = x 0. Se escrevermos x(t) = (x 1 (t),..., x n (t)), e f(t, x) = (f 1 (t, x),..., f n (t, x)), então podemos interpretar a equação diferencial vetorial dx dt = f(t, x) Rn como um sistema de equações diferenciais escalares da forma dx 1 dt dx 2 dt = f 1 (t, x 1 (t),..., x n (t)), = f 2 (t, x 1 (t),..., x n (t)),. dx n dt = f n (t, x 1 (t),..., x n (t)). Uma solução deste sistema é dada por um conjunto de funções deriváveis φ l : I R, l = 1,..., n, tais que, para cada t I dφ l dt = f l (t, φ 1 (t),..., φ n (t)). Equivalentemente, o vetor φ = (φ 1, φ 2,..., φ n ) representa a solução da equação vetorial dx = f(t, x) R n. Desta forma, o teorema de existência e unicidade de soluções dt apresentado acima para uma equação diferencial ordinária vetorial com condição inicial conhecida, também garante a existência e unicidade de solução para um sistema de equações diferenciais escalares em R dada sua condição inicial. Neste caso, uma condição inicial é dada por x 1 (t 0 ) = x 0 1,..., x n (t 0 ) = x 0 n. 6
22 2.3 Métodos Numéricos de Aproximação de Solução de EDO No início do desenvolvimento das equações diferenciais, a principal preocupação dos pesquisadores era a obtenção da solução destas equações de forma explícita. Entretanto, logo se verificou que a quantidade de equações que podiam ser resolvidas analiticamente era muito pequena. Segundo Figueiredo e Neves (2005), com a rigorosidade da análise matemática, surgida mais formalmente no século XIX, passou-se a se considerar primeiro a existência e unicidade de soluções. Além disto, vários estudos vêm sendo realizados ao longo dos anos, em favor de métodos para solução destas equações. Quando a equação é linear, a solução torna-se mais simples, podendo ser resolvida analiticamente. outro lado, quando a equação apresenta algum termo não-linear, em geral, não conseguimos resolvê-la exatamente. Neste contexto, métodos numéricos de aproximação de soluções surgem como ferramentas poderosas, principalmente com o advento e grande progresso da capacidade computacional. No caso de EDO s, dada sua condição inicial, alguns dos métodos mais conhecidos são o método de Euler (e alguns aprimoramentos deste método), o método de Runge-Kutta, o método da colocação e métodos de passos múltiplos. Estes métodos são descritos, por exemplo, em Boyce e DiPrima (2001) e Burden e Faires (1993). Vários softwares matemáticos possuem pacotes com alguns destes métodos implementados, que são de fácil manipulação. A fim de fazermos uma breve introdução dos métodos numéricos citados, vamos considerar o caso mais simples. Seja uma equação diferencial ordinária não-linear da forma dx dt = f(t, x), tais que (t, x) (R R) e com condição inicial x(t 0 ) = x 0. Além disto, vamos supor que f e a derivada parcial com relação à segunda variável, f, sejam contínuas. Daí, pelo teorema 2.1, existe uma única solução do problema dado x num intervalo aberto I, contendo t 0. O primeiro passo para utilizarmos um método numérico consiste em definir N subintervalos do intervalo I. Denotaremos estes N subintervalos por [t n 1, t n ], n = 1,..., N. O Por 7
23 próximo passo é aproximar a solução dentro de cada um destes subintervalos Método de Euler Embora o método de Euler seja raramente utilizado em aplicações reais, sua importância teórica para o entendimento de vários outros métodos mais robustos justifica sua apresentação. Este é um método iterativo, no qual devemos calcular recursivamente a seguinte equação x n = x n 1 + hf(t n 1, x n 1 ), n = 1,..., N, tal que h é o tamanho dos subintervalos [t n 1, t n ], n = 1,..., N. Este é o método de Euler mais simples, pois consideramos aqui o tamanho dos intervalos igual para todos. Desta forma, obtemos uma sequência de pares (t 0, x 0 ),..., (t N, x N ) que aproximam a solução alvo no intervalo I = [t 0, t N ]. Embora a utilização de métodos numéricos seja extremamente útil, esta solução obtida sofre algumas penalizações devido a erros de aproximação. Os principais erros de aproximação nestes casos são os erros de truncamento e de arredondamento. Basicamente, o erro de arredondamento está relacionado à limitação computacional no sentido de que os cálculos são restritos a um número finito de dígitos. Já o erro de truncamento origina-se do fato de que em cada passo somente utilizamos uma fórmula aproximada para obter cada x n, e, além disto, supondo este valor correto, o utilizamos para obter o próximo termo da sequência. Em particular, chamamos de erro de truncamento local o erro efetuado em cada passo. Para uma vasta classe de métodos numéricos é possível obter estimativas destes erros. O método de Euler pode também ser visto como um caso particular dos métodos de série de Taylor, com a presença apenas do termo de primeira ordem (primeira derivada). Na prática, estes métodos com ordens superiores são dificilmente empregados pois requerem as avaliações de derivadas também de ordens superiores, que nem sempre existem, ou ainda, tornam a implementação mais complexa. 8
24 2.3.2 Método de Runge-Kutta O método de Runge-Kutta também pode ser considerado uma generalização do método de Euler. Neste método, podemos variar o tamanho dos subintervalos tomados para a aproximação. Dentre as várias formas deste método, o mais amplamente utilizado é o método de Runge-Kutta de quarta ordem, e sua fórmula de recorrência é dada por ( ) kn 1,1 + 2k n 1,2 + 2k n 1,3 + k n 1,4 x n = x n 1 + h, n = 1,..., N, 6 tais que k n 1,1 = f(t n 1, x n 1 ), k n 1,2 = f(t n h, x n hk n 1,1), k n 1,3 = f(t n h, x n hk n 1,2), k n 1,4 = f(t n h, x n hk n 1,3), para n = 1,..., N. O termo que acompanha h pode ser interpretado como uma média ponderada de valores de f(t, x) em diferentes pontos do intervalo [t n 1, t n ], ou ainda como um coeficiente angular médio. A ampla utilização deste método pode ser explicada devido a ser um método razoavelmente simples de ser manipulado e que retorna uma aproximação satisfatória em muitos casos Métodos de Passos Múltiplos Métodos que utilizam informação de mais de um passo anterior são chamados de métodos de passos múltiplos. Deste modo, os métodos anteriores não fazem parte desta classe de métodos. Entretanto, existem vários métodos que apresentam esta característica. Como ilustração, vejamos um destes métodos, conhecido como método de Adams. A idéia básica deste método é aproximar f(t, x) por um polinômio P k (t) de grau k 1, e usar este polinômio para calcular a integral tn+1 t n f(t, x) dt = x n+1 x n. 9