from IPython.core.display import HTML
css_file = './stylesheets/custom.css'
HTML(open(css_file, "r").read())
Current version is adapted from KAUST Python Tutorials 2015 by Dmitry Kabanov.
Presented on 26 February, 2017.
Prerequisites: Introduction to Python, basic knowledge of numerical methods
Python includes a package for numerical computation called numpy, which is an essential tool for scientific computing. It is based on multidimensional arrays (vectors and matrices) and provides a large library of functions and operators that work efficiently on these objects. In this way, any numerical algorithm is expressed as operations on arrays and in many cases runs almost as quickly as the equivalent C code.
Numpy
serves as a foundation library that makes Python a competitor to MATLAB for scientific computing.
NumPy
is used in:
Numerical Analysis, Linear algebra, Solution of differential equations, Statistical analysis and many more...
Consider the following list and suppose we want to perform some simple algebraic operations:
price = [5.99, 10.25, 2.0, 40.99, 5.60, 63.49]
price*2
[5.99, 10.25, 2.0, 40.99, 5.6, 63.49, 5.99, 10.25, 2.0, 40.99, 5.6, 63.49]
price+0.5
--------------------------------------------------------------------------- TypeError Traceback (most recent call last) <ipython-input-4-ba6ebb983a2c> in <module>() ----> 1 price+0.5 TypeError: can only concatenate list (not "float") to list
Clearly, if we wanted to double the value of each element and add half, this cannot be done easily. Instead, we would have to use a loop:
for i in range(len(price)):
price[i] = price[i]*2
price[i] = price[i]+0.5
print(price)
[12.48, 21.0, 4.5, 82.48, 11.7, 127.48]
However, it is slow (first) and requires too much code (second). It can be done more efficient using NumPy
.
We start, by importing the NumPy module. Import statements like this are the typical way of getting access to functions in Python. The most general way to import a module is the following:
import numpy as np
First example: compute $\sin (3\pi/2)$:
np.sin(3*np.pi/2)
-1.0
Second example: compute $e + \pi$:
np.exp(1) + np.pi
5.8598744820488378
In above examples we compute scalars.
What makes NumPy brilliant is that it allows you make computations on multidimensional quantities (such as vectors and matrices) using thing called ndarray
($N$-dimensional array).
Computer scientists call such a thing a data structure.
We will consider only one-dimensional (vectors) and two-dimensional (matrices) arrays in this tutorial.
Let's create a 1D array $p$ as we did before and compute $2p + 0.5$:
price = np.array([1, 2, 3, 4, 5])
type(price)
numpy.ndarray
price = price*2+0.5
price
array([ 2.5, 4.5, 6.5, 8.5, 10.5])
There are other ways to create one-dimensional arrays:
Function linspace
has three arguments: the value of the left point, the value of the right point and the number of points in total:
x = np.linspace(0, 1, num=5)
x
array([ 0. , 0.25, 0.5 , 0.75, 1. ])
Another function, arange
is very similar but it takes a step between points as the third argument:
y = np.arange(0, 1, 0.25)
y
array([ 0. , 0.25, 0.5 , 0.75])
Notice, that arange
does not include the right point in the result!
Other functions that can be used to create an array:
x = np.zeros(5)
y = np.ones(3)
print(x, '\n', y)
[ 0. 0. 0. 0. 0.] [ 1. 1. 1.]
Now let's compute $\sin 2\pi x$ and see its graph:
x = np.linspace(0, 1, num=21)
y = np.sin(2*np.pi*x)
import matplotlib.pyplot as plt
%matplotlib inline
plt.figure()
plt.plot(x, y, '-')
[<matplotlib.lines.Line2D at 0x1108d5470>]
Change the inputs of the linspace and arange function until you understand exactly what they do. Using both ways, create a vector of starting from -1.5 to 4.5 inclusive with a step of 0.5.
# type your solution here
# %load solutions/session2/1D_array.py
x = np.linspace(-1.5,4.5,13)
y = np.arange(-1.5,5,0.5)
print (x)
print (y)
We can always find out what a function does by typing its name, followed by a question mark:
np.arange?
Resize the help window (to get it out of your way) by dragging the divider. Double click on the divider to close it.
We can create multidimensional arrays, like the following:
A = np.array([
[1, 2.4, -13],
[4.1, 5, 0],
[7.2, 8, 9]
])
print("A = \n", A)
print(A.shape)
A = [[ 1. 2.4 -13. ] [ 4.1 5. 0. ] [ 7.2 8. 9. ]] (3, 3)
To see all different attributes of the array object type its name followed by a fullstop and press tap. Note that some attributes are basically functions and require parentheses (with or without arguments), e.g. transpose and others don't, e.g. ndim:
print(A.transpose())
print(A.ndim)
[[ 1. 4.1 7.2] [ 2.4 5. 8. ] [-13. 0. 9. ]] 2
We can also define the shape of the array:
np.random.uniform(0, 1, size=(5,5)) # 5 x 5 matrix with elements from a uniform distribution
array([[ 0.83140826, 0.57806063, 0.02359354, 0.8937764 , 0.59164341], [ 0.13891259, 0.31376313, 0.82242634, 0.80782011, 0.25030733], [ 0.65847677, 0.95811486, 0.1987836 , 0.7035087 , 0.28968929], [ 0.15904085, 0.97726721, 0.41033834, 0.49373835, 0.79527837], [ 0.27062751, 0.15695301, 0.70117262, 0.25383595, 0.03348545]])
We can use functions zeros
and ones
(that we used before) to create matrices.
Also we can use function eye
to create an identity matrix and function diag
to create a diagonal matrix:
print(np.zeros((3,4)))
print('\n')
print(np.ones((2,3)))
print('\n')
print(np.eye(3))
print('\n')
print(np.diag([1, 2, 3]))
[[ 0. 0. 0. 0.] [ 0. 0. 0. 0.] [ 0. 0. 0. 0.]] [[ 1. 1. 1.] [ 1. 1. 1.]] [[ 1. 0. 0.] [ 0. 1. 0.] [ 0. 0. 1.]] [[1 0 0] [0 2 0] [0 0 3]]
Notice how many parentheses we should use: to construct zero matrix of size $2\times3$, we should pass the size as a tuple: (2,
3)
.
To index an element in an array, we use brackets []
:
x = np.array([1, 1, 2, 3, 5, 8, 13])
x[4]
5
And with matrices we put comma between the index of a row and a column:
A = np.array([
[1, 2, 3],
[4, 5, 6],
[7, 8, 9]
])
print(A)
print(A[0, 0])
print(A[-1, 1])
[[1 2 3] [4 5 6] [7 8 9]] 1 8
Notice, how the negative indices are used to index from the end instead of the beginning.
We can index several elements simultaneously: this is called slicing.
print(A)
[[1 2 3] [4 5 6] [7 8 9]]
A[0, 0:2]
array([1, 2])
Note that the fundamental difference from MATLAB is that the indexing starts from zero and is exclusive of the second index. Hence, in the above example 1:3 refers to the second and third row and 0:2 to the first and second column.
We can index in some clever ways too:
print(A[0, :]) # row selection
print(A[:, 2]) # column selection
[1 2 3] [3 6 9]
print(A[:, 1:]) # all rows, from 2nd column until the end
print(A[-1, :]) # last row, all columns
[[2 3] [5 6] [8 9]] [7 8 9]
We can change array's elements:
print(A)
print()
A[:, 2] = [23, 26, 29]
print(A)
[[1 2 3] [4 5 6] [7 8 9]] [[ 1 2 23] [ 4 5 26] [ 7 8 29]]
There is an easy way to reverse a vector:
a = np.arange(0, 10)
print(a)
print(a[: :-1]) # reverse a
[0 1 2 3 4 5 6 7 8 9] [9 8 7 6 5 4 3 2 1 0]
We can also index an array using boolean arrays:
b = a > 5
print(b)
a[b] = 0 # All elements of 'a' higher than 5 become 0
print(a)
[False False False False False False True True True True] [0 1 2 3 4 5 0 0 0 0]
Now we will use the knowledge we obtained to compute a first derivative to $\sin 2\pi x$ numerically.
We will use the finite-difference formula of the second order of accuracy:
$$ y'(x_i) \approx \frac{y_{i+1} - y_{i-1}}{2\Delta x} $$with formulas of the first order of accuracy at endpoints:
$$ y'(x_0) \approx \frac{y_1 - y_0}{\Delta x}, \quad y'(x_N) \approx \frac{y_N - y_{N-1}}{\Delta x} $$#follow the instructor
We can use array functions min
or max
to find minimum or maximum element in the array:
x = np.array([np.pi, 1, 27, 34.0, 5, 44, 15])
print(x.min())
print(x.max())
1.0 44.0
print(A.sum(axis=0)) # sum of each column
print(A.min(axis=1)) # min of each row
print(A.cumsum(axis=1)) # cumulative sum along each row
[12 15 78] [1 4 7] [[ 1 3 26] [ 4 9 35] [ 7 15 44]]
In Python axis refers to the dimensions of arrays; axis=0 refers to the colums and axis=1 to the rows.
# type your solution here
# %load solutions/session2/array_operations.py
print('Part 1:')
print(10+np.arange(90))
print('\nPart 2:')
print(np.nonzero([1,2,0,0,4,0])[0])
print('\nPart 3:')
chess = np.zeros((8,8))
chess[1::2,::2] = 1
chess[::2,1::2] = 1
print(chess)
print('\nPart 4:')
print(np.fromfunction(lambda i, j: (i +1)* (j+1)**2, (5, 3), dtype=int))
print('\nPart 5:')
normal_dist = np.random.normal(3, 1.5, (4,4))
print(normal_dist)
print(normal_dist.sum(axis=1))
Part 1: [10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99] Part 2: [0 1 4] Part 3: [[ 0. 1. 0. 1. 0. 1. 0. 1.] [ 1. 0. 1. 0. 1. 0. 1. 0.] [ 0. 1. 0. 1. 0. 1. 0. 1.] [ 1. 0. 1. 0. 1. 0. 1. 0.] [ 0. 1. 0. 1. 0. 1. 0. 1.] [ 1. 0. 1. 0. 1. 0. 1. 0.] [ 0. 1. 0. 1. 0. 1. 0. 1.] [ 1. 0. 1. 0. 1. 0. 1. 0.]] Part 4: [[ 1 4 9] [ 2 8 18] [ 3 12 27] [ 4 16 36] [ 5 20 45]] Part 5: [[ 0.96015636 2.55634018 4.17590916 2.91010113] [ 5.64964804 3.27467046 4.3661072 0.84234394] [ 0.36539678 4.7330197 3.91428631 1.39736621] [ 5.38352132 0.61407006 4.02491362 1.86491946]] [ 10.60250683 14.13276963 10.41006901 11.88742445]
Now let's work with some real data. Ensure that the stockholm_temp.txt in the same directory with this notebook. This file contains the daily average temperatures according to observations of Stockholm for the years 1756 - 2012. The following code checks if the file exists, opens it, prints the first lines and stores the data in an array.
filename = 'data/session2/stockholm_temp.txt'
data = np.loadtxt(filename) # alternative we can use getfromtxt command
# Print first three lines of data.
print(data[0:3, :])
[[ 1.75600000e+03 1.00000000e+00 1.00000000e+00 -8.70000000e+00] [ 1.75600000e+03 1.00000000e+00 2.00000000e+00 -9.20000000e+00] [ 1.75600000e+03 1.00000000e+00 3.00000000e+00 -8.60000000e+00]]
Since the data are in the form
we would like to express them in a better way. We define a dtype
, that is the data type for each column: int, int, int, float:
dt = np.dtype([('Year', 'int16'), ('Month', 'int8'), ('Day', 'int8'), ('Temp', 'float64')])
data = np.loadtxt(filename, dtype=dt)
data[:10] # first 10 entries of our data
array([(1756, 1, 1, -8.7), (1756, 1, 2, -9.2), (1756, 1, 3, -8.6), (1756, 1, 4, -7.7), (1756, 1, 5, -7.2), (1756, 1, 6, -1.6), (1756, 1, 7, 0.7), (1756, 1, 8, 1.3), (1756, 1, 9, 2.4), (1756, 1, 10, 0.8)], dtype=[('Year', '<i2'), ('Month', 'i1'), ('Day', 'i1'), ('Temp', '<f8')])
We can now look for the example only on the temperatures:
data['Temp']
array([-8.7, -9.2, -8.6, ..., 0.2, 2.6, 2.2])
and plot temperature in July 2012:
years = data['Year']
temp = data['Temp']
plot_data = data[(data['Year'] == 2012) & (data['Month'] == 7)]
plt.figure()
plt.plot(plot_data['Day'], plot_data['Temp'], '-')
plot_data
array([(2012, 7, 1, 17.7), (2012, 7, 2, 16.3), (2012, 7, 3, 17.7), (2012, 7, 4, 17.4), (2012, 7, 5, 16.4), (2012, 7, 6, 17.2), (2012, 7, 7, 19.7), (2012, 7, 8, 17.9), (2012, 7, 9, 19.3), (2012, 7, 10, 17.2), (2012, 7, 11, 16.5), (2012, 7, 12, 16.9), (2012, 7, 13, 14.9), (2012, 7, 14, 15.2), (2012, 7, 15, 15. ), (2012, 7, 16, 16.1), (2012, 7, 17, 14. ), (2012, 7, 18, 14.6), (2012, 7, 19, 14.9), (2012, 7, 20, 15.5), (2012, 7, 21, 13.6), (2012, 7, 22, 16.3), (2012, 7, 23, 15.9), (2012, 7, 24, 20.3), (2012, 7, 25, 20.6), (2012, 7, 26, 17.9), (2012, 7, 27, 16.4), (2012, 7, 28, 19.2), (2012, 7, 29, 18.6), (2012, 7, 30, 17. ), (2012, 7, 31, 13.9)], dtype=[('Year', '<i2'), ('Month', 'i1'), ('Day', 'i1'), ('Temp', '<f8')])
# type your solution here
# %load solutions/session2/temperatures.py
day = 5
indices = np.where(data['Day'] == day) # find indices of 5th day for every month (over all data)
temp = data['Temp'][indices] # find corresponding temperatures
print(temp.size)
print(temp)
3084 [-7.2 -1.2 -1. ..., 8.5 3.9 -4.2]
Handling ndarrays is really powerful. Of the many things one can do with NumPy's arrays is linear algebra:
x = np.random.random(10)
print(x)
[ 0.51551471 0.13270754 0.55979456 0.24638503 0.86521711 0.27927764 0.49640099 0.81195922 0.71634199 0.21532327]
We can find $L_\infty$-norm of this vector:
np.linalg.norm(x, np.inf) # maximum norm of a vector
0.86521710548285347
We can use this module to solve a linear system of equations
$$ \begin{align} x + 2y &= 12, \\ 3x + 4y &= 29 \end{align} $$A = np.array([[1, 2],[3, 4]])
b = np.array([12, 29])
print(A)
print(b)
[[1 2] [3 4]] [12 29]
print(np.linalg.cond(A))
x = np.linalg.solve(A,b)
print(x)
14.9330343737 [ 5. 3.5]
To check that the solution is correct (that is, that $x$ that we computed satisfies $Ax = b$), we do matrix-vector multiplication. For that, function dot
should be used:
np.dot(A, x)
array([ 12., 29.])
and inner products of two vectors are computed with the same function:
x = np.array([1, 2])
y = np.array([3, 4])
np.dot(x, y)
11
We can also solve eigenvalue problems using function linalg.eig
:
lamda, V = np.linalg.eig(A)
print(lamda)
print(V)
[-0.37228132 5.37228132] [[-0.82456484 -0.41597356] [ 0.56576746 -0.90937671]]
Create a diagonal matrix using the lamda variable and then multiply from the left with V and from the right with the inverse of V (use np.linalg.inv()). What do you find?
# type your solution here
# %load solutions/session2/linear_algebra.py
print(np.dot(np.dot(V,np.diag(lamda)),np.linalg.inv(V)))
NumPy
is complemented by SciPy
; a library that adds more MATLAB-like features to Python.
SciPy
consists of several subpackages that can be used for linear algebra, image and signal processing, statistics, numerical integration of differential equations and more.
Here we'll study only a couple of subpackages to get the idea of how useful SciPy
is.
It is customary to import subpackages of SciPy
in the following manner:
Subpackage scipy.integrate
contains functions to compute quadratures and to solve ordinary differential equations.
The numerical evaluation of the definite integral of a function is called numerical quadrature, or simply quadature. SciPy provides a series of functions for different kind of quadrature, for example the quad, dblquad and tplquad for single, double and triple integrals, respectively.
Let's compute the Gauss integral numerically:
$$ \int_{-\infty}^{+\infty} \exp(-x^2) \, dx = \sqrt \pi $$from scipy import integrate, special
xl = -1000 # the lower limit of x
xu = +1000 # the upper limit of x
def f(x):
return np.exp(-x**2)
val, abserr = integrate.quad(f, xl, xu)
print("integral value =", val, ", absolute error =", abserr, 'analytic value = ', np.sqrt(np.pi))
integral value = 1.7724538509055159 , absolute error = 7.768296244985151e-09 analytic value = 1.77245385091
We can pass extra arguments to integrand function, by using the args
keyword argument.
For example we can use args=(3,)
for the third order Bessel function.
# type your solution here
# %load solutions/session2/bessel_integral.py
def integrand(x, n):
"""
Bessel function of first kind and order n.
"""
return special.jv(n, x)
xl = 0 # the lower limit of x
xu = 10 # the upper limit of x
val, abserr = integrate.quad(integrand, xl, xu, args=(0,))
print(val, abserr)
We will have a look at the Lotka-Volterra model, also known as the predator-prey equations, which is a pair of first order, non-linear, differential equations frequently used to describe the dynamics of biological systems in which two species interact, one a predator and the other its prey. The model was proposed independently by Alfred J. Lotka in 1925 and Vito Volterra in 1926, and can be described by
$$ \begin{align} \frac{du}{dt} &= \phantom{-}au - \phantom{d}buv, \\ \frac{dv}{dt} &= -cv + dbuv, \\ \end{align} $$where
$u(t)$: number of preys (for example, rabbits)
$v(t)$: number of predators (for example, foxes)
$a$, $b$, $c$, $d$ are constant parameters defining the behavior of the population:
$a$ is the natural growing rate of rabbits, when there's no fox
$b$ is the natural dying rate of rabbits, due to predation
$c$ is the natural dying rate of fox, when there's no rabbit
$d$ is the factor describing how many caught rabbits let create a new fox
a = 1.0
b = 0.1
c = 1.5
d = 0.75
def dy_dt(t, y):
"""Compute right-hand side of Lotka--Volterra system."""
u, v = y[0], y[1]
return np.array([
a*u - b*u*v,
-c*v + d*b*u*v
])
# Initial value.
y0 = [10, 5]
t0 = 0.0
dt = 0.1
tfinal = 10
# Two lists that will hold the solution.
times = np.linspace(0, tfinal, num=101)
u = [y0[0]]
v = [y0[1]]
s = integrate.ode(dy_dt)
s.set_integrator('dopri5')
s.set_initial_value(y0, t0)
for t in times[1:]:
s.integrate(t)
u.append(s.y[0])
v.append(s.y[1])
plt.figure()
plt.plot(times, u, '-', label='Rabbits')
plt.plot(times, v, '--', label='Foxes')
plt.legend(loc='best');
The interpolate.interp1d
function accepts array of points $(X, Y)$ and returns a new function that can find interpolated values from this array.
You can use an options for interpoland kinds such as linear, quadratic and cubic spline interpolation.
from scipy import interpolate
x = np.linspace(0, 1, num=11)
y = np.sin(2 * np.pi * x)
f1 = interpolate.interp1d(x, y) # linear interpolation
f3 = interpolate.interp1d(x, y, kind='cubic')
# new points where we want to evaluate the interpoland
xnew = np.linspace(0, 1, num=101)
plt.figure(figsize=(12,4))
plt.plot(x, y,'o', label='Given')
plt.plot(xnew, f1(xnew), '-', label='linear')
plt.plot(xnew, f3(xnew), '--', label='cubic')
plt.legend(loc='best');
When you need to study the spectrum of a signal (function of time), you can use Fourier Transform for this.
The implementation of the Fourier transform, known as Fast Fourier Transform, can be used with the help of subpackage scipy.fftpack
.
Let's apply Fast Fourier Transform to the signal:
$$ y(t) = \sin 2\pi f t + \text{normal noise},$$where $f = 4$.
Suppose, that the signal is sampled with sampling rate 50 times per unit time (time step is 0.02). Hence, we can compute the spectrum with upper bound of frequency being 25.
from scipy import fftpack
f0 = 4
timestep = 0.02
tfinal = 10.0
# Create the oscillator
time = np.arange(0, tfinal, timestep)
y = np.sin(2*np.pi*f0 * time)
y = y + np.random.normal(size=len(y))
# Use an FFT to calculate its spectrum
yhat = fftpack.rfft(y)
# Find the positive frequencies
frequency = fftpack.rfftfreq(y.size, d=timestep)
spectrum = 2*timestep*np.abs(yhat)
# Create a figure
fig, (ax1, ax2) = plt.subplots(2, 1)
# Create x-y plots of the amplitude and transform with labeled axes
ax1.set_xlabel('Time')
ax1.set_ylabel('Signal')
ax1.plot(time, y)
ax2.set_xlabel('Frequency')
ax2.set_ylabel('Power spectrum')
ax2.plot(frequency, spectrum, '-')
fig.tight_layout(pad=0.1)
Sparse matrices are often useful in numerical analysis dealing with large systems. Whenever the problem involves vectors or matrices with mostly zero values then one should take advantage of the sparsity. SciPy has an extensive library for sparse matrices, with basic linear algebra operations (such as equation solving, eigenvalue calculations, etc).
To create a sparse matrix we have to choose the format that will be stored in:
from scipy import sparse
# dense matrix
M = np.array([[10,0,0,0], [0,.3,0,0], [0,2,1.5,0], [-3,0,0,100]])
M
array([[ 10. , 0. , 0. , 0. ], [ 0. , 0.3, 0. , 0. ], [ 0. , 2. , 1.5, 0. ], [ -3. , 0. , 0. , 100. ]])
# convert from dense to sparse
M_sparse = sparse.csr_matrix(M); # or csc_matrix(M)
M_sparse
<4x4 sparse matrix of type '<class 'numpy.float64'>' with 6 stored elements in Compressed Sparse Row format>
We can go back to the dense form:
print(M_sparse)
print(M_sparse.todense())
(0, 0) 10.0 (1, 1) 0.3 (2, 1) 2.0 (2, 2) 1.5 (3, 0) -3.0 (3, 3) 100.0 [[ 10. 0. 0. 0. ] [ 0. 0.3 0. 0. ] [ 0. 2. 1.5 0. ] [ -3. 0. 0. 100. ]]
The easiest way to create a sparse matrix is to use the lil_matrix
function.
It creates an empty matrix and then we allocating non-zero entries by using matrix indexing.
This way we avoid creating large dense matrices.
A = sparse.lil_matrix((4,4)) # empty 4x4 sparse matrix
A[0,0] = 1
A[1,1] = 3
A[2,2] = A[2,1] = 1
A[3,3] = A[3,0] = 1
A
<4x4 sparse matrix of type '<class 'numpy.float64'>' with 6 stored elements in LInked List format>
We can also change the format based on the sparsity:
A = sparse.csr_matrix(A); A # if the rows are sparse
<4x4 sparse matrix of type '<class 'numpy.float64'>' with 6 stored elements in Compressed Sparse Row format>
A = sparse.csc_matrix(A); A # if the columns are sparse
<4x4 sparse matrix of type '<class 'numpy.float64'>' with 6 stored elements in Compressed Sparse Column format>
Construct a $10000 \times 10000$ lil_matrix
$A$ and add values to it by setting all rows and columns $2000$ to $3000$ into random values. Also set random diagonal entries (you can do this without a loop!). Then use the spy
function from matplotlib.pyplot
to "spy" the matrix.
Using a vector $b$ with $10000$ random values solve the system $A x = b$.
To do so use the direcet solver linsolve
from scipy.sparse.linalg.dsolve
.
Change to different sparse formats (CSR and CSC) and measure the time needed to solve the system. You can use the time
package.
# type your solution here
# %load solutions/session2/sparse_matrices.py
"""
Construct a 1000x1000 lil_matrix and add some values to it, convert it
to CSC format and solve A x = b for x with a direct solver.
"""
%pylab inline --no-import-all
from matplotlib import pyplot as plt
import numpy as np
import scipy.sparse as sps
from scipy.sparse.linalg.dsolve import linsolve
import time
rand = np.random.rand
A = sps.lil_matrix((10000, 10000), dtype=np.float64)
A[0::100, 2000:3000] = rand(1000)
A.setdiag(rand(10000))
A = A.tocsc() # Compressed Sparse Column format
# Spy matrix
plt.spy(A, marker='.', markersize=2)
plt.show()
b = rand(10000) # random rhs
t1 = time.time() # initial time
x = linsolve.spsolve(A, b) # solve system
t2 = time.time() # final time
print('residual:', np.linalg.norm(A * x - b))
print(t2-t1)
Populating the interactive namespace from numpy and matplotlib
residual: 1.48698589756e-12 0.010876178741455078
Besides NumPy and Scipy, there are many other useful Python packages for scientific computing. Here is a short list:
Copyright 2017, ACM/SIAM Student Chapter of King Abdullah University of Science and Technology