# Scientific Computing with Python
<a rel="license" href="http://creativecommons.org/licenses/by/4.0/"></a><br />This notebook by Xiaozhou Li is licensed under a <a rel="license" href="http://creativecommons.org/licenses/by/4.0/">Creative Commons Attribution 4.0 International License</a>.  
All code examples are also licensed under the [MIT license](http://opensource.org/licenses/MIT).


```python
%matplotlib inline
import matplotlib.pyplot as plt
```

## Numpy
### Introduction
The numpy package (module) is used in almost all numerical computation using Python. It is a package that provide high-performance vector, matrix and higher-dimensional data structures for Python. It is implemented in C and Fortran so when calculations are vectorized (formulated with vectors and matrices), performance is very good.

### Array Creation and Properties
NumPy arrays allow you to write fast (optimized) code that works on arrays of data.  To do this, there are some restrictions on arrays:

* all elements are of the same data type (e.g. float)
* the size of the array is fixed in memory, and specified when you create the array (e.g., you cannot grow the array like you do with lists)

<div style="background-color: #9999ff; padding: 10px;">Note that the traditional way to import numpy is to rename it np.  This saves on typing and makes your code a little more compact.</div>


```python
import numpy as np
```


```python
a = np.arange(12)
```


```python
a
```




    array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11])




```python
a = np.arange(12).reshape(3,4)
```


```python
a
```




    array([[ 0,  1,  2,  3],
           [ 4,  5,  6,  7],
           [ 8,  9, 10, 11]])



an array is an object of the `ndarray` class, and it has methods that know how to work on the array data.  Here, `.reshape()` is one of the many methods.

A NumPy array has a lot of meta-data associated with it describing its shape, datatype, etc.


```python
print(a.ndim)
print(a.shape)
print(a.size)
print(a.dtype)
print(a.itemsize)
print(type(a))
```

    2
    (3, 4)
    12
    int64
    8
    <class 'numpy.ndarray'>



```python
#help(a)
```

we can create a multi-dimensional array of a specified size initialized all to 0 easily, using the zeros() function.


```python
b = np.zeros((5,4))
b
```




    array([[0., 0., 0., 0.],
           [0., 0., 0., 0.],
           [0., 0., 0., 0.],
           [0., 0., 0., 0.],
           [0., 0., 0., 0.]])



There is also an analogous ones() and empty() array routine.  Note that here we can explicitly set the datatype for the array in this function if we wish. 


```python
c = np.eye(10, dtype=np.float64)
c
```




    array([[1., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
           [0., 1., 0., 0., 0., 0., 0., 0., 0., 0.],
           [0., 0., 1., 0., 0., 0., 0., 0., 0., 0.],
           [0., 0., 0., 1., 0., 0., 0., 0., 0., 0.],
           [0., 0., 0., 0., 1., 0., 0., 0., 0., 0.],
           [0., 0., 0., 0., 0., 1., 0., 0., 0., 0.],
           [0., 0., 0., 0., 0., 0., 1., 0., 0., 0.],
           [0., 0., 0., 0., 0., 0., 0., 1., 0., 0.],
           [0., 0., 0., 0., 0., 0., 0., 0., 1., 0.],
           [0., 0., 0., 0., 0., 0., 0., 0., 0., 1.]])



`linspace` creates an array with evenly space numbers.  The `endpoint` optional argument specifies whether the upper range is in the array or not


```python
d = np.linspace(0, 1, 10, endpoint=False)
print(d)
d = np.linspace(0, 1, 10, endpoint=True)
print(d)
d = np.linspace(0, 1, 11)
print(d)
```

    [0.  0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9]
    [0.         0.11111111 0.22222222 0.33333333 0.44444444 0.55555556
     0.66666667 0.77777778 0.88888889 1.        ]
    [0.  0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1. ]


### Slicing


```python
a = np.arange(9)
a
```




    array([0, 1, 2, 3, 4, 5, 6, 7, 8])



Now look at accessing a single element vs. a range (using slicing)

Giving a single (0-based) index just references a single value


```python
a[2]
```




    2



Giving a range uses the range of the edges to return the values


```python
print(a[2:3])
```

    [2]



```python
a[2:4]
```




    array([2, 3])



The : can be used to specify all of the elements in that dimension


```python
a[:]
```




    array([0, 1, 2, 3, 4, 5, 6, 7, 8])



### Multidimensional Arrays

Multidimensional arrays are stored in a contiguous space in memory -- this means that the columns / rows need to be unraveled (flattened) so that it can be thought of as a single one-dimensional array.  Different programming languages do this via different conventions:

Storage order:

* Python/C use *row-major* storage: rows are stored one after the other
* Fortran/matlab use *column-major* storage: columns are stored one after another

The ordering matters when 

* passing arrays between languages (we'll talk about this later this semester)
* looping over arrays -- you want to access elements that are next to one-another in memory
  * e.g, in Fortran:
  ```
  double precision :: A(M,N)
  do j = 1, N
     do i = 1, M
        A(i,j) = …
     enddo
  enddo
  ```
  
  * in C
  ```
  double A[M][N];
  for (i = 0; i < M; i++) {
     for (j = 0; j < N; j++) {
        A[i][j] = …
     }
  }  
  ```
  

In python, using NumPy, we'll try to avoid explicit loops over elements as much as possible

Let's look at multidimensional arrays:


```python
a = np.arange(15).reshape(3,5)
a
```




    array([[ 0,  1,  2,  3,  4],
           [ 5,  6,  7,  8,  9],
           [10, 11, 12, 13, 14]])



Notice that the output of `a` shows the row-major storage.  The rows are grouped together in the inner `[...]`

Giving a single index (0-based) for each dimension just references a single value in the array


```python
a[1,1]
```




    6



Doing slices will access a range of elements.  Think of the start and stop in the slice as referencing the left-edge of the slots in the array.


```python
a[0:2,0:2]
```




    array([[0, 1],
           [5, 6]])



Access a specific column


```python
a[:,1]
```




    array([ 1,  6, 11])



Sometimes we want a one-dimensional view into the array -- here we see the memory layout (row-major) more explicitly


```python
a.flatten()
```




    array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14])



