ODE15i/s sparse Jacobian pattern trick

Hi,

I am using ODE15i/ODE15s to solve a large-stiff-sparse Differential Algebraic Equation (DAE) system (method of characteristics for a hyperbolic PDE-ODE system) in the order of thousands of equations. Now, without providing the sparsity pattern for the Jacobian ('JPattern',{dFdy, dFdyp}) this can take a very long time to solve. However, deriving the sparsity matrix for this system, particularly for dFdy is tricky. Instead, I have generated a sparsity pattern, based on the numerical Jacobian from the odenumjac function and fed this to odeset. This seems to work well, significantly reducing the solver time (often by an order of magnitude) and giving results consistent with those without providing the sparsity pattern.

Firstly, I was interested in whether this is a sound approach for large systems where the Jacobian sparsity pattern is possibly too difficult to derive analytically?

Secondly, my understanding is that implicit solvers, like ODE15i need to evaluate the Jacobian frequently. Has Matlab missed a trick, or am I missing something? If the above approach is indeed sound, then why does ODE15i/s not generate a sparsity pattern using the method above (if the user has not provided one) and use this to accelerate subsequent calculations?

Thanks

Answer 1

[Torsten]

Your approach is sound if the sparsity pattern does not change with time.

I guess this possible time-dependence is the reason why MATLAB does not generate it automatically at the beginning of the integration.

Best wishes

Torsten.

Answer 2

[Damjan Lasic]

Wanted to add some remarks even if this is an old question.

Indeed MATLAB has to recompute the entire Jacobian, because one can never be sure when and if some of the components may change. This has to do both with time dependence, but in both the sense that the functions may have some time-dependant terms as well as that the values of some functions may be for example zero at some times, which may give a value of 0 at some particular point in the Jacobian, then the value goes >0 as the function value changes.

Your trick is good and i think it can work wonders in many cases, just be careful when using - consider the above effects. I think a safe(r) approach would be to loop through various randomly generated input time and y values, so you can really be sure that the 0's in the Jacobian matrix are due to actual non-dependence rather than due to current input conditions.

[James506]

Thanks for your suggestions. I attempt to account for the time-dependence effects you describe by generating a random non-zero input vector for the odenumjac function, which seems somewhat similar to what you suggest, e.g:

rhs = feval(ode, 0, y0_rand);
Jac = odenumjac(ode,{0 y0_rand}, rhs, dfdy_opts);
Jac_sparse = sparse(Jac~=0.0);