Commits

Randy LeVeque  committed 0f819c0

Added homework solutions

  • Participants
  • Parent commits 9e85b76

Comments (0)

Files changed (19)

File solutions/homework2/hw2a.py

+"""
+Demonstration script for quadratic interpolation.
+Modified for the problem described in Homework 2 problem #1.
+"""
+
+import numpy as np
+import matplotlib.pyplot as plt
+from numpy.linalg import solve
+
+# Set up linear system to interpolate through data points:
+
+# Data points:
+xi = np.array([-1., 1., 2])
+yi = np.array([0., 4., 3.])
+
+# Second row of A changed from demo1.py:
+A = np.array([[1., -1., 1.], [1., 1., 1.], [1., 2., 4.]])
+b = yi
+
+# Solve the system:
+c = solve(A,b)
+
+print "The polynomial coefficients are:"
+print c
+
+# Plot the resulting polynomial:
+x = np.linspace(-2,3,1001)   # points to evaluate polynomial
+y = c[0] + c[1]*x + c[2]*x**2
+
+plt.figure(1)       # open plot figure window
+plt.clf()           # clear figure
+plt.plot(x,y,'b-')  # connect points with a blue line
+
+# Add data points  (polynomial should go through these points!)
+plt.plot(xi,yi,'ro')   # plot as red circles
+plt.ylim(-2,8)         # set limits in y for plot
+
+plt.title("Data points and interpolating polynomial")
+
+plt.savefig('hw2a.png')   # save figure as .png file

File solutions/homework2/hw2b.py