we can also iterate -- this is done over the first axis (rows)


```python
for row in a:
    print(row)
```

    [0 1 2 3 4]
    [5 6 7 8 9]
    [10 11 12 13 14]


or element by element


```python
for e in a.flat:
    print(e)
```

    0
    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14


## Scipy
The SciPy framework builds on top of the low-level NumPy framework for multidimensional arrays, and provides a large number of higher-level scientific algorithms. 

SciPy is a collection of numerical algorithms with python interfaces.  In many cases, these interfaces are wrappers around standard numerical libraries that have been developed in the community and are used with other languages.  Usually detailed references are available to explain the implementation.

There are many subpackages generally, you load the subpackages separately, e.g.

```
from scipy import linalg, optimize
```
then you have access to the methods in those namespaces

## Matplotlib
### Introduction
Matplotlib is an excellent 2D and 3D graphics library for generating scientific figures.

### Reading and writing a panda

**(1) original figure**


```python
plt.figure()
img = plt.imread('../data/panda.jpg')
plt.imshow(img)
plt.imsave("original.jpg",img)
```

**(2) red channel displayed in grey**


```python
plt.figure()
img_red = img[:, :, 0]
plt.imshow(img_red, cmap=plt.cm.gray)
```

**(3) lower resolution (compression)**


```python
plt.figure()
img_tiny = img[::8, ::8]
plt.imshow(img_tiny, interpolation='nearest') 
#plt.savefig("compressed.jpg")
plt.imsave("compressed.jpg",img_tiny)
```

### Mandlebrot Set (Mandelbrot fractal)


```python
def compute_mandelbrot(N_max, some_threshold, nx, ny):
    # A grid of c-values
    x = np.linspace(-2, 1, nx)
    y = np.linspace(-1.5, 1.5, ny)

    c = x[:,np.newaxis] + 1j*y[np.newaxis,:]

    # Mandelbrot iteration

    z = c
    for j in range(N_max):
        z = z**2 + c

    mandelbrot_set = (abs(z) < some_threshold)

    return mandelbrot_set

mandelbrot_set = compute_mandelbrot(50, 50., 601, 401)

plt.imshow(mandelbrot_set.T, extent=[-2, 1, -1.5, 1.5])
plt.gray()
plt.show()
```

### A simple example of 3D plotting
$$ z = \sin(\sqrt{x^2 + y^2}) $$


```python
from mpl_toolkits.mplot3d import Axes3D

fig = plt.figure()
ax = Axes3D(fig)
X = np.arange(-4, 4, 0.25)
Y = np.arange(-4, 4, 0.25)
X, Y = np.meshgrid(X, Y)
R = np.sqrt(X ** 2 + Y ** 2)
Z = np.sin(R)

ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap=plt.cm.hot)
ax.contourf(X, Y, Z, zdir='z', offset=-2, cmap=plt.cm.hot)
ax.set_zlim(-2, 2)

plt.show()
```

### An example displaying the contours of a function
$$ f(x,y) = \left(1 - \frac{x}{2} + x^5 + y^3\right)e^{-x^2-y^2}.$$


```python
def f(x,y):
    return (1 - x / 2 + x**5 + y**3) * np.exp(-x**2 -y**2)

n = 256
x = np.linspace(-3, 3, n)
y = np.linspace(-3, 3, n)
X,Y = np.meshgrid(x, y)

plt.axes([0.025, 0.025, 0.95, 0.95])

plt.contourf(X, Y, f(X, Y), 8, alpha=.75, cmap=plt.cm.hot)
C = plt.contour(X, Y, f(X, Y), 8, colors='black')
plt.clabel(C, inline=1, fontsize=10)

plt.xticks(())
plt.yticks(())
plt.show()
```

## Sympy - Symbolic algebra in Python

### Introduction
There are two notable Computer Algebra Systems (CAS) for Python:

