[Help-gsl] How to avoid exponential increase in ODEs solutions ?

[Ruben Farinelli]

Hi,
I hope this kind of problem has already been faced by users of this
mail-list.

I have a code that solves a relatively simple equation, it is somewhat an
Helmotz equation (with negative parameter) whose solution is combination of
an exponential increasing+ exponential decreasing function
(one can treat first this very simply analytically).

This is the equation

y''[r]+2/r y'[r]-k y[r]=0

analytical solutions are

 y[r]=C1 Exp[- sqrt(k) r] /r+ C2 Exp[ sqrt(k) r] /(Sqrt[k]* r)

If one is interested in only the exponential decreasing solution
can just simply set C2=0.
For a numerical solution, one must instead put initial conditions on the
function and its derivative obtained from the analytical treatment
(which is equivalent to dropping the exponential-increasing solution), so
that the numerical code should of course integrate and obtain only the
exponential decreasing solution (C2=0).

After doing that, depending on the method and accuracy, you get the
desired exponential decay but after a while, (for higher initial K-values)
at some point the exponential increasing solution starts to dominate.
It is of course only matter of numerical integration and not of wrong
initial conditions.

What I'm wondering is whether there is some way to force the code anyway to
avoid the fake exponential increasing, or it must be held by the user
(e.g., stopping when derivative becomes positive).

Presently I use

gsl_odeiv_control *k = gsl_odeiv_control_y_new(eps,eps);

gsl_odeiv_control *g = gsl_odeiv_control_yp_new(eps,eps);

For interested users, I also attach the code.

The command line after standard compilation  is very simple, e.g.,

./test_hmz yi=1 kpar=0.2 meth=bsimp xmax=100 eps=1e-14

(for meth you can set bsimp, rk4, gear1, gear2, rk8pd, rkck),

yi is th initial value of the function (doesn't modify the general solution
behavior)

kpar is the parameter of the function (the higher kpar, the steepest the
exponential decay)

xmax is the right bound of the integration domain (which starts from x=1)

eps is the precision of the GSL controlling function.

You'll see the output written in the two-column file "results_hmz.dat",
whose plot immediately shows the behaviour.

I hope there are some tricks I don't know to get the right behavior.


Thank you in advance

Ruben

test_hmz.c
Description: Text Data