+
+"""
+Demonstration module for quadratic interpolation.
+
+Sample solutions for Homework 2 problems #2 through #7.
+"""
+
+
+import numpy as np
+import matplotlib.pyplot as plt
+from numpy.linalg import solve
+
+def quad_interp(xi,yi):
+    """
+    Quadratic interpolation.  Compute the coefficients of the polynomial
+    interpolating the points (xi[i],yi[i]) for i = 0,1,2.
+    Returns c, an array containing the coefficients of
+      p(x) = c[0] + c[1]*x + c[2]*x**2.
+
+    """
+
+    # check inputs and print error message if not valid:
+
+    error_message = "xi and yi should have type numpy.ndarray"
+    assert (type(xi) is np.ndarray) and (type(yi) is np.ndarray), error_message
+
+    error_message = "xi and yi should have length 3"
+    assert len(xi)==3 and len(yi)==3, error_message
+
+    # Set up linear system to interpolate through data points:
+
+    A = np.vstack([np.ones(3), xi, xi**2]).T
+    c = solve(A,yi)
+
+    return c
+
+def plot_quad(xi, yi):
+    """
+    Perform quadratic interpolation and plot the resulting function along
+    with the data points.
+    """
+
+    # Compute the coefficients:
+    c = quad_interp(xi,yi)
+
+    # Plot the resulting polynomial:
+    x = np.linspace(xi.min() - 1,  xi.max() + 1, 1000)
+    y = c[0] + c[1]*x + c[2]*x**2
+
+    plt.figure(1)       # open plot figure window
+    plt.clf()           # clear figure
+    plt.plot(x,y,'b-')  # connect points with a blue line
+
+    # Add data points  (polynomial should go through these points!)
+    plt.plot(xi,yi,'ro')   # plot as red circles
+    plt.ylim(-2,8)         # set limits in y for plot
+
+    plt.title("Data points and interpolating polynomial")
+
+    plt.savefig('quadratic.png')   # save figure as .png file
+
+
+def cubic_interp(xi,yi):
+    """
+    Cubic interpolation.  Compute the coefficients of the polynomial
+    interpolating the points (xi[i],yi[i]) for i = 0,1,2,3
+    Returns c, an array containing the coefficients of
+      p(x) = c[0] + c[1]*x + c[2]*x**2 + c[3]*x**3.
+
+    """
+
+    # check inputs and print error message if not valid:
+
+    error_message = "xi and yi should have type numpy.ndarray"
+    assert (type(xi) is np.ndarray) and (type(yi) is np.ndarray), error_message
+
+    error_message = "xi and yi should have length 4"
+    assert len(xi)==4 and len(yi)==4, error_message
+
+    # Set up linear system to interpolate through data points:
+
+    A = np.vstack([np.ones(4), xi, xi**2, xi**3]).T
+    c = solve(A,yi)
+
+    return c
+
+
+def plot_cubic(xi, yi):
+    """
+    Perform cubic interpolation and plot the resulting function along
+    with the data points.
+    """
+
+    # Compute the coefficients:
+    c = cubic_interp(xi,yi)
+
+    # Plot the resulting polynomial:
+    x = np.linspace(xi.min() - 1,  xi.max() + 1, 1000)
+    y = c[0] + c[1]*x + c[2]*x**2 + c[3]*x**3
+
+    plt.figure(1)       # open plot figure window
+    plt.clf()           # clear figure
+    plt.plot(x,y,'b-')  # connect points with a blue line
+
+    # Add data points  (polynomial should go through these points!)
+    plt.plot(xi,yi,'ro')   # plot as red circles
+    plt.ylim(-2,8)         # set limits in y for plot
+
+    plt.title("Data points and interpolating polynomial")
+
+    plt.savefig('cubic.png')   # save figure as .png file
+
+
+def poly_interp(xi,yi):
+    """
+    General polynomial interpolation. 
+
+    Compute the coefficients of the polynomial
+    interpolating the points (xi[i],yi[i]) for i = 0,1,2,...,n-1
+    where n = len(xi) = len(yi).
+
+    Returns c, an array containing the coefficients of
+      p(x) = c[0] + c[1]*x + c[2]*x**2 + ... + c[N-1]*x**(N-1).
+
+    """
+
+    # check inputs and print error message if not valid:
+
+    error_message = "xi and yi should have type numpy.ndarray"
+    assert (type(xi) is np.ndarray) and (type(yi) is np.ndarray), error_message
+
+    error_message = "xi and yi should have the same length "
+    assert len(xi)==len(yi), error_message
+
+    # Set up linear system to interpolate through data points:
+    # Uses a list comprehension, see
+    # http://docs.python.org/2/tutorial/datastructures.html#list-comprehensions
+
+    n = len(xi)
+    A = np.vstack([xi**j for j in range(n)]).T
+    c = solve(A,yi)
+
+    return c
+
+def plot_poly(xi, yi):
+    """
+    Perform polynomial interpolation and plot the resulting function along
+    with the data points.
+    """
+
+    # Compute the coefficients:
+    c = poly_interp(xi,yi)
+
+    # Plot the resulting polynomial:
+    x = np.linspace(xi.min() - 1,  xi.max() + 1, 1000)
+
+    # Use Horner's rule:
+    n = len(xi)
+    y = c[n-1]
+    for j in range(n-1, 0, -1):
+        y = y*x + c[j-1]
+
+    plt.figure(1)       # open plot figure window
+    plt.clf()           # clear figure
+    plt.plot(x,y,'b-')  # connect points with a blue line
+
+    # Add data points  (polynomial should go through these points!)
+    plt.plot(xi,yi,'ro')   # plot as red circles
+    plt.ylim(yi.min()-1, yi.max()+1)       # set limits in y for plot
+
+    plt.title("Data points and interpolating polynomial")
+
+    plt.savefig('poly.png')   # save figure as .png file
+
+
+def test_quad1():
+    """
+    Test code, no return value or exception if test runs properly.
+    """
+    xi = np.array([-1.,  0.,  2.])
+    yi = np.array([ 1., -1.,  7.])
+    c = quad_interp(xi,yi)
+    c_true = np.array([-1.,  0.,  2.])
+    print "c =      ", c
+    print "c_true = ", c_true
+    # test that all elements have small error:
+    assert np.allclose(c, c_true), \
+        "Incorrect result, c = %s, Expected: c = %s" % (c,c_true)
+
+    # Also produce plot:
+    plot_quad(xi,yi)
+
+
+def test_quad2():
+    """
+    Test code, no return value or exception if test runs properly.
+    """
+
+    # Generate a test by specifying c_true first:
+    c_true = np.array([7., 2., -3.])
+    # Points to interpolate:
+    xi = np.array([-1.,  0.,  2.])
+    # Function values to interpolate:
+    yi = c_true[0] + c_true[1]*xi + c_true[2]*xi**2
+
+    # Now interpolate and check we get c_true back again.
+    c = quad_interp(xi,yi)
+    print "c =      ", c
+    print "c_true = ", c_true
+    # test that all elements have small error:
+    assert np.allclose(c, c_true), \
+        "Incorrect result, c = %s, Expected: c = %s" % (c,c_true)
+
+    # Also produce plot:
+    plot_quad(xi,yi)
+
+
+def test_cubic1():
+    """
+    Test code, no return value or exception if test runs properly.
+    """
+
+    # Generate a test by specifying c_true first:
+    c_true = np.array([7., -2., -3., 1.])
+    # Points to interpolate:
+    xi = np.array([-1.,  0.,  1., 2.])
+    # Function values to interpolate:
+    yi = c_true[0] + c_true[1]*xi + c_true[2]*xi**2 + c_true[3]*xi**3
+
+    # Now interpolate and check we get c_true back again.
+    c = cubic_interp(xi,yi)
+    print "c =      ", c
+    print "c_true = ", c_true
+    # test that all elements have small error:
+    assert np.allclose(c, c_true), \
+        "Incorrect result, c = %s, Expected: c = %s" % (c,c_true)
+
+    # Also produce plot:
+    plot_cubic(xi,yi)
+
+
+def test_poly1():
+    """
+    Test code, no return value or exception if test runs properly.
+    Same points as test_cubic1.
+    """
+    # Generate a test by specifying c_true first:
+    c_true = np.array([7., -2., -3., 1.])
+    # Points to interpolate:
+    xi = np.array([-1.,  0.,  1., 2.])
+    # Function values to interpolate:
+    # Use Horner's rule:
+    n = len(xi)
+    yi = c_true[n-1]
+    for j in range(n-1, 0, -1):
+        yi = yi*xi + c_true[j-1]
+
+    # Now interpolate and check we get c_true back again.
+    c = poly_interp(xi,yi)
+
+    print "c =      ", c
+    print "c_true = ", c_true
+    # test that all elements have small error:
+    assert np.allclose(c, c_true), \
+        "Incorrect result, c = %s, Expected: c = %s" % (c,c_true)
+
+    # Also produce plot:
+    plot_poly(xi,yi)
+    
+
+def test_poly2():
+    """
+    Test code, no return value or exception if test runs properly.
+    Test with 5 points (quartic interpolating function).
+    """
+    # Generate a test by specifying c_true first:
+    c_true = np.array([0., -6., 11., -6., 1.])
+    # Points to interpolate:
+    xi = np.array([-1.,  0.,  1., 2., 4.])
+    # Function values to interpolate:
+    # Use Horner's rule:
+    n = len(xi)
+    yi = c_true[n-1]
+    for j in range(n-1, 0, -1):
+        yi = yi*xi + c_true[j-1]
+
+    # Now interpolate and check we get c_true back again.
+    c = poly_interp(xi,yi)
+
+    print "c =      ", c
+    print "c_true = ", c_true
+    # test that all elements have small error:
+    assert np.allclose(c, c_true), \
+        "Incorrect result, c = %s, Expected: c = %s" % (c,c_true)
+
+    # Also produce plot:
+    plot_poly(xi,yi)
+    
+
+if __name__=="__main__":
+    print "Running test..."
+    test_quad1()
+    test_quad2()
+    test_cubic1()
+    test_poly1()
+    test_poly2()

