Task solution - numerical integration in fortran

This repository contains the realization of the third task of the fortran programming language course at AGH in Cracov.

Directory structure

Makefile
README.md
run.sh
measure_times.sh
src/
    functions.f90
    quadratures.f90
    main.f90
res/
    1image_results
    5images_results
    times

Makefile

This file contains recipes for building the fortran program and using it to generate files in res/.

run.sh

This scripts is called from Makefile to compute numerical integrals of some example functions using different quadratures and numbers of integration subintervals. The script uses ./integral's output to create res/1image_results and res/5images_results.

measure_times.sh

Another script, used in Makefile to measure computation speed with different numbers of images and write the results to res/times.

src/

This directory contains the fortran sources of the project.

res/

This directory contains result data files with running times of the program and computed integrals with their errors.

Building

As a prerequisite, make, gfortran and opencoarrays are required (I already have my grade, so I didn't do ifort 🙃).
Run

$ make

to build the program executables integrator and integrator_single (compiled with -fcoarray=single instead of lib).

Running

The executables are integrator and integrator_single. The first expected argument is 'analytical' for analytically computed integral or quadrature type ('newton-cotes' or 'gauss') for numerical integral. Second argument is function to be integrated - 'exp', 'sin' or 'poly' (which is some 10 degree polynomial I made up). Third and fourth arguments are start and end of integration interval. Fifth arg is polynomial order (explained later in this README) and and sixth is the number of subintervals. Arguments 5 and 6 are only needed for numerical integrals.
To run a multi-image version:

 $ cafrun -np IMAGES_NUMBER ./integrator ARGUMENTS

For a single image execution:

 $ ./integrator_single ARGUMENTS

Examples:

 $ ./integrator_single analytical exp 0 1 
 1.7182818284590451
 $ cafrun -np 3 ./integrator gauss exp 0 1 2 1000000
 1.71828182845905E+00  2.22044604925031E-16

The small second number of the output is the absolute difference between numerically and analytically computed integrals (numerical error).

Data in res/ can be regenerated by issuing

$ make results

from the main project directory.

Interpretation

This task was very vaguely described. How was I supposed to interpret the "polynomial order" argument in the context of a quadrature? Since Newton-Cotes quadratures come from integration of an interpolating polynomial, I assumed this would be the polynomial, the order of which is passed. As a result, polynomial order of 0 corresponds to rectangles method, order 1 - trapeze method and order 2 - Simpson method. I know this idea lacks consistency (rectangles method is an open quadrature while other 2 are closed), but these 3 were the quadratures required by the task. Applying gaussian quadratures is also equivalent to integrating an interpolating polynomial - just with different distribution of interpolation points. So 1-point quadrature corresponds to interpolating with degree 0 polynomial, 2-points quadrature - degree 1 polynomial, etc. This road of thinking IS extremely unintuitive and misleading, yet, this was the only way I could make use of this extra "polynomial order" parameter. Another idea was to ignore it as at least one other person did... I'm pretty sure what I did here is NOT what author of the task meant. But with such a broken task he provided, nothing better could be done. At least I used function pointers, which I believe we were supposed to do, so no one can say, that I made things easier for myself here...

Another ugly thing is the way I pass the number of subintervals... The required interface of integrating functions was specified in the task without a parameter for the number of subintervals. Perhaps it should be there istead of this stupid "polynomial order"?

In a real project I would obviously do things a different way.

Results analyzing

Running time

I did everything on 2-core CPU. I am surprised there's such a huge time penalty for bigger number of images. Does creation and synchronization of 8 images take over 2.5 seconds, as res/times suggests!? When doing stuff in C with pthreads I easily had this number of threads create and die below 0.01 second! "He did something wrong here", one could say about me. Well, when we did coarrays on a lab class, code compiled under ifort would also take SECONDS to synchronize ;_;

Quadratures accuracy

When summing a lot of numbers, the order is important. If we sum 1 000 000 numbers one by one, we are likely to get higher error related to computer arithmetic, than if we sum numbers by 100 000 and add ten sums together. This suggests, that for high number of subintervals, program ran with more images should give better results (because each image first sums integrals over its corresponding subintervals and first image then adds up all the sums), right? Not in this case, because I employed a smarter adding technique. The results for different numbers of images still vary a little, but the error is on the same level.

I also see some weird things, like 2-point gaussian quadratures being in many cases more accurate, than 3-point ones... I do not exclude the possibility of my own bugs here, but it's also worth noting, that gaussian quadrature using Legendre polynomials is designed for integration of polynomials and the only polynomial among our test functions is a 10-degree one, while the 3-point gaussian quadrature will only exactly integrate polynomials of up to 5 degree.

One could write more about the results here, but I'm already bored enough and already got my grade anyway 😎