differentiation. Reply

[Dimitre Novatchev]

> how about this:
> 
> diff f x = (f (x + epsilon) - f (x - epsilon)) / (2 * epsilon) 
>    where epsilon = 0.00001
> 
> every call to say "diff sin" will cost you about as much as every call to 
> "sin"

Following John Hughes' article "Why Functional Programming Matters", here's one
implementation:

> within :: (Num a, Ord a) => a -> [a] -> a
> within eps (a:b:xs)     | abs(a - b) <= eps = b
>                         | otherwise         = within eps (b:xs)

> relative :: (Num a, Ord a) => a -> [a] -> a
> relative eps (a:b:xs)    | abs(a - b) <= eps * abs(b) = b
>                          | otherwise                  = relative eps (b:xs)

> easydiff :: Fractional a => (a -> a) -> a -> a -> a
> easydiff f x h = (f(x + h) - f x) / h

> halve :: Fractional a => a -> a
> halve x = x /2

> log2 x = log(x) / log(2)

> elimerror :: Floating a => a -> [a] -> [a]
> elimerror n (a:b:xs) = (((b * (2 ** n) - a) / (2 ** n - 1)):elimerror n (b:xs))


> order :: (Num a, Floating b, RealFrac b) => [b] -> a
> order (x:y:z:xs) = fromInteger(round(log2((x-z)/(y-z) - 1)))

> improve xs = elimerror (order xs) xs


> differentiate :: Fractional a => a -> (a -> a) -> a -> [a]
> differentiate h0 f x = map (easydiff f x) (iterate halve h0)

> super :: (Floating a, RealFrac a) => [a] -> [a]
> super xs = map second (myRepeat improve xs)

> second (x:y:ys) = y

> e2 = within 0.001 (super (differentiate 0.1 (^3) 2))


> data Dif_args = Diffargs {eps, h0, xx :: Float, func :: (Float -> Float)}

> myDiff   :: Dif_args -> Float
> myDiff  df = within (eps df) (super (differentiate (h0 df) (func df) (xx df)))

> e3 = myDiff Diffargs {eps=0.001, h0=0.1, func= sin, xx=(3.1415926 / 2) }


I am new to Haskell and this was my first haskell application. Therefore, I'd
appreciate any comments on improving the style of expression or on any defficiencies
of this implementation.

Cheers,
Dimitre Novatchev.

__________________________________________________
Do You Yahoo!?
Send your FREE holiday greetings online!
http://greetings.yahoo.com