[Haskell-cafe] How do I get a long iteration to run in constant space

W M caeshmer at yahoo.com
Mon Oct 4 11:18:17 EDT 2004


Hi, I'm starting to learn haskell, and I was
considering it for a simulation I've been designing,
but I ran into this problem with a simple numerical
calculation.  The attached program is supposed to
solve a stiff differential equation, and you have to
use a very small stepsize and lots of steps to resolve
the part of the solution near a singularity.  It looks
to me like this program should run in constant space
because the iterate2 function is tail-recursive, but
when I run it on hugs or ghc 2.6.1, I get an error
that it runs out of stack space on the third
calculation.  (And the quick fix is to give it a flag
that increases the stack, but if I'm to use haskell
for a potentially large and long simulation, I can't
afford to let it run for three days only to get an
error and figure out that it needs more stack space
than the computer can give.)  I suppose in my ODE
example it's building up expressions somewhere for
lazy evaluation, so I've tried putting seq and $! in
various places.  I have another example where using
"else seq x (iterate2 ...)" helps, but that doesn't
solve my problem here.  I've also tried -O but it
doesn't help.  Can anyone give me some idea as to what
to do?

Thanks,
-wgm

iterate2 f x n = 
    if n == 0 then
      x 
    else
      iterate2 f (f x) (n - 1)

integrate step f xInit xFinal yInit nSteps =
    let
      h = (xFinal - xInit) / fromIntegral nSteps
      next = step f h
    in
      iterate2 next (xInit, yInit) nSteps

rk4Next f h (x, y) =
    let
      a1 = h * f x y
      a2 = h * f (x + h/2) (y + a1/2)
      a3 = h * f (x + h/2) (y + a2/2)
      a4 = h * f (x + h) (y + a3)
    in
      (x + h, y + a1/6 + a2/3 + a3/3 + a4/6)

rk4 = integrate rk4Next

stiff x y = -( ((exp (x - 1)) + (x - 1) * y)
	       / (1 - x - 0.0001 * y) )

main =
    do
      putStr "Begin\n"
      putStr (show (rk4 stiff 0 1 (exp (-1)) 1000))
      putStr "\n"
      putStr (show (rk4 stiff 0 1 (exp (-1)) 10000))
      putStr "\n"
      putStr (show (rk4 stiff 0 1 (exp (-1)) 100000))
      putStr "\n"



		
__________________________________
Do you Yahoo!?
Yahoo! Mail - You care about security. So do we.
http://promotions.yahoo.com/new_mail


More information about the Haskell-Cafe mailing list