mattja / sdeint

Numerical integration of Ito or Stratonovich SDEs
GNU General Public License v3.0
154 stars 25 forks source link

BUG: initial conditions' type are not interpreted correctly #20

Closed arashgmn closed 2 years ago

arashgmn commented 2 years ago

There's a bug that I suppose stems from the type interpretation of initial conditions. The following code which simulates a simple damped harmonic oscillator (with no noise at all) illustrates the error:

import numpy as np
import matplotlib.pyplot as plt
import sdeint

w0 = np.pi  # arbitrary
beta = 0.5  # arbitrary
tspan = np.linspace(0.0, 5.0, 5001)

x0 = 1, 1    # types interpreted as int (wrong solution) 
x1 = 1, 1.   # types interpreted as int (wrong solution)
x2 = 1., 1   # types interpreted as float (correct solution)

def f(x, t):
    return np.array([x[1], -w0*x[0]-beta*x[1]])

def g(x, t):
    return np.diag([0, 0])

result0 = sdeint.itoint(f, g, x0, tspan)
result1 = sdeint.itoint(f, g, x1, tspan)
result2 = sdeint.itoint(f, g, x2, tspan)

# vis the results
for i in range(3):
    l = plt.plot(tspan, eval('result'+str(i))[:,0], label=str(I), linewidth=5-2*i)
    plt.plot(tspan, eval('result'+str(i))[:,1], '--', color=l[0].get_color(), linewidth=5-2*i)
plt.legend()
plt.show()

Might be related to #15 and #16.

mattja commented 2 years ago

This is not a bug. You are supposed to provide initial conditions of the same type as your state space. (e.g. scalar / vector / vector space over complex field / whatever)

To make the software easier to use, I'll add a special case to assume float64 if you provide integers.

arashgmn commented 2 years ago

Indeed, what I mean by a bug is inconsistent behavior of type interpretation. i.e., if the first initial condition (IC) is float, all the rest, too, will be considered as float even if they are not explicitly declared as such (case 3 above). Whereas if the first IC is int, then even an explicitly float-declared second IC will be interpreted as int (case 2).

Overall, instead of a particular treatment for integer-typed IC, I suggest converting all the ICs to float by default to avoid weird type interpretation by python. In most SDEs, int is a bit unnatural anyway.

PS: I also changed the linewidth of the snippet above to avoid full line overlaps.