File solutions/homework3/Makefile

+
+OBJECTS = functions.o newton.o test1.o
+OBJECTS2 = functions.o newton.o intersections.o
+FFLAGS = -g
+
+.PHONY: test1 clean intersections
+
+test1: test1.exe
+	./test1.exe
+
+test1.exe: $(OBJECTS)
+	gfortran $(FFLAGS) $(OBJECTS) -o test1.exe
+
+intersections: intersections.exe
+	./intersections.exe
+
+intersections.exe: $(OBJECTS2)
+	gfortran $(FFLAGS) $(OBJECTS2) -o intersections.exe
+
+
+%.o : %.f90
+	gfortran $(FFLAGS) -c  $< 
+
+clean:
+	rm -f *.o *.exe *.mod
+

File solutions/homework3/functions.f90

+! Modified to add f_int and fprime_int for the intersections problem.
+
+module functions
+
+    implicit none
+    real(kind=8) :: epsilon
+    save
+
+contains
+
+real(kind=8) function f_sqrt(x)
+    implicit none
+    real(kind=8), intent(in) :: x
+
+    f_sqrt = x**2 - 4.d0
+
+end function f_sqrt
+
+
+real(kind=8) function fprime_sqrt(x)
+    implicit none
+    real(kind=8), intent(in) :: x
+    
+    fprime_sqrt = 2.d0 * x
+
+end function fprime_sqrt
+
+real(kind=8) function f_int(x)
+    implicit none
+    real(kind=8), intent(in) :: x
+    real(kind=8) :: pi
+    pi = acos(-1.d0)
+
+    f_int = x*cos(pi*x) - (1.d0 - 0.6d0*x**2)
+
+end function f_int
+
+
+real(kind=8) function fprime_int(x)
+    implicit none
+    real(kind=8), intent(in) :: x
+    real(kind=8) :: pi
+    pi = acos(-1.d0)
+    
+    fprime_int = (cos(pi*x) - x*pi*sin(pi*x)) + 1.2d0*x
+
+end function fprime_int
+
+end module functions

