martenlienen / torchode

A parallel ODE solver for PyTorch
https://torchode.readthedocs.io
MIT License
228 stars 13 forks source link

Problem with integration of simple ode system #33

Closed kschischo closed 9 months ago

kschischo commented 9 months ago

Thanks for developin torchode. I have played around a bit and trie to solve the Lotka-Volterra model, just to get used to torchode. However, the solver doesn't find any solution.

!/usr/bin/env python3

-- coding: utf-8 --

""" Created on Tue Jan 23 19:17:47 2024

@author: kschisch """

import matplotlib.pyplot as plt import torch from torchode import Status, solve_ivp

def vdp(t, y, mu): x, xdot = y[..., 0], y[..., 1] return torch.stack((xdot, mu * (1 - x*2) xdot - x), dim=-1)

def lotka_volterra(t,x, mu): x1, x2 = x[..., 0], x[..., 1] return torch.stack((-x1+x1 x2, x2-x1 x2),dim=-1)

batch_size, mu = 5, 10.0 y0 = torch.randn((batch_size, 2)) t_eval = torch.linspace(0.0, 100.0, steps=50)

sol = solve_ivp(vdp, y0, t_eval, method="tsit5", args=mu)

sol = solve_ivp(lotka_volterra, y0, t_eval, method="tsit5",args=mu)

print(sol.status) # => tensor([0, 0, 0, 0, 0]) assert all(sol.status == Status.SUCCESS.value) print(sol.stats)

martenlienen commented 9 months ago

What do you mean with does not find a solution? Is it a wrong solution? Does it fail?

kschischo commented 9 months ago

Dear Marten, it just didn't finish for some hours. I use pytorch under conda with cuda installed. Many thanks, Maik

martenlienen commented 9 months ago

Please put a print(t) statement in your lotka_volterra to look at the solution progress.

kschischo commented 9 months ago

Ok, I am sending the output below. It seems that the last entry in the sol object is zero.

I ran

python test_torchode.py > lv_out_torchode.tx lv_out_torchode.txt

Below is the code again. I have tried to change the parameters of the ode, just in case. But this is not the problem.

kschischo commented 9 months ago

!/usr/bin/env python3

-- coding: utf-8 --

""" Created on Tue Jan 23 19:17:47 2024

@author: kschisch """

import matplotlib.pyplot as plt import torch from torchode import Status, solve_ivp

def vdp(t, y, mu): x, xdot = y[..., 0], y[..., 1] return torch.stack((xdot, mu * (1 - x*2) xdot - x), dim=-1)

def lotka_volterra(t,x, mu): x1, x2 = x[..., 0], x[..., 1]
print(t) return torch.stack((-3x1+x1 x2, 1.5x2-x1 x2),dim=-1)

batch_size, mu = 5, 10.0 y0 = torch.randn((batch_size, 2)) t_eval = torch.linspace(0.0, 100.0, steps=50)

sol = solve_ivp(vdp, y0, t_eval, method="tsit5", args=mu)

sol = solve_ivp(lotka_volterra, y0, t_eval,args=mu)

print(sol.status) # => tensor([0, 0, 0, 0, 0]) assert all(sol.status == Status.SUCCESS.value) print(sol.stats)

kschischo commented 9 months ago

I think the mistake is on my side. I used randn for the initial conditions, which yields negative values. Then, the solution can go to infinity, which leads to nans in the state. Many thanks again. Sorry for not seeing this trivial mistake.

It works fine with uniformly sampled initial conditions, where only positive values are allowed. Negative rabbits don't make sense, obviously.

martenlienen commented 9 months ago

No worries :)