Lecture 7

	Compute integrals Approximate solution by polynomial Store solution on each element Line / surface / volume integrals Map from reference element and each elem. Change of variables…

Consider an integral

What method should you use to compute this integral?

Recall that for a grid \(x_i = X_L + ih, \ \ i = 0,\ldots,n, \ \ h = \frac{X_R-X_L}{n}\), the composite trapezoidal rule is:

For fun: write a Fortran program that uses the composite trapezoidal rule to approximate any integral for \(n = 2,3,\ldots,N\).

There is ONE situation where trapezoidal rule is the go-to method. Which situation am I talking about?

The trapezoidal rule belongs to a class of quadrature called Newton-Cotes quadrature that approximates integrals using equidistant grids. The order of a composite Newton-Cotes method is typically \(s\) or math:s+1, where \(s\) is the number of points in each panel (2 for the Trapezoidal rule).

Another class of methods is Gauss quadrature. In Gauss quadrature the location of the grid-points (usually referred to as nodes) and weights, \(\omega_i\), are chosen so that the order of the approximation to the weighted integral

is maximized. Here the weight function \(w(z)\) is positive and integrable (for example \(w(z) = 1\).) For a given \(w(z)\) the nodes, \(z_i\) are the zeros of the polynomial \(\tau_n = (z-z_0)(z-z_1)\cdots(z-z_n)\) satisfying

for all polynomials \(q(z)\) of degree less than \(n\).

Nodes, \(z_i\) are the zeros of the polynomial \(\tau_n = (z-z_0)(z-z_1)\cdots(z-z_n)\) satisfying

for all polynomials \(q(z)\) of degree less than \(n\). Or, equivalently, the nodes are the zeros to the degree \(n\) orthogonal polynomial associated with the weight function \(w(z)\). Don’t worry if this sounds complicated, we will only consider the case \(w(z) = 1\) and obtain the weights by computation. For example you can use the repository subroutine lglnodes.f90 (which is a f90 version of Greg von Winckel’s matlab version.)

call lglnodes(x,w,n)
f = exp(cos(pi*pi*x))
Integral_value = sum(f*w)

$ gfortran -c type_defs.f90
$ gfortran -c quat_module.f90
$ gfortran -c apa.f90
$ gfortran -o apa.x apa.o quad_module.o type_defs.o

Once the number of files become moderately large this method does not work and we will use make instead.

1. differentiate_v1.f90
2. Makefile1, $ make  -f Makefile1

1
2
3
4
5

diff.x: differentiate_v1.o
        gfortran differentiate_v1.o -o diff.x

differentiate_v1.o: differentiate_v1.f90
        gfortran -c differentiate_v1.f90

General structure is:

Target: Dependencie(s)
<TAB> Rule 1
<TAB> Rule 2

For example diff.x depends on differentiate_v1.o which in turn depends on differentiate_v1.f90.

1
2
3
4
5

diff.x: differentiate_v1.o
        gfortran differentiate_v1.o -o diff.x

differentiate_v1.o: differentiate_v1.f90
        gfortran -c differentiate_v1.f90

Dependencies are checked by date-stamp.

If the .f90 is newer than .o file make recreates the .o file.

1
2
3
4
5
6
7
8

diff.x: differentiate_v3.o weights.o lglnodes.o
        gfortran differentiate_v3.o weights.o lglnodes.o -o diff.x
differentiate_v3.o: differentiate_v3.f90
        gfortran -c differentiate_v3.f90
weights.o: weights.f
        gfortran -c weights.f
lglnodes.o: lglnodes.f90
        gfortran -c lglnodes.f90

Notes:

1. Now diff.x depends on three files.
2. The files are compiled in the same way
3. We can avoid this by using macros and rules.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11

FC = gfortran
FFLAGS = -O3
F90FLAGS = -O3
run_it: diff.x
        ./diff.x
diff.x: differentiate_v3.o weights.o lglnodes.o
        gfortran differentiate_v3.o weights.o lglnodes.o -o diff.x
%.o : %.f90
        $(FC) $(F90FLAGS) -c $<
%.o : %.f
        $(FC) $(FFLAGS) -c $<

Notes:

1. A single rule for all .f90 and .f files.
2. FC is a macro or makefile variable, use by $(FC).