Método de Picard

Implementação completa do experimento § 2.6 de Equações Diferenciais Ordinárias (Viana & Espinar, IMPA 2025), cobrindo a iteração de Picard para aproximação de soluções de EDOs. O projeto inclui análise de convergência numérica, comparativo entre os integradores do trapézio e do retângulo, otimização de complexidade $\mathcal{O}(NL^2) \to \mathcal{O}(NL)$ via acumulação incremental e extensão para sistemas de dimensão $d = 2$.

O problema de valor inicial e o operador de Picard

Considera-se o problema de valor inicial

onde $F : \mathcal{U} \to \mathbb{R}^d$ é contínua num aberto $\mathcal{U} \subset \mathbb{R} \times \mathbb{R}^d$ e localmente Lipschitziana em $x$. A ideia central do método de Picard é reescrever (1) como equação integral e aplicar o Teorema do Ponto Fixo para Contrações.

Dado $\varepsilon > 0$, o operador de Picard $\mathcal{L} : Y \to Y$ age em funções contínuas $\gamma : (t_0 - \varepsilon, t_0 + \varepsilon) \to \mathbb{R}^d$, com $\gamma(t_0) = x_0$, por

A sequência de iterados de Picard é definida por $\gamma_0 \equiv x_0$ e $\gamma_{n+1} := \mathcal{L}\gamma_n$ para $n \ge 0$.

Existência, unicidade e raio garantido

Dado $(t_0, x_0) \in \mathcal{U}$, fixe $\delta > 0$ tal que $\overline{B}_\delta(t_0) \times \overline{B}_\delta(x_0) \subset \mathcal{U}$. Define-se

e seja $C(\delta)$ a constante de Lipschitz de $F$ em $x$ sobre essa região. Para todo $\varepsilon > 0$ satisfazendo

existe solução única $x : (t_0 - \varepsilon, t_0 + \varepsilon) \to \mathbb{R}^d$ de (1), e os iterados de Picard convergem uniformemente para ela nesse intervalo (Teorema 2.12 de Viana & Espinar).

Cota do erro

O resultado central para toda a análise numérica é o Corolário 2.14:

O fatorial no denominador garante convergência para qualquer valor de $C\varepsilon$ — apenas a taxa varia:

• Se $C(\delta)\varepsilon \le 1$: o fatorial domina a potência desde o primeiro iterado, e a cota decresce monotonamente em $n$.
• Se $C(\delta)\varepsilon > 1$: o termo $(C\varepsilon)^n/n!$ cresce até atingir um máximo em torno de $n \approx C(\delta)\varepsilon$ e só então decai via o fatorial. É o que ocorre no Objetivo 4, com $\varepsilon_{\text{num}} = 0{,}99$ e $C\varepsilon \approx 3{,}96$.

A equação (6) conecta os parâmetros ao número mínimo de iterados necessários para atingir uma tolerância $\tau$:

Parâmetros ótimos para $F(t,x) = x^2$

Para $F(t,x) = x^2$ e $(t_0, x_0) = (0, 1)$, as quantidades de (3) sobre $[-\delta, \delta] \times [1-\delta, 1+\delta]$ são

onde a última igualdade segue de $(1+\delta)^2 > 1$ para $\delta > 0$, que garante que o mínimo em (4) é sempre o segundo argumento. Maximizando $\varepsilon(\delta)$:

Com $C(\delta^*)\varepsilon^* = 1$, a cota (6) reduz-se a $\|\gamma_n - x\|_\infty \le 4/n!$, que para $n = 8$ fornece $4/8! \approx 9{,}9 \times 10^{-5}$. A convergência observada é tipicamente melhor que esta cota superior, uma vez que (6) usa o supremo sobre toda a bola $B_{\delta^*}(x_0)$ e não apenas sobre a trajetória da solução.

Discretização e erro de quadratura

A curva $\gamma_n$ é representada por $L$ valores em pontos igualmente espaçados

Para que $t_0 = 0$ coincida com um índice da malha, é necessário que $L$ seja ímpar; nesse caso $i_0 := \lfloor L/2 \rfloor$ satisfaz $t_{i_0} = 0$.

A integral em (2) é aproximada por duas regras clássicas:

O erro total do $n$-ésimo iterado discreto $\gamma_n^h$ decompõe-se em dois termos com dependências opostas:

com $p = 2$ para o trapézio e $p = 1$ para o retângulo. Aumentar $n$ reduz o primeiro termo via o fatorial; aumentar $L$ reduz o segundo como $L^{-p}$. Para $n$ fixo e suficientemente grande, o retângulo atinge um piso de erro irredutível proporcional a $\varepsilon^2/L$, que não pode ser eliminado aumentando o número de iterados — apenas aumentando $L$ ou diminuindo $\varepsilon$.

