SciML / diffeqpy

Solving differential equations in Python using DifferentialEquations.jl and the SciML Scientific Machine Learning organization
MIT License
531 stars 39 forks source link

Sparse Matrix Jacobian / Mass Matrix support #110

Closed ChrisProv closed 11 months ago

ChrisProv commented 1 year ago

I am attempting to solve a DAE system of the form: M dxdt = F(x) The mass matrix "M" and the Jacobian of "F(X)" are both sparse.

I can compute SciPy CSC matrices for M and J, but these do not appear to be compatible with the "mass_matrix" or "jac_prototype" arguments to ODEFunction. Thus, the only way I can obtain a (slow) solution is to use the dense matrices:

fun = de.ODEFunction(System_F_Julia,mass_matrix=M.toarray(),jac_prototype=J.toarray())

Is there another way I can get around this limitation? Can I convert the SciPy sparse matrices to Julia CSC matrices before using them as input?

Thanks,

ChrisRackauckas commented 1 year ago

Linking to the comments https://discourse.julialang.org/t/diffeqpy-with-sparse-mass-matrix-jacobian/91984

ChrisRackauckas commented 11 months ago

Mass matrices and sparse matrices are now supported from the Python side. Example:

from diffeqpy import de
from juliacall import Main as jl
import numpy as np

def rober(du, u, p, t):
    y1, y2, y3 = u
    k1, k2, k3 = p
    du[0] = -k1 * y1 + k3 * y2 * y3
    du[1] = k1 * y1 - k3 * y2 * y3 - k2 * y2**2
    du[2] = y1 + y2 + y3 - 1
    return

M = np.array([[1.0,0.0,0.0],[0.0,1.0,0.0],[0.0,0.0,0.0]])
f = de.ODEFunction(rober, mass_matrix = jl.convert(jl.Array,M))
prob_mm = de.ODEProblem(f, [1.0, 0.0, 0.0], (0.0, 1e5), (0.04, 3e7, 1e4))
sol = de.solve(prob_mm, de.Rodas5P(), reltol = 1e-8, abstol = 1e-8)

import scipy
spM = scipy.sparse.csr_array(M)
jl.seval("using SparseArrays")

f = de.ODEFunction(rober, mass_matrix = jl.convert(jl.SparseMatrixCSC,M))
prob_mm = de.ODEProblem(f, [1.0, 0.0, 0.0], (0.0, 1e5), (0.04, 3e7, 1e4))
sol = de.solve(prob_mm, de.Rodas5P(), reltol = 1e-8, abstol = 1e-8)

and documented.