5
$\begingroup$

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$!),

Difference between calculated and expected at t=0

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.

enter image description here

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

Improve this question
$\endgroup$

1 Answer 1

Reset to default
1
$\begingroup$

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).

Improve this answer
$\endgroup$
4
  • $\begingroup$ That's a good suggestion, and indeed I was leaning in that direction, with the model of a damped pendulum in mind. The problem I have is this seems to only work for short (albeit longer than when it was absent) times, and can interfere with the dynamics in a major way- that is, I'm taking out more than just that un-physical, fast growing behavior. I am wondering if there is a more clever way to go about this. $\endgroup$
    – Nick P
    Commented Apr 29, 2014 at 4:32
  • 1
    $\begingroup$ I am not an expert on these type of problems, but it sounds like whatever time integrator that you are using is unconditionally unstable. Is it possible to change the time integrator? If not, have you looked at the energy spectrum? Is energy building up in the small scales? Maybe some kind of filtering would help (Gaussian, cut-off, box, etc)? $\endgroup$
    – James
    Commented Apr 29, 2014 at 13:44
  • $\begingroup$ Hi @user2687246, I'm integrating this in time using an implicit-differential algebraic solver, namely, the SUNDIALS IDA suite, which ostensibly has unconditionally stable routines. Energy is in fact building up in the smaller scales, but filtering it out of the solutions explicitly is not ideal, since this would lead to discontinuities between steps, which does not go well with the Jacobian based method the IDA solver employs. This is why I was trying to add this filtering implicitly via a numerical diffusion term. $\endgroup$
    – Nick P
    Commented Apr 29, 2014 at 22:21
  • $\begingroup$ Well if energy is building up at the small scales indefinitely then its not surprising that it becomes unstable. The thing is, by adding a viscosity term you are essentially filtering those high wave numbers out, so adding viscosity is the same as filtering in some ways. The larger your viscosity, the more wavenumbers that are going to be damped. I would suggest try and implement the above viscosity term and look at the spectrum as you increase c. Does the high wavenumber energy get removed more and more as c is increased? Maybe try squaring the above viscosity term (like a hyperviscosity). $\endgroup$
    – James
    Commented Apr 30, 2014 at 13:29

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.