CS357 NUMERICAL METHODS I

Summer 2014

Machine Problem 2

Posted: Friday  July 5^th, 2014

Due Date: 9:00 P.M., Friday, July 18^th, 2014

This MP has three parts.

Handin is available here.

You will have to submit (using svn) the following files:

newton.py

halley.py

newtoni.py

newtonm.py

intervalnewton.py

fractal.py

 parta.txt (answer to question Q1 in Part A)

 Preparation

 Download the tar file  mp2.tar.gz to your Linux directory and open the file using the Linux command  

tar -xzvf mp2.tar.gz

(26 points) Part A: Finding Roots for a function f: R ->R

Note: Due to various forms of roundoff error incurred by writing code in different ways your results may differ slightly from the test results
shown below. For this reason, you must send us your .py files.

Section 1. Newton's method

Write a Python function named newton. This function should use newton's method to find a root of the function y = f(x)  given x = guess. The function named newton must have five input parameters in the order shown below.

Implementation (newton.py file)

 Use the following function header.

 def newton(f, df, x, tol, display=False):

where f is the function, f: R -> R.
            df is the derivative of the function
            x is an approximation of the root and is a scalar value
            tol specifies that the search will iterate until  | f(x_k)| < tol. You may assume that tol > 0.
            display has the value True  or False, True means display the values of the individual iterations, False means don't display.
                         If the user doesn't submit a value for the fifth argument then the default value is display = False. ( This is the purpose of display=False as your last argument in the function definition line.)

                         Use the following code to display values

               if display:
                   print('iteration'
                              '\tx(iteration)'
                              '\t\terror(iteration)')

            and
              print('{0:5d} {1:25.17f} {2:25.17f}'.format(counter, x, error))  <=== where you use your own three variable names and error equals abs(f(x))
          

You must have at least one comment in your program. Copy the function definition line (header) and paste this as a comment immediately below the function definition line.

You must count the number of iterations in your code.

We assume that the x is counted as the first iteration. If the numer of iterations equals twenty then terminate your program and display the message 'no convergence' and newton should return the last value of its iteration.

Tests

At the Linux prompt type:

$ module load python3

$ python3  newton.py

Your code will also be graded by passing the following tests:

root = 0 
 (grader:  1 point)
( result for running
    root = newton(np.sin, np.cos, 0, 1.0e-1)
    print ('root =', root, end='\n\n')   )

root = 0
(grader:  1 point)
( result for running
   root = newton(np.sin, np.cos, 0, 1.0e-1, False)
    print ('root =', root, end='\n\n')  )

iteration      x(iteration)                            error(iteration)
    1       0.00000000000000000       0.00000000000000000
root = 0
(grader:  2 points)
( result for running
     root = newton(np.sin, np.cos, 0, 1.0e-1, True)
    print ('root =', root, end='\n\n'))

iteration      x(iteration)                           error(iteration)
    1       0.50000000000000000       0.47942553860420301
    2      -0.04630248984379048       0.04628594680720108
    3       0.00003311802132639       0.00003311802132034
    4      -0.00000000000001211       0.00000000000001211
root = -1.21079836153e-14
(grader:  2 points)
( result for running
    root = newton(np.sin, np.cos, 0.5, 1.0e-6, True)
    print ('root =', root, end='\n\n'))


iteration       x(iteration)                         error(iteration)
    1       0.50000000000000000       0.46364760900080609
    2      -0.07955951125100758       0.07939228286534084
    3       0.00033530220400547       0.00033530219143974
    4      -0.00000000002513147       0.00000000002513147
    5       0.00000000000000000       0.00000000000000000
root = 0.0

(grader:  2 points)
( result for running
    root = newton(np.arctan, lambda x: 1/(1 + x**2), 0.5, 1.0e-11, True)
    print ('root =', root, end='\n\n'))

root = 0.652922853427
(grader:  2 points)
( result for running
    root = newton(lambda x: np.exp(-x**2) - x, lambda x: -2*x*np.exp(-x**2) -1, 1.5, 1.0e-3)
    print ('root =', root, end='\n\n'))

iteration    x(iteration)        error(iteration)
    1       1.50000000000000000       5.06250000000000000
    2       1.12500000000000000       1.60180664062500000
    3       0.84375000000000000       0.50682163238525391
    4       0.63281250000000000       0.16036153212189674
    5       0.47460937500000000       0.05073939102294389
    6       0.35595703125000000       0.01605426044085334
    7       0.26696777343750000       0.00507966834261375
    8       0.20022583007812500       0.00160723881153013
    9       0.15016937255859375       0.00050854040521071
   10       0.11262702941894531       0.00016090536258620
   11       0.08447027206420898       0.00005091146238079
root = 0.08447027206420898
(grader:  2 points)
( result for running
    root = newton(lambda x: x**4, lambda x: 4*x**3, 1.5, 1.0e-4, 1)
    print ('root =', root, end='\n\n'))

