Re: [boost] Numerical integration library?

From: Frederick Akalin (akalin_at_[hidden])
Date: 2006-07-31 23:38:14

"Numerical integration" is often loosely used to mean quadrature, ODE  
solving, and PDE solving.  I was just following the terminology on  
the Wiki; I'm sorry if I wasn't clear.

I was talking about ODE solving in my original e-mail, but of course  
there might be a need for PDE solving and quadrature libraries, too.   
That PDE library looks interesting, although I don't have much  
experience with solving PDEs.  It looks like a bunch of people have  
quadrature routines, too, so it might be beneficial to integrate that  
into a coherent library.

Anyway, I hacked up an example program demonstrating what I'm talking  
about, available at: http://akalin.cx/boost/ode_solver_example.cpp ,  
which should be compilable by itself.  The interesting stuff is the  
classes at the top.  Some notes:

- The step routines are taken straight out of a numerical analysis  
reference.  Note that they are completely container-agnostic, so they  
can be used with floats, vectors of floats, vectors of particles  
which have their own floats, etc.  Other libraries are more  
restricted.  For example, the GSL library routines only work with  
double arrays: http://www.gnu.org/software/gsl/manual/html_node/ 
Ordinary-Differential-Equations.html#Ordinary-Differential- 
Equations .  Also, the ODE++ library ( http://www.minet.uni-jena.de/ 
www/fakultaet/iam/ode++/ode++_e.html#Problem ) requires adapting your  
code to their framework, including using their vector/matrix classes.  
(I only took a cursory look at the documentation -- correct me if I'm  
misrepresenting anything)
- I haven't done any performance tests, but it seems to me that a  
reasonably smart container/adapter (with expression templates) would  
incur only a small performance hit compared with hand-coded C or  
Fortran. (Tell me if I'm dreaming :) )
- The member functions are templated instead of the classes  
themselves.  That's just for convenience, as more sophisticated  
steppers (ones that store state) probably need to be class- 
templated.  I also omitted a Jacobian parameter for simplicity, as  
these solvers don't need them.
- A real ODE-solver library would of course include higher-level  
drivers for different types of adaptive stepping.  Also, it would  
probably include adaptors for common containers (Blitz++ vectors, etc.)
- I didn't test this example program thoroughly, so let me know if  
you spot any errors.

Has anything like this been done before?  I have a hard time  
believing that such a simple idea would go unnoticed, even though my  
own searches turned up nothing.

Fred

On Jul 31, 2006, at 8:46 AM, Sohail Somani wrote:

>
> I'm unsure what ODE's have to do with integrals.
>
> Anyway, for solving PDEs one should look at
> http://www10.informatik.uni-erlangen.de/~pflaum/expde/public_html/ 
> index.
> html for an example.
>
> I'd be interested to know if anyone has used this library :)
>
> Sohail
> _______________________________________________
> Unsubscribe & other changes: http://lists.boost.org/mailman/ 
> listinfo.cgi/boost