AMath 583 Lab 4 -- Sample solution

Where to find this notebook:

This notebook will be in \$UWHPSC/labs/lab4/lab4b.ipynb if you have cloned the class repository. You may need to "git pull" to update.

You do not need to do the next command when using the IPython notebook on SageMathCloud, but it may be needed if you open this notebook on another system where the defaults are set differently. It does "from pylab import *" so numpy and plotting commands are available, and sets it up for plots to appear in the notebook webpage



In [4]:

    
%pylab inline









    



Populating the interactive namespace from numpy and matplotlib

Newton's method



In [5]:

    
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

New problem:

Find the two points where the circle and line intersect in the plot below. The equation of the circle is $(x-1)^2 + y^2 = 4$ and the curve $y = \cos(x)$. Combine these equations by plugging $\cos(x)$ in for $y$ in the equation for the circle. This gives a single nonlinear equation in $x$ that you can solve using Newton's method.

Try to find both intersections and then replot the circle with these two points as red dots.



In [15]:

    
theta = linspace(0, 2*pi, 1000)
x = 1 + 2*cos(theta)
y = 2*sin(theta)
plot(x,y,'b')

x = linspace(-2,4,1000)
y = cos(x)
plot(x,y,'b')
axis('scaled')  # so circle looks circular









    Out[15]:





(-2.0, 4.0, -2.0, 2.0)

Solution:

Define the function and derivative and plot f(x):



In [16]:

    
def fvals(x):
    f = (x-1)**2 + cos(x)**2 - 4.
    fprime = 2*(x-1) - 2*cos(x)*sin(x)
    return f,fprime



In [17]:

    
x = linspace(-2,4,1000)
fx, fpx = fvals(x)
plot(x,fx)
plot(x,zeros(len(x)), 'k')  # x-axis









    Out[17]:





[<matplotlib.lines.Line2D at 0x3649110>]

Find the two roots via two different starting values:



In [18]:

    
x0 = -1.
xk1 = newton(fvals, x0, 1e-10)









    



   k               xk                      f(xk)
   0       -1.000000000000000       0.291926581726429
   1       -0.905546853890051       0.012113075692203
   2       -0.901281425757293       0.000022473588790
   3       -0.901273482597076       0.000000000077586



In [19]:

    
x0 = 2.
xk2 = newton(fvals, x0, 1e-10)









    



   k               xk                      f(xk)
   0        2.000000000000000      -2.826821810431806
   1        3.025398741927668       1.088799688651651
   2        2.771071614076918       0.005577410491355
   3        2.769749075408306       0.000000459718415
   4        2.769748966379979       0.000000000000004

Plot the solutions found:



In [20]:

    
theta = linspace(0, 2*pi, 1000)
x = 1 + 2*cos(theta)
y = 2*sin(theta)
plot(x,y,'b')

x = linspace(-2,4,1000)
y = cos(x)
plot(x,y,'b')
axis('scaled')  # so circle looks circular

plot([xk1,xk2], [cos(xk1),cos(xk2)],'ro')









    Out[20]:





[<matplotlib.lines.Line2D at 0x34e4150>]



In [20]: