On 8/27/10, Tuomo Keskitalo <[email protected]> wrote: > A quick comment: You need to specify the function for calculating the > Jacobian for sys-henon, the new implicit codes require it. Please look > at documentation in doc/ode-initval.texi. Use of fixed step size is > similar to old ode-initval.
Right. Fixed version (with jacobian, attached below) seems to be working. Actually, it works only for eps_abs ~ round off error. Already for eps_abs = 1e-10 it's easy to see small (linear?) drift in hamiltonian from the constant. Why mix step size control tolerances and root-finding tolerances? Why not just use prefered fixed value of eps_abs for root-finding (e.g.: GSL_DBL_EPSILON), as does attached patch for bug #30540 or other implicit solvers (e.g. http://www.unige.ch/~hairer/prog/stiff/radau5.f), but instead complicate low level solver routines?
bug-ode2.c
Description: Binary data
_______________________________________________ Bug-gsl mailing list [email protected] http://lists.gnu.org/mailman/listinfo/bug-gsl
