Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Discrepancy between lsoda and standard scipy #9

Closed
edumur opened this issue Nov 4, 2021 · 1 comment
Closed

Discrepancy between lsoda and standard scipy #9

edumur opened this issue Nov 4, 2021 · 1 comment

Comments

@edumur
Copy link

edumur commented Nov 4, 2021

Hello,

I got some weird behavior using LSODA.
I compared them to scipy ode and scipy looked right while LSODA is off.
I wrote a test code to demonstrate my point:

# This Python file uses the following encoding: utf-8

import numba as nb
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import ode
from NumbaLSODA import lsoda_sig, lsoda



@nb.cfunc(lsoda_sig)
def rhs(x: float,
        i: tuple,
        di: tuple,
        p: tuple) -> None:

    # Reconstruct the complex number
    i0 = i[0] + 1j*i[1]
    i1 = i[2] + 1j*i[3]
    i2 = i[4] + 1j*i[5]

    kp     = p[0]
    ks     = p[1]
    ki     = p[2]
    ep     = p[3]
    xi     = p[4]

    delta = kp-ks-ki

    di0 = 1j*kp*ep*i1*i2*            np.exp(-1j*delta*x)/4 + 1j*kp*xi*i0/8*(  abs(i0)**2 + 2*abs(i1)**2 + 2*abs(i2)**2)
    di1 = 1j*ks*ep*i0*i2.conjugate()*np.exp( 1j*delta*x)/4 + 1j*ks*xi*i1/8*(2*abs(i0)**2 +   abs(i1)**2 + 2*abs(i2)**2)
    di2 = 1j*ki*ep*i0*i1.conjugate()*np.exp( 1j*delta*x)/4 + 1j*ki*xi*i2/8*(2*abs(i0)**2 + 2*abs(i1)**2 +   abs(i2)**2)


    di[0] = di0.real
    di[1] = di0.imag
    di[2] = di1.real
    di[3] = di1.imag
    di[4] = di2.real
    di[5] = di2.imag
funcptr = rhs.address

def vectorfield(x: float,
                i: tuple,
                p: tuple) -> None:

    # Reconstruct the complex number
    i0 = i[0]
    i1 = i[1]
    i2 = i[2]

    kp     = p[0]
    ks     = p[1]
    ki     = p[2]
    ep     = p[3]
    xi     = p[4]

    delta = kp-ks-ki

    di0 = 1j*kp*ep*i1*i2*            np.exp(-1j*delta*x)/4 + 1j*kp*xi*i0/8*(  abs(i0)**2 + 2*abs(i1)**2 + 2*abs(i2)**2)
    di1 = 1j*ks*ep*i0*i2.conjugate()*np.exp( 1j*delta*x)/4 + 1j*ks*xi*i1/8*(2*abs(i0)**2 +   abs(i1)**2 + 2*abs(i2)**2)
    di2 = 1j*ki*ep*i0*i1.conjugate()*np.exp( 1j*delta*x)/4 + 1j*ki*xi*i2/8*(2*abs(i0)**2 + 2*abs(i1)**2 +   abs(i2)**2)

    r = np.zeros(3, dtype=np.complex128)
    r[0] = di0
    r[1] = di1
    r[2] = di2

    return r

def currents_lsoda():

    i0 = 7e-3 + 0j
    i1 = 7e-5 + 0j
    i2 = 0 + 0j

    # build the initial values array
    w0 = np.array([i0.real, i0.imag,
                   i1.real, i1.imag,
                   i2.real, i2.imag], dtype=np.float64)

    # build the parameters array
    p = np.array([1507,
                  670,
                  670,
                  55,
                  19623], np.float64)

    x = np.linspace(0, 10e-3, 100)

    usol, success = lsoda(funcptr, w0, x, data=p)

    return (usol[:,0] + 1j*usol[:,1],
            usol[:,2] + 1j*usol[:,3],
            usol[:,4] + 1j*usol[:,5])


