ODE15i/s sparse Jacobian pattern trick

8 次查看（过去 30 天）
James506
James506 2016-7-16
评论： James506 2017-10-27
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

请先登录，再回答此问题。

采纳的回答

Torsten
Torsten 2016-7-18
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.

更多回答（1 个）

Damjan Lasic
Damjan Lasic 2017-10-19
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.
  1 个评论
James506
James506 2017-10-27
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);

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息

标签

产品

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by