File solutions/homework3/intersections.f90

+! Main program
+
+program intersections
+
+    use newton, only: solve
+    use functions, only: f_int, fprime_int
+
+    implicit none
+    real(kind=8) :: x, x0, fx
+    real(kind=8) :: x0vals(4)
+    integer :: iters, itest
+	logical :: debug
+
+    print *, "Test routine for computing zero of f"
+    debug = .true.
+
+    ! values to test as x0:
+    x0vals = (/-2.2d0, -1.6d0, -0.8d0, 1.45d0/)
+
+    do itest=1,4
+        x0 = x0vals(itest)
+		print *, ' '  ! blank line
+        call solve(f_int, fprime_int, x0, x, iters, debug)
+
+        print 11, x, iters
+11      format('solver returns x = ', es22.15, ' after', i3, ' iterations')
+
+        fx = f_int(x)
+        print 12, fx
+12      format('the value of f(x) is ', es22.15)
+
+        enddo
+
+end program intersections

File solutions/homework3/intersections.png

Added
New image

File solutions/homework3/intersections.py

+"""
+Script to solve the intersections problem.
+"""
+
+from newton import solve
+from numpy import linspace, pi, cos, sin, abs, where
+import matplotlib.pyplot as plt
+
+def g1vals(x):
+    g1 =  x*cos(pi*x)
+    g1p =  cos(pi*x) - x*pi*sin(pi*x)
+    return g1, g1p
+
+def g2vals(x):
+    g2 = 1. - 0.6*x**2
+    g2p = -1.2*x
+    return g2, g2p
+
+
+def fvals(x):
+    g1,g1p = g1vals(x)
+    g2,g2p = g2vals(x)
+    f = g1 - g2
+    fp = g1p - g2p
+    return f, fp
+
+xx = linspace(-10,10,1000)
+g1xx, g1pxx = g1vals(xx)
+g2xx, g2pxx = g2vals(xx)
+
+plt.figure(1)
+plt.clf()
+plt.plot(xx, g1xx, 'b')
+plt.plot(xx, g2xx, 'r')
+plt.legend(['g1','g2'])
+
+plt.figure(2)
+plt.clf()
+fxx, fpxx = fvals(xx)
+plt.plot(xx, fxx, 'b')
+plt.plot(xx, 0*xx, 'k')  # x-axis
+plt.title("plot of f(x) = g2(x) - g1(x)")
+
+plt.figure(1)
+for x0 in [-2.2, -1.6, -0.8, 1.45]:
+
+    x,iters = solve(fvals, x0)
+    print "\nWith initial guess x0 = %22.15e," % x0
+    print "      solve returns x = %22.15e after %i iterations " % (x,iters)
+    g1x,g1px = g1vals(x)
+    fx,fpx = fvals(x)
+    print "    f(x) = %22.15e" % fx
+    plt.plot(x,g1x,'ko')
+
+plt.axis((-5,5,-4,4))
+plt.savefig('intersections.png')

File solutions/homework3/newton.f90

