Métodos numéricos para equações diferenciais ordinárias - Numerical methods for ordinary differential equations

Os métodos numéricos para equações diferenciais ordinárias são métodos usados ​​para encontrar aproximações numéricas para as soluções de equações diferenciais ordinárias (EDOs). Seu uso também é conhecido como " integração numérica ", embora este termo também possa se referir ao cálculo de integrais .

Muitas equações diferenciais não podem ser resolvidas usando computação simbólica ("análise"). Para fins práticos, no entanto - como em engenharia - uma aproximação numérica da solução costuma ser suficiente. Os algoritmos estudados aqui podem ser usados ​​para calcular tal aproximação. Um método alternativo é usar técnicas de cálculo para obter uma expansão em série da solução.

Equações diferenciais comuns ocorrem em muitas disciplinas científicas, incluindo física , química , biologia e economia . Além disso, alguns métodos em equações diferenciais parciais numéricas convertem a equação diferencial parcial em uma equação diferencial ordinária, que deve então ser resolvida.

O problema

Uma equação diferencial de primeira ordem é um problema de valor inicial (IVP) da forma,

${\ displaystyle y '(t) = f (t, y (t)), \ qquad y (t_ {0}) = y_ {0},}$

( 1 )

onde é uma função e a condição inicial é um determinado vetor. Primeira ordem significa que apenas a primeira derivada de y aparece na equação e as derivadas superiores estão ausentes. ${\ displaystyle f}$${\ displaystyle f: [t_ {0}, \ infty) \ times \ mathbb {R} ^ {d} \ to \ mathbb {R} ^ {d}}$${\ displaystyle y_ {0} \ in \ mathbb {R} ^ {d}}$

Sem perda de generalidade para sistemas de ordem superior, nos restringimos às equações diferenciais de primeira ordem , porque uma ODE de ordem superior pode ser convertida em um sistema maior de equações de primeira ordem pela introdução de variáveis ​​extras. Por exemplo, a equação de segunda ordem y ′ ′ = - y pode ser reescrita como duas equações de primeira ordem: y ′ = z e z ′ = - y .

Nesta seção, descrevemos métodos numéricos para IVPs e observamos que os problemas de valor limite (BVPs) exigem um conjunto diferente de ferramentas. Em um BVP, define-se valores ou componentes da solução y em mais de um ponto. Por causa disso, diferentes métodos precisam ser usados ​​para resolver BVPs. Por exemplo, o método de disparo (e suas variantes) ou métodos globais como diferenças finitas , métodos de Galerkin ou métodos de colocação são apropriados para essa classe de problemas.

O teorema de Picard-Lindelöf afirma que existe uma solução única, desde que f seja contínuo de Lipschitz .

Métodos

Os métodos numéricos para resolver IVPs de primeira ordem geralmente se enquadram em uma de duas grandes categorias: métodos lineares de várias etapas ou métodos de Runge-Kutta . Uma divisão adicional pode ser realizada dividindo os métodos em explícitos e implícitos. Por exemplo, os métodos de várias etapas lineares implícitos incluem métodos de Adams-Moulton e métodos de diferenciação retroativa (BDF), enquanto os métodos Runge – Kutta implícitos incluem Runge – Kutta (DIRK) diagonalmente implícito, Runge – Kutta (SDIRK) implícito diagonalmente implícito, e Gauss– Métodos numéricos de Radau (baseado na quadratura gaussiana ). Exemplos explícitos da família linear de várias etapas incluem os métodos Adams – Bashforth e qualquer método Runge – Kutta com um tableau de Butcher diagonal inferior é explícito . Uma regra geral flexível dita que equações diferenciais rígidas requerem o uso de esquemas implícitos, enquanto problemas não rígidos podem ser resolvidos de forma mais eficiente com esquemas explícitos.

Os chamados métodos lineares gerais (GLMs) são uma generalização das duas grandes classes de métodos acima.

Método de Euler

De qualquer ponto de uma curva, você pode encontrar uma aproximação de um ponto próximo na curva movendo uma curta distância ao longo de uma linha tangente à curva.

Começando com a equação diferencial ( 1 ), substituímos a derivada y ′ pela aproximação de diferença finita

${\ displaystyle y '(t) \ approx {\ frac {y (t + h) -y (t)} {h}},}$

( 2 )

que quando reorganizado produz a seguinte fórmula

${\ displaystyle y (t + h) \ approx y (t) + hy '(t)}$

e usando ( 1 ) dá:

${\ displaystyle y (t + h) \ approx y (t) + hf (t, y (t)).}$

( 3 )

Esta fórmula é geralmente aplicada da seguinte maneira. Escolhemos um passo h , e construímos a sequência . Denotamos por uma estimativa numérica da solução exata . Motivados por ( 3 ), calculamos essas estimativas pelo seguinte esquema recursivo${\ displaystyle t_ {0}, t_ {1} = t_ {0} + h, t_ {2} = t_ {0} + 2h, ...}$${\ displaystyle y_ {n}}$${\ displaystyle y (t_ {n})}$

${\ displaystyle y_ {n + 1} = y_ {n} + hf (t_ {n}, y_ {n}).}$

( 4 )

Este é o método de Euler (ou método de Euler progressivo , em contraste com o método de Euler retroativo , a ser descrito abaixo). O método tem o nome de Leonhard Euler, que o descreveu em 1768.

O método de Euler é um exemplo de método explícito . Isso significa que o novo valor y _{n +1} é definido em termos de coisas que já são conhecidas, como y _n .

Método Backward Euler

Se, em vez de ( 2 ), usarmos a aproximação

${\ displaystyle y '(t) \ approx {\ frac {y (t) -y (th)} {h}},}$

( 5 )

obtemos o método de Euler reverso :

${\ displaystyle y_ {n + 1} = y_ {n} + hf (t_ {n + 1}, y_ {n + 1}).}$

( 6 )

O método de Euler retroativo é um método implícito , o que significa que temos que resolver uma equação para encontrar y _{n +1} . Freqüentemente, usa -se a iteração de ponto fixo ou (alguma modificação do) método de Newton-Raphson para conseguir isso.

Custa mais tempo para resolver essa equação do que métodos explícitos; esse custo deve ser levado em consideração ao selecionar o método a ser usado. A vantagem de métodos implícitos como ( 6 ) é que eles geralmente são mais estáveis ​​para resolver uma equação rígida , o que significa que um tamanho de passo maior h pode ser usado.

Método integrador exponencial de primeira ordem

Os integradores exponenciais descrevem uma grande classe de integradores que recentemente viram muito desenvolvimento. Eles datam de pelo menos 1960.

No lugar de ( 1 ), assumimos que a equação diferencial tem a forma

${\ displaystyle y '(t) = - A \, y + {\ mathcal {N}} (y),}$

( 7 )

ou foi linearizado localmente sobre um estado de fundo para produzir um termo linear e um termo não linear . ${\ displaystyle -Ay}$${\ displaystyle {\ mathcal {N}} (y)}$

Os integradores exponenciais são construídos multiplicando ( 7 ) por e integrando exatamente o resultado ao longo de um intervalo de tempo : ${\ textstyle e ^ {At}}$${\ displaystyle [t_ {n}, t_ {n + 1} = t_ {n} + h]}$

${\ displaystyle y_ {n + 1} = e ^ {- Ah} y_ {n} + \ int _ {0} ^ {h} e ^ {- (h- \ tau) A} {\ mathcal {N}} \ left (y \ left (t_ {n} + \ tau \ right) \ right) \, d \ tau.}$

Esta equação integral é exata, mas não define a integral.

O integrador exponencial de primeira ordem pode ser realizado mantendo constante durante todo o intervalo: ${\ displaystyle {\ mathcal {N}} (y ​​(t_ {n} + \ tau))}$

${\ displaystyle y_ {n + 1} = e ^ {- Ah} y_ {n} + A ^ {- 1} (1-e ^ {- Ah}) {\ mathcal {N}} (y ​​(t_ {n })) \.}$

( 8 )

Generalizações

O método de Euler muitas vezes não é preciso o suficiente. Em termos mais precisos, ele só tem ordem um (o conceito de ordem é explicado a seguir). Isso fez com que os matemáticos procurassem métodos de ordem superior.

Uma possibilidade é usar não apenas o valor previamente calculado y _n para determinar y _{n +1} , mas fazer a solução depender de mais valores anteriores. Isso produz o chamado método de várias etapas . Talvez o mais simples seja o método do salto, que é de segunda ordem e (falando grosso modo) depende de dois valores de tempo.

Quase todos os métodos de várias etapas práticos se enquadram na família dos métodos lineares de várias etapas , que têm a forma

${\ displaystyle {\ begin {alinhados} & {} \ alpha _ {k} y_ {n + k} + \ alpha _ {k-1} y_ {n + k-1} + \ cdots + \ alpha _ {0 } y_ {n} \\ & {} \ quad = h \ left [\ beta _ {k} f (t_ {n + k}, y_ {n + k}) + \ beta _ {k-1} f ( t_ {n + k-1}, y_ {n + k-1}) + \ cdots + \ beta _ {0} f (t_ {n}, y_ {n}) \ right]. \ end {alinhado}} }$

Outra possibilidade é usar mais pontos no intervalo . Isso leva à família dos métodos Runge – Kutta , nomeados em homenagem a Carl Runge e Martin Kutta . Um de seus métodos de quarta ordem é especialmente popular. ${\ displaystyle [t_ {n}, t_ {n + 1}]}$

Características avançadas

Uma boa implementação de um desses métodos para resolver uma ODE envolve mais do que a fórmula de passo de tempo.

Freqüentemente, é ineficiente usar o mesmo tamanho de passo o tempo todo, então foram desenvolvidos métodos de tamanho de passo variável . Normalmente, o tamanho do passo é escolhido de forma que o erro (local) por passo esteja abaixo de algum nível de tolerância. Isso significa que os métodos também devem calcular um indicador de erro , uma estimativa do erro local.

Uma extensão dessa ideia é escolher dinamicamente entre diferentes métodos de diferentes ordens (isso é chamado de método de ordem variável ). Métodos baseados na extrapolação de Richardson , como o algoritmo Bulirsch-Stoer , são freqüentemente usados ​​para construir vários métodos de diferentes ordens.

Outros recursos desejáveis ​​incluem:

• saída densa : aproximações numéricas baratas para todo o intervalo de integração, e não apenas nos pontos t ₀ , t ₁ , t ₂ , ...
• localização do evento : encontrar os horários em que, digamos, uma função específica desaparece. Isso normalmente requer o uso de um algoritmo de localização de raiz .
• suporte para computação paralela .
• quando usado para integração com respeito ao tempo, reversibilidade do tempo

Métodos alternativos

Muitos métodos não se enquadram na estrutura discutida aqui. Algumas classes de métodos alternativos são:

• métodos multiderivativos , que usam não apenas a função f, mas também suas derivadas. Esta classe inclui os métodos Hermite – Obreschkoff e os métodos Fehlberg , bem como métodos como o método Parker – Sochacki ou o método Bychkov – Scherbakov, que calculam os coeficientes da série de Taylor da solução y recursivamente.
• métodos para EDOs de segunda ordem . Dissemos que todos os ODEs de ordem superior podem ser transformados em ODEs de primeira ordem da forma (1). Embora isso certamente seja verdade, pode não ser a melhor maneira de proceder. Em particular, os métodos de Nyström trabalham diretamente com equações de segunda ordem.
• métodos de integração geométrica são especialmente projetados para classes especiais de EDOs (por exemplo, integradores simpléticos para a solução de equações hamiltonianas ). Eles cuidam para que a solução numérica respeite a estrutura ou geometria subjacente dessas classes.
• Os métodos de sistemas de estado quantizados são uma família de métodos de integração ODE baseados na ideia de quantização de estado. Eles são eficientes ao simular sistemas esparsos com descontinuidades frequentes.

Métodos paralelos no tempo

Para aplicativos que requerem computação paralela em supercomputadores , o grau de simultaneidade oferecido por um método numérico torna-se relevante. Em vista dos desafios dos sistemas de computação exascale , métodos numéricos para problemas de valor inicial que podem fornecer simultaneidade na direção temporal estão sendo estudados. Parareal é um exemplo relativamente conhecido de tal método de integração paralelo no tempo , mas as primeiras ideias remontam à década de 1960.

Análise

A análise numérica não é apenas o projeto de métodos numéricos, mas também sua análise. Três conceitos centrais nesta análise são:

• convergência : se o método se aproxima da solução,
• ordem : quão bem ele se aproxima da solução, e
• estabilidade : se os erros são amortecidos.

Convergência

Um método numérico é considerado convergente se a solução numérica se aproxima da solução exata quando o tamanho do passo h vai para 0. Mais precisamente, exigimos que para cada ODE (1) com uma função de Lipschitz f e todo t ^*  > 0,

${\ displaystyle \ lim _ {h \ to 0 ^ {+}} \ max _ {n = 0,1, \ dots, \ lfloor t ^ {*} / h \ rfloor} \ left \ | y_ {n, h } -y (t_ {n}) \ right \ | = 0.}$

Todos os métodos mencionados acima são convergentes.

Consistência e ordem

Suponha que o método numérico seja

${\ displaystyle y_ {n + k} = \ Psi (t_ {n + k}; y_ {n}, y_ {n + 1}, \ dots, y_ {n + k-1}; h). \,}$

O erro local (truncamento) do método é o erro cometido por uma etapa do método. Ou seja, é a diferença entre o resultado dado pelo método, supondo que nenhum erro foi cometido nas etapas anteriores, e a solução exata:

${\ displaystyle \ delta _ {n + k} ^ {h} = \ Psi \ left (t_ {n + k}; y (t_ {n}), y (t_ {n + 1}), \ dots, y (t_ {n + k-1}); h \ right) -y (t_ {n + k}).}$

O método é considerado consistente se

${\ displaystyle \ lim _ {h \ to 0} {\ frac {\ delta _ {n + k} ^ {h}} {h}} = 0.}$

O método tem ordem se ${\ displaystyle p}$

${\ displaystyle \ delta _ {n + k} ^ {h} = O (h ^ {p + 1}) \ quad {\ mbox {as}} h \ to 0.}$

Portanto, um método é consistente se tiver uma ordem maior do que 0. O método de Euler (progressivo) (4) e o método de Euler retroativo (6) introduzidos acima têm ambos ordem 1, portanto, são consistentes. A maioria dos métodos usados ​​na prática atinge uma ordem superior. A consistência é uma condição necessária para a convergência, mas não suficiente; para um método ser convergente, ele deve ser consistente e estável em zero .

Um conceito relacionado é o erro global (truncamento) , o erro sustentado em todas as etapas necessárias para atingir um tempo fixo . Explicitamente, o erro global no momento é onde . O erro global de um método de uma etapa da ordem é ; em particular, tal método é convergente. Essa afirmação não é necessariamente verdadeira para métodos de várias etapas. ${\ displaystyle t}$${\ displaystyle t}$${\ displaystyle y_ {N} -y (t)}$${\ displaystyle N = (t-t_ {0}) / h}$${\ displaystyle p}$${\ displaystyle O (h ^ {p})}$

Estabilidade e rigidez

Para algumas equações diferenciais, a aplicação de métodos padrão - como o método de Euler, métodos explícitos de Runge-Kutta ou métodos de várias etapas (por exemplo, métodos de Adams-Bashforth) - exibe instabilidade nas soluções, embora outros métodos possam produzir soluções estáveis. Esse "comportamento difícil" na equação (que pode não ser necessariamente complexo em si) é descrito como rigidez e geralmente é causado pela presença de escalas de tempo diferentes no problema subjacente. Por exemplo, uma colisão em um sistema mecânico como em um oscilador de impacto normalmente ocorre em uma escala de tempo muito menor do que o tempo para o movimento dos objetos; esta discrepância cria "curvas fechadas" nas curvas dos parâmetros de estado.

Problemas rígidos são onipresentes em cinética química , teoria de controle , mecânica de sólidos , previsão do tempo , biologia , física de plasma e eletrônica . Uma forma de superar a rigidez é estender a noção de equação diferencial para a de inclusão diferencial , que permite e modela a não suavidade.

História

Abaixo está uma linha do tempo de alguns desenvolvimentos importantes neste campo.

• 1768 - Leonhard Euler publica seu método.
• 1824 - Augustin Louis Cauchy prova a convergência do método de Euler. Nessa prova, Cauchy usa o método implícito de Euler.
• 1855 - Primeira menção dos métodos de várias etapas de John Couch Adams em uma carta escrita por Francis Bashforth .
• 1895 - Carl Runge publica o primeiro método Runge – Kutta .
• 1901 - Martin Kutta descreve o popular método Runge-Kutta de quarta ordem .
• 1910 - Lewis Fry Richardson anuncia seu método de extrapolação , extrapolação de Richardson .
• 1952 - Charles F. Curtiss e Joseph Oakland Hirschfelder cunham o termo equações rígidas .
• 1963 - Germund Dahlquist introduz a estabilidade A dos métodos de integração.

Soluções numéricas para problemas de valor limite unidimensional de segunda ordem

Os problemas de valor limite (BVPs) são geralmente resolvidos numericamente resolvendo um problema de matriz aproximadamente equivalente obtido pela discretização do BVP original. O método mais comumente usado para resolver numericamente BVPs em uma dimensão é chamado de Método de Diferença Finita . Este método aproveita as combinações lineares de valores de pontos para construir coeficientes de diferença finita que descrevem as derivadas da função. Por exemplo, a aproximação da diferença central de segunda ordem para a primeira derivada é dada por:

${\ displaystyle {\ frac {u_ {i + 1} -u_ {i-1}} {2h}} = u '(x_ {i}) + {\ mathcal {O}} (h ^ {2}), }$

e a diferença central de segunda ordem para a segunda derivada é dada por:

${\ displaystyle {\ frac {u_ {i + 1} -2u_ {i} + u_ {i-1}} {h ^ {2}}} = u '' (x_ {i}) + {\ mathcal {O }} (h ^ {2}).}$

Em ambas as fórmulas, é a distância entre os valores x vizinhos no domínio discretizado. Em seguida, constrói-se um sistema linear que pode então ser resolvido por métodos de matriz padrão . Por exemplo, suponha que a equação a ser resolvida seja: ${\ displaystyle h = x_ {i} -x_ {i-1}}$

${\ displaystyle {\ begin {alinhado} & {} {\ frac {d ^ {2} u} {dx ^ {2}}} - u = 0, \\ & {} u (0) = 0, \\ & {} u (1) = 1. \ end {alinhado}}}$

A próxima etapa seria discretizar o problema e usar aproximações derivadas lineares, como

${\ displaystyle u '' _ {i} = {\ frac {u_ {i + 1} -2u_ {i} + u_ {i-1}} {h ^ {2}}}}$

e resolver o sistema resultante de equações lineares. Isso levaria a equações como:

${\ displaystyle {\ frac {u_ {i + 1} -2u_ {i} + u_ {i-1}} {h ^ {2}}} - u_ {i} = 0, \ quad \ forall i = {1 , 2,3, ..., n-1}.}$

À primeira vista, este sistema de equações parece ter dificuldade em ser associado ao fato de que a equação não envolve termos que não sejam multiplicados por variáveis, mas na verdade isso é falso. No i = 1 e n - 1 não é um termo que envolve os valores de limite e e uma vez que estes dois valores são conhecidos, pode-se simplesmente substituí-los para esta equação e como resultado, têm um sistema linear não homogénea de equações que tem não- soluções triviais. ${\ displaystyle u (0) = u_ {0}}$${\ displaystyle u (1) = u_ {n}}$

Veja também

• Condição de Courant-Friedrichs-Lewy
• Deriva de energia
• Métodos lineares gerais
• Lista de tópicos de análise numérica # Métodos numéricos para equações diferenciais ordinárias
• Algoritmo de propagação do sistema de referência reversível
• Software Modelica Language e OpenModelica

Notas

Referências

• Bradie, Brian (2006). Uma introdução amigável à análise numérica . Upper Saddle River, Nova Jersey: Pearson Prentice Hall. ISBN 978-0-13-013054-9.
• JC Butcher , métodos numéricos para equações diferenciais ordinárias , ISBN  0-471-96758-0
• Ernst Hairer, Syvert Paul Nørsett e Gerhard Wanner, Resolvendo equações diferenciais ordinárias I: Nonstiff problems, segunda edição, Springer Verlag, Berlin, 1993. ISBN  3-540-56670-8 .
• Ernst Hairer e Gerhard Wanner, Resolvendo equações diferenciais ordinárias II: Stiff and diferencial-algebraic problems, segunda edição, Springer Verlag, Berlin, 1996. ISBN  3-540-60452-9 .
  (Esta monografia de dois volumes cobre sistematicamente todos os aspectos do campo.)
• Hochbruck, Marlis ; Ostermann, Alexander (maio de 2010). "Integradores exponenciais". Acta Numerica . 19 : 209–286. Bibcode : 2010AcNum..19..209H . CiteSeerX  10.1.1.187.6794 . doi : 10.1017 / S0962492910000048 .
• Arieh Iserles, A First Course in the Numerical Analysis of Differential Equations, Cambridge University Press, 1996. ISBN  0-521-55376-8 (capa dura), ISBN  0-521-55655-4 (brochura).
  (Livro didático, voltado para alunos avançados de graduação e pós-graduação em matemática, que também discute equações diferenciais parciais numéricas .)
• John Denholm Lambert, Numerical Methods for Ordinary Differential Systems, John Wiley & Sons, Chichester, 1991. ISBN  0-471-92990-5 .
  (Livro didático, um pouco mais exigente do que o livro de Iserles.)

links externos

• Joseph W. Rudmin, Aplicação do Método Parker-Sochacki à Mecânica Celestial , 1998.
• Dominique Tournès, L'intégration approchée des équations différentielles ordinaires (1671-1914) , thèse de doctorat de l'université Paris 7 - Denis Diderot, junho de 1996. Réimp. Villeneuve d'Ascq: Presses universitaires du Septentrion, 1997, 468 p. (Amplo material online sobre história de análise numérica de ODE, para material em inglês sobre a história da análise numérica de ODE, veja, por exemplo, os livros de papel de Chabert e Goldstine citados por ele.)
• Pchelintsev, AN (2020). "Um método e algoritmo numérico preciso para a construção de soluções de sistemas caóticos" (PDF) . Journal of Applied Nonlinear Dynamics . 9 (2): 207–221. doi : 10.5890 / JAND.2020.06.004 .
• kv no GitHub ( biblioteca C ++ com solucionadores ODE rigorosos)
• INTLAB (uma biblioteca feita por MATLAB / GNU Octave que inclui solucionadores ODE rigorosos)