Getting started with Scientific Programming


This article was written by me some time back and was first published in Linux For You. Sharing it here with community members

Scientific programming, or in broader terms, scientific computing, deals with solving scientific problems with the help of computers, so as to obtain results more quickly and accurately. Computers have long been used for solving complex scientific problems — however, advancements in computer science and hardware technologies over the years have also allowed students and academicians to play around with robust scientific computation tools.

Although tools like Mathematica and Matlab remain commercial, the open source community has also developed some equally powerful computational tools, which can be easily used by students and independent researchers. In fact, these tools are so robust that they are now also used at educational institutions and research labs across the globe.

So, let’s move on to setting up a scientific environment.

Setting up the environment

Most UNIX system/Linux distributions have Python installed by default. We will use Python 2.6.6 for the purposes of this article. It’s recommended to install IPython, as it offers enhanced introspection, additional shell syntax, syntax highlighting and tab-completion.

Next, we’ll install the two most basic scientific computational packages for Python: NumPy and SciPy. The former is the fundamental package needed for scientific computing with Python. It contains a powerful N-dimensional array object, sophisticated functions, tools for integrating C/C++, and Fortran code with useful linear algebra, Fourier transforms, and random-number capabilities. The SciPy library is built to work with NumPy arrays, and provides many user-friendly and efficient numerical routines for numerical integration and optimisation.

Open the Synaptic Package Manager and install the python-numpy and python-scipy packages. Now that we have NumPy and SciPy installed, let’s get our hands dirty with some mathematical functions and equations!


Numerical computations with NumPy, SciPy and Maxima

NumPy offers efficient array computations with fixed-size, homogeneous, multi-dimensional array types, and a plethora of functions to perform various array operations. Array-programming languages like NumPy generalise operations in scalars to apply transparently to vectors, matrices and other higher-dimensional arrays. Python does not have a default array data type, and processing data with Python lists and for loops is dramatically slower compared to corresponding operations in compiled languages like FORTRAN, C and C++. NumPy comes to the rescue, with its dynamically typed environment for array computation, similar to basic Matlab. You can create a simple array with the array function in NumPy:

In[1]: import numpy as np

In[2]: a = np.array([1,2,3,4,5])

In[2]: b = np.array([6,7,8,9,10])

In[3]: type(b) #check datatype

Out[3]: type numpy.ndarray #array

In[4]: a+b

Out[4]: array([7,9,11,13,15])

You can also convert a simple array to a matrix array using the shape attribute.

In[1]: import numpy as np

In[5]: c = np.array([1,4,5,7,2,6])

In[6]: c.shape = (2,3)

In[7]: c

Out[7]: array([1,4,5],[7,2,6]) // converted to a 2 column matrix

Matrix operations

Now let us take a look at some simple matrix operations. The following matrix can be simply defined as:


Defining a matrix and matrix multiplication

In[1]: import numpy as np

In[2]: x = np.array([[1,2,3],[4,5,6],[7,8,9]])

In[3]: y = np.array([[1,4,5],[2,6,5],[6,8,3]]) #another matrix

In[4]: z =,y) #matrix multiplication using dot attribute

In[5]: z

Out[5]: z = ([[23, 40,24], [50,94,63],[77,148,102]])

You can also create matrices in NumPy using the matrix class. However, it’s preferable to use arrays, since most NumPy functions return arrays, and not matrices. Moreover, matrix objects have a maximum of Rank-2. To hold Rank-3 data, you need an array. Also, arrays are closer in semantics to tensor algebra, compared to matrix objects.

The following example shows how to transpose a matrix and define a diagonal matrix:

In[7]: import numpy as np

In[8]: x = np.array([[1,2,3],[4,5,6],[7,8,9]])

In[7]: xT = np.transpose(x) #take transpose of the matrix

In[8]: xT

Out[8]:xT = ([[1,4,7],[2,5,8],[3,6,9]])

In[9]:n = diag(range(1,4)) #defining a diagnol matrix

In[10]: n

Out[10]:n = ([[1,0,0],[0,2,0],[0,0,3]])


The scipy.integrate package provides several integration techniques, which can be used to solve simple and complex integrations. The package provides various methods to integrate functions. We will be discussing a few of them here. Let us first understand how to integrate the following functions:


Simple Integration of x^2

In[1]: from scipy.integrate import quad

In[2]: import scipy as sp

In[3]: sp.integrate.quad(lambda x: x**2,0,3)

Out[3]: (9.0, 9.9922072216264089e-14)

Integration of 2^sqrt(x)/sqrt(x)

In[4]: sp.integrate.quad(lambda x: 2**sqrt(x)/sqrt(x),1,3)

Out[4]: (3.8144772785946079, 4.2349205016052412e-14)

The first argument to quad is a “callable” Python object, i.e., a function, method, or class instance. We have used a lambda function as the argument in this case. A lambda function is one that takes any number of arguments — including optional arguments — and returns the value of a single expression. The next two arguments are the limits of integration. The return value is a tuple, with the first element holding the estimated value of the integral, and the second element holding an upper bound on the error.


We can get a derivative at a point via automatic differentiation, supported by FuncDesigner and OpenOpt, which are scientific packages based on SciPy. Note that automatic differentiation is different from symbolic and numerical differentiation. In symbolic differentiation, the function is differentiated as an expression, and is then evaluated at a point. Numerical differentiation makes use of the method of finite differences.