+! Modified to add module parameter tol
+
+module newton
+
+    ! module parameters:
+    implicit none
+    integer, parameter :: maxiter = 40
+    real(kind=8), parameter :: tol = 1.d-14
+    save
+
+contains
+
+subroutine solve(f, fp, x0, x, iters, debug)
+
+    ! Estimate the zero of f(x) using Newton's method. 
+    ! Input:
+    !   f:  the function to find a root of
+    !   fp: function returning the derivative f'
+    !   x0: the initial guess
+    !   debug: logical, prints iterations if debug=.true.
+    ! Returns:
+    !   the estimate x satisfying f(x)=0 (assumes Newton converged!) 
+    !   the number of iterations iters
+     
+    implicit none
+    real(kind=8), intent(in) :: x0
+    real(kind=8), external :: f, fp
+    logical, intent(in) :: debug
+    real(kind=8), intent(out) :: x
+    integer, intent(out) :: iters
+
+    ! Declare any local variables:
+    real(kind=8) :: deltax, fx, fxprime
+    integer :: k
+
+
+    ! initial guess
+    x = x0
+
+    if (debug) then
+        print 11, x
+ 11     format('Initial guess: x = ', es22.15)
+	    endif
+
+    ! Newton iteration to find a zero of f(x) 
+
+    do k=1,maxiter
+
+        ! evaluate function and its derivative:
+        fx = f(x)
+        fxprime = fp(x)
+
+        if (abs(fx) < tol) then
+            exit
+            endif
+
+        ! compute Newton increment x:
+        deltax = fx/fxprime
+
+        ! update x:
+        x = x - deltax
+
+        if (debug) then
+            print 12, k,x
+ 12         format('After', i3, ' iterations, x = ', es22.15)
+	        endif
+
+        enddo
+
+    if (k > maxiter) then
+        ! might not have converged
+
+        fx = f(x)
+        if (abs(fx) > tol) then
+            print *, '*** Warning: has not yet converged'
+            endif
+        endif 
+
+    ! number of iterations taken:
+    iters = k-1
+
+
+end subroutine solve
+
+end module newton

File solutions/homework3/newton.py

+
+import numpy as np
+
+def solve(fvals, x0, debug=False):
+    tol = 1e-14
+    maxiter = 20 
+    x = x0
+    if debug:
+        print "Initial guess: x = %22.15e" % x0
+    for k in range(maxiter):
+        f,fp = fvals(x)
+        if abs(f) < tol:
+            break
+        delta_x = f / fp
+        x = x - delta_x
+        if debug:
+            print "After %s iterations, x = %22.15e" % (k+1,x)
+
+    if k==maxiter-1:
+        f,fp = fvals(x)
+        if abs(f) > tol:
+            print "*** warning: has not yet converged"
+        return x, k+1
+    else:
+        return x, k
+
+def fvals_sqrt(x):
+    """
+    Return f(x) and f'(x) for applying Newton to find a square root.
+    """
+    f = x**2 - 4.
+    fp = 2.*x
+    return f, fp
+
+def test1(debug_solve=False):
+    """
+    Test Newton iteration for the square root with different initial
+    conditions.
+    """
+    from numpy import sqrt
+    for x0 in [1., 2., 100.]:
+        print " "  # blank line
+        x,iters = solve(fvals_sqrt, x0, debug=debug_solve)
+        print "solve returns x = %22.15e after %i iterations " % (x,iters)
+        fx,fpx = fvals_sqrt(x)
+        print "the value of f(x) is %22.15e" % fx
+        assert abs(x-2.) < 1e-14, "*** Unexpected result: x = %22.15e"  % x
+

File solutions/homework3/problem7/Makefile

+
+OBJECTS = functions.o newton.o test1.o
+OBJECTS3 = functions.o newton.o test_quartic.o
+FFLAGS = -g
+
+.PHONY: test1 clean test_quartic
+
+test1: test1.exe
+	./test1.exe
+
+test1.exe: $(OBJECTS)
+	gfortran $(FFLAGS) $(OBJECTS) -o test1.exe
+
+test_quartic: test_quartic.exe
+	./test_quartic.exe
+
+test_quartic.exe: $(OBJECTS3)
+	gfortran $(FFLAGS) $(OBJECTS3) -o test_quartic.exe
+
+%.o : %.f90
+	gfortran $(FFLAGS) -c  $< 
+
+clean:
+	rm -f *.o *.exe *.mod
+

File solutions/homework3/problem7/functions.f90