O experimento está organizado em cinco objetivos progressivos, correspondentes às instruções da Seção 2.6:

1. Escrever o método de Picard em código computacional.
2. Considerar $x' = x^2$ com $x(0) = 1$, calcular os primeiros iterados num domínio $[0, \varepsilon]$ adequado, comparar com a solução exata $x(t) = 1/(1-t)$ e estimar o erro do $n$-ésimo iterado.
3. Calcular soluções no domínio $(-\varepsilon, \varepsilon)$, observando o comportamento das duas regras de quadratura no extremo esquerdo da curva.
4. Investigar a convergência dos iterados em domínios maiores que $(-\varepsilon, \varepsilon)$, inclusive em todo o intervalo $(-\infty, 1)$ de definição da solução.
5. Repetir os passos anteriores para $x'' = -x + 2\operatorname{sen} t$ com $x(0) = 1$, $x'(0) = 0$, cuja solução exata é $x(t) = (1-t)\cos t + \operatorname{sen} t$. Neste caso o problema tem dimensão $d = 2$.

Todos os códigos usam Python com NumPy e Matplotlib.

O domínio $[0, \varepsilon^*]$ com $\varepsilon^* = 0{,}25$ é escolhido de modo a maximizar o raio garantido para $F(t,x) = x^2$. Com $C(\delta^*)\varepsilon^* = 1$, a cota $\|\gamma_n - x\|_\infty \le 4/n!$ decresce monotonamente desde o primeiro iterado.

import numpy as np

eps = 0.25;  L = 500;  N = 8;  x0 = 1.0
t = np.linspace(0, eps, L)

def F(t, x):
  return x**2

y = np.zeros((N+1, L))
y[0, :] = x0

for n in range(N):
  for k in range(L):
      integrando = F(t[:k+1], y[n, :k+1])
      y[n+1, k]  = x0 + np.trapezoid(integrando, t[:k+1])

A matriz y[n, k] armazena $\gamma_n(t_k)$. O laço interno discretiza a recorrência do operador (2): para cada $t_k$, integra $F(s, \gamma_n(s))$ de $t_0 = 0$ a $t_k$ pela regra do trapézio sobre os $k+1$ pontos já calculados.

O Objetivo 2 estende a análise para $\varepsilon = 0{,}5 > \varepsilon^* = 0{,}25$, saindo da região de garantia do teorema, e introduz a comparação entre os dois integradores.

Embora $\varepsilon = 0{,}5$ ultrapasse a garantia teórica, a convergência numérica ocorre em todo o intervalo $[0,\, 0{,}5]$. Isso decorre de um fato estrutural: os iterados de Picard para $F(t,x) = x^2$ com $x(0) = 1$ coincidem exatamente com as somas parciais da série geométrica

cujo raio de convergência é $|t| < 1$. Portanto, para qualquer $\varepsilon < 1$, os iterados convergem para a solução exata, ainda que o teorema só garanta isso para $\varepsilon \le \varepsilon^* = 0{,}25$. Este é um dos pontos matematicamente mais elegantes do experimento: o método de Picard constrói, iterado a iterado, os coeficientes da expansão em série de Taylor da solução.

h = t[1] - t[0]

for n in range(N):
  for k in range(L):
      integrando   = F(t[:k+1], y_tp[n, :k+1])
      y_tp[n+1, k] = x0 + np.trapezoid(integrando, t[:k+1])

for n in range(N):
  for k in range(L):
      integrando   = F(t[:k+1], y_rt[n, :k+1])
      y_rt[n+1, k] = x0 + np.sum(integrando[:-1]) * h

Na regra do retângulo, integrando[:-1] seleciona os $k$ pontos à esquerda de cada subintervalo, implementando (11).

O trapézio atinge piso numérico da ordem de $10^{-7}$; o retângulo estabiliza em $\approx 6{,}93 \times 10^{-4}$ — um piso que não diminui aumentando $N$.

O Objetivo 3 estende o domínio para $(-\varepsilon, \varepsilon)$, exigindo integração bidirecional a partir de $t_0 = 0$. Essa extensão expõe um problema assimétrico entre os dois integradores: enquanto o trapézio mantém convergência em $L$, o retângulo acumula um piso de erro irredutível que não pode ser reduzido aumentando $N$.

Para $t < t_0 = 0$, a equação integral requer

i0 = L // 2      # índice de t = 0

