Table of Contents

• 1  Exercises for Numpy, Scipy and Matplotlib
  • 1.1  Notebook examples
  • 1.2  Data-analysis-1
  • 1.3  Data-analysis-2
  • 1.4  Mandelbrot Set in Python
  • 1.5  Solving differential equation
  • 1.6  Data fitting - half time of radioactive source
  • 1.7  Polynom fit -- which polynom degree is needed?
  • 1.8  Extra: Mathematical model of an infectious disease

Exercises for Numpy, Scipy and Matplotlib

Notebook examples

• Check in detail the examples presented
• Test variations (data creation, plot options, ...)

Software Carpentry Tutorial

• Go through the 1-day Python introduction tutorial, which uses numpy and matplotlib for data analysis
• A useful complement and rehearsal of this course

Data-analysis-1

The file rohr1.dat contains a list of measurements of wire-positions of drift tubes used in the ATLAS Muon-Chamber System Read the numbers with:

data = np.loadtxt('rohr1.dat')
# or
data = np.loadtxt('http://www-static.etp.physik.uni-muenchen.de/kurs/Computing/python/nb/data/rohr1.dat')

1. Determine mean and standard-deviation (Hint: numpy-Functions)

2. Fill the values in a histogram and plot it.

In a similar way read the (x,y) coordinates of file rohr2.dat using

x,y = numpy.loadtxt('rohr2.dat',unpack=True)

Determine for both x and y mean and standard-deviation as well as the correlation.

Visualize the data:

1. 1D histogram of both x and y
2. (x,y) point plot
3. 2D histogram

Data-analysis-2

The file faithful.csv contains measurement data of the Old Faithful Geysir , i.e. duration of the eruption and time since the last eruption.

Download the file ( wget http://...) and read the data with numpy:

data=np.loadtxt('faithful.csv',delimiter=',',skiprows=1)
# or directly
data=np.loadtxt('https://people.sc.fsu.edu/~jburkardt/data/csv/faithful.csv',delimiter=',',skiprows=1)

(Why the options in in loadtxt(...)?)

1. Determine again mean and standard-deviation for duration and wait-time
2. Fill histograms for both values and plot it. Is it compatible with Gaussian distribution?
3. Are there correlations between wait-time and duration or the duration of sub-sequent eruptions? (Make (x,y) plots for both)

Mandelbrot Set in Python

The notebook Mandelbrotmenge.ipynb contains a showcase how to develop visualisation of the Mandelbrot set.

Solving differential equation

(a) Extend the example for free fall and introduce a term for air-friction:

$$F_R = \frac{1}{2} c_W \rho_{Luft} A v^2 $$

This results in a non-linear DE which can only be solved numerically:

$$ y'' = - g + \frac{1}{2} c_W \rho_{Luft} A / m_{fb} \times y'^{ 2}$$

Take numbers for a football:

    # Luftwiderstandskraft Fußball
    rhol = 1.184 # Luftdichte (kg/m^3)
    rfb  =  0.11 # Fußballradius (m)
    cwfb =  0.45 # Kugel CW-Wert
    mfb  =  0.43 # Masse Fußball (kg)
    g = 9.81 # Erdbeschleunigung
    v = y[1] # Geschwindigkeit
    fR = 0.5 * cwfb * rhol * math.pi * rfb**2 * v**2

(b) The next step would be the extension to a trajectory of a thrown object with air resistance ("Schiefer Wurf mit Luftreibung"), the combined movement in horizontal and vertical direction.

How far will a football fly which is shot under 45 degree angle?

(The aerodynamic of a football is rather complex, see an interesting discussion at https://www.weltderphysik.de/thema/hinter-den-dingen/coca-cola-formel/)

Data fitting - half time of radioactive source

In an experiment you measure the activity of a radioactive source with a counter. The counter takes readings every 30 seconds and shows the integrated number of counts during that period. The list of numbers below is the result of such a measurement. What is the lifetime of the radioactive source? Hint: Fit an exponential

counts = [728, 612, 581, 460, 381, 324, 305, 211, 203, 166, 139, 110, 120, 92, 85, 77, 48, 53, 47, 38, 30, 20, 22, 17, 17, 16, 18, 11, 9]

Polynom fit -- which polynom degree is needed?

Below is notebook with test-data and polynom fit, check which degree of polynom is sensible

%matplotlib inline
import numpy as np
import scipy as sp
import matplotlib.pyplot as plt

# test data
x = [ -3.5, -2.5, -1.5, -0.5, 0.5, 1.5,  2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, 10.5, 11.5]
y = [ 219., 155., 89., 11., 54., 80., 48., 5., -47., -64., -101., -125., -124., 13., 239., 525. ]
ey = [ 65., 41, 30., 26., 26., 24., 18., 4., 17., 33., 45., 49., 38., 6., 74., 160.]
 
ey = np.array(ey)
w=1./ey

n=1 #polynom degree

pars=sp.polyfit(x,y,n,w=w) # do the fit

yr=sp.polyval(pars,x) # compute residuals based on fitted line
chi2=sp.sum((yr-y)**2/ey**2) #compute Chi2

print (pars)
print ('Chi2, n, chi2/ndf:', chi2, n, chi2/(len(x)-n))

plt.title('Polynom Regression Example')
plt.errorbar(x,y,ey,fmt='ob')
#plt.plot(x,yr,'r--')

xr = np.linspace(x[0],x[-1],200)
yr2 = sp.polyval(pars,xr) 
plt.plot(xr,yr2,'r--')

plt.legend(['regression','values']);
plt.grid(True);

Extra: Mathematical model of an infectious disease

Consider the SIR model characterized through the following coupled systems of ordinary differential equations:

• ${\frac {dS}{dt}}=-c{\frac {S}{N}}I$
• ${\frac {dI}{dt}}=c{\frac {S}{N}}I-wI$
• ${\frac {dR}{dt}}=wI$

with the boundary condition that the total number of individuals N is the sum of S, I, and R and constant

• $N=S+I+R$

(Obviously, the boundary condition is trivially fulfilled as the sum of the derivatives is zero, i.e. no births or deaths are included in this version of the model.)

Herein are:

• S: number of susceptible individuals
• I: number of infected individuals
• R: number of recovered individuals
• c: probability of disease transmission
• w: rate of recovery
• t: time (independent variable)

Your tasks:

• Solve this system of equations with suitable methods from those presented above and plot S, I and R as a function of time.
• Change the parameters c and w and study their impact.
• Make the model more dynamic by introducing a birth and a death rate.

The SIR model is one of the simplest of a group of compartmental models in theoretical biology.

Proposed solution: SIRmodel.py and SIRmodel_dyn.py