[

Lab 1: Numerical Integration

]

The goal of the lab is to develop for Fortran code to perform numerical integration. This code will be developed further next week as part of an assignment.

PROGRAMMING ASSIGNMENTS WILL ACCOUNT FOR 15% OF THE EXAMINATION MARK FOR THIS SUBJECT

O

Some things you'll need to know:

Trapezoid Rule:

$J = J_{T,n} + \mathcal{E}_{T,n}$ where

$$J_{T,n} = \frac{h}{2} \left(f_0 + 2 \sum_{i=1}^{n-1} f_i + f_n \right),$$

and

$$\mathcal{E}_{T,n} = - \frac{h^2}{12} (b-a) f''(\tau), \qquad \tau \in [a,b].$$

Simpson's Rule:

(With $n=2m$ subintervals of equal length ) $J = J_{S,n} + \mathcal{E}_{S,n}$ where

$$J_{S,n} = \frac{h}{3} \left(f_0 + 4 \sum_{i=1}^{n/2} f_{2i-1} + 2 \sum_{i=1}^{n/2-1} f_{2i} + f_n \right),$$

and

$$\mathcal{E}_{S,n} = - \frac{h^4}{180} (b-a) f^{(iv)}(\tau), \qquad \tau \in [a,b].$$

O

The following program implements the Trapezoid Rule to approximate

$$\int_a^b f(x) \mathrm{d}x= \int_a^b \frac{x}{1 + x^2} \mathrm{d}x.$$

The user inputs $a$, $b$, and $n$ (the number of subintervals used)

program numint
  implicit none  
  integer n
  real a, b, j
  
  print*, 'a, b, n, = ?'
  read*, a, b, n
      
  call trap(a, b, n, j)
  print*, 'Int f(x) dx =', j
  
end program numint
!     **************


real function f(x)
  real x
  f = x/(1.0 + x*x)
end function f
!     ***************


subroutine trap(a, b, n, j)
  implicit none
  real, intent(in) :: a,b
  integer, intent(in) :: n
  real, intent(out) :: j
  real h, x, f
  integer i
  
  h = (b - a)/n
  j = f(a)
  x = a
  loopi:   do i = 1,n-1
     x = x + h
     j = j + 2.0*f(x)
  end do loopi
  
  j = j + f(b)
  j = j * h/2.0
end

Q1

Compile and run the above program. Test it by computing $\int_0^2 \frac{x}{1 + x^2} \mathrm{d}x$ with $n=4$, $n=8$, and $n=16$. Ensure that the results are converging to the true solution, i.e.,

$$\frac{1}{2} \ln(5) =0.8047189560$$

Q2

For $n = 4, 8, 16$, compute approximations $J_{T,n}$ to

$$\int_0^2 x \mathrm{d}x, \quad \int_0^2 x^2 \mathrm{d}x, \quad \int_0^2 x^3 \mathrm{d}x, \quad \int_0^2 \exp(x) \mathrm{d}x.$$

Complete the following table:

	$n=4$	$n=8$	$n=16$
$\Bigl\vert \int_0^2 x \mathrm{d}x- J_{T,n} \Bigr\vert$
$\Bigl\vert \int_0^2 x^2 \mathrm{d}x- J_{T,n} \Bigr\vert$
$\Bigl\vert \int_0^2 x^3 \mathrm{d}x- J_{T,n} \Bigr\vert$
$\Bigl\vert \int_0^2 \exp(x) \mathrm{d}x- J_{T,n} \Bigr\vert$

Q3

For $\int_0^2 \exp(x) \mathrm{d}x$, what value of $n$ must you choose so that $\vert J - J_{T,n}\vert \leq 10^{-5}$?

O

Next Week: you will have to

1. Extend this program to also approximate integrals using Simpson's Rule.
2. Implement ``error estimation by halving $n$''.