+! Modified to add f_quartic and fprime_quartic for quartic function,
+! and to add epsilon as a module variable.
+
+module functions
+
+    implicit none
+    real(kind=8) :: epsilon
+    save
+
+contains
+
+real(kind=8) function f_sqrt(x)
+    implicit none
+    real(kind=8), intent(in) :: x
+
+    f_sqrt = x**2 - 4.d0
+
+end function f_sqrt
+
+
+real(kind=8) function fprime_sqrt(x)
+    implicit none
+    real(kind=8), intent(in) :: x
+    
+    fprime_sqrt = 2.d0 * x
+
+end function fprime_sqrt
+
+real(kind=8) function f_quartic(x)
+    implicit none
+    real(kind=8), intent(in) :: x
+
+    f_quartic = (x - 1.d0)**4 - epsilon
+
+end function f_quartic
+
+
+real(kind=8) function fprime_quartic(x)
+    implicit none
+    real(kind=8), intent(in) :: x
+    
+    fprime_quartic = 4.d0 * (x - 1.d0)**3
+
+end function fprime_quartic
+
+end module functions

File solutions/homework3/problem7/newton.f90

+! Modified to make tol a module variable
+
+module newton
+
+    ! module parameters:
+    implicit none
+    integer, parameter :: maxiter = 40
+    real(kind=8) :: tol
+    save
+
+contains
+
+subroutine solve(f, fp, x0, x, iters, debug)
+
+    ! Estimate the zero of f(x) using Newton's method. 
+    ! Input:
+    !   f:  the function to find a root of
+    !   fp: function returning the derivative f'
+    !   x0: the initial guess
+    !   debug: logical, prints iterations if debug=.true.
+    ! Returns:
+    !   the estimate x satisfying f(x)=0 (assumes Newton converged!) 
+    !   the number of iterations iters
+     
+    implicit none
+    real(kind=8), intent(in) :: x0
+    real(kind=8), external :: f, fp
+    logical, intent(in) :: debug
+    real(kind=8), intent(out) :: x
+    integer, intent(out) :: iters
+
+    ! Declare any local variables:
+    real(kind=8) :: deltax, fx, fxprime
+    integer :: k
+
+
+    ! initial guess
+    x = x0
+
+    if (debug) then
+        print 11, x
+ 11     format('Initial guess: x = ', es22.15)
+	    endif
+
+    ! Newton iteration to find a zero of f(x) 
+
+    do k=1,maxiter
+
+        ! evaluate function and its derivative:
+        fx = f(x)
+        fxprime = fp(x)
+
+        if (abs(fx) < tol) then
+            exit
+            endif
+
+        ! compute Newton increment x:
+        deltax = fx/fxprime
+
+        ! update x:
+        x = x - deltax
+
+        if (debug) then
+            print 12, k,x
+ 12         format('After', i3, ' iterations, x = ', es22.15)
+	        endif
+
+        enddo
+
+    if (k > maxiter) then
+        ! might not have converged
+
+        fx = f(x)
+        if (abs(fx) > tol) then
+            print *, '*** Warning: has not yet converged'
+            endif
+        endif 
+
+    ! number of iterations taken:
+    iters = k-1
+
+
+end subroutine solve
+
+end module newton

File solutions/homework3/problem7/test_quartic.f90

+! Main program for quartic function
+
+program test_quartic
+
+    use newton, only: solve, tol
+    use functions, only: f_quartic, fprime_quartic, epsilon
+
+    implicit none
+    real(kind=8) :: x, x0, fx, xstar
+    real(kind=8) :: tolvals(3), epsvals(3)
+    integer :: iters, itest1, itest2
+	logical :: debug
+
+    x0 = 4.d0
+    print 10 ,x0
+10  format("Starting with initial guess ", es22.15,/)
+    debug = .false.
+
+    ! values to test as epsilon:
+    epsvals = (/1.d-4, 1.d-8, 1.d-12/)
+
+    ! values to test as tol:
+    tolvals = (/1.d-5, 1.d-10, 1.d-14/)
+
+    print *, '    epsilon        tol    iters          x                 f(x)        x-xstar'
+
+    do itest1=1,3
+        epsilon = epsvals(itest1)
+        xstar = 1.d0 + epsilon**(0.25d0)
+        do itest2=1,3
+
+        tol = tolvals(itest2)
+
+        call solve(f_quartic, fprime_quartic, x0, x, iters, debug)
+
+        fx = f_quartic(x)
+        print 11, epsilon, tol, iters, x, fx, x-xstar
+11      format(2es13.3, i4, es24.15, 2es13.3)
+
+        enddo
+		print *, ' '  ! blank line
+        enddo
+
+end program test_quartic

