Computational Physics with Python : Mathematical Operations in NumPy

Mathematical Operations in Numpy

Numpy has a lot of built-in flexibility when it comes to doing maths. The most obvious thing you might want to do it to treat numpy arrays as matrices or vectors, depending on the circumstances. The recommended way to do this is to use the functions that specify vector multiplication:

import numpy as np


A = np.linspace(1,9,9)
B = A.reshape((3,3))
print(B)
print(B*B) # here we do the element wise multiplication

print(np.dot(B,B)) # here we do the matrix multiplication

Out: [[ 30. 36. 42.]
[ 66. 81. 96.]
[102. 126. 150.]]

print(np.dot(A,A)) # vector dot product

Out: 285.0

print(np.dot(np.array((1,2,3)), B)) # left matrix multiplication

print(np.dot(B, np.array((1,2,3)))) # right multiplication

Out: [30. 36. 42.], [14. 32. 50.]

np.dot(C, np.linalg.inv(C))

Out: array([[ 1.00000000e+00, -3.33066907e-16, 6.66133815e-16],
[-5.32907052e-15, 1.00000000e+00, 1.33226763e-15],
[-1.77635684e-15, -1.11022302e-15, 1.00000000e+00]])

Random Numbers

Generating random numbers is an extremely useful tool in scientific computing. We will see some applications in a second, but let’s just get the syntax.

The first thing you might want is a random number drawn from a uniform probability between 0 and 1

np.random.rand()

Out: 0.9364363627444482

This function actually will give you a bunch of random numbers in any shape.

np.random.rand(10)

Out: array([0.66503376, 0.18556123, 0.66901048, 0.11084154, 0.8679818 ,
0.94987479, 0.18319166, 0.55419005, 0.88653254, 0.81304748])

np.random.rand(2,2,2)

Out: array([[[0.81066867, 0.51840196],
[0.58857379, 0.28347615]],

[[0.9803307 , 0.28988763],
[0.32737379, 0.8205811 ]]])

Given such a distribution, we can make our own function to cover any uniform distribution

def my_uniform(low, high, number):
  out = np.random.rand(number)
  out*=(high-low)
  out+=low
  return out

my_uniform(10,11,20)

my_uniform(-102.3,99.2,20)

Naturally numpy has their own version for the same function

np.random.uniform(-102.3,99.2,20)

It is important to realize that once you have a source of random numbers, you can mold it to do other things you want.

The other common kind of random number you want is one drawn from a normal distribution. This means that the probability density of drawing the number x is

The number mu is called the mean and sigma is the variance.

Python has a nice way of generating numbers with mu=0 and sigma = 1

np.random.randn(10)

Out: array([ 2.05664013, 0.97541087, 0.43781921, -0.27368105, 1.00171526,
0.16865747, 0.75375289, 1.40993338, 0.00240816, -0.75028019])

Let’s make a histogram to see what this is doing

import matplotlib.pyplot as plt
plt.hist(np.random.randn(100000), bins= np.linspace(-3,3,1000))
plt.show()

number_of_flips=10
number_of_tails=1000

results=np.zeros(number_of_tails)

for i in range(number_of_tails):
  results[i] = np.random.randint(2, size=number_of_flips).sum()/float(number_of_flips)

plt.hist(results,bins=np.linspace(0,1,10))
plt.show()

Demonstration: Calculating the value of π

To get a quick sense of why this is useful, we can demonstrate how we can calculate π using what we have just learned. This will actually serve as an important example of a much broader and more powerful set of techniques later in the course, but for now it is just to give you a taste of how these things can work for you.

The idea is the following: image you pick a point randomly inside a square of length 2 on each side. The area of this square is 4. A circle placed within the square with radius 1 has area π.

If a random draw numbers that land inside the square with a uniform probability, then π/4 of them (on average) should also be inside the circle. Said different, after picking a point in the square, we ask if it is also in the circle and keep track. At the end we take number in circle / total we have calculated π/4

def pi_calculator(N):
    x=np.random.uniform(-1,1,N)  # make a list of N random numbers of x-axis of box
    y=np.random.uniform(-1,1,N) # make a list of N random numbers of y-axis of box
    z=(x**2+y**2<1) # make a list of every time x^2 + y^2 < 1 (inside the cicle)
    return z.sum()/float(N)*4 # add all the points in the circle up and return 4*cicle/N


pi_calculator(10**7)

Out: 3.1417844

x=np.random.uniform(-1,1,5)
y=np.random.uniform(-1,1,5)
z=(x**2+y**2<1)

print(x,y)
print(x**2+y**2)
print(z)

To see, how this is working, let’s write a slower version to see the steps.

def pi_slow(N):
    circle=0
    for i in range(N):
        x=np.random.uniform(-1,1,1)  # pick a x coordinate in the box
        y=np.random.uniform(-1,1,1) # pick a y coordinate in the box
        if x**2+y**2<1: # make a list of every time x^2 + y^2 < 1 (inside the cicle)
            circle+=1 
    return 4*circle/N # add all the points in the circle up and return 4*cicle/N



pi_slow(10**6)

Out: 3.14154

The slow implementation makes it clear what we are doing, but it clearly takes much longer.

import time
t1 = time.time()
pi_slow(1000000)
print(time.time() - t1)

Out: 13.578989744186401

t1 = time.time()
pi_calculator(1000000)
print(time.time() - t1)

Out: 0.03968214988708496

Let's see how our error in the measurement of π scales with the number of random points we pick

short=int(10**5)
medium=int(10**6)
Long=int(10**(7))
trials=100
pi_list_short=np.zeros(trials)
pi_list_medium=np.zeros(trials)
pi_list_Long=np.zeros(trials)
for i in range(trials):
    pi_list_short[i]=pi_calculator(short)
    pi_list_medium[i]=pi_calculator(medium)
    pi_list_Long[i]=pi_calculator(Long)

fig, ax = plt.subplots()
p1=ax.hist(pi_list_short,bins=np.linspace(3.13,3.15,100),label='10^5')
p2=ax.hist(pi_list_medium,bins=np.linspace(3.13,3.15,100),label='10^6')
p3=ax.hist(pi_list_Long,bins=np.linspace(3.13,3.15,100),label='10^7')
ax.plot([np.pi,np.pi],[0,40])
plt.ylim(0,40)
plt.xlim(3.1375,3.145)
leg = ax.legend()
plt.show()

By eye, it looks like the blue is a approximately 10 times wider than the green. This would make sense if the error on the value of π decreased by 1/sqrt(N) where N is the number of random points used in the calculation. This is indeed what is going on and is a much more general fact about random numbers.

Summary

Numpy has a number of basic math operations we will make a lot of use off. Random numbers are a particularly valuable tool that is employed in all areas of science and engineering. E.g. simulating the behavior of any measurement involves adding random numbers to you signal. We can always use the output of the random number library to create the type of noise we want for a given application.