def currents_scipy():

    i0 = 7e-3 + 0j
    i1 = 7e-5 + 0j
    i2 = 0 + 0j

    # build the initial values array
    w0 = np.array([i0,
                   i1,
                   i2], dtype=np.complex128)

    # build the parameters array
    p = np.array([1507,
                  670,
                  670,
                  55,
                  19623], np.float64)

    x = np.linspace(0, 75e-3, 100)
    ys = np.array([])
    xs = np.array([])

    r = ode(vectorfield).set_integrator('zvode').set_f_params(p).set_initial_value(w0, t=x[0])

    dx = x[1]-x[0]
    while r.successful() and r.t-dx <=x[-1]:
        r.integrate(r.t+dx)
        xs = np.append(xs, r.t)
        ys = np.append(ys, r.y)
    ys = np.reshape(ys, (int(len(ys)/len(w0)), len(w0))).T

    return xs, ys



i0_l, i1_l, i2_l = currents_lsoda()
xs, (i0_s, i1_s, i2_s) = currents_scipy()


x = np.linspace(0, 75e-3, len(i0_l))


fig, ax = plt.subplots(3, 2)

ax[0][0].set_title('lsoda')
ax[0][0].plot(x*1e3, np.abs(i0_l)*1e6)
ax[0][0].set_ylabel('i0')

ax[1][0].plot(x*1e3, np.abs(i1_l)*1e3)
ax[1][0].set_ylabel('i1')

ax[2][0].plot(x*1e3, np.abs(i2_l)*1e3)
ax[2][0].set_ylabel('i2')

ax[-1][0].set_xlabel('distance (mm)')


ax[0][1].set_title('scipy')
ax[0][1].plot(xs*1e3, np.abs(i0_s)*1e6)
ax[0][1].set_ylabel('i0')

ax[1][1].plot(xs*1e3, np.abs(i1_s)*1e3)
ax[1][1].set_ylabel('i1')

ax[2][1].plot(xs*1e3, np.abs(i2_s)*1e3)
ax[2][1].set_ylabel('i2')

ax[-1][1].set_xlabel('distance (mm)')

for i in fig.get_axes():
    i.grid('both')
fig.tight_layout()
plt.show()

Which gives:

image

The right column is the expected results.

Is there something wrong in the way I am using LSODA?
@Nicholaswogan
Copy link
Owner

Everything looks ok. The problem is that you are not comparing NumbaLSODA to Scipy's implementation of LSODA. You are comparing NumbaLSODA to the zvode integrator. Different integrators are better suited for different problems, and might also have different default error tolerances.

For this problem i compared NumbaLSODA to scipy.integrate.solve_ivp for method="LSODA". They got the same answer (see code below).

So the discrepancy isn't a problem with NumbaLSODA, but has to do with the best integrator and error tolerances for this problem. Or I suppose there could be a typo in rhs or vectorfield causing them to be different.

import numba as nb
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import ode, solve_ivp
from NumbaLSODA import lsoda_sig, lsoda

@nb.njit
def rhs_numba(x,i,p):
    i0 = i[0] + 1j*i[1]
    i1 = i[2] + 1j*i[3]
    i2 = i[4] + 1j*i[5]

    kp     = p[0]
    ks     = p[1]
    ki     = p[2]
    ep     = p[3]
    xi     = p[4]

    delta = kp-ks-ki

    di0 = 1j*kp*ep*i1*i2*            np.exp(-1j*delta*x)/4 + 1j*kp*xi*i0/8*(  abs(i0)**2 + 2*abs(i1)**2 + 2*abs(i2)**2)
    di1 = 1j*ks*ep*i0*i2.conjugate()*np.exp( 1j*delta*x)/4 + 1j*ks*xi*i1/8*(2*abs(i0)**2 +   abs(i1)**2 + 2*abs(i2)**2)
    di2 = 1j*ki*ep*i0*i1.conjugate()*np.exp( 1j*delta*x)/4 + 1j*ki*xi*i2/8*(2*abs(i0)**2 + 2*abs(i1)**2 +   abs(i2)**2)
    di = np.empty(6)
    di[0] = di0.real
    di[1] = di0.imag
    di[2] = di1.real
    di[3] = di1.imag
    di[4] = di2.real
    di[5] = di2.imag
    
    return di

@nb.cfunc(lsoda_sig)
def rhs(x, i, di, p):

    dii = rhs_numba(x,i,p)
    for i in range(len(dii)):
        di[i] = dii[i]
        
