Matemática Computacional

Resumo da matéria [Versão 3.1]

[Capítulo 1] [Capítulo 2] [Capítulo 3] [Capítulo 4] [Capítulo 5]

Capítulo V

Métodos Numéricos para EDO's

Começamos por apresentar uma breve introdução

Equações Diferenciais Ordinárias (EDO's)

O estudo de equações diferenciais é iniciado com Euler, sendo o caso mais simples, de equações diferenciais ordinárias, sistematizado por Cauchy e também por Picard. É no final do séc. XIX, com Runge (ver também Kutta) que o estudo de métodos numéricos é desenvolvido, levando aos métodos de Runge-Kutta e depois já no início do séc.XX aos métodos de Adams Adams.
Encontrar métodos eficazes para a solução de grandes sistemas de EDO's é ainda um problema actual, com dificuldades ao nível da estabilidade do métodos.

O estudo de EDO's pode ser concentrado no caso de sistemas de EDO's de 1ª ordem, no problema de Cauchy:

╭
│	x'(t) = f(t, x(t))
┤
│	x(t₀) = x₀
╰

em que x e f são funções escalares ou vectoriais.
A incógnita é a função x que será diferenciável, e está definida num instante inicial t₀ onde se conhece o seu valor x₀. O objectivo será determinar a função x(t) para t > t₀.

A primeira observação que fazemos, é que esta atribuição no instante inicial permite conhecer também a derivada nesse mesmo instante inicial t₀, pois

x'(t₀) = f(t₀,x(t₀)) = f(t₀,x₀).

Sendo possível obter x''(t), x'''(t), ... a partir da expressão de x'(t) é então intuitivo que será possível obter todas as derivadas em t₀.
Isto significa que, enquanto a solução x for uma função analítica em torno de t₀, é possível obter o seu desenvolvimento de Taylor, a partir do cálculo das derivadas nesse ponto inicial. Este é o argumento inicial de Cauchy para a existência e unicidade de solução (local).

Exemplo:

╭
│	x'(t) = x(t)²
┤
│	x(0) = a
╰

Neste caso f(t,x) = x², e t₀=0, x₀=a.
Reparamos que

x''(t) = 2x'(t)x(t)= 2x(t)³; ...; x^(k)(t) = k! x(t)^k+1

ou seja, a expansão de Taylor de x em t₀=0, dá

x(t) =

∞
∑
k=0

x^(k)(0) t^k/k! =

∞
∑
k=0

x(0)^k+1 t^k = a

∞
∑
k=0

(a t)^k =

1 − a t

quando t < 1/a. É fácil confirmar que é solução, e a sua construção pode mostrar que é única.

Mais importante é notar que esta solução Explode quando t=1/a.
Ou seja, ainda que f seja analítica, a solução x pode estar definida apenas num pequeno intervalo, próximo do instante inicial.
Por outro lado, a função f pode não ser analítica, e por isso é conveniente estabelecer um resultado mais geral de existência e unicidade local:

Teorema: Seja (t₀,x₀)∈ D, em que D é um domínio onde f é contínua e Lipschitziana na segunda variável, ou seja, existe L>0:

|f(t,x) − f(t,y)| ≤ L|x − y|, ∀ (t,x),(t,y)∈ D.

Então existe uma vizinhança de t₀ onde a solução x, do problema de Cauchy, existe e é única.

dem: Resulta do teorema do ponto de fixo em espaços funcionais (espaços de Banach, eg. [1])
Convém apenas referir que a condição Lipschitz pode ser simplificada, no caso de funções diferenciáveis, para

|∂₂ f(t,x)| ≤ L, ∀ (t,x)∈ D

Observações: (dependência numa só das variáveis)
Há dois casos de EDO's que podem ser tratados como problemas de integração:
(1) Trivialmente, quando f(t,x) = f(t), ou seja, a função só depende de t, então a solução reduz-se ao integral

x'(t) = f(t) ⇒ x(t) = ∫_[t₀,t] f(t) dt

e basta aplicar a primitivação ou regras de quadratura.
(2) Quando f(t,x)=f(x), a equação diferencial x'(t) = f(x(t)) diz-se autónoma, e a solução também pode ser obtida por primitivação.
Com efeito, sendo G=∫ 1/f temos

x'(s) = f(x(s)) ⇔ x'(s)G'(x(s)) = 1 ⇒ ∫_[t₀,t] x'(s)G'(x(s)) ds = t − t₀

e como (G(x(s)))' = x'(s)G'(x(s)), obtemos

G(x(t)) − G(x₀) = t − t₀ ,

uma equação não linear, que define implicitamente o valor de x(t).

Exemplo: Retomando o exemplo anterior, x'(t) = x(t)², temos

G(x)= ∫ x^-2 dx = -x^-1

e assim

-x(t)^-1 + a^-1 = t − 0 ⇒ x(t) =

1 − a t

obtendo por outra forma a solução encontrada pela expansão em série de Taylor.

Observações: (EDO's de Ordem Superior e Sistemas)
Quando as equações diferenciais ordinárias são de ordem superior à primeira, é possível reduzir o problema a um sistema de EDO's de primeira ordem.
Se tivermos

u^(p) = F ( t, u(t), u'(t), ⋯, u^(p-1)(t) )

definimos uma variável vectorial

x = (x₀, x₁, ⋯, x_p-1) := (u, u', ⋯, u^(p-1))

que tem os valores atribuídos inicialmente pelos dados de Cauchy

x(t₀) = ( u(t₀), u'(t₀), ⋯, u^(p-1)(t₀) )

A EDO é então reescrita em termos de um sistema usando essa mudança.

x₀'(t) = x₁(t)

x₁'(t) = x₂(t)

⋮

x_p-2'(t) = x_p-1(t)

x_p-1'(t) = F ( t, x₀(t), x₁(t), ⋯, x_p-1(t) )

ou ainda, abreviadamente

x'(t) = F( t, x(t) ).

Exemplo: quando os coeficientes são constantes, a equação

u'''(t) = a u(t) + b u'(t) + c u''(t)

pode escrever-se num sistema

x₀'(t) = x₁(t)

x₁'(t) = x₂(t)

x₂'(t) = a x₀(t) + b x₁(t) + c x₂(t)

que neste caso (simples) é mesmo um sistema linear

x'(t) =

┌	x₀'(t)	┐
│	x₁'(t)	│
└	x₂'(t)	┘

┌	0	1	0	┐
│	0	0	1	│
└	a	b	c	┘

┌	x₀(t)	┐
│	x₁(t)	│
└	x₂(t)	┘

= M x(t)

e a sua solução é dada através de uma exponencial matricial

x(t) = x(0) e^Mt.

A exponencial matricial pode ser definida pelos valores próprios...
Nesse caso isso é fácil de obter pois a matriz M tem na última linha os coeficientes que definem o polinómio característico.
Os valores próprios serão as (3) raízes complexas λ_k desse polinómio

λ³ = a + bλ + cλ².

Assim, quando a matriz é diagonalizável, podemos mesmo obter

u(t) = C₁ e^tλ₁ + C₂ e^tλ₂ + C₃ e^tλ₃

em que os coeficientes C_k são obtidos a partir dos dados iniciais de Cauchy.

Métodos de Taylor e Runge-Kutta

No que se segue iremos formular os métodos em termos escalares, mas a sua aplicação a sistemas é exactamente igual, encarando x e f como vectores.

Iremos começar por abordar métodos que podem ser deduzidos por expansões em série de Taylor − por isso são denominados Métodos de Taylor.
A partir dos métodos de Taylor deduzimos uma variante − os Métodos de Runge-Kutta − que evitam o cálculo de derivadas para ordens p>1.

Método de Euler

Começamos por introduzir o método numérico mais simples, que pode ser deduzido pela aproximação de Taylor de 1ª ordem (x ∈ C²)

x(t+h) = x(t) + h x'(t) + x''(ξ)/2 h² = x(t) + h f(t,x(t)) + O(h²).

Isto significa que é possível prever a evolução num instante posterior, se desprezarmos o termo O(h²), o que é aceitável para valores h pequenos.

Dessa forma, quando queremos avaliar a solução num tempo T podemos subdividir a distância usando um espaçamento constante

h=(T-t₀)/N,

para um certo N, definindo pontos intermédios de cálculo

t₀, t₁ = t₀ + h, ⋯, t_N=t₀ + N h = T.

Sendo conhecido o valor x₀=x(t₀), podemos aplicar sucessivamente a aproximação de 1ª ordem,

x(t_m+1) ≈ x(t_m) + h f(t_m, x(t_m))

e assim definir o Método de Euler

Dado x₀, iteramos:

x_m+1 = x_m + h f_m

abreviando-se f_m = f(t_m, x_m).
Cada valor x_m é uma aproximação de x(t_m), sendo claro que x_N será a aproximação final de x(T).

Exemplo: Consideremos a EDO x'(t)=x(t) com x(0)=1.
É bem conhecido que x(t)=e^t é a sua solução única.
Supondo que pretendemos aproximar x(T)=e^T num instante fixo T=N h, usando o Método de Euler, obtemos o seguinte:

x₀=1; x_m+1 = x_m + x_m h ⇒ x_m=(1+h)^m

Portanto, como h=T/N reencontramos uma aproximação clássica:

x_N = (1 + T/N)^N → e^T (quando N→∞, ou h → 0).

Outros Métodos de Taylor

O Método de Euler é um caso particular de métodos que se baseiam na expansão em série de Taylor.
Considerando uma expansão de 2ª ordem (x ∈ C³)

x(t+h) = x(t) + h x'(t) + x''(t)/2 h² + O(h³).

Podemos desprezar o termo O(h³), obtendo uma aproximação melhor, de segunda ordem.
No entanto, se a substituição x'(t)=f(t,x(t)) é imediata, o cálculo de x''(t) envolve já uma derivação adicional.
Isso pode ser feito recorrendo directamente à derivação da expressão x'(t), e procedendo à eventual substituição de x'(t) e depois dos valores t por t_m e x(t) por x_m.
Abreviadamente, o Método de Taylor de ordem 2 fica

Dado x₀, iteramos:

x_m+1 = x_m + h f_m + ¹/₂ h² [x''(t)]|_{t=t_m; x(t)=x_m} .

Exemplo: Consideremos a aplicação do Método de Taylor de ordem 2 à equação

x'(t)= -t x(t), com x(0)=1.

Como x''(t)= -x(t) − t x'(t) = (t² − 1)x(t), obtemos

x_m+1 = x_m + h (-t_mx_m) + ¹/₂ h² (t_m² − 1)x_m ,

fazendo as substituições indicadas.

O cálculo de x''(t) pode ainda ser obtido pela regra em cadeia, recorrendo às derivadas parciais

x''(t) = d/dt f(t,x(t)) = ∂₁f(t,x(t)) + f(t,x(t)) ∂₂f(t,x(t)).

e no ponto t_m, abreviamos a notação

[x''(t)]_m = (∂₁f)_m + f_m (∂₂f)_m.

No entanto, este termo requer o cálculo prévio de derivadas, que pode ser dispensado usando uma técnica que leva aos métodos de Runge-Kutta.

Métodos de Runge-Kutta

Os métodos de Runge-Kutta de ordem p partem de um método de Taylor de ordem p, procurando substituir as derivadas parciais por outras avaliações pontuais de f.
No caso mais simples, o método de Euler, com p=1 não há derivação a realizar, e por isso também pode ser encarado como um Runge-Kutta trivial, de ordem 1.

Podemos escrever os métodos de Taylor de ordem 2 na forma

x_m+1 = x_m + h Φ_m ,
com Φ_m = f_m + ^h/₂ ( (∂₁f)_m + f_m (∂₂f)_m )

e a ideia será substituir as derivadas no cálculo de Φ_m.
Consideramos a expansão de Taylor de f em duas variáveis,

f(t+H₁, x+H₂) = f(t,x) + ∂₁f(t,x)H₁ + ∂₂f(t,x)H₂ +O(H₁H₂).

Podendo desprezar o termo O(H₁H₂), obtemos a aproximação

f(t_m+H₁, x_m+H₂) ≈ f_m + (∂₁f)_mH₁ + (∂₂f)_mH₂.

Quando H₁= ^h/₂ , H₂= f_m ^h/₂ , temos

f(t_m+ ^h/₂ , x_m+f_m ^h/₂ ) ≈ Φ_m ,

com um erro de O(h²).
Substituindo Φ_m obtemos assim um Método de Runge-Kutta de ordem 2 (do ponto-médio)

Dado x₀, iteramos:

x_m+1 = x_m + h f ( t_m+ ^h/₂ , x_m + f_m ^h/₂ ) .

Métodos de Runge-Kutta de ordem 2 (caso geral)
Variando um parâmetro θ≠ 0, e combinando com f_m, obtemos
(1-θ) f_m+ θ f(t_m+H₁, x_m+H₂) ≈ f_m + θ ( (∂₁f)_mH₁ + (∂₂f)_mH₂ ) ,
resultando noutras aproximações para Φ_m:
(1-θ) f_m+ θ f(t_m+H₁, x_m+H₂) ≈ Φ_m ,

quando θ H₁= ^h/₂ , θ H₂= f_m ^h/₂ .
Obtemos assim uma forma geral

x_m+1 = x_m + (1-θ) h f_m + θ f ( t_m + ^h_{/_2θ} , x_m + f_m ^h_{/_2θ} )
que podemos concretizar:
- O caso θ = 1 já foi visto (ponto-médio);
- O caso θ = 1/2 também é conhecido como método de Heun:
  x_m+1 = x_m + ¹/₂ h f_m + ¹/₂ f ( t_m + h, x_m + f_mh )
- Nota: estes métodos também são por vezes designados como métodos de Euler modificados.
  Cada passo necessita de 2 cálculos da função f.
Métodos de Runge-Kutta de ordem 3
Tal como no caso anterior são obtidos a partir do Método de Taylor de ordem 3, e cada passo necessita de 3 cálculos da função f. Não são tão utilizados como os de ordem 2 ou 4.

Métodos de Runge-Kutta de ordem 4
Tal como no caso anterior são obtidos a partir do Método de Taylor de ordem 4, e exigem 4 cálculos da função f.
O exemplo mais conhecido é

x_m+1 = x_m + (F_h1 + 2F_h2 + 2F_h3 + F_h4)/6,

definindo recursivamente:

F_h1 = h f_m

F_h2 = h f ( t_m + ^h/₂ , x_m + F_h1/2 )

F_h3 = h f ( t_m + ^h/₂ , x_m + F_h2/2 )

F_h4 = h f ( t_m + h , x_m + F_h3 )

Métodos de Runge-Kutta de ordem ≥ 5
Raramente são utilizados sem adaptações. Foi demonstrado por J. C. Butcher (1960's) que para p ≥ 5 o número de avaliações é superior à ordem p. Isso torna estes métodos claramente menos eficazes, a menos que se introduzam modificações (e.g. Runge-Kutta-Fehlberg, 1969).

Erros: Consistência, Estabilidade e Convergência

Os métodos numéricos apresentados consistem em aproximações locais, em cada subintervalo de comprimento h, cujo erro local acumula sucessivamente até à aproximação no tempo final T > t₀. Esse valor final envolve já um erro global, que poderá ser encarado com soma dos erros locais, se o método for estável.
Assim há 3 noções a considerar:
(i) ordem de consistência ( ~ ordem do erro da aproximação local)
(ii) estabilidade do método ( ~ a influência local é somada globalmente)
(iii) ordem de convergência ( ~ ordem do erro global)

De um modo geral iremos definir os métodos (unipasso, explícitos) na forma

x_m+1 = x_m + h Φ_h(t_m, x_m)

onde Φ_h é uma função que define o método:
- no caso de Euler: Φ_h = f;
- no caso do RK-2-ponto-médio: Φ_h(t,x) = f(t + h/2, x + h f(t,x)/2);
Iremos assumir que tal como f, a função Φ é também contínua e Lipschitziana na 2ª variável.

Ordem de Consistência

A ordem de consistência é definida localmente, através do erro no cálculo de x_m+1, admitindo que os valores anteriores x_m são exactos.
Assim, recuperando a expansão de Taylor, para o método de Euler temos

E(h) = x(t_m+1) − x_m+1 = ( x(t_m) + h x'(t_m) + x''(ξ_m)/2 h² ) − ( x_m − h f_m ) = x''(ξ_m)/2 h²,

assumindo que x(t_m)=x_m.
Há assim um erro local de ordem 2 − designado erro de truncatura local:

E(h) = x''(ξ_m)/2 h² = O(h²)

e a ordem de consistência é definida um valor abaixo, ou seja, neste caso a ordem de consistência é 1.

De forma semelhante, para outros métodos de Taylor de ordem p, obtemos

E(h) = x^(p+1)(ξ_m)/(p+1)! h^p+1 = O(h^p+1)

havendo um erro de truncatura local de ordem p+1 e uma correspondente ordem de consistência p.

Observação: No caso dos métodos de Runge-Kutta de ordem p haverá também um erro de truncatura O(h^p+1) e uma correspondente ordem de consistência p.
Vejamos para o caso do Método R-K de ordem 2 do ponto-médio.
Designando por X_m+1 o valor dado pelo método de Taylor de ordem 2, e por x_m+1 o valor dado pelo R-K, recordamos que

X_m+1 − x_m+1 = h O(H₁H₂) = h O(h²) = O(h³).

Portanto, como o método de Taylor verifica x(t_m+1) − X_m+1 = O(h³) temos

x(t_m+1) − x_m+1 = (x(t_m+1) − X_m+1) + (X_m+1 − x_m+1) = O(h³).

Estabilidade e Convergência

Há várias noções de estabilidade.
Interessa-nos aqui a noção de estabilidade-zero associada à limitação da propagação de erro que resulta da função Φ, que define o método, ser Lipschitziana.
No contexto dos métodos unipasso, que vimos, isso é suficiente para assegurar a estabilidade. Iremos ver que isso não é suficiente no caso de métodos multipasso (onde o cálculo de x_m+1 não depende apenas de x_m mas também de x_m-1, x_m-2, ...).

A ordem de convergência é dada pela ordem p do erro global

e_N = x(t_N) − x_N = O(h^p)

resultando da acumulação dos erros locais durante os N passos.
É intuitivo que, quando o erro de truncatura local é O(h^p+1), temos

e_N = x(t_N) − x_N ≈ N E(h) = N O(h^p+1) = Nh O(h^p) = O(h^p)

(pois N h = T − t₀ será constante).
Ou seja, a ordem de consistência p implicará convergência da mesma ordem p, quando garantida a estabilidade do método.

Teorema: (convergência de métodos unipasso explícitos)
Se um método tem consistência de ordem p, sendo estável, terá também ordem de convergência p.

dem:
Se admitirmos que os valores são correctos em t_m obtemos

X_m+1 = x(t_m) + h Φ_h(t_m,x(t_m)) ,

e a ordem de consistência p dá-nos

|E(h)| = |x(t_m+1) − X_m+1| ≤ C h^p+1

Por isso, decompondo o erro global

e_m+1 = x(t_m+1) − x_m+1 = (x(t_m+1) − X_m+1) + (X_m+1 − x_m+1)

basta analisar a parte

|X_m+1 − x_m+1| = | ( x(t_m) + h Φ_h(t_m,x(t_m)) ) − ( x_m + h Φ_h(t_m,x_m) ) |

|X_m+1 − x_m+1| ≤ |x(t_m) − x_m| + h |Φ_h(t_m,x(t_m)) − Φ_h(t_m,x_m) ) | ≤ |e_m| + h L_Φ|e_m|

concluindo-se

|e_m+1| ≤ C h^p+1 + (1 + h L_Φ)|e_m|

escrevendo K=1 + h L_Φ resulta recursivamente

|e_N| ≤ C (1 + K + ⋯ + K^N-1) h^p+1 + K^N|e₀|

Usando a progressão aritmética, e exp(L_Φh) ≥ K obtemos

|e_N| ≤ C

K^N − 1

K − 1

h^p+1 + K^N|e₀| = C

exp(L_ΦNh) − 1

L_Φ

h^p + exp(L_ΦNh)|e₀|

Como o erro inicial é nulo, conclui-se a fórmula

|e_N| ≤ C

exp(L_Φ(T-t₀)) − 1

L_Φ

h^p = O(h^p)

Corolário: No Método de Euler, temos a estimativa

|e_N| ≤ C

exp(L(T-t₀)) − 1

L_Φ

h + exp(L(T-t₀))|e₀|

em que L é a constante de Lipschitz,

L = max_D |∂₂ f|, C = max_[t₀,T] |x''|/2.

dem: É resultado da demonstração anterior, notando que no caso do método de Euler

Φ(t,x)=f(t,x)

e a consistência é p=1 pois

|E(h)| = |x(t_m+1) − X_m+1| = |x''(ξ_m)| h² /2 ≤ C h²

com C = max |x''|/2.

Métodos Implícitos

Começamos por notar que há alguns casos em que apesar da convergência estar demonstrada quando h → 0, o método de Euler pode apresentar instabilidades para valores de h que não sejam suficientemente pequenos.

Exemplo: Consideremos, por exemplo, o caso simples em que (α > 0)

x'(t) = − α x(t), com x(0) = 1.

A solução é x(t)=e^{-α t} e da aplicação do método de Euler

x_m = (1 − α h)^m.

Apesar de x(t)→ 0, quando t→∞, se o h não for suficientemente pequeno podemos ter |x_m|→∞.
Basta reparar que para h = 3/α obtemos

|x_m|=|(-2)^m| → ∞

Neste caso há uma relação directa com a constante de Lipschitz L=α e vemos assim que para valores de L mais elevados, somos forçados a considerar h muito pequenos, e consequentemente muitos passos.
A estabilidade do método está condicionada pelo valor de h e é designada por A-estabilidade.
Os métodos explícitos, que vimos, têm uma zona de A-estabilidade reduzida, condicionando a validade da aproximação apenas a valores de h suficientemente pequenos.

Para contornar este tipo de problema consideram-se métodos implícitos.
Por exemplo, o Método de Euler implícito

x_m+1 = x_m + h f(t_m+1, x_m+1)

em que o valor x_m+1 está definido implicitamente, sendo necessário resolver uma equação para o obter.

Exemplo: No exemplo apresentado antes, essa equação é fácil de resolver:

x_m+1 = x_m − h α x_m+1 ⇔ x_m+1 = x_m/(1 + α h)

ou seja,

x_m = (1 + α h)^-m.

Verificamos que independentemente do valor de α ou h a sucessão converge para zero.

Conclui-se que a aplicação de métodos implícitos pode resolver o problema de A-estabilidade.

Observação: Euler Implícito − Aplicação do Método do Ponto Fixo
Para resolver a equação que define implicitamente x_m+1 no Método de Euler implícito, podemos usar directamente o método do ponto fixo:

x_m+1 = g(x_m+1) = x_m + h f(t_m+1, x_m+1)

em que temos assim definida a função iteradora

g(z) = x_m + h f(t_m+1, z)

e uma condição para garantir convergência local será

|g'(z)|< 1 ⇔ h |∂₂f|< 1

ou ainda L h < 1, em termos da constante de Lipschitz.

A aplicação do método do ponto fixo leva à iteração:

y₀ = x_m

y₁ = x_m + h f(t_m+1,y₀) = x_m + h f(t_m+1,x_m)

y₂ = x_m + h f(t_m+1,y₁) = x_m + h f(t_m+1, x_m + h f(t_m+1,x_m))

⋱

reparando que a cada iteração isso envolve cálculos adicionais de f que podem não acrescentar nada à precisão (que é apenas de 1ª ordem) e comprometem a eficácia do método − face a considerar metade do passo h no método explícito.
Assim, normalmente, basta considerar duas ou três iterações no método do ponto fixo, usando o valor y₂ (ou y₃,...) como aproximação para x_m+1.
Para além disso, podemos começar o método implícito com uma iterada inicial melhor do que x_m, usando um método explícito. A combinação entre estes dois aspectos é conhecida como técnica preditor-corrector.

Preditor-Corrector. Métodos Multipasso.

Um outro processo de definir métodos para EDO's consiste na utilização da formulação integral equivalente. Em particular, integrando a equação diferencial obtemos

x(t_m+1) − x(t_m) = ∫_{[t_m,t_m+1]} x'(t) dt = ∫_{[t_m,t_m+1]} f(t,x(t)) dt

Sendo Q uma regra de integração aplicada ao cálculo de

I(F) = ∫_{[t_m,t_m+1]} F(t) dt

podemos fazer a aproximação de

x(t_m+1) − x(t_m) = I(f(⋅,x(⋅)))

através da regra Q considerando

x_m+1 = x_m + Q(f(⋅,x(⋅))).

Por exemplo considerando a regra dos trapézios

T(f(⋅,x(⋅))) = ( f(t_m,x(t_m)) + f(t_m+1,x(t_m+1)) ) h/2

obtemos o Método dos Trapézios

x_m+1 = x_m + ( f(t_m,x_m) + f(t_m+1,x_m+1 ) h/2.

que é ainda um método implícito.

Para a resolução da equação não linear

y = x_m + ( f(t_m,x_m) + f(t_m+1,y) ) h/2.

podemos considerar de novo o método do ponto fixo

y₀ = x_m

y₁ = x_m + (f_m + f(t_m+1,y₀)) h/2.

⋱

y_k+1 = x_m + (f_m + f(t_m+1,y_k)) h/2.

mas existe uma outra estratégia mais eficaz − a técnica Preditor-Corrector.

Técnica Preditor-Corrector: Consiste num par de métodos (explícito, implícito), em que o método explícito serve para inicializar o método implícito.
Por exemplo, considerando o par (Método de Euler, Método dos Trapézios), obtemos

y₀ = x_m + h f_m

y₁ = x_m + (f_m + f(t_m+1,y₀)) h/2 = x_m + (f_m + f(t_m+1,x_m + h f_m)) h/2.

que é a expressão de um método de Runge-Kutta de ordem 2 (de Heun).

Observação: Reparamos ainda que se considerássemos a regra do ponto-médio para a integração, teríamos

Q(f(⋅,x(⋅))) = h f(t_m+h/2, x(t_m+h/2) )

e usando a aproximação de Euler

x(t_m+h/2)≈ x_m + h/2 f_m

obtemos o método de Runge-Kutta do ponto-médio

x_m+1 = x_m + h f(t_m+h/2, x_m + h/2 f_m ).

Observação: Finalmente notamos que a aplicação de duas iteradas do método de Euler com metade do passo =h/2 conduz a algo diferente

x_m+1/2 = x_m + h/2 f_m

x_m+1 = x_m+1/2 + h/2 f(t_m+h/2, x_m+1/2) = x_m + h/2 f_m + h/2 f(t_m+h/2, x_m + h/2 f_m)

Não se trata de nenhum método Runge-Kutta, e de facto o erro local mantém-se na ordem O(h²) (apesar de decrescer 1/4), pelo que há conveniência em considerar os métodos de Runge-Kutta que permitem aumentar a ordem de convergência.

Métodos Multipasso

Até agora considerámos métodos em que o valor de x_m+1 era dado com base em x_m, independentemente de ser por forma explícita ou implícita.
Estes métodos são designados unipasso por oposição aos métodos multipasso em que o valor x_m+1 é definido a partir dos p valores anteriores:

x_m, x_m-1, ⋯, x_m-p+1

correspondente a um método a p passos. Nesta classe de métodos encontram-se os métodos de Adams

Os métodos de Adams com p-passos são deduzidos a partir da formulação integral

x(t_m+1) = x(t_m) + ∫_{[t_m,t_m+1]} f(t,x(t)) dt

considerando uma regra de integração baseada nos nós

x_m, x_m-1, ⋯, x_m-p+1 : Métodos Adams-Bashforth (explícitos)

x_m+1, x_m, x_m-1, ⋯, x_m-p+1 : Métodos Adams-Moulton (implícitos)

A regra de integração é construída pelo método dos coeficientes indeterminados, procurando que I(F)=Q(F) para F(t) monomiais 1, t, t², ⋯ . Dessa forma terá grau p no caso dos métodos explícitos, e grau p+1 no caso dos métodos implícitos.

Exemplo: (Regra de Adams-Bashforth de ordem 2)
Sendo

I(F) = ∫_{[t_m,t_m+1]} F(t) dt ; Q(F)= w₀ F(t_m) + w₁ F(t{m-1}) ;

procuramos que I(F)=Q(F) para F(t) = 1, t :

I(1) = Q(1) ⇔ h = w₀ + w₁

I(t) = Q(t) ⇔ (t_m+1)²/2 − (t_m)²/2 = w₀ t_m + w₁ t_m-1

o que dá o sistema linear

h = w₀ + w₁

h t_m + h²/2 = (w₀ + w₁) t_m − w₁ h

obtendo-se w₁= − h/2 e w₀ = 3h/2.
Ficamos assim com a regra de quadratura

Q(F) = ( 3F(t_m) − F(t_m-1) ) h/2

que origina a Regra de Adams-Bashforth de ordem 2

x_m+1 = x_m + ( 3f_m − f_m-1 ) h/2

[Capítulo 1] [Capítulo 2] [Capítulo 3] [Capítulo 4] [Capítulo 5]

C J S Alves (2008, 2009)

Matemática Computacional Resumo da matéria [Versão 3.1]

Capítulo V Métodos Numéricos para EDO's

Equações Diferenciais Ordinárias (EDO's)

Métodos de Taylor e Runge-Kutta

Método de Euler

Outros Métodos de Taylor

Métodos de Runge-Kutta

Erros: Consistência, Estabilidade e Convergência

Ordem de Consistência

Estabilidade e Convergência

Métodos Implícitos

Preditor-Corrector. Métodos Multipasso.

Métodos Multipasso

Matemática Computacional

Resumo da matéria [Versão 3.1]

Capítulo V

Métodos Numéricos para EDO's