fitting a variable in a system of ODE as a function of time

[pabloacera]

Hi,

Thanks for developing this awesome package! I want to solve an inverse problem modelled as a system of ODE evolving through time. I want to calculate several parameters that also depend on time. My problem is similar to the Lorenz system example but some of the parameters I want to estimate are as a function of time, so instead of one value they will have a different value per time point, and I have no idea the shape or nature of that function. Please, Is there a way to do this with your package? if so, what would be the best way of doing this? Thanks a lot for your time! Cheers, Pablo.

[shuheng-liu]

Hi Pablo, can you provide a minimal example of what you're trying to solve so that we can have a more specific context?

[shuheng-liu]

If I understand correctly, you are trying to solve ODEs w.r.t. u(t) involving a time-dependent coefficient (say θ(t)), where both the solution u(t) and the coefficient θ(t) are unknown. This is essentially the same as solving ODEs w.r.t. two unknown functions (u(t), θ(t)). If this is the case, the Lotka-Volterra equation in the README is a good example.

(If you don't want to enforce any initial condition for θ, you can put a neurodiffeq.conditions.NoCondition() instance as a placeholder.

[pabloacera]

Hi Shuheng-liu,

Thanks for the quick response. I can provide the systems of ODE that I am trying to solve.
I have a dataset with u and s over time and I was wondering if I could use this package to find the other variables through time. Thanks a lot! any insights would be very valuable.

[shuheng-liu]

Hi Pablo,

Assuming your data is free of noise, you can recover u(t) and s(t) by interpolating discrete data points (ti, ui) and (ti, si). You can also recover the derivatives du/dt and ds/dt if your interpolant is differentiable (say Lagrangian polynomial).

This way, the problem is essentially solving an ODE system w.r.t. three unkonwn functions alpha(t), beta(t), and gamma(t).

from neurodiffeq import diff
from neurodiffeq.solvers import Solver1D
from neurodiffeq.conditions import NoCondition
import matplotlib.pyplot as plt

u_interp = lambda t: ... # some interpolant function (note: t is a torch.Tensor)
s_interp = lambda t: ... # some interpolant function

def eq(alpha, beta, gamma, t):
    u, s = u_interp(t), s_interp(t)
    dudt, dsdt = diff(u, t), diff(s, t)

    return [
        alpha - beta * u - dudt,
        beta * u - gamma * s - dsdt,
    ]

conditions = [
    NoCondition(),  # no constraint for alpha
    NoCondition(),  # no constraint for beta
    NoCondition(),  # no constraint for gamma
]

T_MIN, T_MAX = 0, 1

solver = Solver1D(
    eq, 
    conditions,
    t_min=T_MIN,
    t_max=T_MAX,
    n_batches_valid=1,  # [Optional] to speed up computation
)

solver.fit(max_epochs=10)

solution = solver.get_solution(best=True, copy=True)

t = np.linspace(T_MIN, T_MAX, 100)
alpha, beta, gamma = solution(t, to_numpy=True)

plt.plot(t, alpha, label='alpha')
plt.plot(t, beta, label='beta')
plt.plot(t, gamma, label='gamma')
plt.legend()

If you have constraints for alpha, beta, or gamma, you can change the corresponding entry in conditions to neurodiffeq.conditions.IVP(..., ...) or neurodiffeq.conditions.DirichletBVP(..., ...). See documentation for details.