File solutions/homework3/test1.f90

+! $UWHPSC/codes/fortran/newton/test1.f90
+
+program test1
+
+    use newton, only: solve
+    use functions, only: f_sqrt, fprime_sqrt
+
+    implicit none
+    real(kind=8) :: x, x0, fx
+    real(kind=8) :: x0vals(3)
+    integer :: iters, itest
+	logical :: debug
+
+    print *, "Test routine for computing zero of f"
+    debug = .false.
+
+    ! values to test as x0:
+    x0vals = (/1.d0, 2.d0, 100.d0 /)
+
+    do itest=1,3
+        x0 = x0vals(itest)
+		print *, ' '  ! blank line
+        call solve(f_sqrt, fprime_sqrt, x0, x, iters, debug)
+
+        print 11, x, iters
+11      format('solver returns x = ', es22.15, ' after', i3, ' iterations')
+
+        fx = f_sqrt(x)
+        print 12, fx
+12      format('the value of f(x) is ', es22.15)
+
+        if (abs(x-2.d0) > 1d-14) then
+            print 13, x
+13          format('*** Unexpected result: x = ', es22.15)
+            endif
+        enddo
+
+end program test1

File solutions/homework4/quadrature.f90

+
+module quadrature
+
+contains
+
+real(kind=8) function trapezoid(f, a, b, n)
+
+    ! Estimate the integral of f(x) from a to b using the
+    ! Trapezoid Rule with n points.
+
+    ! Input:
+    !   f:  the function to integrate
+    !   a:  left endpoint
+    !   b:  right endpoint
+    !   n:  number of points to use
+    ! Returns:
+    !   the estimate of the integral
+     
+    implicit none
+    real(kind=8), intent(in) :: a,b
+    real(kind=8), external :: f
+    integer, intent(in) :: n
+
+    ! Local variables:
+    integer :: j
+    real(kind=8) :: h, trap_sum, xj
+
+    h = (b-a)/(n-1)
+    trap_sum = 0.5d0*(f(a) + f(b))  ! endpoint contributions
+    
+    do j=2,n-1
+        xj = a + (j-1)*h
+        trap_sum = trap_sum + f(xj)
+        enddo
+
+    trapezoid = h * trap_sum
+
+end function trapezoid
+
+subroutine error_table(f,a,b,nvals,int_true)
+    implicit none
+    real(kind=8), intent(in) :: a,b, int_true
+    real(kind=8), external :: f
+    integer, dimension(:), intent(in) :: nvals
+
+    ! Local variables:
+    integer :: j, n
+    real(kind=8) :: sum, ratio, last_error, error, int_trap
+
+    print *, "    n         trapezoid            error       ratio"
+    last_error = 0.d0 
+    do j=1,size(nvals)
+        n = nvals(j)
+        int_trap = trapezoid(f,a,b,n)
+        error = abs(int_trap - int_true)
+        ratio = last_error / error
+        last_error = error  ! for next n
+
+        print 11, n, int_trap, error, ratio
+ 11     format(i8, es22.14, es13.3, es13.3)
+        enddo
+
+end subroutine error_table
+
+
+end module quadrature
+

File solutions/homework4/quadrature_omp.f90

