Scheme to alleviate (numerical?) instability in system of coupled nonlinear ODEs

I'm solving a system of nonlinear ODEs that take the form

$Q_{nm} \ddot{y}_m + S_{nkl}\dot{y}_k\dot{y}_l +V_n = 0$

where Einstein summation is assumed, $y_i$ are the dependent (complex) variables, $Q_{nm}$,$S_{nkl}$ and $V_n$ are all low order polynomials in the dependent variables, $n =1,...,N$, with $N$ the $resolution$ of the model. These equations, together with the initial conditions $y_i(0)$ and $\dot{y}_i(0)$ completely specify the problem.

The scheme is solved using (the Fortran interface to) SUNDIALS. In particular, because of the nontrivial coefficient $Q$, an implicit differential algrebraic solver is used (IDA) to solve the equations.

I have successfully implemented these equations and test cases yield coherent results. The computation time is constrained by $N$, the resolution of our model.

Now, no matter how precisely I specify my initial conditions for certain physically relevant scenarios, I find a spurious growth in my system, that leads to the equations becoming stiff, and eventually Q is no longer invertible and the integration fails. The problem is (temporarily) alleviated by increasing $N$, but no matter how high I can reasonably make $N$, the instability still appears. So I would like to see if anyone has any recommendations on a smoothing operation, or perhaps adding some type of viscosity, to attenuate this instability, while keeping the results physically relevant.

As an example of this, I put IC that correspond to equations with known solution. The physical space solutions are related to a summation of the fourier series, with coefficients $y_i(t)$. At time t=0, I see that the difference between the numerical solution to my system, and what I expect is good (see fig 1; note this is for $N=16$!),

and is order $10^{-8}$, which is the order of my tolerances in the integration. Note, the solutions are order $10^{-1}$. Now, at time t=0.1, the difference between what I expect and the calculated solution is no longer so negligible, namely we see that the difference is 5 order of magnitude larger! This increases at this exponential rate until the system eventually becomes stiff and stops integrating. I am not sure if this is some kind of instability due to the physics, or the numerics, but in either case I want to get ride of it.

Naively, it seems like some kind of smoothing between integration steps could be useful, but this is has major drawbacks. The IDA uses an internal Jacobian, so it would not be a good idea to impulsively change the solutions between integration steps.

Probably the most attractive technique would be to add some kind of $viscosity$ to the system, to dampen or get rid of these growing instabilities. The part of the solution that is unstable is growing fastest, so adding something to these equations that can pick this out would be ideal. Exactly how to do this is unclear to me, so any suggestions or references to similar problems would be greatly appreciated!

Thanks,

Nick

If you are trying to add some viscosity where the solution grows fastest then maybe something like $\nu = c\left|\frac{dy}{dx}\right|$ would work where c is some constant that you choose (probably by trial and error).