Derivatives Lab

Derivatives Lab#

This lab allows students to experiment with different differentiation methods, including symbolic, finite difference, complex step differentiation, and automatic differentiation.

Symbolic Differentiation with SymPy#

We first try out symbolic differentiation using scipy. This is very similar to the class notes, except we will take derivatives of the \(\text{sinc}(x) = \sin(x) / x\).

# HANDSON: make sure that sympy is installed

!pip install sympy

import sympy as sp

Collecting sympy
  Downloading sympy-1.13.3-py3-none-any.whl.metadata (12 kB)

Collecting mpmath<1.4,>=1.1.0 (from sympy)
  Downloading mpmath-1.3.0-py3-none-any.whl.metadata (8.6 kB)
Downloading sympy-1.13.3-py3-none-any.whl (6.2 MB)
?25l   ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 0.0/6.2 MB ? eta -:--:--

   ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 6.2/6.2 MB 88.7 MB/s eta 0:00:00
?25hDownloading mpmath-1.3.0-py3-none-any.whl (536 kB)
?25l   ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 0.0/536.2 kB ? eta -:--:--
   ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 536.2/536.2 kB 45.1 MB/s eta 0:00:00
?25h

Installing collected packages: mpmath, sympy

Successfully installed mpmath-1.3.0 sympy-1.13.3

# HANDSON: use sympy to define the equation f(x) = sin(x) / x ...

x = ...
f = ...

# HANDSON: ... and get its derivative

f_prime = ...
f_prime_simplified = ...

# Display f'(x)

f_prime

Ellipsis

# Display simplified f'(x)

f_prime_simplified

Ellipsis

# We didn't do this in the notes,
# but it is possible to obtain numerical values of functions from sympy.

display(f.evalf(subs={x:1}))
display(f_prime_simplified.evalf(subs={x:1}))

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Cell In[6], line 4
      1 # We didn't do this in the notes,
      2 # but it is possible to obtain numerical values of functions from sympy.
----> 4 display(f.evalf(subs={x:1}))
      5 display(f_prime_simplified.evalf(subs={x:1}))

AttributeError: 'ellipsis' object has no attribute 'evalf'

# HANDSON: modify function `f(x)` and confirm that `sympy` is able to
# compute its derivatives

pass

To help visualize the results of derivative, let’s copy the plotting function from the class notes here:

import numpy as np
from matplotlib import pyplot as plt
from typing import Callable

def mkplot(g, X, Fp):
    if isinstance(g, Callable):
        f = g
    else:
        f = lambda x: g.evalf(subs={'x': x})
        
    Xd = np.linspace(min(X), max(X), num=1001)
    Fd = [f(x) for x in Xd]
    
    plt.plot(Xd, Fd, lw=5, alpha=0.25)
    for (x, fp) in zip(X, Fp):
        y = f(x)
        plt.plot(
            [x-0.05,    x+0.05],
            [y-0.05*fp, y+0.05*fp],
        )

X       = range(10)
F_prime = [f_prime.evalf(subs={'x':x}) for x in X]

mkplot(f, X, F_prime)

Enhancing Dual Number Autodiff#

In the class, we implemented a Dual-number based autodiff scheme. It supports many basic operators except division and power.

Try implementing these extra operators and test them out with our visualization tool.

def V(x):
    """Select the value from a dual number.

    Work for both python built-in numbers (often used in function) and dual numbers.
    """
    if isinstance(x, Dual):
        return x[0]
    else:
        return x

def D(x):
    """Select the derivative from a dual number.

    Work for both python built-in numbers (often used in function) and dual numbers.
    """
    if isinstance(x, Dual):
        return x[1]
    else:
        return 0

class Dual(tuple):
    """Dual number for implementing autodiff in pure python"""

    def __new__(self, v, d=1): # tuple is immutable so we cannot use __init__()
        return tuple.__new__(Dual, (v, d))

    def __add__(self, r):
        return Dual(
            V(self) + V(r),
            D(self) + D(r),
        )
    def __radd__(self, l):
        return self + l # addition commutes

    def __sub__(self, r):
        return Dual(
            V(self) - V(r),
            D(self) - D(r),
        )
    def __rsub__(self, l):
        return Dual(
            V(l) - V(self),
            D(l) - D(self),
        )

    def __mul__(self, r):
        return Dual(
            V(self) * V(r),
            D(self) * V(r) + V(self) * D(r),
        )
    def __rmul__(self, l):
        return self * l # multiplication commutes

    def __truediv__(self, r):
        return Dual(
            V(self) / V(r),
            ..., # HANDSON: implement chain-rule for division
        )
    def __rtruediv__(self, l):
        return Dual(
            V(l) / V(self),
            ..., # HANDSON: implement chain-rule for division
        )

    def __pow__(self, r): # assume r is constant
        if r == 0:
            return ... # HANDSON: implement chain-rule for power
        elif r == 1:
            return ... # HANDSON: implement chain-rule for power
        else:
            return Dual(
                V(self)**r,
                ..., # HANDSON: implement chain-rule for power
            )

def sin(x):
    return Dual(
        np.sin(V(x)),
        np.cos(V(x)) * D(x)  # chain rule: d/dx sin(x) = cos(x) * x'
    )

def f(x):
    return (x**2 + 3) / x
    # return sin(x) / x

X     = np.linspace(1,2,num=11)
F, Fp = f(Dual(X))

mkplot(f, X, Fp)
Contents