Solutions for homework assignment #2

Gabe's solutions in a Jupyter notebook [hw2_sol.ipynb]

Problem1

The purpose of this assignment was to investigate the scaling of the error as a function of the number of intervals used. Below is a link to an F90 program that integrates the function using both integration formulas in the same run. It uses a trick of using a previously calculated integral when doubling the number of intervals (no need to do that but saves about 50% computation time)..

F90 Solution program [openint.f90]

A) Non-singular integrand

These are results along with lines with slope 2 and 3. The expected error scaling for small h is confirmed. Note that there are small deviations from the linear behavior for the largest h, showing the presence of higher-order corrections to the scaling.

B) Singular integrand

The exponent can be extracted by plottng the logarithm of the error versus log(h) [or log(N) or n].

Using the line fitting program from the previous assignment (with all errors set to 1 since they are irrelevant but should be the same for every point to give them equal weight in the fit) gave the following exponents (slopes): 
   alpha=1/2, 1st order: slope = 0.49983
              2nd order: slope = 0.50003
   alpha=3/4, 1st order: slope = 0.24996
              2nd order: slope = 0.25001
This is in very good agreement with the expectation; slope=1-alpha.

C) Almost singular integrand

The figure below shows results using alpha=1/2 and epsilon=0.00001.

The lines show the expected slopes 1/2 (for h much larger than epsilon) 2 (for h much less than epsilon with 1st-order formula) and 3 (for h much less than epsilon with second-order formula).

Problem 2

F90 solution program [inertia.f90]

Running the program with 1 million sampling points per bin and nbi bins gave the following results and error bars for Iz and Iy 
  nbi =   50:  Iz = 0.00469302 +/- 0.00000081,  Ix = 0.00494854 +/- 0.00000091
  nbi =  500:  Iz = 0.00469148 +/- 0.00000027,  Ix = 0.00494760 +/- 0.00000031    
  nbi = 5000:  Iz = 0.00469192 +/- 0.00000009,  Ix = 0.00494771 +/- 0.00000010
Note that when increasing the number of bins by a factor of 10, the error bar is reduced approximately by sqrt(10), as expected statistically. The nbi=5000 run gave bin averages with the following distributions (using 50 histogram bins within the range where data appeared):
Gaussians have been fitted to the data sets, demonstrating normal distributions as expected by the central limit theorem when the number of points per bin is large.