AMath 583 Lab 5

Where to find this notebook:

This notebook will be in \$UWHPSC/labs/lab5/lab5a.ipynb if you have cloned the class repository. You may need to "git pull" to update.
Or you can bring it down to your computer (or SageMathCloud project) by typing this in a terminal shell:
- wget http://faculty.washington.edu/rjl/classes/am583s2014/lab5a.ipynb

Announcements:

Please sit at assigned tables.
If you're having problems with notebooks in SMC, note that you can also use IPython from a terminal.
Homework 2 will be posted by Thursday. Part of the homework will be to read this paper: "Best Practices for Scientific Computing, by G. Wilson, D. A. Aruliah", by C. T. Brown, et. al. that can be found at http://www.plosbiology.org/article/info%3Adoi%2F10.1371%2Fjournal.pbio.1001745



In [ ]:

    
%pylab inline

SymPy

SymPy is a package of Python tools for symbolic calculations in algebra and calculus. This is notebook has a very brief demo of some features, but you can learn much more about it from various on-line resources such as:

You can do "from sympy import *" to get all of sympy in your namespace, but note that this will overwrite standard numpy variables such as pi and functions such as sin, cos, exp, sqrt, with sympy versions. Since we may want to mix symbolic and numerical computing, we will be explicit about what's coming from SymPy:



In [ ]:

    
import sympy as S  # shorthand



In [ ]:

    
S.init_printing()   # so symbolic math is printed nicely

Define a polynomial and factor it:



In [ ]:

    
x = S.symbols('x')  # defines a symbol
f = x**3 - 3*x**2 - 18*x + 40   # a new symbolic expression
f                   # print it nicely



In [ ]:

    
S.factor(f)

We can also differentiate:



In [ ]:

    
S.diff(f,x)  # differentiate f with respect to x

A messier function and its derivative:



In [ ]:

    
f = (x**3 * S.exp(5*x)*S.cos(S.pi*x)) / (1 + S.sqrt(x))
f



In [ ]:

    
g = S.diff(f,x,n=3)
g

Note that if you "print g" it does not come out so pretty, but this might be more useful if you want to cut and paste this into a Fortran program, for example:



In [ ]:

    
print g

Evaluate this derivative at the point $x = 0.2$:



In [ ]:

    
g2 = g.subs(x, 0.2)
g2

Note this substituted 0.2 for $x$, but the special symbol $\pi$ has not be converted to a number. The SymPy function N converts special constants like $\pi$ and $e$ to numerical values:



In [ ]:

    
S.N(g2)

You can specify how many digits to use when it does this substitution, e.g.



In [ ]:

    
S.N(S.pi,n=50)

Symbolic differentiation for a Newton iteration.

In the last two labs you implemented Newton's method for different functions and probably had to compute the derivatives by hand. Here's a way to write and fvals function using SymPy:



In [ ]:

    
def fvals(x, debug=False):
    from sympy import symbols,diff,sqrt
    # First specify f symbolically and differentiate using SymPy:
    xs = symbols('xs')
    fs = xs**2 - 2.
    fprimes = diff(fs, xs)
    
    # Now evaluate numerically at the value x passed in to this function:
    f = fs.subs(xs, x)
    fprime = fprimes.subs(xs, x)
    
    # The next lines are just for illustrating that this is working:
    if debug:
        print "fs = ",fs
        print "fprimes = ",fprimes
        print "x = ",x
        
    return f, fprime

# Try it out:
fv = fvals(3., debug=True)
print "fvals returns: ", fv

Try this out with Newton's method:



In [ ]:

    
def newton(fvals, x0, tol):
    xk = x0
    kmax = 30   # maximum number of iterations
    print "   k               xk                      f(xk)"
    for k in range(kmax):
        
        fxk, fpxk = fvals(xk)  # evaluate f(x) and f'(x)
        print "%4i   %22.15f  %22.15f"  % (k, xk, fxk)

        if abs(fxk) < tol:
            break  #leave the loop if tolerance satisfied
            
        xk = xk - fxk / fpxk  # update xk using Newton's method

    return xk



In [ ]:

    
newton(fvals, 2, 1e-10)

Exercises:

Use Newton's method to find a root of $f(x) = (x^2 - 2)\exp(-0.1 x^2) + 0.5$ from Lab 4.
Let $f(x) = \sqrt{\cos(2\pi x) e^x}$. Compute $f''(0.1)$, the second derivative evaluated at $x=0.1$. This value is needed for the Lab 5 Quiz.

Remember you can get help in IPython by adding a ? to the end of an object and running the cell...



In [ ]:

    
S.integrate?



In [ ]:

    
S.diff?