However, automatic differentiation is the decomposition of differentials provided by the chain rule. A complete understanding of automatic differentiation is beyond the scope of this article, so I’d recommend that interested readers refer to Wikipedia. Automatic differentiation works by decomposing the vector function into elementary sequences, which are then differentiated by a simple table lookup.

Unfortunately, a deeper understanding of automatic differentiation is required to make full use of the scientific packages provided in Python. Hence, in this article, we’ll focus on symbolic differentiation, which is easier to understand and implement. We’ll be using a powerful computer algebra system known as Maxima for symbolic differentiation.

Maxima is a version of the MIT-developed MACSYMA system, modified to run under CLISP. Written in Lisp, it allows differentiation, integration, solutions for linear or polynomial equations, factoring of polynomials, expansion of functions in the Laurent or Taylor series, computation of the Poisson series, matrix and tensor manipulations, and two- and three-dimensional graphics.

Open the Synaptic Package Manager and install the maxima package. Once installed, you can run it by executing the maxima command in the terminal. We’ll be differentiating the following simple functions with the help of Maxima:

d / dx(x4) d / dx(sin x + tan x) d / dx(1 / log x)


You have to simply define the function in diff() and maxima will calculate the derivative for you.

(%i1) diff(x^4)

(%o1) 4x^3 del(x)

(%i2) diff(sin(x) + tan(x))

(%o2) (sec^2(x) + cos(x))del(x)

(%i3) diff(1/log(x))

(%o3) - del(x)/x log^2(x)-

The command diff(expr,var,num) will differentiate the expression in Slot 1 with respect to the variable entered in Slot 2 a number of times, determined by a positive integer in Slot 3. Unless a dependency has been established, all parameters and variables in the expression are treated as constants when taking the derivative. Similarly, you can also calculate higher order differentials with Maxima.

Ordinary differential equations

Maxima can also be used to solve ODEs. We’ll dive straight into some examples to understand how to solve ODEs with Maxima. Consider the following differential equations:

dx/dt = e-t + x d2x / dt2 – 4x = 0



Let’s rewrite our example ordinary differential equations using the noun form diff, which uses a single quote. Then use ode2, and call the general solution gsoln.

The function ode2 solves an ordinary differential equation of the first or second order. This takes three arguments: an ODE given by eqn, the dependent variable dvar, and the independent variable ivar. When successful, it returns either an explicit or implicit solution for the dependent variable. %c is used to represent the integration constant in the case of first-order equations, and %k1 and %k2 the constants for second-order equations.

We can also find the solution at predefined points using ic1 and call this particular solution, psoln. Consider the following non-linear first order differential equation:

(x2y)dy / dx = xy +x3 – 1

Let’s first define the equation, and then solve it with ode2. Further, let us find the particular solution at points x=1 and y=1 using ic1.

We can also solve ODEs with NumPy and SciPy using the FuncDesigner and OpenOpt packages. However, both these packages make use of automatic differentiation to solve ODEs. Hence, Maxima was chosen over these packages. ODEs can also be solved using the scipy.integrate.odeint package. We will later use this package for mathematical modelling.

Curve plotting with MatPlotLib

It’s said that a picture is worth a thousand words, and there’s no denying the fact that it’s much more convenient to make sense of a scientific experiment by looking at the plots as compared to looking just at the raw data.

In this article, we’ll be focusing on MatPlotLib, which is a Python package for 2D plotting that produces production-quality graphs. Matlab is customisable and extensible, and is integrated with LaTeX markup, which is really useful when writing scientific papers. Let us make a simple plot with the help of MatPlotLib:

Simple Plot with MatPlotLib


import matplotlib.pyplot as plt

x = range(10)

plt.plot(x, [xi**3 for xi in x])


You can create various types of plots using MatPlotLib. Let us take a look at Pie Plot and Scatter Plot.

import matplotlib.pyplot as plt


plt.title('Distribution of Dark Energy and Dark Matter in the Universe')

x = [74.0,22.0,3.6,0.4]

labels = ['Dark Energy', 'Dark Matter', 'Intergalatic gas', 'Stars,etc']

plt.pie(x, labels=labels, autopct='%1.1f%%');



import matplotlib.pyplot as plt

import numpy as np

plt.title('Scatter Plot')

x = np.random.randn(200)

y = np.random.randn(200)

plt.xlabel('X axis')

plt.ylabel('Y axis')


In this article, we have covered some of the most basic operations in scientific computing. However, we can also model and simulate more complex problems with NumPy and SciPy. These tools are now actively used for research in quantum physics, cosmology, astronomy, applied mathematics, finance and various other fields. With this basic understanding of scientific programming, you’re now ready to explore deeper realms of this exciting world!

The author of this article is Adit. Here is Adit's bio:

Founder & CEO at Function Space. I'm all about design, programming and astrophysics. In other words, you can call me a valence electron.

  • 0.00
    Amitabh Posted 2013-07-12 14:14:55

    Wow! a greatly compressed and info packed article. I'm trying to learn python myself and this is article helped a lot. You might want to move the images of Maxima solving ODE's to come after the text in proper sequence (kinda confused me for a while).

    P.S.: Here's the guide I'm using for learning python: Link

  • 3.90
    Goutham Posted 2013-12-21 12:45:45

    Hi,can you help me to install scikit-image on os x ?

Sign up for free to share your opinion. Already have an account? Sign in to post opinion