Python para Equações Diferenciais

Neste capítulo, o Python será usado para aproximar soluções de problemas de valor inicial. Os métodos de Euler, Runge-Kutta e Adams-Bashforth repetem atualizações a partir de um ponto inicial. Como funções, laços e listas já foram apresentados na introdução ao Python, aqui o foco será organizar a malha de cálculo, armazenar a trajetória da solução, encapsular passos numéricos e comparar resultados com solve_ivp.

Definindo a EDO

Uma equação diferencial ordinária de primeira ordem pode ser escrita como:

Em Python, representamos \(f(x,y)\) como uma função que recebe o ponto atual e devolve a inclinação da solução naquele ponto:

def f(x, y):
    return 2*x - y + 1

print(f(0.2, 1))

Valores iniciais e malha

Um problema de valor inicial informa o ponto inicial, o valor inicial da solução e o passo de cálculo:

import numpy as np

x0 = 0
y0 = 1
h = 0.2
xf = 1

n_passos = int(round((xf - x0) / h))
malha_x = np.linspace(x0, xf, n_passos + 1)

Usar uma malha com quantidade definida de passos evita problemas de comparação entre números decimais que podem ocorrer quando controlamos o laço apenas com while x < xf.

Método de Euler em forma de função

A fórmula do método de Euler é:

Em Python, podemos escrever um passo do método como uma função independente. Isso facilita trocar Euler por Runge-Kutta mantendo a mesma estrutura de repetição:

def euler(x, y, h):
    return y + h*f(x, y)

print(euler(0.2, 1, 0.2))

Repetindo os passos com histórico

Para avançar de \(x_0\) até \(x_f\), guardamos a trajetória calculada. Em vez de depender de uma comparação contínua com \(x_f\), percorremos a quantidade de passos definida pela malha:

valores_x = [x0]
valores_y = [y0]

x = x0
y = y0

for passo in range(n_passos):
    y = euler(x, y, h)
    x = x + h

    valores_x.append(x)
    valores_y.append(y)

Tabela de resultados

Para conferir o método, podemos reunir os valores calculados em pares \((x_i,y_i)\). Essa tabela será útil para comparar Euler, Runge-Kutta e Adams-Bashforth:

tabela = list(zip(valores_x, valores_y))

for xi, yi in tabela:
    print(xi, yi)

Preparando métodos de ordem maior

Métodos como Runge-Kutta calculam inclinações intermediárias antes de atualizar \(y\). A estrutura abaixo mostra como guardar esses coeficientes em variáveis separadas:

def rk2(x, y, h):
    k1 = f(x, y)
    k2 = f(x + h, y + h*k1)

    return y + h/2 * (k1 + k2)

Com essa organização, a função que executa os passos pode receber o método como argumento:

def resolver_pvi(metodo, x0, y0, h, n_passos):
    valores_x = [x0]
    valores_y = [y0]

    x = x0
    y = y0

    for passo in range(n_passos):
        y = metodo(x, y, h)
        x = x + h
        valores_x.append(x)
        valores_y.append(y)

    return valores_x, valores_y

Usando SciPy

Em algumas situações, podemos comparar os métodos estudados com uma função pronta da biblioteca scipy. A função solve_ivp recebe o intervalo, o valor inicial e os pontos nos quais desejamos avaliar a solução:

from scipy.integrate import solve_ivp

solucao = solve_ivp(f, [x0, xf], [y0], t_eval=malha_x)

print(solucao.t)
print(solucao.y[0])

Atividade

Para \(f(x,y)=x+y\), defina x0 = 0, y0 = 1, h = 0.1 e xf = 0.5. Calcule n_passos, gere a malha e use uma função resolver_pvi para guardar todos os pares \((x_i,y_i)\) obtidos pelo método de Euler.