While NumPy arrays are useful in their own right for a wide range of computations, they also serve as the lingua franca for a number of other libraries that implement other sorts of algorithms. Foremost among these is SciPy, an open-source library for scientific computing in Python. It depends on NumPy, but its scope is much broader. SciPy comprises modules for linear algebra, optimization, quadrature, interpolation, Fourier analysis, and signal and image processing, as well as ODE solvers, special mathematical functions, and other common tools needed in science and engineering. Many of the SciPy functions are Python wrappers around widely-used, time-tested routines developed over the years within the numerical analysis community, and written in compiled languages such as Fortran and C. Rather than crafting such numerical routines yourself, you can take advantage of both the expertise of the community and a convenient calling interface within the SciPy wrapper.

SciPy simplifies the process of calling these compiled routines, as well as allowing users to pass Python functions to compiled library routines. Those Python functions might represent the instantaneous time derivatives of a set of ODEs being integrated, an objective function being optimized, or a function whose roots are being sought. NumPy arrays are often used to represent data produced by these routines, such as the dynamical trajectory returned by an ODE integrator or the set of optimal parameters found by an optimization function. With this convenience of operation comes some impact on performance, as described below.

Consider the integration of ordinary differential equations (ODEs), supported by the function solve_ivp in the scipy.integrate module. (More precisely, solve_ivp is used to numerically solve initial value problems, where an ODE is integrated forward in time from a specific initial condition.) SciPy enables users to integrate arbitrary ODEs of the form dy/dt = f(y,t,parameters) for a vector-valued state variable y, by writing the right-hand-side function f(y,t,parameters), which returns instantaneous time derivatives of y. In the simplest case, the function f(y,t,parameters) is written in Python and passed as an argument to solve_ivp. (And since functions are first-class objects in Python, they can be passed directly to other functions as arguments, without the need to wrap up a function in a handle as one might need to do in other languages.) Consider the following code that uses SciPy to integrate the Lorenz equations, a famous model that demonstrates low-dimensional chaos:

The function lorenz is written in Python, but is passed to solve_ivp, which is a wrapper around a variety of compiled routines that implement different integration methods. Every time the right-hand-side function lorenz needs to be evaluated as part of the integration, control passes from the compiled library up to the Python interpreter for function evaluation, and then back down to the integrator. This situation is depicted schematically in the figure below, in what is essentially the inverse of what was depicted previously in the discussion of hybrid codes. This process of repeatedly popping up to the interpreter does slow down the time integration as compared to that seen in a fully compiled routine, and for longer integrations, the SciPy version can be perhaps 10-30x slower than the equivalent compiled version.

This type of performance hit might not be a bottleneck for your particular application, so do not immediately assume that you need to find some alternative to calling solve_ivp or some other SciPy routine with a Python function as an argument. If, however, you decide that you need to speed up such a routine, there are ways to circumvent the back-and-forth between the interpreter and the compiled library. For example, you could write the lorenz function in C, and pass a pointer to that function to the odeint integrator. Describing such a process is beyond the scope of this topic (and is not necessarily for the faint of heart), but there are some tools that support this sort of approach, although any particular application is likely to differ in some of the details. An online tutorial describes the use of code generation and compilation techniques in SymPy and Cython to generate compiled code for use with odeint. In addition, scipy has introduced the notion of a LowLevelCallable to provide an mechanism for calling compiled functions from a handful of scipy routines (although not as yet from solve_ivp).