NanoPy #5

Contents

SciPy / scipy.integrate.odeint / scipy.optimize.fsolve / scipy.optimize.fmin*

SciPy

 SciPy (www.scipy.org) is a suite of Python modules that are broadly useful for scientific computing, and most are interfaces to compiled numerical and scientific codes, typically written in either Fortran or C. SciPy functions use numpy arrays to communicate between the Python layer and the lower-level compiled routines. SciPy provides support for numerical integration, function minimization, root finding, fitting, fft's, interpolation, plotting, etc. The SciPy tutorial is a good place to start to learn about what routines are available and how to use them.

Unlike most packages, scipy does a full import of its package hierarchy with a single import statement; i.e., import scipy imports the entire scipy package, although individual modules still need to be accessed in their appropriate namespaces. scipy also replicates the entire numpy namespace within its own namespace. 

import scipy
x = scipy.arange(0., 10., 0.1)            # scipy.arange is numpy.arange, and so on
y = scipy.integrate.odeint(dydt, y0, ts)  # use scipy.integrate.odeint function
z = scipy.optimize.fsolve(f, z0)          # use scipy.integrate.fsolve function
Individual package imports are also possible, however, using the usual syntax, e.g.: 
import scipy.arange, scipy.integrate.odeint, scipy.optimize.fsolve
# or
from scipy import arange
from scipy.integrate import odeint
from scipy.optimize import fsolve

scipy.integrate.odeint

The odeint function in the scipy.integrate package is a wrapper around lsoda, an ODE integrator developed at Lawrence Livermore National Lab, as part of the odepack package. lsoda/odeint automatically switches between stiff and non-stiff integration routines, depending on the characteristics of the solution, and does adaptive time-stepping to achieve a desired level of solution accuracy.

The full calling sequence for odeint is: 

odeint(func, y0, t, args=(), Dfun=None, col_deriv=0, 
   full_output=0, ml=None, mu=None, rtol=None, atol=None, 
   tcrit=None, h0=0.0, hmax=0.0, hmin=0.0, ixpr=0, 
   mxstep=0, mxhnil=0, mxordn=12, mxords=5, printmessg=0)
Run pydoc scipy.integrate.odeint from the command line, or help(scipy.integrate.odeint) from within the python interpreter after having imported scipy, for a description of all the function arguments.

The most important arguments to the odeint function are the first three, which must be supplied. 

y_trajectory = odeint(func, y0, t_array)
will integrate the ODE dy/dt = func(y, t), starting at y=y0, and returning the solution y at the timepoints t specified in the t_array. If the state y has multiple components (e.g., (theta, thetaDot) for the pendulum problem, or (theta, thetaDot, phi, phiDot) for the simple walker), then those components are assembled into an array object. The returned trajectory is an (NT X NC) array, where NT is the number of timepoints in the t_array, and NC is the number of components in the state vector y. Typically one wants to reset y after the integration has been done, and this is done by extracting the last element of the y_trajectory array: 
y = y_trajectory[-1]   # equivalent to y = y_trajectory[len(y_trajectory)-1]
If the state y has only one component, it need not be packed into an array (although odeint will convert it to a one-element array); the requested timepoints in t_array must be passed as an array, rather than a Python list, since odeint does not apparently do that conversion.

The ODE function func is assumed, at a minimum, to be a function of y and t, and must be written as such, even if the function does not explicitly depend on both variables. It is assumed that the first argument is y, and the second is t: 

def func(y, t):
         # return array of same shape as y describing the 
         # time rate of change of each element in y
If the function depends on other parameters, those are passed in through the optional args argument. args is to consist of a tuple of arguments, beyond the required y and t. For example, if func depends on two additional parameters alpha and beta, then the following is required: 
def func(y, t, alpha, beta):
         # return array of same shape as y describing the 
         # time rate of change of each element in y

# set values of y0, t_array, alpha and beta
y_trajectory = odeint(func, y0, t_array, args=(alpha, beta))

scipy.optimize.fsolve

The fsolve function in the scipy.optimize package is a wrapper around the hybrd and hybrj algorithms in the MINPACK library developed at Argonne National Lab. fsolve finds a solution x to the equation func(x)=0, for a scalar function func of multiple variables x. Its full calling sequence is: 
fsolve(func, x0, args=(), fprime=None, full_output=0, col_deriv=0, 
   xtol=1.49012e-08, maxfev=0, band=None, epsfcn=0.0, factor=100, diag=None)
The core arguments passed to fsolve are similar in spirit to those in odeint:

scipy.optimize.fmin, fmin_powell, fmin_cg, fmin_bfgs, etc.

These various functions seek to minimize a function func of multiple parameters x, using different underlying algorithms (Nelder-Mead, Powell's method, Conjugate Gradient, etc.). They all have calling conventions similar to what we've seen previously with odeint and fsolve. Consider fmin as a specific example: 
fmin(func, x0, args=(), xtol=0.0001, ftol=0.0001, maxiter=None, maxfun=None, 
   full_output=0, disp=1, retall=0)
The syntax of fmin is similar in spirit to that of odeint and fsolve: