ODE15i/s sparse Jacobian pattern trick
4 views (last 30 days)
Show older comments
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
0 Comments
Accepted Answer
Torsten
on 18 Jul 2016
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.
0 Comments
More Answers (1)
Damjan Lasic
on 19 Oct 2017
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.
See Also
Categories
Find more on Geometry and Mesh in Help Center and File Exchange
Products
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!