class to hold the problem data

[ghost]

Hi, thank you for your work on this project.

Did sundae support to use class to hold the problem data, like scikits.odes did, such as the example showed here?
https://github.com/bmcage/odes/blob/master/ipython_examples/Double%20Pendulum%20as%20DAE%20with%20roots.ipynb,
actually what I want to do is, pass a python class (which include some data and defined function ) to the residual and jac function to calculate the res and jac results.
thanks!

[c-randall]

Hi,

Yes, scikit-sundae supports classes. There are a variety of ways users can define and input their problem. The official API is built around users providing functions, but it also accepts methods for callables. For example, if you just wanted to define your userdata in a class and still use typical functions to define the problem, you could use:

import numpy as np
import sksundae as sun
import matplotlib.pyplot as plt


def resfn(t, y, yp, res, userdata):
    
    p0 = userdata.p0
    p1 = userdata.p1
    p2 = userdata.p2 
    
    res[0] = yp[0] + p0*y[0] - p1*y[1]*y[2]
    res[1] = yp[1] - p0*y[0] + p1*y[1]*y[2] + p2*y[1]**2
    res[2] = y[0] + y[1] + y[2] - 1
    
    
class UserData:
    p0 = 0.04
    p1 = 1e4
    p2 = 3e7
    

userdata = UserData()

tspan = np.logspace(-6, 6, 50)
y0 = np.array([1, 0, 0])
yp0 = np.array([-0.04, 0.04, 0])

solver = sun.ida.IDA(resfn, atol=1e-8, algebraic_idx=[2], userdata=userdata)
soln = solver.solve(tspan, y0, yp0)
print(soln)

soln.y[:,1] *= 1e4  # scale the y1 values for plotting

plt.semilogx(soln.t, soln.y)
plt.legend(['y0', 'y1', 'y2'])
plt.xlabel(r"$t$ [s]")
plt.ylabel("concentration")

If you wanted to define your entire problem as a class, where the resfn and jacfn are given as methods, you could instead do:

import numpy as np
import sksundae as sun
import matplotlib.pyplot as plt


class DAEProblem:
    p0 = 0.04
    p1 = 1e4
    p2 = 3e7

    def resfn(self, t, y, yp, res):
        
        p0 = self.p0
        p1 = self.p1
        p2 = self.p2 
        
        res[0] = yp[0] + p0*y[0] - p1*y[1]*y[2]
        res[1] = yp[1] - p0*y[0] + p1*y[1]*y[2] + p2*y[1]**2
        res[2] = y[0] + y[1] + y[2] - 1
        
    def jacfn(self, t, y, yp, res, cj, JJ):
        pass
        
    def solve(self, tspan, y0, yp0, **options):
        
        # could also add jacfn to **options if you want
        # options = {**options, 'jacfn': self.jacfn}
        
        solver = sun.ida.IDA(self.resfn, **options)
        soln = solver.solve(tspan, y0, yp0)
        
        return soln
    

tspan = np.logspace(-6, 6, 50)
y0 = np.array([1, 0, 0])
yp0 = np.array([-0.04, 0.04, 0])
    
dae = DAEProblem()
soln = dae.solve(tspan, y0, yp0, atol=1e-8, algebraic_idx=[2])
print(soln)

soln.y[:,1] *= 1e4  # scale the y1 values for plotting

plt.semilogx(soln.t, soln.y)
plt.legend(['y0', 'y1', 'y2'])
plt.xlabel(r"$t$ [s]")
plt.ylabel("concentration")

I'm sure there are also plenty of other organizational and design decisions that are compatible with the package. However, while the API is similar to scikits.odes, we do not provide base classes like ida.IDA_RhsFunction or ida.IDA_JacRhsFunction. Hope this helps.

[ghost]

thanks for your feedback, by the way, is this possble to support python 3.8, I try to install it in python 3.8 envirement, and failed, then I noticed this need >=3.9, thanks!

[c-randall]

No, 3.8 is not supported. Each release only covers the versions that are actively maintained at the time of the release. The first release was late last year and covered 3.9-3.13. The next release will drop 3.9 and add 3.14.

You can see the schedule for which versions have reached end-of-life and future scheduled plans to stop support for older versions on the official Python website: https://devguide.python.org/versions/.