SUNDIALS-CLN：An Arbitrary Precision Extension of the CVODE Library

原始代码作者网页

CVODE是美国livermore劳伦斯国家实验室应用计算中心推出的早期DVODE（fortran版本）常微分方程求解器的C语言改写版。其应用极为广泛。CVODE现在也是该实验室的一些新成员开发维护的SUNDIALS (实际是一系列英语单词首字母缩写，组合在一起的新单词直译：日晷）软件包的一个组件。

Mathematica的NDSolve求微分方程数值解的部分功能直接引用了Sundials的代码（相当于为这个包提供了一个Mathematica下直接和更方便调用的界面）。——但是我在使用过程中发现，NDSolve中调用CVODE部分的代码，无法像Mathematica中其它计算一样设置为支持复数或任意精度的计算。这说明，Wolfram的确是老老实实为sundials提供了一个mathematica界面，并没有改写底层代码为mathematica风格。

偶然从Sundials用户群的邮件来往中发现了一个Sundials 2.2版本的任意精度支持的改写版，转发过来，希望Wolfram开发者也看到了类似的改写工作，将来依法炮制。——此外，我个人更喜欢的是Eigen C++ template library + mpir+mpfr+gmp+ mpfr C++ wrapper形式的改进，从个人开发者的角度看似乎更合理。

CVODE is a widely used C library for solving ordinary differential equations. It is part of the SUNDIALS suite of libraries. CLN is a C++ library that implements arbitrary precision numbers. I have converted the CVODE library to use the CLN floating point class as its representation of real numbers. This means CVODE can solve ordinary differential equations with as much numerical precision as desired.

The development of this software was partially funded by the National Science Foundation (Grant DMS-0514468).

Download

Copy one of these files to your home directory: sundials_cln.tar.bz2 (8.9M), sundials_cln.tar.gz (9.9M).

More information about the software, including how to install it, is included in the file README-WW.

This software is still in the experimental stage. If you have any questions about it, please contact me. Also let me know if you download and install it. Even if you have no problems, I am interested in knowing how this software is used.

Example

We solve the following system of differential equations:

x’ = (-y + x - x3/3)/ε

y’ = x-a

The following plot shows two solutions. Both were computed with 100 digits of precision. The parameter values were a=-1.5 and ε=0.005. The initial x coordinate for both solutions is x(0)=-2. For the dashed curve, the initial y coordinate is

y(0)=-0.699903492086728166143978401208788824227818610553271339023675758554576281195181698

and for the solid curve it is

y(0)=-0.699903492086728166143978401208788824227818610553271339023675758554576281195181809

The difference in the two initial y values is less than 2e-69. An initial value problem solver that used standard double precision floating point numbers could not compute these solutions.