Week 6 : Numpy arrays & exploring chaotic phenomena¶
Outcome: Students will learn more about numpy, extending last lesson’s knowledge to 2D arrays. Students will also be made acquainted with famous chaotic systems, and how code enables the exploration and study of these systems.
- What we will do:
Setup check
2D arrays in numpy
Exercise: Sum of all elements in a 2D array
Chaotic systems
Exercise: Lorenz system
Exercise: Mandelbrot set
Setup check¶
For this week’s class, you will need to be using a local Python installation, just like last week. In fact, you will need it for the rest of this course! Given that we will need to import libraries every week, online REPLs will not be able to support our work any further.
If you can run the following code chunk without errors, you are good to go for today!
import numpy as np
import matplotlib.pyplot as plt
plt.plot([1, 2, 3, 4], [1, 2, 3, 4])
plt.show()
2D arrays in numpy¶
Previously, we learned on how to initialize 1D arrays like below
arr_a = np.array([57, 82, 147, 71, 111, 115, 101, 89, 260, 287])
We live in a 3D world, so obviously not everything is in 1D! We should be able to have arrays that have more than one dimension. This is what we’ll be learning today.
2D arrays?¶
Before we get into that, we should understand what 2D arrays are and what kind of information we can put into them. Let’s say for example in GENIUSpintar’s senior batch there are 49 students and 7 classes in total. If we had a 1D array, we would only be able to keep all the student’s names and nothing else. However, if we had a 2D array, we can keep all types of information such as the student’s names, id number, class number, and so on.In short, imagine a 2D array like Microsoft Excel where you can store data in different columns and rows.
Initializing 2D arrays¶
This is how we initialize a 2D array using np.array. It is fairly similar to what we did with 1D arrays with a small difference. The following code sets up a 2D array containing a few rows of data, where the first column is person number, and the second column is age. Verify that by printing this array yourself:
arr_a = np.array([[1, 24], [2, 19], [3, 21]])
Last week, we also used np.zeros and np.ones to initialize arrays that are all zeros or all ones, in 1D. Below is how we initialize a 2D array with 10 rows and 4 columns:
np.zeros(shape=(10, 4))
np.ones(shape=(10, 4))
And the following is how we use np.zeros and np.ones to initialize a 3D array with 9 rows, 8 columns, and a third dimension with length 2.
np.zeros(shape=(9, 8, 2))
np.ones(shape=(9, 8, 2))
Do you see the format here? Basically, the shape of an array is specified by the length of each dimension in sequence. Typically, instead of the word “dimension”, we use the word “axis”, or in plural, “axes”.
Array indexing on 2D axes¶
We have seen last week that indexing a 1D array is similar to indexing for lists. For arrays with two dimensions, the syntax is similar, just that we need to place the indices in the corresponding axes. We will use the very first array initialized above as the example array to index.
- Following are a few examples:
accessing the first value in the second axis
accessing the third value in the first axis
accessing the last value in the last axis
accessing all values in the first axis
arr_a[0, 1]
arr_a[2, 0]
arr_a[-1, -1]
arr_a[:, -1]
- In simpler words:
Row 1 Column 2, thus indexed by [0, 1]
Row 3 Column 1, thus indexed by [2, 0]
Row -1 Column -1, thus indexed by [-1, -1]
All rows Column 1, thus indexed by [:, -1]
What’s going on in the last example?
Note
When indexing arrays from the first element to the last, you can specify :, as a shortcut to represent 0:-1.
Array functions in 2D¶
Last week, we learned how to use the functions sum, mean, max and argmax for 1D arrays. Just like array indexing, the syntax is similar, just that we need to be mindful of the axes.
This is how to find the sum of the values along each row, i.e. sum by columns:
arr_a.sum(axis=0)
This is how to calculate the sum of the values along each column, i.e. sum by rows:
arr_a.sum(axis=1)
Array conditionals in 2D¶
Array conditionals still work for 2D arrays as well. However, the format of the output is very different; numpy will first reshape the array to one dimension, then run the conditional. Try the following code, where we tell numpy to only return us elements in arr_a that are divisible by two. There are two elements that fulfill this condition, at arr_a[0, 1] and arr_a[1, 0] respectively.
print(arr_a[arr_a % 2 ==0])
# [24 2]
Notice that the output is in one dimension.
Array arithmetic in 2D¶
Element-wise operations in 2D work similarly as in 1D. Arithmetic between 2D arrays and numbers work as expected thanks to array broadcasting:
np.zeros(shape=(5, 3)) + 2
However, any arithmetic between arrays will require them to have the same shape. You can see this for yourself below:
# These lines of code will return errors!
np.ones(shape=(5, 3)) + np.ones(shape=(5, ))
np.ones(shape=(5, 3)) + np.ones(shape=(3, ))
# Only this one will not!
np.ones(shape=(5, 3)) + np.ones(shape=(5, 3))
With 2D arrays, we can also perform operations on them if they have the same shape.
Matrix operations¶
As 2D arrays are equivalent to matrices in mathematics, it is simple to perform matrix operations on 2D arrays. For instance, you can perform a dot product using np.dot(array1, array2):
arr_b = np.array([[ 0, 1, 2],
[ 3, 4, 5],
[ 6, 7, 8],
[ 9, 10, 11]])
arr_c = np.array([[ 0, 1, 2, 3],
[ 4, 5, 6, 7],
[ 8, 9, 10, 11]])
print(np.dot(arr_b, arr_c))
Exercise: Sum of all elements in a 2D array¶
The following code below creates an array and fills it with random numbers. Find the total sum of the numbers in the array.
big_array = np.random.random(size=(10000, 10000))
# YOUR CODE BELOW
Chaotic systems¶
Chaotic systems are systems that appear to be random but are actually deterministic, and are highly sensitive to initial conditions.
Think of it like this, remember Conway’s Game of Life and how it is important to have rules and its initial state defined? We then get seemingly random but deterministic patterns such as spaceships and oscillators just by specifying these initial conditions.
The butterfly effect in pop culture knowledge is a good example of this.
Note
- “A butterfly flaps its wings in the Amazonian jungle, and subsequently a storm ravages half of Europe”
-Terry Pratchett and Neil Gaiman”
The butterfly effect is defined as the sensitive dependence on initial conditions in which a small change in one state of a deterministic nonlinear system can result in large differences in a later state. Just by comparing the differences of the definition between Chaotic Systems/Chaos Theory and the Butterfly Effect, we can see that they are very similar. The concept “Butterfly Effect” can be imagined like this.
The founder of chaos theory Edward Lorenz states he used the Butterfly as an example as it had symbolic representation of unknowable quantity.
Note
“For want of a nail the shoe was lost,
For want of a shoe the horse was lost,
For want of a horse the rider was lost,
For want of a rider the battle was lost,
For want of a battle the kingdom was lost,
- And all for the want of a horseshoe nail.”
-Benjamin Franklin
Let’s go back to foundations of chaos theory by looking at Edward Lorenz did. Lorenz was a meteorologist and mathematician who combined both fields to become what we know today as chaos theory. In an experiment to model weather predictions, he decided to enter an initial condition of 0.506 instead of 0.506127. He thought that the weather prediction wouldn’t change by that much, but he got a completely different result. He deduced that a miniscule change in initial conditions had long-term implications. in 1963, he wrote a paper entitled Deterministic Nonperiodic Flow which he theorized that “weather prediction models are inaccurate because knowing the precise starting conditions is impossible, and a tiny change can throw off the results.”
It shows science to be less accurate than we assume, as we have no means of making accurate predictions due to the exponential growth of errors.
Another example of chaos theory is the 2x pendulum. The 2x pendulum is set up in a way that the 2nd pendulum is attached to edge of the 1st pendulum. This is so we can see how the 2nd pendulum moves about when we release at a specific angle. Now, we set 3 different double pendulum with the exact same specifications (same height, same weight). Then, we release the pendulums at the same angle. Just by doing this we can observe that the 3 pendulums do not follow each other. Simple miniscule difference in initial conditions change the result significantly.
Exercise: Lorenz System¶
Let’s go into more detail on the work of Edward Lorenz. He was working on a set of equations to model atmospheric convection in 1963, which describe how a two-dimensional layer of fluid would be warmed from below and cooled from above. Below is a simplified version of the Lorenz equations:
where x, y and z represent the dimensions, and a, b and c represent constants. This is a trio of ordinary differential equations, where the rate of change of x, y and z are dependent on each other over time. (Speaking of differential equations, we will cover it later in the course!)
Lorenz eventually found that a certain combination of constants resulted in chaotic behaviour, where given slightly different initial conditions, the equations produce wildly different output. We will take a look for ourselves in this exercise.
Part A¶
- Implement a function called lorenz that has the following arguments:
a, b and c, representing the constants of the Lorenz equations
init_state, a list of three numbers that represents initial conditions
- Inside the function, initialize two variables below:
dt to represent time step size. Set this to 0.01
duration to represent total length of time to simulate. Set this to 40.
Initialize the variable state as an array of all zeros, with shape (int(duration/dt), 3).
Set the first row of state to be equal to init_state
- Loop for i from 0 to len(state)-1:
Set x, y and z to be equal to state[i, 0], state[i, 1], state[i, 2]
Implement the Lorenz equations:
dx = (a*(y-x)) * dt dy = (x*(c-z) - y) * dt dz = (x*y - b*z) * dt state[i+1, 0] = x + dx state[i+1, 1] = y + dy state[i+1, 2] = z + dz
Return state after the loop.
Part B¶
What we have done in Part A is to use numerical methods to approximate the solutions of the Lorenz equations over time. We will also explore this particular topic later in the course!
1. Copy paste this function below.
def show_ts(state):
xs, ys, zs = state[:, 0], state[:, 1], state[:, 2]
plt.plot(range(len(xs)), xs)
plt.show()
plt.plot(range(len(xs)), ys)
plt.show()
plt.plot(range(len(xs)), zs)
plt.show()
Use it to explore the outputs of different runs of the lorenz function like so:
state = lorenz(a=10, b=8/3, c=5, init_state=[1, 1, 1])
show_ts(state)
Use constants from the example above, change the values in init_state a few times. What do you observe? Print the final values to compare how they differ from each other.
Repeat (2) using c equals 12, then 15, then 28. What do you observe?
Part C¶
- Modify show_ts to create the function show_2d, that:
plots xs vs ys, then xs vs zs, then ys vs zs
labels the x and y labels appropriately
Re-run Steps 2 and 3 from Part B. What do you observe?
Exercise: Mandelbrot Set¶
Take a complex number c. We square it, then add itself to the result. When we iterate over this process, we will get certain numbers of c where the outcome will either rise to infinity or drop to negative infinity. We will also get certain numbers of c where the absolute value will stay bounded, i.e. stay finite over multiple iterations. The latter numbers form the definition of the Mandelbrot set.
Note
The reason complex numbers are considered, is because the aforementioned relation is trivial for real numbers. On one hand, at a certain negative number and below, the square of itself will always exceed adding itself, causing the value to tend to infinity over iterations. On the other hand, at a certain positive number and above, the same happens.
Complex numbers however follow a more interesting squaring rule. A complex number is represented as \(a+bi\), where a and b are constants, while i represents the unit imaginary number. The square of a complex number would be:
From the above equation, you can see that many combinations of a and b can lead to different outcomes!
The Mandelbrot set is famous for demonstrating fractal behaviour, i.e. self-similarity, at different scales. As you zoom in, you will see more and more detail. In this module, we will look at how certain regions of the Mandelbrot set exhibits chaotic behaviour, and also visualize the Mandelbrot set while we’re at it.
Part A¶
Implement a function, mandelbrot, to determine if a complex number is part of the Mandelbrot set.
The function takes two arguments, c for the complex number, and iterations for number of times to repeat the equation.
Inside the function, define a variable z that is equal to c. Create a variable zs that is a list, and store z in it.
- Run a for-loop to loop over the range of numbers from 1 to iterations. In the loop:
Assign z to be equal to the square of itself, plus c.
Append z to zs
Return zs
Modify the function definition from def mandelbrot(c, iterations): to def mandelbrot(c, iterations=10):. This assigns a default value of 10 if we don’t specify iterations.
Note
Assign a value to the argument in a function definition to assign a default value. The proper name for this is optional arguments. Example:
def add_nums(num1, num2=10):
print(num1 + num2)
print_number(5, 5) # Returns 10
print_number(5) # Returns 15
6. Copy-paste the following code chunk. Modify the complex number passed to the function. What do you observe?
zs = mandelbrot(complex(1.21, 1.5))
print(zs[-1])
plt.plot(zs)
plt.show()
Note
Complex numbers in Python are initialized by the complex function. Example:
a = complex(1, 2)
print(a)
# (1+2j)
j is used to represent i, as i is already frequently used in for-loops!
You can also initialize complex numbers by specify it directly as a + bj, where a and b are numbers. Example:
print(1 + 5j)
# (1+5j)
Part B¶
Given that we’re in a coding class, we can do better than checking different complex numbers by hand. Why not use code to check a range of numbers for us? We can then visualize them on a plot, using the x-axis for real numbers, and y-axis for imaginary numbers.
- Write the function is_bounded to check if a given complex number is bounded, i.e. is part of the Mandelbrot set.
The function takes two arguments, real and imaginary. Instead of passing a complex number directly, we pass both parts separately because we will be running this function on a 2D grid of numbers! We will be using real numbers along the column (i.e. y-axis) to represent the imaginary components.
The function will also take another two arguments, iterations and thresh that have default values 50 and 10.
Inside the function, create the variable c as the complex number, using both arguments real and imaginary.
Initialize the variable z to have the value of c
- Create a for-loop over i, for i in the range of 0 to iterations:
Do z = z ** 2 + c
Check if abs(z) is large than thresh. If yes, return 0.
Outside the for-loop, add code to return 1.
2. Copy-paste the function definition below to run is_bounded on a sequence of xs and ys:
def calc_mandelbrot(xs, ys):
xs_count = len(xs)
ys_count = len(ys)
# Make array to store results
# First dim is how many rows, second dim is how many cols
res = np.zeros(shape=(ys_count, xs_count))
# For each col
for i in range(ys_count):
# For each row
for j in range(xs_count):
# Assign col, row results of checking (real, imaginary) which is (row, col)
res[i, j] = is_bounded(xs[j], ys[i])
return res
3. Generate xs and ys, inputs to calc_mandelbrot using np.arange:
xs = np.arange(-3, 2, 0.1)
ys = np.arange(-2, 2, 0.1)
Note
np.arange(start, end, step) allows specifying a 1D array where values from start to end are created, spaced by step.
4. Use the following code-chunk to generate a low resolution plot of the Mandelbrot set!
res = calc_mandelbrot(xs, ys)
plt.figure(figsize=(18, 10))
plt.imshow(res, origin="lower", extent=[xs[0], xs[-1], ys[0], ys[-1]])
plt.show()
Note
plt.imshow can be used to visualize 2D arrays as an image!
Enhance the resolution by reducing 0.1 in np.arange to 5e-3. Note that this might take a while to run depending on your computer. You will see a plot of the Mandelbrot set, where all numbers in the set are visualized in the blob.
Using the plot, pick complex numbers to verify if they are part of the set or not, by repeating (6) in Part A.
Part C¶
Setting a single threshold might not capture the full variation in the Mandelbrot set, as the points that are unbounded increase/decrease at different speeds!
Modify is_bounded to instead return the number of iterations required to exceed the threshold.
Re-run Part B Q4. You should see a lot more detail in the plot.
We want to study further the boundary between numbers in the Mandelbrot set and numbers that are not. We will create a function to simplify plotting the Mandelbrot set. Instead of specifying the bounds of the x-axis and y-axis, we will instead specify a midpoint, and the bounds of the image away from it, specified by imgsize.
Complete the function below:
def show_mandelbrot(x, y, imgsize, precision, iterations=50, thresh=10):
# YOUR CODE BELOW
# Use x, y and imgsize to find the following
xstart =
xend =
ystart =
yend =
xs = np.arange(xstart, xend, precision)
ys = np.arange(ystart, yend, precision)
res = calc_mandelbrot(xs, ys, iterations=iterations, thresh=thresh)
plt.figure(figsize=(18, 10))
plt.imshow(res, extent=[xs[0], xs[-1], ys[0], ys[-1]], origin="lower", cmap="magma")
plt.axis("image")
plt.colorbar()
plt.show()
4. We will zoom in to see the boundaries of the Mandelbrot set better. Run the following code line by line. Was it simple to determine the boundaries?
show_mandelbrot(-0.5, 0, 2, precision=5e-3)
show_mandelbrot(-0.75, 0.2, 0.5, precision=1e-3)
show_mandelbrot(-0.76, 0.1, 0.05, precision=5e-4)
show_mandelbrot(-0.77, 0.12, 0.03, precision=3e-4)
show_mandelbrot(-0.773, 0.125, 5e-3, precision=5e-5)
show_mandelbrot(-0.77315, 0.127, 1e-3, precision=1e-5)
show_mandelbrot(-0.77295, 0.127, 3e-4, precision=1e-6)
show_mandelbrot(-0.772975, 0.12704, 1e-5, precision=1e-7)
show_mandelbrot(-0.772972, 0.12704, 3e-6, precision=1e-8)
show_mandelbrot(-0.772972, 0.1270404, 5e-7, precision=5e-9)
show_mandelbrot(-0.7729719, 0.12704049, 3e-8, precision=1e-10)
show_mandelbrot(-0.772971895, 0.127040495, 3e-9, precision=5e-12)
Increase iterations to 200 and re-run. What has changed?
iterations = 200
show_mandelbrot(-0.75, 0.2, 0.5, precision=1e-3, iterations=200)
show_mandelbrot(-0.76, 0.1, 0.05, precision=3e-4, iterations=200)
show_mandelbrot(-0.77, 0.12, 0.03, precision=1e-4, iterations=200)
show_mandelbrot(-0.773, 0.125, 5e-3, precision=3e-5, iterations=200)
show_mandelbrot(-0.77315, 0.127, 1e-3, precision=5e-6, iterations=200)
show_mandelbrot(-0.77295, 0.127, 3e-4, precision=1e-6, iterations=200)
show_mandelbrot(-0.772975, 0.12704, 1e-5, precision=5e-8, iterations=200)
show_mandelbrot(-0.772972, 0.12704, 3e-6, precision=1e-8, iterations=200)
show_mandelbrot(-0.772972, 0.1270404, 5e-7, precision=5e-9, iterations=200)
show_mandelbrot(-0.7729719, 0.12704049, 3e-8, precision=1e-10, iterations=200)
show_mandelbrot(-0.772971895, 0.127040495, 3e-9, precision=5e-12, iterations=200)
Conclusion¶
- Take-away message for this week:
We learnt about working with arrays in two dimensions
We learnt about the concept of chaotic behaviour and inspected a few classic chaotic systems