# CME 193 - Lecture 5 - SciPy

SciPy is a library built on top of NumPy that contains functionality common in scientific computing. Here is a list of subpackages it contains:

- cluster:	Clustering algorithms
- constants:	Physical and mathematical constants
- fftpack:	Fast Fourier Transform routines
- **integrate:	Integration and ordinary differential equation solvers**
- interpolate:	Interpolation and smoothing splines
- io:	Input and Output
- **linalg:	Linear algebra**
- ndimage:	N-dimensional image processing
- odr:	Orthogonal distance regression
- **optimize:	Optimization and root-finding routines**
- signal:	Signal processing
- sparse:	Sparse matrices and associated routines
- spatial:	Spatial data structures and algorithms
- special:	Special functions
- stats:	Statistical distributions and functions

Today we'll "tour" the ones in **bold**.

# Solving linear systems of equations
By far the most common operation in scientific computing is solving a system of equations, specified by $Ax = b$. For example:

$$ 
\begin{align}
3x + 2y + 5z &=  2 \\
x - y + z &= 4 \\
6y  + z &= -1
\end{align}
$$


```python
import numpy as np
import scipy.linalg as linalg
```


```python
a = np.array([[3, 2, 5], [1, -1, 1], [0, 5, 1]])
```


```python
b = np.array([2, 4, -1])
```


```python
x = linalg.solve(a, b)
x
```


```python
a @ x
```


```python
# A matrix with more rows than columns (more equations than unknowns)
a = np.array([[1, 2, 1],
              [1, 1, 2],
              [2, 1, 1],
              [1, 1, 1]])
b = np.array([4, 3, 5, 4])
```


```python
linalg.solve(a, b)
```


```python
linalg.lstsq(a, b)
```


```python
x = linalg.lstsq(a, b)[0]
a @ x
```


```python
linalg.norm(a @ x - b)**2
```

# Optimization using `scipy.optimize`

Often in scientific computing you'll want to find the minimum of a function. For example, consider the following function.


```python
import matplotlib.pyplot as plt

def f(x):
    return x**4 - x**2

x = np.linspace(-1.2, 1.2, 100)
y = f(x)

plt.plot(x, y)
plt.grid()
plt.xlabel("$x$", fontsize=16)
plt.ylabel("$x^4 - x^2$", fontsize=16)
plt.show()
```


```python
import scipy.optimize as opt
```

The simplest way to minimize a function is using [`scipy.optimize.minimize`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html).


```python
result = opt.minimize(f, x0=0.5)
result
```


```python
plt.plot(x, y)
plt.plot(result.x, f(result.x), "ro")

plt.grid()
plt.xlabel("$x$", fontsize=16)
plt.ylabel("$x^4 - x^2$", fontsize=16)
plt.show()
```


```python
def g(z):
    x, y = z
    return np.sin(x + 2) * np.sin(y + 2)

result = opt.minimize(g, x0=np.array([0, 1]))
result
```

## Linear programming
A very common case is linear programming (LP). These are optimization problems that can be written in the form

$$
\begin{equation}
\begin{split}
\text{minimize} \;\; & c^{T}x  \\
\text{subject to} \;\; & A_{ub}x \leq b_{ub} \\
& A_{eq}x = b_{eq}
\end{split}
\end{equation}
$$

Here, we are finding the vector $x$ that minimizes the dot product $c^T x$, where $c$ is some fixed vector, out of all $x$ that satisfy $A_{ub}x \leq b_{ub}$ and $A_{eq}x = b_{eq}$, where $A_{ub}$ and $A_{eq}$ are matrices and $b_{ub}$ and $b_{eq}$ are vectors.

## Exercise 6
Using Google and reading documentation are important parts of programming. `scipy.optimize` comes with specialized functions for solving linear programming problems.

Figure out how to solve LPs using `scipy.optimize`, and solve the following LP:
$$
\begin{equation}
\begin{split}
\text{minimize} \;\; & x_1 + 2 x_2  \\
\text{subject to} \;\; & x_1 \leq 1 \\
& 5 x_1 + x_2 \geq 0
\end{split}
\end{equation}
$$

Note that the problem is equivalent to

$$
\begin{equation}
\begin{split}
\text{minimize} \;\; & c^T x  \\
\text{subject to} \;\; & A_{ub}x \leq b_{ub} \\
& A_{eq}x = b_{eq}
\end{split}
\end{equation}
$$
where
$$
c = \begin{pmatrix} 1 \\ 2 \end{pmatrix} \qquad 
A_{ub} = \begin{pmatrix} 1 & 0 \\ -5 & -1 \end{pmatrix}, \qquad
b_{ub} = \begin{pmatrix} 1 \\ 0 \end{pmatrix}, \qquad
A_{eq} = 0, \qquad
b_{eq} = 0.
$$

We can see this because 
$$
 \begin{pmatrix} 1 & 0 \\ -5 & -1 \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \end{pmatrix} = \begin{pmatrix} x_1 \\ -5x_1 - x_2 \end{pmatrix} 
$$


```python
# YOUR CODE HERE

```

## Curve fitting


```python
def f(x, a, b):
    return a * np.sin(b * x)

x = np.linspace(-2, 2, 30)
y = f(x, 1.3, 3) + 0.3 * np.random.randn(30)
plt.plot(x, y, 'ro')
plt.show()
```


```python
(a, b), cov = opt.curve_fit(f, x, y)
a, b
```


```python
x_fine = np.linspace(-2, 2, 200)
plt.plot(x, f(x, a, b))
plt.plot(x, y, 'ro')
plt.show()
```

# Numerical integration using `scipy.integrate`


```python
import scipy.integrate as integrate
```

Here's how to compute $$
\int_0^\pi \sin(x) \,e^{-x} \,dx 
$$


```python
def f(x):
    return np.sin(x) * np.exp(-x)
```


```python
integrate.quad(f, 0, np.pi)
```


```python
def f_antiderivative(x):
    return -np.exp(-x) * (np.sin(x) + np.cos(x)) / 2
```


```python
f_antiderivative(np.pi) - f_antiderivative(0)
```


```python
integrate.quad(f, 0, np.inf)
```


```python
f_antiderivative(np.inf) - f_antiderivative(0)
```

## Physical simulation
SciPy can solve ODEs (ordinary differential equations). Let's simulate a throwing a ball in two dimensions. The ball is described as a function of time by four functions: $x(t), y(t), v_x(t), v_y(t)$, which are governed by the ODE system:

$$
    \frac{dx}{dt} = v_x, \qquad
    \frac{dy}{dt} = v_y, \qquad
    \frac{dv_x}{dt} = 0, \qquad
    \frac{dv_y}{dt} = -g.
$$

We can think of this in vectorized form as:

$$
    \frac{d}{dt} \begin{pmatrix} x \\ y \\ v_x \\ v_y \end{pmatrix} = \begin{pmatrix} v_x \\ v_y \\ 0 \\ -g \end{pmatrix}
$$


```python
def dzdt(z, t):
    x, y, vx, vy = z
    g = 1
    return np.array([vx, vy, 0, -g])
```


```python
x, y = 0, 0
vx, vy = 1, 20
z = np.array([x, y, vx, vy])
t = np.linspace(0, 50, 20)

result = integrate.odeint(dzdt, z, t)
print(result.shape)
plt.scatter(result[:, 0], result[:, 1])
plt.show()
```

## Exercise 7
We can add drag with the following slight modification to the ODE, which adds a force with direction opposing the current velocity, and with magnitude proportional to the velocity squared.
    
$$
    \frac{d}{dt} \begin{pmatrix} x \\ y \\ v_x \\ v_y \end{pmatrix} = \begin{pmatrix} v_x \\ v_y \\ -\alpha v_x \sqrt{v_x^2 + v_y^2} \\ -g -\alpha v_y \sqrt{v_x^2 + v_y^2} \end{pmatrix}
$$

Implement this with $\alpha = 0.001$, plot the resulting trajectory, and compare to the dragless case.


```python
# YOUR CODE HERE
```


```python

```

5b343a43e8a31a333ec93daa1b1a447465d22839
ipynb
Jupyter Notebook
nb/2019_spring/Lecture_5.ipynb
samuelcheang0419/cme193
609e4655544292a28dbb9ca0301637b006970af2
MIT
2016-02-17T06:03:51.000Z
2021-11-30T03:47:27.000Z
nb/2019_spring/Lecture_5.ipynb
samuelcheang0419/cme193
609e4655544292a28dbb9ca0301637b006970af2
MIT
nb/2019_spring/Lecture_5.ipynb
samuelcheang0419/cme193
609e4655544292a28dbb9ca0301637b006970af2
MIT
2016-01-19T18:23:46.000Z
2020-12-23T03:08:23.000Z
Qwen/Qwen-72B
1. YES
2. YES