Hi,

I'm after some general advice on selecting an appropriate ODE solver for a single nonlinear ODE, say dy/dx=f(y,x,B). the ODE solution is inside a Monte Carlo loop and is solved a zillion times with different parameters B in the RHS (that is fairly simple), so efficiency is quite important. Since the ODE itself is approximate, I would be happy with accuracy of the order of <5%.

Runge-Kutta schemes was my first guess, until I considered the possibility of stiffness due to RHS terms such as (-B2*y^B3), so an implicit or semi-implicit solver seemed more appropriate. However, since {B} is a parameter, the stiffness may vary depending on the current trial B (and may not even occur). But since the ODE solver is inside the loop, there would seem to be a need to switch the solver at runtime. Code-wise this is no problem if I knew when to switch solvers. Is there some cheap way to determine when would an explicit Runge-Kutta start to fall over itself due to stiffness?

I had 2 ideas:

1. Look at the stepsizes - if they are too small change to implicit. This requires setting some threshold time step size.

2. Run with an implicit code at all times: this has the obvious robustness advantage but likely far slower for nonstiff cases, plus need for Jacobians (which are analytically available).

I would appreciate if somebody could share some insights into such problem.

thank you in advance,
Dmitri

Hi Dmitri,

As NR explains at the beginning of Section 16.6, stiffness is a property of systems of equations. It should not be a problem for a single equation. For 5% accuracy, I would imagine 2nd order RK would be fine.

Saul,

That was my first reaction, since I did see the statement in NR and my experience with finite elements suggests stiffness increases with number of equations. However, how does that square with the model ODE dy/dx=ky which imposes stability restrictions on the forward Euler but obvioulsy not on the implicit Euler? I realise the truncation error of the two is almost identical, but the explicit scheme does seem less robust in general, especialy for exponential decay. Presumably explicit RK schemes have better stability but still at some point would suffer stepsize limitations. So my feeling was that stiffness is far more common in ODE systems, but still may appear in single ODEs. Am I missing something trivial here?

your insight would be appreciated,
Dmitri

Hello everbody,

I want to ask a question about the folowing problem:
the equation is : d(d(y)/dt)/dt=G(t,d(y)/dt,y)
this equation reduces to
the system
d(y)/dt=z
d(z)/dt=G(t,z,y)
and the boundary conditions are:
t=0, y=0
t=1, y=1
My question is if this problem is solvable with realaxation method or not? Because I tried to do but I get the message 'singular matrix, row all 0 in pinvs'!!
Thank you