[vox-tech] ODE solvers in C

[Peter Jay Salzman]

begin Matt Holland <mdholland@ucdavis.edu> 
> I must admit I know very little about numerical methods, but I have 
> written little Euler method solvers on occasion (albeit in more 
> forgiving languages than C).  However, my advisor tells me that I'll 
> need to use a method that's good for stiff ODEs.  This doesn't mean much 
> to me as my ODE training is woefully inadequate, but he said I should 
> use the Gear method.

ok.  in a nutshell...

stiff ODE's arise when you can't rescale the problem in a rationalized
unit system.  for instance, the following equation:

   10^999 y' + 2 = 0

looks awful from a numerical standpoint, but you can play around with
the scaling of y to make it into something pleasant:

   y' + 2 = 0

however, when you have (for example) coupled equations, you may not be
able to find a scale that makes the problem amenable for computer
solving.   another thing that can happen is that the constants can
"conspire" against you in the sense that they can be both large and
small, making it impossible to scale the problem.   this is called a
stiff equation.

stiff equations are infamous for runaway solutions (error feedback in
your algorithm) and dramatic loss of precision (this should make sense
since you're adding/subtracting large and small numbers with each
other).

i don't know much about gear's method, but i do know it's a way to deal
with stiff equations.

> I was just wondering what the C programmers out 
> there use, since I was having trouble finding stuff.  Leave it to Pete 
> to write his own solvers...
 
heh.  thanks for the compliment (i think!).    :-)

> Thanks to Shawn for the suggestions.  I'll look into linking Fortran 
> code and I'll take a look at NIST.
 
on C and fortran...

before going down this route, be aware that there are some serious
issues with integrating fortran and C object files.  serious enough that
when i was in your position a few years ago, i opted to NOT do the
fortran/C thing.  for a good discussion on the many pitfalls of
integrating C and fortran, check out "computational physics" by landau
and paez (wiley 1997) chapter 15.  you'll hate what you see.

on using canned code...

i didn't mention anything because i want to discourage using other
peoples' code (especially when you have no idea what a "good" answer
will look like ahead of time).  but if you have your heart set on it, some
standard libraries that i know are:

netlib (everything under the sun)
imsl
essl
nag (i can give you a copy.  shhh...)
slatec
lapack (linear algebra)
cern
blas (linear algebra)

some of them have C as well as fortran.  what you'll find is a bloody
mess.  something like netlib has many hundreds of programs which you'll
have to pour through to see if they do what you want.  it took me a few
days of hunting through netlib to find some routines that did what i
want.

then figuring out how to compile them, what arguments they take, what
domain of problem they're good for, what form the output is in (for
instance, some of them express an answer g(x) as G(x) = g(x)/x.  you'll
see this alot for ODE solvers that have a singularity, like the phi and
theta equations in the separated ODE for laplace's equation in polar
coordinates.   the whole thing is one big headache.  this is why i
decided to roll my own PDE solvers.  i couldn't deal with these issues.
btw, the book i cited above is VERY GOOD.  it discusses some of these
libraries and how to use some of the components.



what i would do...

1. look into the canned libraries and be disheartened (you're on your
   way to completing this).

2. look up the problem in 'numerical recipes in C'.  don't take their
   code seriously; it sucks majorly.  however, their theory is right on
   the money.  use their ideas but ignore their code.

3. pick up a couple of good books on the subject.  go to the bookstore
   (shield library/phys sci library are all but useless, but you can
   try) and get any book that has a large section on "gear's method" and
   "stiff equations".

4. do a web search.  there's quite alot out there for numerical work.

5. seek out a math professor who teaches this kind of stuff.  angela
   cheer was very competant, but kind of on the unfriendly side.  she
   knows her stuff though.  john hunter knows his stuff and is supremely
   friendly.  he was so friendly, in fact, that he was a no-brainer for
   my qualifying exam.  make use of your friends.

6. spend a week reading what you got and stewing over it in your brain.

7. don't get distracted by linux.   ;)

disclaimer...

take what i say with a grain of salt.

i had to solve a stiff 2-dimensional complex valued non-linear PDE with
singularities of a complex variable.  it took me 2 years to cross my t's
and dot my i's.  using canned solvers wasn't an option for me; there
simply wasn't anything out there.

there might be for your problem.  there might even be C code that you
have a chance of understanding the source code.

but if you find yourself getting frustrated, consider looking at
numerical recipes.  it ain't all that bad.

pete (sorry for the long email, but i'm finding myself seriously
identifying with your situation)