# ODE forth-order very stiff equation with large errors

I’m using Mathematica home edition software to numerically solve a specific inflation equation in cosmology. The ODE equation is forth- order, non-linear, stiff. I was using the stiffness switching method and with Working Accuracy of 60, 100, 200, but I’m not getting a stable result. Plugging in the result back into the equations, I’m getting a large deviation, especially around the area that I wish to extract the info. It is a computation issue, and I wonder if someone can recommend different method, or a different software to try. I'm sure there are plenty of researchers that use to deal with similar cosmology inflation equations all the time. thx Ezra

• Can you provide the ODE and some numerical results? – Anton Menshov Jun 4 '20 at 22:22
• Can i send you the Mathematica notebook file? it has the equation and the numerical results – user258923 Jun 4 '20 at 23:10
• no, with restricting to Mathematica I would suggest posting at Mathematica SE. And sending things directly is not how SE works. – Anton Menshov Jun 4 '20 at 23:12
• understood., i will provide the information in a day or so, thank u – user258923 Jun 4 '20 at 23:15

## 1 Answer

This is not an answer to your question, but more of an observation: More often than not, an ODE is "stiff" or a linear system is "nearly singular" because of a mistake either in deriving the equation to be solved, or in implementing it. Trying to find a way to solve what you have is then just a way to paper over the problem. If you were to find a way to stably solve for what you're interested in, you'd get a number that is, nevertheless, meaningless.

As a consequence, I would recommend to go back to the derivation of where your model comes from and to really understand why it is stiff, and if that is something that you expect for physical reasons. Only when you are convinced that the equation should indeed be stiff does it make sense to think about what method you should use to solve it.

Empirically, in 90% of the cases I've seen, the problem is then in fact no longer stiff/illconditioned (as this was an artifact/bug) and in the other 10% thinking long enough about the origin of the problem suggests how one should approach solving the problem.

• Thank u. I did study the equation over a month, checking the poles, got the asymptotic analytic solution and i think that i know what to expect as a solution. – user258923 Jun 4 '20 at 23:13
• Well done, then! :-) – Wolfgang Bangerth Jun 4 '20 at 23:20
• 99% of the time where a user is like "this equation is hard!" what really happened was the model was implemented incorrectly, and the integrator was correctly exploding to infinity on the wrong model. So it's definitely the first thing to try. The best way to diagnose this is just to try 3 different solvers, preferably different languages and different authors. If all fail, do you think all methods ever are the problem, or the model? Could be the former, but it's usually the latter. If the model is correct and the issue is due to a singularity, you might need a different integrator. – Chris Rackauckas Jun 5 '20 at 2:20
• Most of the integration techniques that are used are not robust to having a singularity in the middle of a solution. – Chris Rackauckas Jun 5 '20 at 2:20
• you are correct, and the singularity is the middle of the equation, which solvers do you recommend for that ? – user258923 Jun 5 '20 at 14:31