Section 2. Halley's method(as in Halley's comet)

 Let's try a higher order method. Remember from Newton's method we started with a Taylor Series,

f(x_k+1)  ₌  f( x_k) + f'( x_k)* ( x_{k+1  -}  x_k ) + higher order terms

and we discarded the higher order terms and set  f(x_k+1) = 0 and then solved for x_k+1  to obtain Newton's formula,

x_{k+1   =} x_k  -   f(x_k ) / f'(x_k )

In Halley's Method we will keep one more term of the Taylor Series approximation,

f(x_k+1)  ₌  f( x_k) + f'( x_k)* ( x_{k+1  -}  x_k ) + 1/2 *f''( x_k)* ( x_{k+1  -}  x_k )²   + higher order terms

next we will discard the higher order terms, set f(x_k+1) = 0  and then factor (x_{k+1  -}  x_k ) to obtain,

0 ₌  f( x_k) + (f'( x_k) + 1/2 *f''( x_k)* ( x_{k+1  -}  x_k )) * ( x_{k+1  -}  x_k )

We next approximate the inner term ( x_{k+1  -}  x_k )  by using Newton's Method  x_{k+1   =} x_k  -   f(x_k ) / f'(x_k ) or x_k+1  -  x_k = - f(x_k ) / f'(x_k ) and obtain,

0 ₌  f(x_k) + (f'( x_k)  - 1/2 *f''( x_k)*  f(x_k ) / f'(x_k )  ) * ( x_{k+1  -}  x_k )

and now, solving the above for x_k+1  gives Halley's formula,

x_{k+1   =} x_k  -   2*f(x_k )_* f'(x_k ) /  ( 2* (f'(x_k ))² - f''(x_k )*f(x_k ) )

Implementation (halley.py file)
Assuming that you correctly coded the newton function in Section 1, copy and paste the body of your function to the provided file halley.py and make the necessary changes to change the program to implement Halley's Method.

  Note that the function named halley has one more input parameter ddf  than your newton function, which will contain the second derivative of f.

 def halley(f, df, ddf, x, tol, display=False):

Don't forget to modify the comment and make it agree with the new function halley.

 
Tests

At the Linux prompt type:

$ module load python3

$ python3  halley.py

Your code will also be graded by passing the following tests:

 =====================================================IMPORTANT============================================================

 Q1 (grader:  2 points) Indicate the order of convergence for Halley's Method for roots f(r) = 0 where f'(r) is not equal to 0. Do NOT simply state the order, you must show by example why the order is the value you state. Hint: Try to find a function that takes several iterations to converge. Write your answer in a file named parta.txt. You will handin parta.txt along with your other files. You can use gedit to create this file. Make sure that this file is in the same directory as your .py files when you handin.

Section 3. Newton's Method using an approximation to the derivative.

The secant method approximates the derivative. One problem with the secant method is that we must compute the slope of the secant line which may involve subtracting terms that are almost equal. We can overcome this by using the following trick. Assuming that j represents the imaginary number, the square root of -1, every complex number can be written as  a+b*j (or just a+bj) where a and b are real numbers. Python implements complext numbers. Note that Python uses 'j' to represent the square root of -1, and so we have 'j' here, as opposed to the more commonly used 'i'.  For example,

>>> z = complex(2.0, 3.0)
>>> z
(2+3j)
>>> z.imag
3.0
>>> z.real
2.0
>>>

Consider the Taylor Series of the expression  f(x+hj) and forget that hj is complex. We write,