for n in range(N):
  for k in range(L):
      if k > i0:      # t > 0
          intg         = F(t[i0:k+1], y_tp[n, i0:k+1])
          y_tp[n+1, k] = x0 + np.trapezoid(intg, t[i0:k+1])
      elif k < i0:    # t < 0: sinal negativo
          intg         = F(t[k:i0+1], y_tp[n, k:i0+1])
          y_tp[n+1, k] = x0 - np.trapezoid(intg, t[k:i0+1])
      else:
          y_tp[n+1, k] = x0

O erro de quadratura global da regra do retângulo ao integrar sobre $[t, 0]$ para $t < 0$ é

Para o trapézio, o erro local por subintervalo é $\mathcal{O}(h^3)$, resultando em erro global $\mathcal{O}(\varepsilon^3/L^2)$. O trapézio converge duas ordens a mais em $L$, tornando-o sempre preferível do ponto de vista da convergência.

Os decaimentos observados nos gráficos de erro pontual são pontos onde a curva iterada toca ou se aproxima muito da solução exata. A assimetria entre os lados $t > 0$ e $t < 0$ deve-se ao acúmulo de erro pelo uso de pontos à esquerda.

O Objetivo 4 concentra as contribuições computacionais e teóricas mais relevantes do trabalho. Do lado algorítmico, a implementação direta de complexidade $\mathcal{O}(NL^2)$ é substituída por uma versão incremental de complexidade $\mathcal{O}(NL)$, tornando viável o uso de $L = 20001$ pontos.

Redução de complexidade

Nos Objetivos 1–3, o cálculo de $\gamma_{n+1}(t_k)$ aplica a quadratura no subintervalo $[t_{i_0}, t_k]$, o que exige o recálculo de toda a integral desde o ponto inicial a cada passo, com custo $\Theta(k - i_0)$ operações. Somando sobre todos os pontos à direita de $i_0$:

e de forma análoga no lado negativo, totalizando $\mathcal{O}(L^2)$ por iterado e $\mathcal{O}(NL^2)$ no total. Para $N = 30$ e $L = 20001$, isso representa $\approx 3 \times 10^9$ operações.

A redução a $\mathcal{O}(NL)$ baseia-se na propriedade de acumulação incremental: a regra do trapézio satisfaz

e a do retângulo, $I_k^+ = I_{k-1}^+ + f_{k-1}\,h$. Cada $I_k^+$ custa $\Theta(1)$ operações, reduzindo o custo por iterado a $\mathcal{O}(L)$.

for n in range(N):
  f_tp = F(t, y_tp[n])        # avalia F uma vez em todos os pontos

  I_tp_pos = np.zeros(L)
  for k in range(i0 + 1, L):          # direcao positiva
      I_tp_pos[k] = I_tp_pos[k-1] + (f_tp[k-1] + f_tp[k]) / 2.0 * h

  I_tp_neg = np.zeros(L)
  for k in range(i0 - 1, -1, -1):     # direcao negativa
      I_tp_neg[k] = I_tp_neg[k+1] + (f_tp[k] + f_tp[k+1]) / 2.0 * h

  y_tp[n+1, i0+1:] = x0 + I_tp_pos[i0+1:]
  y_tp[n+1, :i0]   = x0 - I_tp_neg[:i0]
  y_tp[n+1, i0]    = x0

Iterados necessários por tolerância

O valor $N_{\max}(\tau)$ é calculado via a recorrência

que evita overflow numérico ao não calcular as potências e os fatoriais separadamente.

Convergência fora do raio garantido

O código usa $\varepsilon_{\text{num}} = 0{,}99 \gg \varepsilon^* = 0{,}25$. O teorema não garante convergência nesse regime; contudo, a série (15) converge para $|t| < 1$, de modo que os iterados convergem numericamente em $|t| \le 0{,}99$. Para $t \in (0{,}99,\, 1)$ a solução cresce rapidamente em direção à singularidade em $t = 1$, e os iterados precisam de $N \gg N_{\max}(\varepsilon^*)$ iterações para acompanhá-la — consequência de $C\varepsilon_{\text{num}} = 4 \times 0{,}99 \approx 3{,}96 > 1$, que faz os primeiros iterados piorarem antes de convergir.

A equação $x'' = -x + 2\operatorname{sen} t$ com $x(0) = 1$, $x'(0) = 0$ é reescrita como sistema de primeira ordem:

A solução exata é

Ao contrário do Objetivo 4, $\mathbf{x}$ é uma função inteira em $\mathbb{R}$ (sem singularidade real), implicando que os iterados convergem em qualquer domínio compacto.

Adaptações para d = 2

A adaptação é mínima estruturalmente: o vetor de iterados passa a ter forma $(N+1, L, 2)$, os acumuladores $I^+$ e $I^-$ têm forma $(L, 2)$, e a recorrência (17) é aplicada elemento a elemento pelo NumPy sem alteração de código.

