Suspected bad Jacobian for a stiff system of ODEs with ode15s

Dear All,

I have a quick question about ODE15s and solving a system of stiff ODEs that are particularly problematic.

First of all, I have a system of DAEs where algebraic variables can be eliminated. So inside the right hand side function, I compute the algebraic variables as a function of state variables and use the results to evaluate the right hand side outputs, i.e., the time derivatives. This strategy has been working well, especially with the case where we are imposing a constraint on the system of DAEs and computing the solution for the algebraic variables.

However, there is a case where we want to use an empirical formula to compute the algebraic variables. The empirical relationship again requires us to know the states at time t, which is always given to us as input to the right hand side function at time t, and some known parameters. The problem is that, without reducing the magnitude of a parameter by , the numerical integration takes an exorbitant amount of time; of course by tuning down the parameter, we are not simulating the original system anyways.

Here are some statistics from ode15s:

• Numerical Integration Summary (no tuning down of the parameter)
• *Solver : ode15s
• *Number of Successful Steps : 3271
• *Number of Failed Steps : 3119
• *Number of Function Evaluations : 2067372
• *Number of Jacobian Evaluations : 1671
• *Number of LU Decompositions : 4081
• *Number of Linear Solves : 11593
• Numerical Integration Summary (with a tuned down parameter)
• *Solver : ode15s
• *Number of Successful Steps : 152
• *Number of Failed Steps : 14
• *Number of Function Evaluations : 11857
• *Number of Jacobian Evaluations : 9
• *Number of LU Decompositions : 45
• *Number of Linear Solves : 338

One key note here is that, when we do not tune down the paremter, we still get the right physical behavior form the numerical integration of the system of ODEs. However, as evidenced by the numer of function evaulations, the numerical integration takes forever... Furthermore, it is noteworthy to see that, without tuning down the parameter, not only we are seeing many more number of failed steps, but also the proportion of the failed steps is significantly larger. And we also see that due to the increased number of the Jacobian approximation via finite difference, we are seeing an increase in the number of function evaluation as well.

I am stuck as to what might be causing this behavior, which is caused by simply turning down the magnitude of the parameter used in computing the algebraic variables, inside the right hand side function, as a function of state variables. I think the relative and absolute tolerance could be playing a role here somehow. Also, for the expression for the algebraic variable, it is a simple linear function of the state variables where we have a constant of proportionality and a difference between some state variables: i.e., , where x's are the algebraic variables and y's are the state variables. The problem is when we are computing the difference term, the state variables can be very similar in its magnitude, thereby causing some catastrophic cancellation problem... And this may cause problem with the given tolerance values in ode15s that we specify... Tightening the tolerance did not help in this case...

Do you have any suggestions as to what we can do to see why ode15s is having a much harder time without tuning down the parameter?

Thank you very much for your thoughts and support!

Sincerely,

Tae

p.s. Also, I realized that specifying a pattern for the Jacobian matrix can save a significant amount of computation effort. However, without doing the actual analysis, we are not certain which entries in the Jacobian matrix is indeed zero. We could use odenumjac.m to get the Jacobian matrix pattern at the beginning of the numerical integration. However, the pattern can change over time and I've been running to issues in the ode solver because of the wrongly specified Jacobian pattern causing the integrator to stall...

[Sam Chak]

Hi Kim

Can you show the mathematical model of the stiff system in question?

[Taehun Kim]

Dear Sam,

Thanks a lot for your interest in helping. It is quite complicated to provide the entire mathematical model. But in a nutshell, it is a discretized system of convective-diffusion-reaction equation (Convection–diffusion equation - Wikipedia). I am solving the coupled system of PDEs using the method of lines and I am using the first order finitve volume discretization scheme. The algebraic variables are related to the convective flow terms and the state variables corresponds to the dependent variables in the PDE model.

Sincerely,

Tae

[Sam Chak]

No worries, Tae. I understand the some systems are quite complicated to show. You probably need advices from users who are very good at math, physics, and MATLAB. I know a few of them: @Torsten, @Bruno Luong, @Matt J, @David Goodmanson, @John D'Errico.

[Taehun Kim]

Thank you so much Sam!

[Torsten]

@Taehun Kim

I don't think that a reason for this behaviour can reliably be diagnozed without your code.

Even with your code, it might be impossible to find the real reason.

[Taehun Kim]

I see... This sounds like doing a research! Thanks for your supports always!

[Torsten]

Or try another integrator, e.g. RADAU available under

https://www.mathworks.com/matlabcentral/mlc-downloads/downloads/submissions/56162/versions/1/previews/RADAU/radau.m/index.html

[Taehun Kim]

This turned out to be a modeling problem, rather than a numerical problem. By specifying the right set of parameters, I was able to resolve the numerical issues... Thanks for your support in this!