f(x+hj) = f(x)  + f'(x)(hj) + (f''(x)/2!)(hj)² + (f'''(x)/3!)(hj)³ + (f''''(x)/4!)(hj)⁴ + (f''''(x)/5!)(hj)⁵ +...

and using the fact that j² = -1, we can re-write the above simplifying the powers of j,            

f(x+hj) = f(x)  + (f'(x) h ) j   -  (f''(x)/2!)(h)² -  ((f'''(x)/3!)(h)³ ) j+ (f''''(x)/4!)(h)⁴ + ((f''''(x)/5!)(h)⁵ ) j+...

and collecting real and imaginary parts,

f(x+hj) =  f(x) -  (f''(x)/2!)(h)² + (f''''(x)/4!)(h)⁴ - ...                 (real numbers)
                + (f'(x) h -  ((f'''(x)/3!)(h)³  + (f''''(x)/5!)(h)⁵ + ... ) * j

so we can say that

(f(x+hj)/h).imag  = (f'(x) h -  ((f'''(x)/3!)(h)³  + (f''''(x)/5!)(h)⁵ + ... )  / h =   f'(x) -  ((f'''(x)/3!)(h)²  + (f''''(x)/5!)(h)⁴ + ...

and thus

(f(x+ih)/h).imag approx equals f'(x)




 Implementation (newtoni.py file)

Assuming that you correctly coded the newton function in Section 1, copy and paste the function body into the provided file newtoni.py and make the necessary changes to change the program to implement this method.

def newtoni(f, x, tol,  display=False, h = 1.0e-16):

Don't forget to modify the comment and make it agree with the new function newtoni.

Use a default value of h = 1.0e-16 (user can specify a different value by passing h while calling the function) for computing the derivative using (f(x+ih)/h).imag.
Use the Python function named complex to compute   x+hj ,e.g.  complex(x,h).
If z = a + bj is a complex number then use the Python expression z.imag to extract the imaginary part  b of the number.
The abs function works for both real and complex numbers. The abs(z) = sqt(a**2 + b**2)where z = a + bj

Tests

At the Linux prompt type:

$ module load python3

$ python3  newtoni.py

(amended: 

At the Linux prompt type:

$ python /usr/bin/easy_install --user pyinterval

)

Your code will also be graded by passing the following tests:

Section 4. Newton's Method using an approximation to the derivative and improving the order of convergence when f'(root) = 0.

We showed in lecture that Newton's method converges at quadratic order when f'(root) is not equal to zero. But if  f'(root) = 0 then as seen in the last example above in the Section 3 we didn't get quadratic order. It can be shown that given y = f(x) and f(root) = 0 , f'(root) = 0, then if we define g(x) = f(x)/f'(x) we have g(root) = 0 (in the limit since f(root)/f'(root) = 0/0) and g'(root) does not equal zero.

Let's modify the results from Section 3 and compute the roots of g(x) = f(x)/f'(x). To make this example simple we will always use g(x) to compute
a root whether or not f'(root) = 0.

Implementation (newtonm.py file)

Assuming that you correctly coded the newtoni  function in Section 3, copy and paste the body of your function into newtonm.py and make the necessary changes to change the program to implement this method.

  Name the function newtonm .

 def newtonm(f, df, x, tol, display=False, h=1.0e-16):

Don't forget to modify the comment and make it agree with the new function newtonm .

Use a fixed value of h = 1.0e-16.
Do not bother to trap the error when f'(x_k) = 0.

Create an auxiliary function with the following code.

g = lambda x: f(x)/df(x)



Important: You will need to make several changes in your code since we are evaluating the function g and not f.

Tests

At the Linux prompt type:

$ module load python3

$ python3  newtonm.py

Your code will also be graded by passing the following tests:

iteration    x(iteration)        error(iteration)
    1       1.50000000000000000       0.37500000000000000
    2       0.00000000000000022       0.00000000000000006
root = 2.220446049250313e-16
(grader:  2 points)
( result for running
  root = newtonm(lambda x: x**4, lambda x: 4*x**3, 1.5, 1.0e-4, 1)
    print ('root =', root, end='\n\n'))

iteration    x(iteration)        error(iteration)
    1       5.00000000000000000    1350.34732126256380980
    2       4.82355247086866257     513.35142443708184601
    3       4.63469110445832122     196.11807400937146895
    4       4.43039809285584951      75.40760621598593616
    5       4.20631358998010896      29.24268852070842684
    6       3.95591163446007688      11.46776359049319538
    7       3.66919413964655483       4.55891144393153080
    8       3.33120185308977623       1.83412663090670769
    9       2.92410558774785478       0.72987134867967063
   10       2.45598204476744719       0.25342468810375218
   11       2.08072509932272753       0.04049434791362784
   12       2.00052320660784932       0.00026160333972116
   13       2.00000000014326451      -0.00000000000000000
root = 2.00000000014
(grader:  2 points)
( result for running
  root = newtonm(lambda x: 1 - np.exp(-(x-2)**2),  lambda x: -2*(x-2)*np.exp(-(x-2)**2),  5, 1.0e-4, 1)
    print ('root =', root, end='\n\n'))

def gasket(x, y, tol):

# function to draw the Sierpinski Gasket

# x is a vector of the left base, top, and right base x coordinates of the biggest triangle

# y is a vector of the left base, top, and right base y coordinates of the biggest triangle

# tolerance is a scalar tolerance.  if the base plus the height is smaller than this, quit

    fill(x,y,'r')                     # color in the big triangle in red

    show(block=False)                 # show the figure

    remove_triangle(x,y, tol)         # make the call to our recursive removal function

# This is a function to recursively remove smaller and smaller upside-down triangles

# if triangle too small just quit

if max(np.diff(x)+np.diff(y)) < tol:

    return

 # calculate the points of the upside-down triangle

    xr = # your code goes here

    yr = # your code goes here

    fill(xr,yr,'w')   # remove our new triangle by coloring it white

    draw()            # update the figure to show newly removed triangle

    time.sleep(0.1)   # added to slow things down so we can see what's going on...

    #make three recursive calls to remove each of our new (smaller) red triangles

    #  3 statements need to be written here, each calls remove_triangle