* [SymPy](http://sympy.org/en/index.html) - A python module that can be used in any Python program, or in an IPython session, that provides powerful CAS features. 
* [Sage](http://www.sagemath.org/) - Sage is a full-featured and very powerful CAS enviroment that aims to provide an open source system that competes with Mathematica and Maple. Sage is not a regular Python module, but rather a CAS environment that uses Python as its programming language.

Sage is in some aspects more powerful than SymPy, but both offer very comprehensive CAS functionality. The advantage of SymPy is that it is a regular Python module and integrates well with the Jupyter notebook. 


```python
from sympy import *

init_printing()
```

### Expand, factor and simplify


```python
x, y = symbols('x y')

(x+1)*(x+2)*(x+3)*(x+4)*(x+5)
```


```python
expand((x+1)*(x+2)*(x+3)*(x+4)*(x+5))
```


```python
sin(x+y)
```


```python
expand(sin(x+y), trig=True)
```


```python
expand((x+y)**8)
```


```python
x**3 + 6 * x**2 + 11*x + 6
```


```python
factor(x**3 + 6 * x**2 + 11*x + 6)
```


```python
sin(x)**2 + cos(x)**2
```


```python
simplify(sin(x)**2 + cos(x)**2)
```


```python
cos(x)/sin(x)
```


```python
simplify(cos(x)/sin(x))
```

### Calculus
**(1) differentiation and integration**

$f(x) = (x+1)^2$


```python
f = (x+1)**2
f
```

Computing $\frac{d f}{dx}$, $\frac{d f^2}{dx}$


```python
diff(f,x)
```


```python
diff(f**2,x)
```

Computing $\frac{d \sin(f)}{dx}$, $\frac{d^2 \sin(f)}{dx^2}$


```python
diff(sin(f),x)
```


```python
diff(sin(f),x,2)
```


```python
diff(sin(f),x,4)
```

$$ f(x,y) = \sin(xy) + \cos(xy),$$
computing
$$ \frac{\partial^3 f}{\partial x \partial y^2},\quad \int f(x,y)\,dx,\quad \int_{-1}^{1}f(x,y)\,dx$$


```python
f = sin(x*y) + cos(y*x)
f
```


```python
diff(f, x, 1, y, 2)
```


```python
integrate(f, x)
```


```python
integrate(f, (x, -1, 1))
```

Computing $\int_{-\infty}^\infty e^{-x^2}\,dx$


```python
integrate(exp(-x**2), (x, -oo, oo))
```

**(2) limits**
$$ \lim\limits_{x\rightarrow 0}\frac{\sin(x)}{x},\quad \lim\limits_{x\rightarrow 0^{+}}\frac{1}{x},\quad \lim\limits_{x\rightarrow 0^{-}}\frac{1}{x}$$


```python
limit(sin(x)/x, x, 0)
```


```python
limit(1/x, x, 0, dir="+")
```


```python
limit(1/x, x, 0, dir="-")
```

**(3) series**


```python
exp(x)
```


```python
series(exp(x), x)
```


```python
series(exp(x), x, 1)
```


```python
series(sin(x), x, 0, 12)
```


```python
series(sin(x)*cos(x), x, 0, 8)
```


```python
series(sin(x)*cos(x)*exp(x), x, 0, 12)
```

### Linear algebra: Matrices


```python
m11, m12, m21, m22 = symbols("m11, m12, m21, m22")
b1, b2 = symbols("b1, b2")
```


```python
A = Matrix([[m11, m12],[m21, m22]])
A
```


```python
b = Matrix([[b1], [b2]])
b
```


```python
A**2
```


```python
A**3
```


```python
A * b
```


```python
A.det()
```


```python
A.inv()
```

### Solving equations
Solving 
$$ x^2 - 1 = 0,\quad x^4 - x^2 - 1 = 0$$


```python
solve(x**2 - 1, x)
```


```python
solve(x**4 + x**3 - x**2 - 1, x)
```

Solving systems:
$$ x + y - 1 = 0,\quad x - y - 1 = 0,$$
and
$$ x + y - a = 0,\quad x - y - b = 0.$$


```python
solve([x + y - 1, x - y - 1], [x,y])
```


```python
a, b = symbols('a, b')
solve([x + y - a, x - y - b], [x,y])
```


```python

```

eb142285f6d9c5d85329dc799ad4bbaa73cf0cdf
ipynb
Jupyter Notebook
Numerical_PDEs/Scientific_Python.ipynb
xiaozhouli/Jupyter
68d5a384dd939b3e8079da4470d6401d11b63a4c
MIT
2020-02-27T13:09:06.000Z
2021-11-14T09:50:30.000Z
Numerical_PDEs/Scientific_Python.ipynb
xiaozhouli/Jupyter
68d5a384dd939b3e8079da4470d6401d11b63a4c
MIT
Numerical_PDEs/Scientific_Python.ipynb
xiaozhouli/Jupyter
68d5a384dd939b3e8079da4470d6401d11b63a4c
MIT
2018-10-18T10:20:56.000Z
2021-09-24T08:09:27.000Z
Qwen/Qwen-72B
1. YES
2. YES