-
Lists are not the best datastructure for this type of code (with lots of (++), and (last)). You lose a lot of time constucting and deconstructing lists. I’d use Data.Sequence or arrays, as in C versions.
-
There is no chance for thunks of makeu0 to be garbage-collected, since you need to retain all of them (well, all of the results of “diffuse”, to be exact) all the way till the end of computation in order to be able to do “reverse” in applyBC. Which is very expensive thing, considering that you only need two items from the tail of the list for your “zeroflux”.
Here is fast hack of you code that tries to achieve better list fusion and does less list (de)constructing:
module Euler1D
( stepEuler
) where
-- impose zero flux condition
zeroflux mu (boundary:inner:xs) = boundary+mu*2*(inner-boundary)
-- one step of integration
stepEuler mu n = (applyBC . (diffused mu)) $ makeu0 n
where
diffused mu (left:x:[]) = [] -- ignore outer points
diffused mu (left:x:right:xs) = -- integrate inner points
let y = (x+mu*(left+right-2*x))
in y `seq` y : diffused mu (x:right:xs)
applyBC inner = lbc + sum inner + rbc -- boundary conditions
where
lbc = zeroflux mu ((f 0 n):inner) -- left boundary
rbc = zeroflux mu ((f n n):(take 2 $ reverse inner)) -- right boundary
-- initial condition
makeu0 n = [ f x n | x <- [0..n]]
f x n = ((^2) . sin . (pi*) . xi) x
where xi x = fromIntegral x / fromIntegral n
For 200000 points, it completes in 0.8 seconds vs 3.8 seconds for initial version