+
+module quadrature_omp
+
+    ! Sample solution for Homework 4.
+
+    use omp_lib
+
+contains
+
+real(kind=8) function trapezoid(f, a, b, n)
+
+    ! Estimate the integral of f(x) from a to b using the
+    ! Trapezoid Rule with n points.
+
+    ! Input:
+    !   f:  the function to integrate
+    !   a:  left endpoint
+    !   b:  right endpoint
+    !   n:  number of points to use
+    ! Returns:
+    !   the estimate of the integral
+     
+    implicit none
+    real(kind=8), intent(in) :: a,b
+    real(kind=8), external :: f
+    integer, intent(in) :: n
+
+    ! Local variables:
+    integer :: j
+    real(kind=8) :: h, trap_sum, xj
+
+    h = (b-a)/(n-1)
+    trap_sum = 0.5d0*(f(a) + f(b))  ! endpoint contributions
+    
+    !$omp parallel do private(xj) reduction(+ : trap_sum) 
+    do j=2,n-1
+        xj = a + (j-1)*h
+        trap_sum = trap_sum + f(xj)
+        enddo
+
+    trapezoid = h * trap_sum
+
+end function trapezoid
+
+subroutine error_table(f,a,b,nvals,int_true)
+    implicit none
+    real(kind=8), intent(in) :: a,b, int_true
+    real(kind=8), external :: f
+    integer, dimension(:), intent(in) :: nvals
+
+    ! Local variables:
+    integer :: j, n
+    real(kind=8) :: sum, ratio, last_error, error, int_trap
+
+    print *, "      n         trapezoid            error       ratio"
+    last_error = 0.d0 
+    do j=1,size(nvals)
+        n = nvals(j)
+        int_trap = trapezoid(f,a,b,n)
+        error = abs(int_trap - int_true)
+        ratio = last_error / error
+        last_error = error  ! for next n
+
+        print 11, n, int_trap, error, ratio
+ 11     format(i8, es22.14, es13.3, es13.3)
+        enddo
+
+end subroutine error_table
+
+
+end module quadrature_omp
+

File solutions/homework4/test1.f90

+
+program test1
+
+    use quadrature, only: trapezoid, error_table
+
+    implicit none
+    real(kind=8) :: a,b,int_true
+    integer :: nvals(7), i
+
+    a = 0.d0
+    b = 2.d0
+    int_true = (b-a) + (b**4 - a**4) / 4.d0
+
+    print 10, int_true
+ 10 format("true integral: ", es22.14)
+    print *, " "  ! blank line
+
+    ! values of n to test:
+    do i=1,7
+        nvals(i) = 5 * 2**(i-1)
+        enddo
+
+    call error_table(f, a, b, nvals, int_true)
+
+contains
+
+    real(kind=8) function f(x)
+        implicit none
+        real(kind=8), intent(in) :: x 
+        
+        f = 1.d0 + x**3
+    end function f
+
+end program test1

File solutions/homework4/test2.f90

+
+program test2
+
+    ! Sample solution for Homework 4.
+
+    use quadrature, only: trapezoid, error_table
+
+    implicit none
+    real(kind=8) :: a,b,int_true,k
+    integer :: nvals(12)
+    integer :: i
+
+    k = 1.d3
+    a = 0.d0
+    b = 2.d0
+    int_true = (b-a) + (b**4 - a**4) / 4.d0 - (1.d0/k) * (cos(k*b) - cos(k*a))
+
+    print 10, int_true
+ 10 format("true integral: ", es22.14)
+    print *, " "  ! blank line
+
+    ! values of n to test:
+    do i=1,12
+        nvals(i) = 5 * 2**(i-1)
+        enddo
+
+    call error_table(f, a, b, nvals, int_true)
+
+contains
+
+    real(kind=8) function f(x)
+        implicit none
+        real(kind=8), intent(in) :: x 
+        
+        f = 1.d0 + x**3 + sin(k*x)
+    end function f
+
+end program test2

File solutions/homework4/test2_omp.f90

+
+program test2_omp
+
+    ! Sample solution for Homework 4.
+
+    use omp_lib
+
+    use quadrature_omp, only: trapezoid, error_table
+
+    implicit none
+    real(kind=8) :: a,b,int_true,k
+    integer :: nvals(12)
+    integer :: i
+
+
+    !$ call omp_set_num_threads(2)
+
+
+    k = 1.d3
+    a = 0.d0
+    b = 2.d0
+    int_true = (b-a) + (b**4 - a**4) / 4.d0 - (1.d0/k) * (cos(k*b) - cos(k*a))
+
+    print 10, int_true
+ 10 format("true integral: ", es22.14)
+    print *, " "  ! blank line
+
+    ! values of n to test:
+    do i=1,12
+        nvals(i) = 5 * 2**(i-1)
+        enddo
+
+    call error_table(f, a, b, nvals, int_true)
+
+contains
+
+    real(kind=8) function f(x)
+        implicit none
+        real(kind=8), intent(in) :: x 
+        
+        f = 1.d0 + x**3 + sin(k*x)
+    end function f
+
+end program test2_omp