Jacobian, W, and inverse W calculation for ODE systems about modelingtoolkit.jl HOT 3 CLOSED

The Jacobian for the ODE system sys is the matrix J[i,j] = Derivative(sys.dvs[j])*sys.eqs[i], i.e. the derivative of the ith equation by the jth dependent variable. The Rosenbrock W is the matrix

M - γJ
# Transformed
M/γ - J

those are two different versions, and some ODE solvers use the former while others use the latter. γ is just some input constant and M is the mass matrix. The mass matrix is the matrix for Mu' = f(t,u) in vector notation, so normally an ODE just has M = I unless explicitly stated.

That's an intermediate step to the calculation of W^(-1) which is what the solvers actually want. We can also looking into just solving for the components of qrfact(W). If our Jacobian is just represented by a matrix of operations, then generic routines should directly give us this results, and then we would want to run simplification routines on it (Julia's compiler automatically does some if we @fastmath so that's an option to add as well).

For now we should make their computation explicit, i.e. add in extra fields and let the user do calculate_jacobian!(sys), etc., and then add a caclulate_performance_functions!(sys) that does them all. Then finally after we have all of that together, we can do an operation check on the Jacobian and calculate W^(-1) only when it's appropriate, i.e. when J is sufficiently small or sufficiently sparse. I am not sure what a good heuristic is there (ParameterizedFunctions.jl just does "try it and see if it fails"), so it'll take some testing which is why I am saying that should be done last.

This is the last issue for becoming feature-complete with ParameterizedFunctions.jl

from modelingtoolkit.jl.

Jacobians are working. W is working. Inversion is all that's left.

from modelingtoolkit.jl.

Inversion is working. We need better simplification routines but it does work.

from modelingtoolkit.jl.