--
-- Explicit Euler method for 1D linear diffusion equation
-- (using lists)
--
module Euler1D ( stepEuler, evalEuler, makeu0 ) where
-- impose zero flux condition
zeroflux mu (boundary:inner:xs) = boundary+mu*2*(inner-boundary)
-- one step of integration
stepEuler mu u0 = (applyBC . (diffused mu)) $ u0
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] ++ inner ++ [rbc] -- boundary conditions
where
lbc = zeroflux mu ((head u0):(take 1 inner)) -- left boundary
rbc = zeroflux mu ((last u0):(take 1 $ 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
-- integrate until time t with time step dt
evalEuler :: (Floating a, Ord a) => a -> a -> a -> [a] -> [a]
evalEuler dt t mu u
| t <= 0.0 = u
| otherwise =
let prevstep = stepEuler mu u
in evalEuler (min dt (t-dt)) (t-dt) mu $! prevstep