I’ve encountered a strange behavior when calculating numerical derivatives, where I don’t know what the problem is. Maybe someone could help. I’ve tried the following code in the Haskell playground. Every ghc version shows the same results. If I write two different versions of the derivative function (one with forward and one with backward difference) and try different combinations when calculating the second derivative, I get correct results.
derivative :: (Double -> Double) -> Double -> Double
derivative f x = cases where
h = 0.000001
delta_f_plus = f (x + h) - f x -- forward difference
delta_f_minus = f x - f (x - h) -- backward difference
cases
| abs delta_f_plus <= abs delta_f_minus = delta_f_plus/h
| otherwise = delta_f_minus/h
derivativeF :: (Double -> Double) -> Double -> Double
derivativeF f x = cases where
h = 0.000001
delta_f_plus = f (x + h) - f x
cases = delta_f_plus/h
derivativeB :: (Double -> Double) -> Double -> Double
derivativeB f x = cases where
h = 0.000001
delta_f_minus = f x - f (x - h)
cases = delta_f_minus/h
main :: IO ()
main = do
let f x = (cos x)^2
-- Calculate f'' 0, must be approx -2:
-- This works:
let f' x = derivativeF f x
let f'' x = derivativeF f' x
putStrLn $ show (f'' 0)
-- This works:
let f' x = derivativeF f x
let f'' x = derivativeB f' x
putStrLn $ show (f'' 0)
-- This works:
let f' x = derivativeB f x
let f'' x = derivativeF f' x
putStrLn $ show (f'' 0)
-- This works:
let f' x = derivativeB f x
let f'' x = derivativeB f' x
putStrLn $ show (f'' 0)
-- This works not:
let f' x = derivative f x
let f'' x = derivative f' x
putStrLn $ show (f'' 0)