bugGNU Scientific Library - Bugs: bug #30540, please, add convergence checks in...


You are not allowed to post comments on this tracker with your current authentication level.

bug #30540: please, add convergence checks in rk4imp/rk2imp

Submitter:  Sergey B Kirpichev <skirpichev>
Submitted:  Sat 24 Jul 2010 12:15:06 PM UTC
Category:  Accuracy problem Severity:  3 - Normal
Operating System:  Status:  Ready For Test
Assigned to:  None Open/Closed:  Open

Thu 26 Aug 2010 12:05:28 PM UTC, comment #2: 

I can't figure out how to use this code for fixed step size.  imprk2/imprk4 seems to be not working at all for this case:
$ ./bug
Segmentation fault

Other steppers works (e.g. rk2) just fine.  See bug-ode2.c.

Second, Newton-type iterations could be inefficient in general.  See: Ernst Hairer, Christian Lubich, Gerhard Wanner, Geometric Numerical Integration.  Structure-Preserving Algorithms for Ordinary Differential Equations. Springer Series in Comput. Mathematics, Vol. 31, Springer-Verlag 2002. pp. 330-332.

(file #21314)

Sergey B Kirpichev <skirpichev>
Fri 06 Aug 2010 01:09:11 PM UTC, comment #1: 

Thanks for the bug report.

I'm merging a completely new and upgraded version of the ode solvers "ode-initval2" written by Tuomo Keskitalo in the next release.  It includes the new solver ode-initval2/imprk4.c for implicit RK4.

You can get the branch with

bzr checkout http://bzr.savannah.gnu.org/r/gsl/ode-initval2/

To compile it use

  autoreconf -i -f; ./configure --enable-maintainer-mode ; make

Please try it and let us know if it solves this problem.

-Deleted Account- <bjg>
Sat 24 Jul 2010 12:15:06 PM UTC, original submission:  

Quote from http://lists.gnu.org/archive/html/help-gsl/2009-12/msg00047.html

It seems, for example, the 2 stages Gauss points are only iterate 3 steps, without checking whether it is converged or not.

The reason I care about this is, Implicit Gauss RK4 is symplectic, which for Hamiltonian Dynamics simulation is a very important property, it has better conservation of the first integral.

Attached example code (bug.c) for Henon-Heiles problem shows the linear error grow in the Hamiltonian.  Compare with correct behavior for patched solvers (rk2imp.png/rk4imp.png).

Please, review attached patch, targeted to solve the problem.

Sergey B Kirpichev <skirpichev>


(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attached Files
file #21314:  bug-ode2.c added by skirpichev (2KiB - application/octet-stream)
file #21050:  bug.c added by skirpichev (2KiB - application/octet-stream)
file #21051:  rk2imp.png added by skirpichev (16KiB - image/png)
file #21052:  rk4imp.png added by skirpichev (13KiB - image/png)


Depends on the following items: None found

Items that depend on this one: None found


Carbon-Copy List
  • -email is unavailable- added by skirpichev (Submitted the item)

    Follow 6 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2010-08-26 skirpichev Attached File- Added bug-ode2.c, #21314
    2010-08-06 bjg StatusNone Ready For Test
    2010-07-24 skirpichev Attached File- Added bug.c, #21050
        Attached File- Added rk2imp.png, #21051
        Attached File- Added rk4imp.png, #21052
        Attached File- Added 0001-Add-convergence-checks-for-rk-imp.patch, #21053

    Back to the top

    Powered by Savane 3.13-8484.
    Corresponding source code