def F(t, x):                # x tem estrutura (L, 2)
  dx1 = x[:, 1]
  dx2 = -x[:, 0] + 2.0 * np.sin(t)
  return np.column_stack([dx1, dx2])

y_tp = np.zeros((N+1, L, 2))
y_tp[0, :, :] = x0      # x0 = array([1., 0.])

I_tp_pos = np.zeros((L, 2))
for k in range(i0 + 1, L):
  I_tp_pos[k] = I_tp_pos[k-1] + (f_tp[k-1] + f_tp[k]) / 2.0 * h

Constante de Lipschitz exata

Para $F(t, \mathbf{x}) = (x_2,\; -x_1 + 2\operatorname{sen} t)$:

Portanto $C(\delta) = 1$ para todo $\delta > 0$. O campo $F$ é afim em $\mathbf{x}$ com parte linear $A\mathbf{x}$, onde

é uma rotação e, portanto, uma isometria — sua norma espectral é $\|A\|_2 = 1$. Já

com $M(\delta) \to 1$ quando $\delta \to 0^+$.

Com $C(\delta) = 1$ e $\varepsilon = 1$, o produto $C\varepsilon = 1$, e a cota (6) reduz-se a $M(\delta^*)/n!$, que decresce monotonamente desde o primeiro iterado — sem o período inicial de piora observado no Objetivo 4. O sistema bidimensional é, portanto, intrinsecamente mais favorável ao método de Picard.

O desempenho do método de Picard depende diretamente da estrutura algébrica do campo $F$. A constante de Lipschitz $C(\delta)$ determina tanto o tamanho do domínio acessível quanto a velocidade de convergência, via a cota $M(\delta)(C(\delta)\varepsilon)^n/n!$. Quanto mais próximo de linear for $F$ em $x$, mais favorável é o método.

1. O trapézio é sempre preferível ao retângulo

O trapézio tem ordem $\mathcal{O}(h^2)$ frente ao $\mathcal{O}(h)$ do retângulo, resultando num piso de erro irredutível mais baixo que não pode ser eliminado aumentando $N$. No Objetivo 2, o piso do retângulo estabiliza em $\approx 6{,}93 \times 10^{-4}$ enquanto o trapézio atinge $\sim 10^{-7}$ — quatro ordens de diferença para os mesmos parâmetros.

2. Convergência além do raio teórico para F = x²

O teorema garante convergência apenas para $\varepsilon \le \varepsilon^* = 0{,}25$, mas os iterados convergem numericamente em todo $|t| < 1$, porque coincidem com as somas parciais da série geométrica $\sum t^k$. Para $|t| > 1$ a divergência é verificável executando o código do Objetivo 4 com $\varepsilon > 1$.

3. Otimização computacional — fator 2500×

A substituição da integração direta pela acumulação incremental reduz a complexidade de $\mathcal{O}(NL^2)$ para $\mathcal{O}(NL)$, tornando viável trabalhar com $L = 20001$ pontos. Para $N = 30$: de $\approx 3 \times 10^9$ para $\approx 1{,}2 \times 10^6$ operações.

4. O sistema 2D é intrinsecamente mais favorável

Com $C(\delta) = 1$ constante — reflexo de a parte linear ser uma isometria — o produto $C\varepsilon = 1$ garante que a cota decai monotonamente desde o primeiro iterado, sem o período de piora inicial do Objetivo 4, onde $C\varepsilon \approx 3{,}96$. Além disso, a solução é inteira em $\mathbb{R}$: não há limitação de domínio análoga ao teto $|t| < 1$.

5. Equilíbrio entre N e L

Os dois termos do erro total têm dependências opostas: aumentar $N$ reduz o erro de convergência exponencialmente via o fatorial; aumentar $L$ reduz o erro de quadratura como $L^{-p}$. O valor $L$ mínimo efetivo para equilíbrio, dado $N$ fixo, satisfaz

Na prática, para $N = 14$ e trapézio, $L \sim 2000$ já equilibra os dois termos.

Comparativo final

Em geral, problemas cujo campo $F$ admite uma decomposição $F(t,\mathbf{x}) = A\mathbf{x} + G(t)$ com $\|A\|_2$ pequena e $G$ uniformemente limitada são intrinsecamente mais favoráveis ao método: a cota decai mais rapidamente e domínios maiores são acessíveis sem aumentar $N$. Há, ainda assim, uma limitação numérica ao aumentar o domínio em ambos os casos, pelo acúmulo de erro do método numérico — erro que pode ser controlado com a otimização dos parâmetros.

Sumário dos cinco objetivos

Referência

M. Viana e J. Espinar, Equações Diferenciais Ordinárias, IMPA, Rio de Janeiro, 2025.