[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