funcptr = rhs.address

def vectorfield(x: float,
                i: tuple,
                p: tuple) -> None:

    # Reconstruct the complex number
    i0 = i[0]
    i1 = i[1]
    i2 = i[2]

    kp     = p[0]
    ks     = p[1]
    ki     = p[2]
    ep     = p[3]
    xi     = p[4]

    delta = kp-ks-ki

    di0 = 1j*kp*ep*i1*i2*            np.exp(-1j*delta*x)/4 + 1j*kp*xi*i0/8*(  abs(i0)**2 + 2*abs(i1)**2 + 2*abs(i2)**2)
    di1 = 1j*ks*ep*i0*i2.conjugate()*np.exp( 1j*delta*x)/4 + 1j*ks*xi*i1/8*(2*abs(i0)**2 +   abs(i1)**2 + 2*abs(i2)**2)
    di2 = 1j*ki*ep*i0*i1.conjugate()*np.exp( 1j*delta*x)/4 + 1j*ki*xi*i2/8*(2*abs(i0)**2 + 2*abs(i1)**2 +   abs(i2)**2)

    r = np.zeros(3, dtype=np.complex128)
    r[0] = di0
    r[1] = di1
    r[2] = di2

    return r

@nb.njit
def currents_lsoda():

    i0 = 7e-3 + 0j
    i1 = 7e-5 + 0j
    i2 = 0 + 0j

    # build the initial values array
    w0 = np.array([i0.real, i0.imag,
                   i1.real, i1.imag,
                   i2.real, i2.imag], dtype=np.float64)

    # build the parameters array
    p = np.array([1507,
                  670,
                  670,
                  55,
                  19623], np.float64)

    x = np.linspace(0, 10e-3, 100)

    usol, success = lsoda(funcptr, w0, x, data=p, rtol=1e-6,atol=1e-20)

    return (usol[:,0] + 1j*usol[:,1],
            usol[:,2] + 1j*usol[:,3],
            usol[:,4] + 1j*usol[:,5])


def currents_scipy_lsoda():

    i0 = 7e-3 + 0j
    i1 = 7e-5 + 0j
    i2 = 0 + 0j

    # build the initial values array
    w0 = np.array([i0.real, i0.imag,
                   i1.real, i1.imag,
                   i2.real, i2.imag], dtype=np.float64)

    # build the parameters array
    p = np.array([1507,
                  670,
                  670,
                  55,
                  19623], np.float64)

    x = np.linspace(0, 10e-3, 100)

    sol = solve_ivp(rhs_numba, [min(x),max(x)], w0, t_eval=x, args=(p,), rtol=1e-6,atol=1e-20,method='LSODA')
    usol = sol.y.T

    return (usol[:,0] + 1j*usol[:,1],
            usol[:,2] + 1j*usol[:,3],
            usol[:,4] + 1j*usol[:,5])

def currents_scipy():

    i0 = 7e-3 + 0j
    i1 = 7e-5 + 0j
    i2 = 0 + 0j

    # build the initial values array
    w0 = np.array([i0,
                   i1,
                   i2], dtype=np.complex128)

    # build the parameters array
    p = np.array([1507,
                  670,
                  670,
                  55,
                  19623], np.float64)

    x = np.linspace(0, 75e-3, 100)
    ys = np.array([])
    xs = np.array([])

    r = ode(vectorfield).set_integrator('zvode').set_f_params(p).set_initial_value(w0, t=x[0])
    
    dx = x[1]-x[0]
    while r.successful() and r.t-dx <=x[-1]:
        r.integrate(r.t+dx)
        xs = np.append(xs, r.t)
        ys = np.append(ys, r.y)
    ys = np.reshape(ys, (int(len(ys)/len(w0)), len(w0))).T

    return xs, ys



i0_ll, i1_ll, i2_ll = currents_lsoda()
xs, (i0_s, i1_s, i2_s) = currents_scipy()
i0_l, i1_l, i2_l = currents_scipy_lsoda()

print(np.isclose(i0_ll.real,i0_l.real))

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Assignees
No one assigned
Labels
None yet
Projects
None yet
Milestone
No milestone
Development

No branches or pull requests
2 participants
@edumur @Nicholaswogan