Jupyter notebook coursenotes/Topic8/Topic8.ipynb
Topic 8: The Python
scientific modules - numpy
and scipy
scipy
Scipy
is another Python
module that offers a large number of pre-written routines for data manipulation including integration, solving differential equations, Fourier analysis and transforms, statistics, image processing, etc.
There is a lot of information about scipy at www.scipy.org
In PHY235 we will limit our use of scipy
to differentiation, integration and Fourier Transforms.
Integration
General integration is done with the quad
function, part of the scipy.integrate
sub-package. To access this function it is necessary to use the following form of the import
command
from scipy.integrate import quad
Calls to the quad
function are of the general form
result = quad(func, limlo, limho)
where func here is the function to be integrated, it should be a callable Python
object (i.e. a function, method or class). For the purposes of this course we will only use a function here. quad
returns a tuple containing 2 values, the first is the value of the integration, the second is an upper bound on the error. Here is an example of a short script to perform the integration
Double, triple and higher level integration can be done with the dblquad
, tplquad
and nqua
d fucntions from the scipy.integrate
package.
Fourier transforms
scipy
offers a way of performing Fourier transforms on data using routines from the fftpack
package. The following example performs a fast Fourier transform (FFT) on the function f where
this is a very similar problem to a question in the PHY250 tutorial questions on Fourier Techniques
** Notes **
in this example the
pyplot
plotting package has been used to visualise the data. The entry
%matplotlib inline
causes the plots to be drawn inline in the output cell of the Jupyter notebook. (Otherwise a separate window is created).
Solving ordinary differential equations
The scipy.integrate
package provides a function called odeint which can be used to solve differential equations.
odeint has many options. In its simplest form it can be executed via a function call such as
y = odeint(func, y0, t)
where
y is the output (solution) from odeint
func is a
Python
version of the mathematical function to be solvedy0 is the value of y(0)
t is a
numpy
array of values for t
The following is a simple example of solving a first order differential equation
** Notes **
This function must have 2 input variables even if the output is only dependent on one of those
Solving second order differential equations
It is possible to use scipy
to solved second order differential equations. In this case the first step is always to transform any nth-order ODE into a system of n first order ODEs. In addition we also need n initial conditions, one for each variable represented in the first order ODEs. This is best illustrated with an example.
For a damped harmonic oscillator of mass m then
where b is the damping coefficient and k is the spring constant
This can be rewritten as
To solve using odeint this needs to be written as two first order ODEs as follows
Taking m=l, k=1.2, b=0.3 and the initial conditions as x=1 and dx/dt=v=0 respectively it is then possible to set up a Python
script to solve these two equations simultaneously, as follows
The same technique can be used to solve 2 simultaneous first order ODEs.
Interpolation
The following example illustrates the interp1d
function (note that this is interp-1-d not interp-l-d).
A sparsely sampled series of data points (just 13 points over the range -2.5$\pi\pi$) of the cosine function is interpolated using the interp1d
function using both Linear and Cubic interpolation. The resulting functions, InterpLin and InterpCub can be used to predict a value of the function for any value of x. For the cubic interpolation the agreement between the function and the actual cosine function is excellent.
It is of interest to re-run the Python
script with varying values of the number of samples (13 in this example) to observe how the accuracy of the interpolation varies. This example also illustrates use of the subplot
command.
Numerical integration
In the case where we have an array of y values without knowing the explicit functional dependence of y and x, it is still possible to perform a numerical integration using scipy's
simps function (within the integrate
package). scipy.integrate.simps
uses Simpson's rule to perform the numerical integration. In the following the data from the Interpolation example are compared in 3 ways:
Integrate the function from the linear interpolation
Integrate the cosine function
Numerically integrate the array of points
A high level of agreement between the three methods is observed.