# Stiffness emerges as number of ODEs increases

I want to solve a system of ordinary differential equations with Matlab. I need this to solve a mechanical engineering related problem. If $n$ is the number of degrees of freedom of my mechanical system, then the number of variables (and also of equations, obviously) is $n^2$. I can easily numerically integrate the ODE for example for $n=6$. I can tell that the problem is stiff because ode45 takes a very long time, while ode15s takes seconds and gives accurate results. For a 13 degrees of freedom system (i.e. 169 variables) instead the computations are extremely slow and that's because the step size is tiny, at the beginning is $\approx 10^{-18}$, after one minute it becomes $\approx 10^{-15}$ and it gets bigger and bigger but at some point I have to stop the computation and the integration hasn't reached $10^{-10}$ (whereas with $n=6$ I easily integrate from $0$ to $100$). I was wondering if this is something I could expect because of the 169 variables. What can I do in this cases? I tried all the available solvers for stiff problems but it's not any better.

EDIT: I understand that the 6 degrees of system and the 13 degrees of freedom system are just different models so I can't make a comparison, but from the mechanical point of view the bigger one is not much more complex. I'm very naive for what concerns numerical solution of ODEs and I expected that doubling the size of the problem (thus solving 4x equations) would have led to (roughly) 4x computation time. Of course it cannot be as simple as this, but I'm surprised by the difference between the two models. The computation time was "seconds" for the small system and it has become "neverending" for the big one!

EDIT: I eventually discovered what my problem was! As I said, the ODE was related to a mechanical engineering problem. The issue was that the variables represented actual physical quantities the were in different units and very different orders of magnitude! I just scaled the problem by means of appropriate scalar factors and the computation time decreased significantly! Obviously, you have to scale the solution back to the right units, but it's not a big deal. For additional information I recommend this paper http://dx.doi.org/10.1137/S0895479803434914

• You will substantially increase the likelihood of getting a useful reply if you include the ODEs in question. There must be some compact way of representing it which does not require you to write out 169 equations in 169 variables. In particular, consider adding a few words about why you need $n^2$ equations/unknowns for a system with $n$ degrees of freedom. Are you differentiating the flow of one ODE with respect to the initial condition? – Carl Christian Jun 3 '16 at 19:22
• Did you try increasing the values of the accuracy tolerances for ode15s, AbsTol and RelTol? Perhaps ode15s is trying to obtain a more accurate solution than you need. – Bill Greene Jun 3 '16 at 20:31
• @CarlChristian I haven't mentioned the ODEs because I can't ask the other people to read a paper to understand my problem, but since you're asking the ODE that I want to solve is for example eq. (23) in dx.doi.org/10.1137/0727062 I need $n^2$ equations because they are the entries of an orthogonal matrix. – Roberto Belotti Jun 6 '16 at 9:56
• @BillGreene I tried to do what you suggested but with little improvement. As you can see from my other comment, I cannot lose much accuracy because the variables of my ODEs form a matrix which I need to be orthogonal. Lowering the values of AbsTol and RelTol actually speeds up the computation, but not dramatically, though. – Roberto Belotti Jun 6 '16 at 10:04
• The computed approximation need not be orthogonal unless you force it. You are trying to compute a trajectory which lives in a manifold, i.e. the set of orthogonal matrices. Truncation and rounding errors will cause you to slide away off this manifold unless you compensate. This is why think you need to treat your problem as a DAE. Time pr. Iteration is often proportional to the number of unknowns, but the iteration count can easily grow as the problem becomes more ill-conditioned/harder. This may be the issue here. This is a guess on my part, and you should get a second opinion. – Carl Christian Jun 7 '16 at 11:53