STLSQ na Prática

import numpy as np
from scipy.integrate import solve_ivp

a, b, c, d = 1.1, 0.4, 0.4, 0.1
def lotka_volterra(t, z):
    x, y = z
    return [a*x - b*x*y, -c*y + d*x*y]

t_eval = np.linspace(0, 60, 3000)
sol = solve_ivp(lotka_volterra, (0,60), [10,5], t_eval=t_eval, rtol=1e-10, atol=1e-10)
X = sol.y.T
dX = np.gradient(X, t_eval[1]-t_eval[0], axis=0)

x, y = X[:,0], X[:,1]
Theta = np.column_stack([np.ones_like(x), x, y, x**2, x*y, y**2])
nomes = ["1","x","y","x^2","xy","y^2"]

def stlsq(Theta, dX, thresh=0.05, n_iter=10):
    Xi = np.linalg.lstsq(Theta, dX, rcond=None)[0]
    for _ in range(n_iter):
        small = np.abs(Xi) < thresh
        Xi[small] = 0
        for c in range(dX.shape[1]):
            big = ~small[:,c]
            Xi[big,c] = np.linalg.lstsq(Theta[:,big], dX[:,c], rcond=None)[0]
    return Xi

Xi = stlsq(Theta, dX)

for nome, linha in zip(nomes, Xi):
    print(f"{nome:>4}: dx/dt={linha[0]:+.3f}  dy/dt={linha[1]:+.3f}")

print(f"verdadeiro: a={a} b={-b} c={-c} d={d}")

import matplotlib.pyplot as plt

def modelo(t, z):
    x, y = z
    th = np.array([1, x, y, x**2, x*y, y**2])
    return [th @ Xi[:,0], th @ Xi[:,1]]

sol2 = solve_ivp(modelo, (0,60), [10,5], t_eval=t_eval, rtol=1e-8, atol=1e-8)
fig, ax = plt.subplots()
_ = ax.plot(sol.y[0], sol.y[1], "-", lw=2, label="verdadeiro")
_ = ax.plot(sol2.y[0], sol2.y[1], "--", lw=2, label="descoberto")
_ = ax.legend()
plt.tight_layout(); plt.show()