0

I have defined a function in Python 2.7.15 which takes two inputs (b,v), but I have noticed that f(50,0.1) produces very different results to those of f(50.0,0.1). This is the function:

def f(b,v):
    h=b*v/math.sqrt(1-v**2)
    def dJ_dp(J,p):
        return [J[1],-J[0]+3.0/2*J[0]**2+1/(2*h**2)]
    J0 = [0.0000001,1/b]
    ps = np.linspace(0,15,50)
    Js = odeint(dJ_dp, J0, ps)
    us = Js[:,0]
    return (ps,1/us)

I've needed to define dJ_dp inside f(b,v) because it needs the value h. Why are the outputs so different? Why are they different at all?

My hypotheses was that something went wrong when defining h but that doesn't seem to be the case.

Jonathan Hall
  • 75,165
  • 16
  • 143
  • 189
tr416
  • 123
  • 2
  • Python 2.7 defines division differently for reals and integers. [Take a look](https://stackoverflow.com/questions/21316968/division-in-python-2-7-and-3-3) – Ma0 Mar 20 '19 at 08:32
  • the second entry of ```J0 = [0.0000001,1/b]``` is an integer division if b is an int, i.e. it's result will be 0 instead of 0.02 – mrks Mar 20 '19 at 08:32

1 Answers1

1

The problem is probably here: J0 = [0.0000001,1/b].

If b is the int 50, 1/b will be done using integer division, and result in 0. If b is the floating point 50.0 it will be done with floating point division and will result in 0.02.

You could use 1.0 instead of 1 to force floating point arithmetic:

J0 = [0.0000001, 1.0/b]
# Here -----------^
Mureinik
  • 297,002
  • 52
  • 306
  • 350