### Numerical differentiation : Classes # In[1]: # Imports import math import cmath # ### Classes # # Classes are a tool to make the programs modular and hence reusable. Also they lend us a bit of abstraction which makes coding a bit easier. To give a simple example, one can define a class `Quadratic` which codes quadratic polynomials. Calling something like `q = Quadratic(a, b, c)` would define a function `q` such that `q(x) = ax^x + bx + c`. Classes can be used for more complicated things like coding and drawing geometric objects and so on. # # Classes have their own set of variables, which we call *attribute*s and functions called *method*s. We have already encountered *method*s in our course, figure out which. # # **All the material in this lecture are in the 7th chapter of Langtangen's *A Primer on Scientific Programming with Python***. # # First let us code the `Quadratic` class to gain some experience. Note that, though it is not technically necessary, all the classes are conventionally named using a capitalized word. # In[2]: class Quadratic : """This class creates quadratic equations. Attributes : a, b, c : Coefficients of ax^2 + bx + c Methods : valueat(x) : Value of the quadrating at x. roots() : Returns a tuple containing the two roots of ax^2 + bx + c = 0 display() : Print the quadratic equation. """ # The doc string is similar to that of functions. # In a class, we need a function to initialize the class. def __init__(self, a, b, c) : """Constructor function for the class Quadratic.""" # I'll explain the variable self below. self.a = a self.b = b self.c = c def valueat(self, x) : """This computes the value of the quadratic at x.""" # In the following line note how we refer to the elements # a, b,c of the quadratic as self.a, self.b etc. self # is the instance of the current class. return (self.a)*x*x + (self.b)*x + (self.c) def roots(self) : """Returns the roots of the quadratic.""" a = self.a b = self.b c = self.c if a != 0 : disc = b*b - 4*a*c if disc >= 0 : r1 = float(-b + math.sqrt(disc))/2*a r2 = float(-b - math.sqrt(disc))/2*a else : r1 = (-b + cmath.sqrt(disc))/2*a r2 = (-b - cmath.sqrt(disc))/2*a retval = (r1, r2) else : if b != 0 : retval = - float(c)/b else : print "No solutions." retval = None return retval def display(self) : """Returns a string printing the quadratic.""" str = "%g x^2 + %g x + %g" % (self.a, self.b, self.c) return str # Now we use the class to do some computations. For fun we define another function to do the computation. # In[3]: def solve_quads(a, b, c) : q = Quadratic(a, b, c) # We call q an instance of the class Quadratic. print "q = ", q.display(), "has", q.roots(), "as roots. q(", print math.pi,") = ", q.valueat(math.pi) return None solve_quads(1, 2, 1) solve_quads(1, 0, 1) solve_quads(1, 1, 0) # Out[3]: # q = 1 x^2 + 2 x + 1 has (-1.0, -1.0) as roots. q( 3.14159265359 ) = 17.1527897083 # q = 1 x^2 + 0 x + 1 has (1j, -1j) as roots. q( 3.14159265359 ) = 10.8696044011 # q = 1 x^2 + 1 x + 0 has (0.0, -1.0) as roots. q( 3.14159265359 ) = 13.0111970547 # # ### Making classes callable. # We eventually want to diffentiate and integrate functions. We'll implement the as classes. What we want to do is to make the instances behave like functions. We do that using the `__call__` method. As an example, we just copy the `Quadratic` class to `Quadratic2` renaming the `valueat` method to be `__call__`. # In[4]: class Quadratic2 : """This class creates quadratic equations. Attributes : a, b, c : Coefficients of ax^2 + bx + c Methods : __call__(x) : Value of the quadrating at x. roots() : Returns a tuple containing the two roots of ax^2 + bx + c = 0 display() : Print the quadratic equation. """ # The doc string is similar to that of functions. # In a class, we need a function to initialize the class. def __init__(self, a, b, c) : """Constructor function for the class Quadratic.""" # I'll explain the variable self below. self.a = a self.b = b self.c = c def __call__(self, x) : """This computes the value of the quadratic at x.""" # In the following line note how we refer to the elements # a, b,c of the quadratic as self.a, self.b etc. self # is the instance of the current class. return (self.a)*x*x + (self.b)*x + (self.c) def roots(self) : """Returns the roots of the quadratic.""" a = self.a b = self.b c = self.c if a != 0 : disc = b*b - 4*a*c if disc >= 0 : r1 = float(-b + math.sqrt(disc))/2*a r2 = float(-b - math.sqrt(disc))/2*a else : r1 = (-b + cmath.sqrt(disc))/2*a r2 = (-b - cmath.sqrt(disc))/2*a retval = (r1, r2) else : if b != 0 : retval = - float(c)/b else : print "No solutions." retval = None return retval def display(self) : """Returns a string printing the quadratic.""" str = "%g x^2 + %g x + %g" % (self.a, self.b, self.c) return str # In[5]: q = Quadratic2(1, -3 , 2) print q(10) # Out[5]: # 72 # # ### Numeric Differentiation # # We know that the mathematical definition of the derivative $f'$ of a function $f$ at $x$ is # # \begin{equation} # \lim_{h \to 0} \frac{f(x + h) - f(x)}{h} # \end{equation} # # Instead of finding limit, we can code the derivative as just the fraction $(f(x+h) - f(x))/h$ for small $h$. *Assuming* that as we reduce $h$, the value of the function will be closer and closer to $f'$, we can also write a check if reducing $h$, say by half, doesn't change the value of the fraction beyond some acceptible error. Then we can use that value as an approximate value for the derivative. # In[6]: class SimpleDerivative : """Implements a callable class for a naive numerical derivative.""" def __init__(self, f, h = 1.0E-3) : self.f = f self.h = float(h) def __call__(self, x) : """Return the value of the derivative at x.""" f = self.f h = self.h der = (f(x + h) - f(x)) / h return der def set_deviation(self, val) : """Set the deviation.""" self.h = val # In[7]: def f(x) : return x*x df = SimpleDerivative(f,.01) print df(2) df.set_deviation(df.h/100) print df(2) # Out[7]: # 4.01 # 4.00010000001 # # Now to end, we modify the class so that when we call, it automatically keeps on halving the deviation till two consecutive derivatives differ by a very small amount. Again as in Newton's method we have to keep track of how many iterations we are doing. # In[8]: class NaiveDerivative : """This uses a loop to compute derivatives with smaller smaller deviations till the difference between successive computations is less than a predetermined error Attributes : f : function whose derivative we seek, h : deviaion err : error N : max number of interations allowed Methods : single_der(x) : returns (f(x+h) - f(x)) / h set_dev(d) : sets the deviation to d set_err(e) : sets the error allowed to e allow_iter(M) : allow M iterations __init__ : Constructor __call__ : the actual code """ def __init__(self, f, h=0.01, err=1E-5, N=10000) : self.f = f self.h = float(h) self.err = err self.N = N def set_dev(self, d) : self.h = d return None def set_err(self, e) : self.err = e return None def allow_iter(self, M) : self.N = M return None def single_der(self, x) : f = self.f h = self.h return (f(x + h) - f(x))/h def __call__(self, x) : N = self.N err = self.err dfold = self.single_der(x) self.set_dev(self.h / 2.0) dfnew = self.single_der(x) no_iter = 0 while no_iter < N and abs(dfold - dfnew) >= err : dfold = dfnew self.set_dev(self.h / 2.0) dfnew = self.single_der(x) no_iter += 1 if abs(dfold - dfnew) >= err : print "Could not converge after %d iterations." % no_iter retval = None else : retval = dfnew return retval # In[9]: def myfn(x) : return x*x*x + x dmyfn = NaiveDerivative(myfn) dmyfn.set_dev(1) dmyfn.set_err(.1) print dmyfn(1) dmyfn.set_err(1E-5) print dmyfn(1) # Out[9]: # 4.0947265625 # 4.00000572205 #