Calculating Derivatives

May 5, 2017

This is easy in a language that supports higher-order functions, and either difficult or impossible in one that doesn’t. Here’s a Python solution:

Python 2.7.10 (default, Jun  1 2015, 18:05:38)
[GCC 4.9.2] on cygwin
Type "help", "copyright", "credits" or "license" for more information.
>>> def deriv(f):
...     dx = 0.0001
...     return lambda x : (f(x+dx) - f(x)) / dx
...
>>> def cube(x):
...     return x*x*x
...
>>> cube_deriv = deriv(cube)
>>> print(cube_deriv(2))
12.00060001
>>> print(cube_deriv(3))
27.0009000101
>>> print(cube_deriv(4))
48.0012000099

You can run the program in Scheme or Python.

The real point of this exercise is to get you to look at Joe Marshall’s paean to Scheme. If your language doesn’t make this as easy as Scheme, maybe you need a new language.

Pages: 1 2

Posted by programmingpraxis
Filed in Exercises
8 Comments »

8 Responses to “Calculating Derivatives”

  1. Rutger said
    May 5, 2017 at 9:01 AM 
    import math

def derivative(f,x, dx = 10**-8):
	return (f(x+dx) - f(x))/dx

def area_circle(r):
	return math.pi*r**2

for i in range(1,10):
	d, p = derivative(area_circle,i), math.pi*2*i
	print(i, d,p,d-p)

# output
# 1 6.283185260969049 6.283185307179586 -4.6210537618662784e-08
# 2 12.566370521938097 12.566370614359172 -9.242107523732557e-08
# 3 18.849555871724988 18.84955592153876 -4.981377088597583e-08
# 4 25.132741399147562 25.132741228718345 1.7042921740539896e-07
# 5 31.41592657129877 31.41592653589793 3.540083781672365e-08
# 6 37.69911245399271 37.69911184307752 6.109151939881485e-07
# 7 43.98229691560118 43.982297150257104 -2.3465592136062696e-07
# 8 50.26548706155154 50.26548245743669 4.604114849371399e-06
# 9 56.54867436533095 56.548667764616276 6.600714677063024e-06
  2. john said
    May 5, 2017 at 9:03 AM

    Using C11 and function pointers:


    #include <stdio.h>
    #include <stdlib.h>

    double deriv(double (*f)(double), double x);

    double cube(double x);

    int main(int argc, char **argv) {
      printf("%f\n", deriv(cube, 2));
      printf("%f\n", deriv(cube, 3));
      printf("%f\n", deriv(cube, 4));
      
      exit(0);
    }

    double deriv(double (*f)(double), double x) {
      const double dx = 0.0000001;
      const double dy = f(x + dx) - f(x);
      return dy / dx;
    }

    double cube(double x) {
      return x * x * x;
    }

  3. Paul said
    May 5, 2017 at 9:44 AM

    Using a slightly different formula.

    import sys
import math

EPS = sys.float_info.epsilon
EPS3 = EPS ** (1/3)

def der(f, x):
    h = EPS3 * abs(x) if x != 0 else EPS3
    xm, xp = x - h, x + h
    return (f(xp) - f(xm)) / (xp - xm)

    Results using the same example as Rutger:

    1  6.283185307  6.283185307 -2.06173e-11
2 12.566370614 12.566370614 -4.12346e-11
3 18.849555921 18.849555922 -4.61107e-11
4 25.132741229 25.132741229 -8.24691e-11
5 31.415926536 31.415926536 -6.01545e-11
6 37.699111843 37.699111843 -9.22213e-11
7 43.982297150 43.982297150 -1.91349e-10
8 50.265482457 50.265482457 -1.64938e-10
9 56.548667764 56.548667765 -2.31651e-10
  4. Josef said
    May 5, 2017 at 2:51 PM

    Haskell

    dx = 0.0000001

derive f = let f' x = (f (x + dx) - f (x)) / dx
           in  f'

cube x = x * x * x

test = do
  print (derive cube 2)
  print (derive cube 3)
  print (derive cube 4)
  5. bavier said
    May 5, 2017 at 4:19 PM

    Gforth

    : cube ( F: f -- f^3)  fdup fdup f* f* ;
: f' ( xt --) ( F: f -- f'-of-x)
   dup fdup  1e-7 f+ execute  fswap execute 
   f-  1e-7 f/ ; 

: test  10 0 DO
     ['] cube  I s>d d>f fdup f. ." => "  f' f.	CR
   LOOP ;
test  bye

    Prints

    0. => 0.00000000000001 
1. => 3.00000030151182 
2. => 12.0000005843224 
3. => 27.0000008484317 
4. => 48.000001413584 
5. => 75.0000016580543 
6. => 108.00000211475 
7. => 147.000002925779 
8. => 192.000001106862 
9. => 243.000000637039
  6. Globules said
    May 6, 2017 at 7:15 PM

    Here’s another Haskell version, similar to Josef’s. This example shows that you can use it with various numeric types. For example, using the “constructive reals” (CReal) from the “numbers” package we can specify an arbitrary precision. (Also, note that ^3 is shorthand for the cube function, \x -> x^3.)

    import Data.Number.CReal

deriv :: Fractional a => a -> (a -> a) -> a -> a
deriv dx f x = (f (x + dx) - f x) / dx

main :: IO ()
main = do
  print $ deriv 0.0000001 (^3) (2 :: Double)
  print $ deriv 0.0000001 (^3) (2 :: Rational)
  print $ deriv 0.000000000000000000000000000001 (^3) (2 :: CReal)

    $ ./derivlim 
12.000000584322379
1200000060000001 % 100000000000000
12.000000000000000000000000000006
  7. Jussi Piitulainen said
    May 7, 2017 at 12:31 PM 
    julia> deriv(f; dx = 0.0000001) = x -> (f(x + dx) - f(x)) / dx
deriv (generic function with 1 method)

julia> "$(deriv(x -> x^3 + 1)(4)) ≈ $(3*4^2)"
"48.00000141358396 ≈ 48"

julia> "$(deriv(sin)(3.1)) ≈ $(cos(3.1))"
"-0.999135150725472 ≈ -0.9991351502732795"

julia> "$(deriv(sin; dx = 0.001)(3.1)) ≈ $(cos(3.1))"
"-0.9991557740801349 ≈ -0.9991351502732795"
  8. Steve said
    May 8, 2017 at 8:20 PM

    MUMPS V1:

    EXTR ; Test of calculating derivatives
    N I
    F I=0:1:9 W !,I,” –> “,$$DERIV(“CUBE”,I,.0000001)
    Q
    ;
    CUBE (N) ; Cube a number
    Q N*N*N
    ;
    DERIV (FUNC,VAL,DX) ;
    Q @(“$$”_FUNC_”(VAL+DX)”)-@(“$$”_FUNC_”(VAL)”)/DX

    MCL> D ^EXTR

    0 –> 0
    1 –> 3.0000003
    2 –> 12.0000006
    3 –> 27.0000009
    4 –> 48.0000012
    5 –> 75.0000015
    6 –> 108.0000018
    7 –> 147.0000021
    8 –> 192.0000024
    9 –> 243.0000027

Leave a comment