I've been learning Haskell recently and was talking to a friend who is working through SICP. We were curious to compare Common Lisp and Scheme and so I decided as an exercise to try to translate exercise 1.29 into Haskell.
This exercise uses a function sigma which represents the mathematical Summation function Sigma. This function takes a function f to apply to each term, a lower bound, a function to apply to each term to get the next term, and an upper bound. It returns the sum of f applied to each term.
simpsonIntegral is supposed to use Simpson's rule to approximate the integral of the function f over the range [a, b] using an "accuracy" n. I'm having trouble getting this function to work because there seems to be something I don't understand about the types involved.
This code will compile with version 6.12.1 of ghc but simpsonIntegral will be given a type context (Integral a, Fractional a) which doesn't make any sense and the function blows up as soon as you call it. I've got this working at one point but what I did was so obviously a hack that I wanted to ask here how this would be handled idiomatically.
How does one idiomatically handle the Integral -> Fractional/Real conversion needed in h? I read a number of things but nothing seemed obvious and clean.
sigma :: (Ord a, Num b) => (a -> b) -> a -> (a -> a) -> a -> b
sigma f a next b = iter a 0
where
iter current acc | current > b = acc
| otherwise = iter (next current) (acc + f current)
simpsonIntegral f a b n = 1.0 * (h / 3) * (sigma simTerm 0 (1+) n)
where
h = (b - a) / n
simTerm k = (yk k) * term
where
yk k = f (a + h * k)
term =
case k of
0 -> 1
1 -> 1
otherwise -> if odd k then 4 else 2
fromIntegral :: (Integral a, Num b) => a -> b
r = fromIntegral i
To follow up on Justice's answer: if you're curious about where to put the fromIntegral
s, the following compiles:
simpsonIntegral :: (Integral a, Fractional b) => (b -> b) -> a -> a -> a -> b
simpsonIntegral f a b n = 1.0 * (h / 3) * (sigma simTerm 0 (1+) n)
where
h = fromIntegral (b - a) / fromIntegral n
simTerm k = (yk k) * term
where
yk k = f (fromIntegral a + h * fromIntegral k)
term =
case k of
0 -> 1
1 -> 1
otherwise -> if odd k then 4 else 2
And seems to work:
*Main> simpsonIntegral (^3) 0 1 100
0.2533333233333334
*Main> simpsonIntegral (^3) 0 1 1000
0.2503333333323334
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With