Homework 2 (Due 6 pm, Friday, 10/20/2017)ΒΆ

Please submit your homework to your git repo by 6 pm, Friday, 10/20/2017.

In this homework, you don’t need to submit Problem 1 but only Problem 2 and Problem 3 by the due date. Although Problem 1 is not required for submission, you’re strongly encouraged to do it for your own educational interests. The Fortran part of Problem 1 will not be graded, but note that there is a bonus problem for creating a shell script.

  1. You’re going to solve all 16 exercise problems in Part 1 of the Prof. Garaud’s notes on Fortran. For Problem 1, you don’t need to check in your solutions to your git repo, but you are still welcome to check them in (recommended). If this is the case, please do the following:

    • Make a directory called hw2 and make 16 subdirectories (e.g., exercise01, ..., exercise16) under hw2.
      • NOTE: Too many subdirectories to generate by hand? You may want to come up with a smart way of generating N many of them (directories, files, etc.) by writing a short shell script on your own. This is definately worth a try, and there will be some bonus points if you implement a short shell script to do this task.
    • Include the corresponding Fortran routines in each exercise** directory.
    • In each problem please write a short note on your answers, findings and any issues you encounter.
    • In all problems, you are going to use double precisions for all real variables. Please do this by using the proper gfortran compilers as we discussed in the class. Indicate what flag options you used.
  2. Write a Fortran 90 program(s) to approximate \pi using the sum of the series

    \pi = \sum_{n=0}^\infty 16^{-n} \Big( \frac{4}{8n+1}-\frac{2}{8n+4}-\frac{1}{8n+5}-\frac{1}{8n+6}\Big).

    Since you can only add up to a finite number of terms in a practical run, say, summing up to the N th term, you are going to need a stopping criterion on N to obtain an accurate enough numerical approximation. For instance, you can stop the summation if the absolute value of the difference between your approximation \pi_{appx} and the reference (or true) solution of \pi_{true} defined as \pi_{true} = acos(-1.d0) becomes smaller than a threshold value, e.g.:

    threshold = 1.e-8
    diff = abs(pi_appx - pi_true)
    if diff > threshold
      continue summing up

    Please report numbers of terms N and the difference diff for four different values of threshold, 1.e-4, 1.e-8, 1.e-12, and 1.e-16. Please make sure you use double precisions for all real variables in your implementation using a proper Fortran flag(s).

    Note that you also need to provide a clear instruction on how to compile your code(s), including compilation command(s), compiler flag(s), etc. so that your code result can be reproducible by TA and the intructor.

  3. Write a Fortran 90 algorithm for approximating integral quantities using the trapezoidal rule. For those who are not familiar with the trapezoidal approximation, study the numerical algorithm in the above Wikipedia link. In your code, the main driver routine will call the trapezoidal approximation implemented in two different forms and an exact calculation:

    • function implementation,
    • subroutine implementation,
    • exact implementation,

    all of which are in a Fortran module.

    For example, a pseudo-code may look like:

    function f(x)
      (implement the specific form of your funciton f(x) = x^2 + 1)
    end function
    function blablah...
      (you may want another function implementation for
      exact integration...)
    ! module trapezoidApprox includes two implementations,
    ! one in function and the other one in subroutine
    module trapezoidApprox
      ! -------------------------!
      ! [1] FUNCTION             !
      ! -------------------------!
      function trapezoidFunc(...)
              (your function implementation here)
      end function trapezoidFunc
      ! -------------------------!
      ! [2] SUBROUTINE           !
      ! -------------------------!
      subroutine trapezoidSub(...)
              (your subroutine implementation here)
      end subroutine trapezoidSub
      ! -------------------------!
      ! [3] EXACT                !
      ! -------------------------!
      subroutine trapezoidExact(...)
              (your subroutine implementation here, assuming that you know
              the function to integrate and its exact integral value)
      end subroutine trapezoidExact
    end module trapezoidApprox
    ! -------------------------!
    ! Driver routine           !
    ! -------------------------!
    program compute_integration
        (this is your driver routine that calls
          [1] trapezoidFunc,
          [2] trapezoidSub,
          [3] trapezoidExact.)
        print*, "[1] Trapezoidal in function   =", trapezoidFunc result
        print*, "[2] Trapezoidal in subroutine =", trapezoidSub result
        print*, "[3] Exact integration         =", trapezoidExact result
        print*, "[4] Error in function         =", trapezoidExact - trapezoidFunc result
        print*, "[5] Error in subroutine       =", trapezoidExact - trapezoidSub result
    end program compute_integration

    Test your code using a function

    f(x) = x^2 + 1

    to approximate its integration over [a,b]=[0,1]

    I =\int^b_a f(x) dx= \int^1_0 (x^2+1) dx

    Your input and output are:

    • INPUT: the endpoints a, b of the interval and the number N of subintervals

      a=x_0 \le x_1 \le x_2 \le \cdots \le x_{N-1} \le x_N=b

    • OUTPUT: approximate values from your function, subroutine, and exact implementations.

    Make sure you use proper sets of Fortran 90 flags for double precision handlings of floating point operations. It is strongly recommended to test your code using valuable debugging flags too.

    Test your code using N=25, 50, 100, 200 and 400, and observe how the errors change as a function of different values of N.