text:
```python
from IPython.core.display import HTML
def css_styling():
    styles = open("./styles/custom.css", "r").read()
    return HTML(styles)
css_styling()

```

### BEFORE YOU DO ANYTHING...
In the terminal:
1. Navigate to __inside__ your ILAS_Python repository.
2. __COMMIT__ any un-commited work on your personal computer.
3. __PULL__ any changes *you* have made using another computer.
4. __PULL__ textbook updates (including homework answers).

1. __Open Jupyter notebook:__   Start >> Programs (すべてのプログラム) >> Programming >> Anaconda3 >> JupyterNotebook
1. __Navigate to the ILAS_Python folder__.
1. __Open today's seminar__  by clicking on 7_Numerical_computation_with_Numpy.

Numerical Computation with Numpy

Lesson Goal

To make solutions to mathematical problems involving large data sets, multiple repeated computations and simultaneous solutions more efficient using programming.

To acheive this we will cover an overview of some widely used tools from the Python Numpy package and related packages such as Matplotlib (plotting), Scipy (scientific computation) and Sympy (symbolic).

## Lesson Structure
- The `array` data structure for numerical computation.
- Subpackages
- Array manipulation (independent study)
- Introduction to plotting.
- Mathematics with arrays
 - Mathematics with vectors
 - Mathematics with matrices
 - Useful matrix operations (independent study)
 - Simultaneous equations
 - Rates of change
- Vectorised functions
- Broadcasting
- Reshaping and resizing

## Why are we studying this?
Numerical computation is central to almost all scientific and engineering problems.

There are programming languages specifically designed for numerical computation:
- Fortran
- MATLAB

There are libraries dedicated to efficient numerical computations:
- Numpy
- Scipy
- Sympy ...

NumPy (http://www.numpy.org/)
 - The most widely used Python library for numerical computations.
 - Large, extensive library of data structures and functions for numerical computation.
 - Useful for perfoming operation you will learn on mathematics-based courses.

Scipy (https://www.scipy.org/)
- Builds on Numpy, additional functionality
- More specialised data structures and functions over NumPy.

Matplotlib (https://matplotlib.org/)
- A library for generating plots.
- Complementary to numerical computation libraries.

If you are familiar with MATLAB, NumPy and SciPy provide similar functionality.

Last week we covered an introduction to some basic functions of Numpy.

NumPy is a very extensive library.

This seminar will:
- Introduce some useful functions
- Briefly discuss how to search for additional functions you may need.

Use online resources to search for functions you need e.g.  http://stackoverflow.com/.

## Importing the NumPy module

To make NumPy functions and variables available to use in our program in our programs, we need to __import__ it using.

`import numpy`

We typically import all modules at the start of a program or notebook.

```python
import numpy as np
```

The shortened name `np` is often used for numpy.

All Numpy functions can be called using `np.function()`.

## Data Structure: The Numpy `array`

### Why do we need another data structure?

Python lists hold 'arrays' of data.

Lists are very flexible. e.g. holding mixed data type.

There is a trade off between flexibility and performance e.g. speed.

Science engineering and mathematics problems often involve large amounts of data and numerous operations.

We therefore use specialised functions and data structures for numerical computation.

```python

```

## Numpy array

A numpy array is a grid of values, *all of the same type*.

To create an array we use the Numpy `np.array()` function.

We can create an array in a number of ways.

Let's start with something that is already familiar to you...

We can give a data structure (list, tuple) as an *argument* to convert it to a numpy array:

```python
a = (4.0,)

b = np.array(a)

print(type(a))
print(type(b))
print(b.dtype)
```

    <class 'tuple'>
    <class 'numpy.ndarray'>
    float64

The method `dtype` tells us the type of the data contained in the array.

__Note:__The data type can be optionally set by the user when creating the array. This can be useful when you need more control over how your data is stored in memory and on disk. Especially in cases where you’re working with large data, it’s good that you know to control the storage type.

```python
c = [4.0, 5, 6.0]

d = np.array(c)

print(type(c))
print(type(d))
print(d.dtype)
```

    <class 'list'>
    <class 'numpy.ndarray'>
    float64

## Multi-dimensional arrays.

Unlike the data types we have studied so far, arrays can have multiple dimensions.

__`shape`:__ a *tuple* of *integers* giving the *size* of the array along each *dimension*.

We define the dimensions of an array using square brackets

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

# 2-dimensional array
b = np.array([[1, 2, 3], [4, 5, 6]])

b = np.array([[1, 2, 3],
              [4, 5, 6]])

print(a.shape)
print(b.shape)

```

    (3,)
    (2, 3)

```python
# 2-dimensional array
c = np.array([[1, 2, 3]])

# 2-dimensional array
d = np.array([[1],
              [4]])

print(c.shape)
print(d.shape)
```

    (1, 3)
    (2, 1)

```python
# 3-dimensional array

c = np.array(
    [[[1, 1],
      [1, 1]],

     [[1, 1],
      [1, 1]]])

print(c.shape)

c = np.array(
    [[[1, 1],
      [1, 1]],

     [[1, 1],
      [1, 1]],

     [[1, 1],
      [1, 1]]])

print(c.shape)
```

    (2, 2, 2)
    (3, 2, 2)

```python
# 3-dimensional array

c = np.array(
    [[[1, 1],
      [1, 1]],

     [[1, 1],
      [1, 1]]])

# 4-dimensional array
d = np.array(
    [[[[1, 1],
       [1, 1]],

      [[1, 1],
       [1, 1]]],

      [[[1, 1],
       [1, 1]],

      [[1, 1],
       [1, 1]]]])

print(c.shape)
print(d.shape)
```

    (2, 2, 2)
    (2, 2, 2, 2)

## Creating a numpy array.

We don't always have to manually create the individual elements of an array.

There are several other ways to do this.

For example, if you don’t know what data you want to put in your array you can initialise it with placeholders and load the data you want to use later.

```python
# Create an empty matrix
# The empty() function argument is the shape.
# Shape: tuple of integers giving the size along each dimension.
import numpy as np
x = np.empty((4))
print(x)

print()

x = np.empty((4,4))
print(x)
```

    [  6.95093668e-310   4.66954888e-310   6.95093690e-310   6.95093051e-310]

    [[  6.95093668e-310   4.66954920e-310   0.00000000e+000   0.00000000e+000]
     [  0.00000000e+000   0.00000000e+000   0.00000000e+000   0.00000000e+000]
     [  0.00000000e+000   0.00000000e+000   0.00000000e+000   0.00000000e+000]
     [  0.00000000e+000   0.00000000e+000   0.00000000e+000   6.95093668e-310]]

```python
# Create an array of elements with the same value
# The full() function arguments are
# 1) Shape: tuple of integers giving the size along each dimension.
# 2) The constant value

y = np.full((1,1), 3)
print(y)
print(y.shape)

print()

y = np.full((2,2), 4)
print(y)
```

    [[3]]
    (1, 1)

    [[4 4]
     [4 4]]

```python
# Create a 1D array of evenly spaced values
# The arange() function arguments are the same as the range() function.
# Shape: tuple of integers giving the size along each dimension.

z = np.arange(5,10)
print(z)

print()

z = np.arange(5, 10, 2)
print(z)
```

    [5 6 7 8 9]

    [5 7 9]

```python
# Create a 1D array of evenly spaced values
# The linspace() function arguments are
# The lower limit of the range of values
# The upper limit of the range of values (inclusive)
# The desired number of equally spaced values

z = np.linspace(-4, 4, 5)
print(z)
```

    [-4. -2.  0.  2.  4.]

```python
# Create an array of all zeros
# The zeros() function argument is the shape.
# Shape: tuple of integers giving the size along each dimension.

a = np.zeros(5)
print(a)

print()

a = np.zeros((2,2))
print(a)
```

    [ 0.  0.  0.  0.  0.]

    [[ 0.  0.]
     [ 0.  0.]]

```python
# Create an array of all ones

b = np.ones(5)
print(b)

print()

b = np.ones((1, 4))
print(b)

```

    [ 1.  1.  1.  1.  1.]

    [[ 1.  1.  1.  1.]]

```python
# Create a constant array
# The second function argument is the constant value

c = np.full(6, 8)
print(c)

print()

c = np.full((2,2,2), 7)
print(c)

```

    [8 8 8 8 8 8]

    [[[7 7]
      [7 7]]

     [[7 7]
      [7 7]]]

## Subpackages
Packages can also have subpackages.

The `numpy` package has a subpackage called `random`.

It contains functions to deal with random variables.

If the `numpy` package is imported with `import numpy as np`, functions in the `random` subpackage can be called using `np.random.function()`.

```python
# Create an array filled with random values in the range 0.0 to 1.0

e = np.random.rand(1)
print(e)
print()

e = np.random.rand(3,2,1)
print(e)
print()

e = np.random.random((2,2))
print(e)
```

    [ 0.13974223]

    [[[ 0.3813017 ]
      [ 0.81091848]]

     [[ 0.3293824 ]
      [ 0.06379465]]

     [[ 0.61971733]
      [ 0.50822455]]]

    [[ 0.73459092  0.79476614]
     [ 0.64264     0.40050299]]

```python
# Create an array filled with random integer values
```

```python
# integer values between 4 and 16
e = np.random.randint(4, 16, size=(4,4))
print(e)

print()

# integer values between 1 and 8
e = np.random.randint(1, 8, size=(2, 2, 2))
print(e)
```

    [[14 12 13  7]
     [ 6  5  4 15]
     [ 5  8  8 12]
     [15  6 13  9]]

    [[[4 4]
      [4 7]]

     [[3 4]
      [7 2]]]

## Indexing into multi-dimensional arrays.

We can index into an array exactly the same way as the other data structures we have studied.

```python
x = np.array([1, 2, 3, 4, 5])

# Select a single element
print(x[4])

# Select elements from 2 to the end
print(x[2:])
```

    5
    [3 4 5]

For an n-dimensional (nD) matrix we need n index values to address an element or range of elements.

Example: The index of a 2D array is specified with two values:
- first the row index
- then the column index.

Note the order in which dimensions are addressed.

```python
# 2 dimensional array

y = np.array([[1, 2, 3],
              [4, 5, 6]])

# Select a single element
print(y[1,2])

# Select elements that are both in rows 1 to the end AND columns 0 to 2
print(y[1:, 0:2])
```

    6
    [[4 5]]

We can address elements by selecting a range with a step:

For example the index:

`z[0, 0:]`

selects every element of row 0 in array, `z`

The index:

`z[0, 0::2]`

selects every *other* element of row 0 in array, `z`

```python
# 2 dimensional array

z = np.zeros((4,8))

# Change every element of row 0
z[0, 0:] = 10

# Change every other element of row 1
z[1, 0::2] = 10

print(z)
```

    [[ 10.  10.  10.  10.  10.  10.  10.  10.]
     [ 10.   0.  10.   0.  10.   0.  10.   0.]
     [  0.   0.   0.   0.   0.   0.   0.   0.]
     [  0.   0.   0.   0.   0.   0.   0.   0.]]

```python
z = np.zeros((4,8))

# Change the last 4 elements of row 2, in negative direction
# You MUST include a step to count in the negative direction
z[2, -1:-5:-1] = 10

# Change every other element of the last 6 elements of row 3
# in negative direction
z[3, -2:-7:-2] = 10

print(z)
```

    [[  0.   0.   0.   0.   0.   0.   0.   0.]
     [  0.   0.   0.   0.   0.   0.   0.   0.]
     [  0.   0.   0.   0.  10.  10.  10.  10.]
     [  0.   0.  10.   0.  10.   0.  10.   0.]]

```python
# 3-dimensional array

c = np.array(
    [[[2, 1, 4],
      [2, 6, 8]],

     [[0, 1, 5],
      [7, 8, 9]]])

print(c[0, 1, 2])

```

    8

Where we want to select all elements in one dimension we can use :

__Exception__: If it is the last element , we can omit it.

```python
print(c[0, 1])

print(c[0, :, 1])
```

    [2 6 8]
    [1 6]

### Boolean array indexing

Recall that we can use *conditional operators* to check the value of a single variable against a condition.

The value returned is a Boolean True or False value.

```python
a = 4
print('a < 2:', a < 2)
print('a > 2:', a > 2)
```

    a < 2: False
    a > 2: True

If we instead use *conditional operators* to check the value of an array against a condition.

The value returned is an *array* of Boolean True or False values.

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

idx = a > 2

print(idx)

```

    [[False False]
     [ True  True]
     [ True  True]]

A particular elements of an array can be are specified by using a boolean array as an index.

Only the values of the array where the boolean array is `True` are selected.

The varaible `idx` can therefore now be used as the index to select all elements greater than 2.

```python
print(a[idx])
```

    [3 4 5 6]

To do the whole process in a single step

```python
print(a[a > 2])
```

    [3 4 5 6]

To apply multiple conditions, use () parentheses to sperate different conditions.

Use `&` for elementwise `and`.

Use `|` for elementwise `or`.

```python
x = np.array([[4, 2, 3, 1],
              [2, 4, 2, 8],
              [2, 3, 3, 27],
              [4, 1, 4, 64]])

# elements of x that are greater then 2 AND less than 10
print(x[(2 < x) & (x < 10)])

# elements of x that are less then 2 OR greater than 10
print(x[(2 < x) & (x < 10)])
```

    [4 3 4 8 3 3 4 4]
    [4 3 4 8 3 3 4 4]

Multiple conditions can also be applied to a subsection of an array.

For example to select elements $>2$ and $<4$ in the first row of `x` only (`x[0]`):

```python
x = np.array([[4, 2, 3, 1],
              [2, 4, 2, 8],
              [2, 3, 3, 27],
              [4, 1, 4, 64]])

print(x[0][(2 < x[0]) & (x[0] < 4)])

```

    [3]

## Iterating over multi-dimensional arrays.
We can iterate over a 1D array in the same way as the data structures we have previously studied.

```python
A = np.array([1, 2, 3, 4, 5])
```

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

    1
    2
    3
    4
    5

To loop through individual elements of a multi-dimensional array, we use a nested loop for each dimension of the array.

```python
B = np.array([[1, 2, 3],
              [4, 5, 6]])

for row in B:
    print("-----")
    for col in row:
        print(col)
```

    -----
    1
    2
    3
    -----
    4
    5
    6

## Manipulating arrays
We can use many of the same operations to manipulate arrays as we use for lists.

However, it is important to note a few subtle differences in how array manipulations behave.

```python
# Length of an array

a = np.array([1, 3, 4, 17, 3, 21, 2, 12])

b = np.array([[1, 3, 4, 17],
              [3, 21, 2, 12]])

print(len(a))
print(len(b))

```

    8
    2

Note the length is the length of the first dimension (e.g. indexing).

```python
# Sort an array

a = np.array([1, 3, 4, 17, 3, 21, 2, 12])

b = np.array([[1, 3, 4, 17],
              [3, 21, 2, 12]])

# The function sorted applies to 1D data structures only
print(sorted(a))
print(sorted(b[1]))

# The method sort() applies to arrays of any size
a.sort()
b.sort()

print(a)
print(b)
```

    [1, 2, 3, 3, 4, 12, 17, 21]
    [2, 3, 12, 21]
    [ 1  2  3  3  4 12 17 21]
    [[ 1  3  4 17]
     [ 2  3 12 21]]

Arrays are *immutable* (unchangeable).

Technically you cannot add or delete items of an array.

However, you can make a *new* array (which may have the same name as the original array), with the values ammended as required:

#### Appending Arrays
Appending connects array-like (integer, list....) value  to the *end* of the original array.

By default, 2D arrays are appended as if joining lists.
The new array is a 1D array

```python
# 2D array
a = np.array([[0], [1], [2]])
print(a)
print()

# 2D array
b = np.array([[3], [4]])
print(b)
print()

# 1D array
c = np.array([3, 4])
print(b)
print()

# integer
d = 1

print(f"original 2D array shapes: a = {a.shape}, b = {b.shape}")
print()

X = np.append(a, b)
print(X)
print(f"new array shape: {a.shape}")
print()

X = np.append(b, d)
print(X)
print(f"new array shape: {a.shape}")
print()

X = np.append(c, d)
print(X)
print(f"new array shape: {a.shape}")
print()
```

    [[0]
     [1]
     [2]]

    [[3]
     [4]]

    [[3]
     [4]]

    original 2D array shapes: a = (3, 1), b = (2, 1)

    [0 1 2 3 4]
    new array shape: (3, 1)

    [3 4 1]
    new array shape: (3, 1)

    [3 4 1]
    new array shape: (3, 1)

The axis on which to append an array can be optionally specified.

e.g. 2D array:
 - 0: columns
 - 1: rows

The arrays must have the same shape, except in the dimension corresponding to the specified axis

```python
# 2D array
a = np.array([[0], [1], [2]])
print(a)
print()

# 2D array
b = np.array([[3], [4]])
print(b)
print()

new2d = np.append(a, b, axis=0)
print(new2d)
print(f"new array shape: {new2d.shape}")
```

    [[0]
     [1]
     [2]]

    [[3]
     [4]]

    [[0]
     [1]
     [2]
     [3]
     [4]]
    new array shape: (5, 1)

For example, in the cell above, if you change `axis=0` to `axis=1`,

you are trying to connect the side of `a` with length=3 to the side of `b` with length=2.

There are dedicated functions to simplify joining or merging arrays.

If you are interested to expeirment further with joiing arrays you can try out the following functions:
 - `np.concatenate()` : Joins a sequence of arrays along an existing axis.
 - `np.vstack()` or `np.r_[]`: Stacks arrays row-wise
 - `np.hstack()` : Stacks arrays horizontally
 - `np.column_stack()` or `np.c_[]` : Stacks arrays column-wise
Refer to last week's seminar for how to inpterpret the function documentation.

It can also be useful to remove individual (single or multiple) elements.

For example, the following expand the locations within the array that you can change beyond the location at the *end* of the array.

#### Adding elements to an array

```python
# Add items to an array
# The insert() function arguments are
# 1) The array to insert to
# 2) The index of the inserted element
# 3) The value of the inserted element

a = ([1, 2, 3])
a = np.insert(a, 1, 4)
print(a)
```

    [1 4 2 3]

Notice that, again, the output is a 1D aray by default

```python
# Add items to an array

b = np.array([[1, 1],
              [2, 2],
              [3, 3]])

print(f"original array shape: {b.shape}")

b = np.insert(b, 1, [4, 4])

print(b)

print(f"new array shape: {b.shape}")
```

    original array shape: (3, 2)
    [1 4 4 1 2 2 3 3]
    new array shape: (8,)

To preserve the multi-dimensional structure of an array, we can specify the axis on which to insert an element or range of elements.

 In the example below, a column is inserted at element 1 of axis 1.

```python
# Add items to an array

b = np.array([[1, 1],
              [2, 2],
              [3, 3]])

b = np.insert(b, 1, [3, 2, 1], axis=1)
print(b)
```

    [[1 3 1]
     [2 2 2]
     [3 1 3]]

Notice what happens when we insert a *single* value on a specified axis

```python
b = np.insert(b, 1, 4, axis=1)
print(b)
```

    [[1 4 3 1]
     [2 4 2 2]
     [3 4 1 3]]

This behaviour is due to a very useful property called *broadcasting*.

We will study the rules governing broadcasting later in this seminar.

#### Deleting items from an array

```python
# Items are deleted from their position in a 1D array by default

z = np.array([1, 3, 4, 5, 6, 7, 8, 9])

z = np.delete(z, 3)
print(z)

z = np.delete(z, [0, 1, 2])
print(z)

```

    [1 3 4 6 7 8 9]
    [6 7 8 9]

```python
# Again, axes to delete can be optionally specified:

z = np.array([[1, 3, 4, 5], [6, 7, 8, 9]])
print(z)
print()

z = np.delete(z, 3, axis=1)
print(z)
print()

z = np.delete(z, [0, 1, 2], axis=1)
print(z)
print()
```

    [[1 3 4 5]
     [6 7 8 9]]

    [[1 3 4]
     [6 7 8]]

    []

#### Changing items in an array

```python
c = np.array([1, 2, 3])
c[1] = 4
print(c)
```

    [1 4 3]

## Introduction to Plotting
It is often useful to represent data visually.

Matplotlib is a library for creating beautiful graphics.

First we will import the plotting part of `matplotlib`, renaming it  `plt`.

Second we will use a command which tells Python to show any graphs inside the Notebook(not in a separate window).

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

### A simple line plot
It is very easy to generate simlpe plots from arrays or other data structures.

A single parameter is automatically plotted against the index of each element.

```python
x = np.array([1,7,2,3,3,4,5])
plt.plot(x)
```

To plot one parameter against another parameter:

```python
x = [1, 2, 5, 7, 12]
y = [3, 2, 0.7, 0.6, 0.1]
plt.plot(x, y)
```

### A scatter plot
To plot a scatter plot of points instead of a connected line add a `formatstring`:

```python
x = [1, 2, 5, 7, 12]
y = [3, 2, 0.7, 0.6, 0.1]
plt.plot(x, y, 'o')
```

## Mathematics with arrays.

Unlike lists, NumPy arrays support common arithmetic operations, such as addition of two arrays.

```python
# To add the elements of two lists we need the Numpy function: add
a = [1, 2, 3]
b = [4, 5, 6]

c = a + b
print(c)

c = np.add(a, b)
print(c)
```

    [1, 2, 3, 4, 5, 6]
    [5 7 9]

To add the elements of two arrays we can just use regular arithmetic operators.

The two arrays should have the same dimensions.

```python
a = np.array([1, 2, 3])
b = np.ones((1,3))

c = a + b
print(c)

```

    [[ 2.  3.  4.]]

Algebraic operations are appled *elementwise* to an array.

This means the function is applied individually to each element in the list.

```python
a = np.array([1.0, 0.2, 1.2])
b = np.array([2.0, 0.1, 2.1])

print(a - b)

print(np.subtract(a, b))
```

    [-1.   0.1 -0.9]
    [-1.   0.1 -0.9]

```python
a = np.array([[1.0, 0.2, 1],
              [1.0, 0.2, 1]])

b = np.array([[2.0, 0.1, 2],
              [1.0, 0.2, 1]])

# Elementwise multiplication of a and b
print(a * b)
print()
print(np.multiply(a, b))

print(),print()

# Elementwise division of a and b
print(a / b)
print()
print(np.divide(a, b))
```

    [[ 2.    0.02  2.  ]
     [ 1.    0.04  1.  ]]

    [[ 2.    0.02  2.  ]
     [ 1.    0.04  1.  ]]

    [[ 0.5  2.   0.5]
     [ 1.   1.   1. ]]

    [[ 0.5  2.   0.5]
     [ 1.   1.   1. ]]

Apply a mathematical function to a range of values.

Example: Evaluate $y = x^4 - 16$ between $x=-10$ and $x=10$

```python
x = np.linspace(-10, 10, 100)
y = x**4 - 16
plt.plot(x, y)

```

```python
x = np.linspace(-2*np.pi, 2*np.pi, 100)
y = np.cos(x)
plt.plot(x, y)
```

## Mathematics with Vectors (1D arrays)
1D arrays are a useful ways to represent vectors.

A vector is a quantity with a direction and a magnitude.

Let's look at a  previous example for computing the dot product of two vectors.

The dot product of two $n$-length-vectors:

 $ \mathbf{A} = [A_1, A_2, ... A_n]$

 $ \mathbf{B} = [B_1, B_2, ... B_n]$

\begin{align}
\mathbf{A} \cdot \mathbf{B} = \sum_{i=1}^n A_i B_i.
\end{align}

We learnt to solve this very easily using a Python `for` loop.

With each iteration of the loop we increase the value of `dot_product` (initial value = 0.0) by the product of `a` and `b`.

```python
A = [1.0, 3.0, -5.0]
B = [4.0, -2.0, -1.0]

# Create a variable called dot_product with value, 0.
dot_product = 0.0

for a, b in zip(A, B):
    dot_product += a * b

print(dot_product)
```

Numpy makes solving the dot product even easier.

We can use the Numpy function `dot()`.

`dot()` can also takes lists as inputs.

We can also use `@`

`@` does not work on lists or in Python 2

Example: Computing the dot product of two vectors.

```python
import numpy as np
A = [9, 2, 7]
B = [4, 8, 10]
print(np.dot(A,B))

A = np.array([9, 2, 7])
B = np.array([4, 8, 10])
print(np.dot(A,B))
print(A @ B)
```

    122
    122
    122

__Try it yourself__

In the cell below write a user-defined function that takes two lists and returns the dot product.

*Hint: use the code from Seminar 4: Data Structures (shown above).*

Use the magic function `%timeit` to compare the speed of the for loop with the Numpy `dot()` function for solving the dot product.

```python
# Write a function for the dot product of two vectors expressed as lists
# Compare the speed of your function to the Numpy function
```

## Mathematics with Matrices (2D arrays)
If you have previously studied matrices, the operations in this section will be familiar.

If you have not yet studied matrices, you may want to refer back to this section once matrices have been covered in your mathematics classes.

Matrix operations will not be  included in the exam.

2D arrays are a convenient way to represents matrices.

For example, the  2x3 matrix

$$
A =
\begin{bmatrix}
3 & 5 & 7\\
2 & 4 & 6
\end{bmatrix}
$$

can be represented as a 2D array.

```python
A = np.array([[3, 5, 7],
              [2, 4, 6]])
print(A)
```

    [[3 5 7]
     [2 4 6]]

We can use shape to return the matrix dimensions.

```python
print(A.shape)
print(f"Number of rows is {A.shape[0]}, number of columns is {A.shape[1]}")
print(f"A is an {A.shape[0]} by {A.shape[1]} matrix")
```

    (2, 3)
    Number of rows is 2, number of columns is 3
    A is an 2 by 3 matrix

#### Matrix multiplication rules.
If the number of __columns in A__

is the same as number of __rows in B__,

we can find the matrix product of $\mathbf{A}$ and $\mathbf{B}$.

 $\mathbf{C} = \mathbf{A} \cdot \mathbf{B}$

For example:

$\mathbf{A}$ has 3 rows and __3 columns__

$\mathbf{B}$ has __3 rows__ and 1 column

($\mathbf{B}$ is a vector represented as a matrix)

\begin{equation*}
\underbrace{
\begin{bmatrix}
1 & 2 & 3 \\
4 & 5 & 6 \\
7 & 8 & 9 \\
\end{bmatrix}
}_{\mathbf{A} \text{ 3 rows} \text{ 3 columns}}
\cdot
\underbrace{
\begin{bmatrix}
10 \\
20 \\
30 \\
\end{bmatrix}
}_{\mathbf{B} \text{  3 rows} \text{  1 column}}
\end{equation*}

So we can multiply them...

In matrix $\mathbf{C}$, the element in __row $i$__, __column $j$__

is equal to the dot product of the $i$th __row__ of $\mathbf{A}$, $j$th __column__ of $\mathbf{B}$.m

\begin{equation*}
\underbrace{
\begin{bmatrix}
\color{red}1 & \color{red}2 & \color{red}3 \\
4 & 5 & 6 \\
7 & 8 & 9 \\
\end{bmatrix}
}_{\mathbf{A} \text{ 3 rows} \text{ 3 columns}}
\cdot
\underbrace{
\begin{bmatrix}
\color{red}{10} \\
\color{red}{20} \\
\color{red}{30} \\
\end{bmatrix}
}_{\mathbf{B} \text{  3 rows} \text{  1 column}}
=\underbrace{
\begin{bmatrix}
\color{red}{1 \cdot 10 \quad + \quad 2 \cdot 20 \quad + \quad 3 \cdot 30} \\
4 \cdot 10 \quad + \quad 5 \cdot 20 \quad + \quad 6 \cdot 30 \\
7 \cdot 10 \quad + \quad 8 \cdot 20 \quad + \quad 9 \cdot 30 \\
\end{bmatrix}
}_{\mathbf{C} \text{  3 rows} \text{  1 column}}
=\underbrace{
\begin{bmatrix}
\color{red}{140} \\
320 \\
500 \\
\end{bmatrix}
}_{\mathbf{C} \text{  3 rows} \text{  1 column1}}
\end{equation*}

\begin{equation*}
\underbrace{
\begin{bmatrix}
1 & 2 & 3 \\
\color{red}4 & \color{red}5 & \color{red}6 \\
7 & 8 & 9 \\
\end{bmatrix}
}_{\mathbf{A} \text{ 3 rows} \text{ 3 columns}}
\cdot
\underbrace{
\begin{bmatrix}
\color{red}{10} \\
\color{red}{20} \\
\color{red}{30} \\
\end{bmatrix}
}_{\mathbf{B} \text{  3 rows} \text{  1 column}}
=\underbrace{
\begin{bmatrix}
1 \cdot 10 \quad + \quad 2 \cdot 20 \quad + \quad 3 \cdot 30 \\
\color{red}{4 \cdot 10 \quad + \quad 5 \cdot 20 \quad + \quad 6 \cdot 30} \\
7 \cdot 10 \quad + \quad 8 \cdot 20 \quad + \quad 9 \cdot 30 \\
\end{bmatrix}
}_{\mathbf{C} \text{  3 rows} \text{  1 column}}
=\underbrace{
\begin{bmatrix}
140 \\
\color{red}{320} \\
500 \\
\end{bmatrix}
}_{\mathbf{C} \text{  3 rows} \text{  1 column1}}
\end{equation*}

\begin{equation*}
\underbrace{
\begin{bmatrix}
1 & 2 & 3 \\
4 & 5 & 6 \\
\color{red}7 & \color{red}8 & \color{red}9 \\
\end{bmatrix}
}_{\mathbf{A} \text{ 3 rows} \text{ 3 columns}}
\cdot
\underbrace{
\begin{bmatrix}
\color{red}{10} \\
\color{red}{20} \\
\color{red}{30} \\
\end{bmatrix}
}_{\mathbf{B} \text{  3 rows} \text{  1 column}}
=\underbrace{
\begin{bmatrix}
1 \cdot 10 \quad + \quad 2 \cdot 20 \quad + \quad 3 \cdot 30 \\
4 \cdot 10 \quad + \quad 5 \cdot 20 \quad + \quad 6 \cdot 30 \\
\color{red}{7 \cdot 10 \quad + \quad 8 \cdot 20 \quad + \quad 9 \cdot 30} \\
\end{bmatrix}
}_{\mathbf{C} \text{  3 rows} \text{  1 column}}
=\underbrace{
\begin{bmatrix}
140 \\
320 \\
\color{red}{500} \\
\end{bmatrix}
}_{\mathbf{C} \text{  3 rows} \text{  1 column1}}
\end{equation*}

\begin{equation*}
\underbrace{
\begin{bmatrix}
1 & 2 & 3 \\
4 & 5 & 6 \\
7 & 8 & 9 \\
\end{bmatrix}
}_{\mathbf{A} \text{ 3 rows} \text{ 3 columns}}
\cdot
\underbrace{
\begin{bmatrix}
10 \\
20 \\
30 \\
\end{bmatrix}
}_{\mathbf{B} \text{  3 rows} \text{  1 column}}
=\underbrace{
\begin{bmatrix}
1 \cdot 10 \quad + \quad 2 \cdot 20 \quad + \quad 3 \cdot 30 \\
4 \cdot 10 \quad + \quad 5 \cdot 20 \quad + \quad 6 \cdot 30 \\
7 \cdot 10 \quad + \quad 8 \cdot 20 \quad + \quad 9 \cdot 30 \\
\end{bmatrix}
}_{\mathbf{C} \text{  3 rows} \text{  1 column}}
=\underbrace{
\begin{bmatrix}
140 \\
320 \\
500 \\
\end{bmatrix}
}_{\mathbf{C} \text{  3 rows} \text{  1 column1}}
\end{equation*}

Matrix $\mathbf{C}$ therefore has:
- the same number of __rows__ as $\mathbf{A}$,
- the same number of __columns__ as $\mathbf{B}$.

```python
# In the equation above, vector B must be represented as a column vector
A = np.array([[1, 2, 3],
              [4, 5, 6],
              [7, 8, 9]])

# In Python, 1D arrays are ALWAYS represented horizontally
# This does not define the array as a row vector
B = np.array([10, 20, 30])

# For example, C is represented horizontally
C = np.dot(A,B)
print(C)

```

    [140 320 500]

As an example, if $\mathbf{B}$ were a row vector:

\begin{equation*}
\underbrace{
\begin{bmatrix}
1 & 2 & 3 \\
4 & 5 & 6 \\
7 & 8 & 9 \\
\end{bmatrix}
}_{\mathbf{A} \text{ 3 rows} \text{ 3 columns}}
\cdot
\underbrace{
\begin{bmatrix}
10 & 20 & 30 \\
\end{bmatrix}
}_{\mathbf{B} \text{  1 row} \text{  3 columns}}
\end{equation*}

We *cannot* find the dot product $\mathbf{B}\cdot\mathbf{A}$.

The number of columns in $\mathbf{A}$ __is not__ the same as number of rows in $\mathbf{B}$.

We can swap the order of $\mathbf{A}$ and $\mathbf{B}$.

The multiplication is now possible.

However, the outcome is different.

\begin{equation*}
\underbrace{
\begin{bmatrix}
10 & 20 & 30 \\
\end{bmatrix}
}_{\mathbf{B} \text{ 1 row} \text{ 3 columns}}
\cdot
\underbrace{
\begin{bmatrix}
\color{red}1 & \color{blue}2 & \color{green}3 \\
\color{red}4 & \color{blue}5 & \color{green}6 \\
\color{red}7 & \color{blue}8 & \color{green}9 \\
\end{bmatrix}
}_{\mathbf{A} \text{  3 rows} \text{  3 columns}}
=\underbrace{
\begin{bmatrix}
\color{red}{10 \cdot 1 + 20 \cdot 4 + 30 \cdot 7} &
\color{blue}{4 \cdot 10 + 5 \cdot 20 + 6 \cdot 30} &
\color{green}{7 \cdot 10 + 8 \cdot 20 + 9 \cdot 30} \\
\end{bmatrix}
}_{\mathbf{C} \text{  1 row} \text{  3 columns}}
=\underbrace{
\begin{bmatrix}
\color{red}{140} &
\color{blue}{320} &
\color{green}{500} \\
\end{bmatrix}
}_{\mathbf{C} \text{  3 rows} \text{  1 column1}}
\end{equation*}

In Python, normal matrix multiplication rules apply to 2D arrays.

This holds even if the length of one of the dimensions of the 2D array is equal to 1.

```python
A = np.array([[1, 2, 3],
              [4, 5, 6],
              [7, 8, 9]])

# 2D array
X = np.array([[10, 20, 30]])

# 2D array
Y = np.array([[10],
              [20],
              [30]])

print(np.dot(X, A)) #, print(np.dot(A, X))
print(np.dot(A, Y)) #, print(np.dot(Y, A))
```

    [[300 360 420]]
    [[140]
     [320]
     [500]]

However, the orientation with which 1D arrays are shown (always horizontal) does not impact their allowbale placement in an expression.

Python will automatially treat the 1D as a column where appropriate.

```python
A = np.array([[1, 2, 3],
              [4, 5, 6],
              [7, 8, 9]])

# 1D array
Z = np.array([10, 20, 30])

print(np.dot(Z, A))
print(np.dot(A, Z))
```

    [300 360 420]
    [140 320 500]

### Useful Matrix Operations
#### Inverse of a square matrix

```python
A = np.array([[1,2],
              [3, 4]])

Ainv = np.linalg.inv(A)

print(f"A = \n {A}")
print(f"Inverse of A = \n {Ainv}")
```

    A =
     [[1 2]
     [3 4]]
    Inverse of A =
     [[-2.   1. ]
     [ 1.5 -0.5]]

#### Determinant of a square matrix

```python
A = np.array([[1,2],
              [3, 4]])

Adet = np.linalg.det(A)

print(f"A = \n {A}")
print(f"Determinant of A = {round(Adet, 2)}")
```

    A =
     [[1 2]
     [3 4]]
    Determinant of A = -2.0

#### Transpose of a matrix

- The columns of the transpose matrix are the rows of the original matrix.
- The rows of the transopse matrix are the columns of the original matrix.

```python
a = np.zeros((2,4))
print(a)
print()

print(a.T)
print()

#or

print(np.transpose(a))
```

    [[ 0.  0.  0.  0.]
     [ 0.  0.  0.  0.]]

    [[ 0.  0.]
     [ 0.  0.]
     [ 0.  0.]
     [ 0.  0.]]

    [[ 0.  0.]
     [ 0.  0.]
     [ 0.  0.]
     [ 0.  0.]]

#### Generate Identity Matrix

```python
I = np.eye(2)
print(I)

print()

I = np.eye(4)
print(I)
```

    [[ 1.  0.]
     [ 0.  1.]]

    [[ 1.  0.  0.  0.]
     [ 0.  1.  0.  0.]
     [ 0.  0.  1.  0.]
     [ 0.  0.  0.  1.]]

### Useful Matrix Applications: Solving Simultaneous Equations

A mathematical problem that arises in a wide range of engineering applications; the simultaneous solution to a set of n linear equations with n unknown values.

Matrices are a widely used and popular technique for solving this type of problem.

Solving problems of this __computationally__ is particularly useful where:
- we need to solve a large number of equations simultaneously e.g.
 - the force distribution in a truss structure
 - the electrical current flowing in and out of each node in a circuit.

- we need to solve the equations over multiple steps
e.g. observing the progression of a system over time.

### Systems of Equations Example: An Electrical Circuit

In electrical systems, network theory is the study of how to solve circuit problems.

The voltage and current in a closed circuit obey two simple rules;
- Kirchoff's current law (first law)
- Kirchhoff's voltage law (second law)

#### Kichoff's Current Law

The algebraic sum of ALL the currents entering and leaving a junction must be equal to zero.

$\sum i_{in} - \sum i_{out} = 0 $

#### Kirchhoff's Voltage Law
For a closed loop series path the algebraic sum of all the *voltages* and *voltage drops* around any closed loop in a circuit is equal to zero.

$\sum E - \sum V = 0 $

#### Electrical Elements Obey Ohm's Law
The current through a conductor (I, units amps) is the voltage measured across the conductor (V, units volts) divided by the resistance (R, units Ohms).

$$V = IR$$

A Wheatstone bridge circuit is a very widely used circuit in scientific equipment.

It is found in strain gauges, heart rate monitors, wind tunnels....etc.

Fundamentally, it is used to measure resistance.

A wheatstone bridge circuit has mutiple loops and nodes.

It has one voltage source, $E$.

It has five resisters $R_1, R_2, R_3, R_4, R_5, R_6$.

*Let's say we want to find the current in each branch $i_{1-6}$...*

The circuit has:
 - 4 loops (A-D)
 - a nodes (a-d)

We can use Kirchhoff's laws and Ohm's law with $E$ and $R_{1-6}$.

$\sum i_{in} - \sum i_{out} = 0 $

$\sum E - \sum V = 0 $

$V = IR$

6 unknown values ($i_{1-6}$), 6 simultaneous equations.

__Loop A:__   $R_1 i_1 + R_2 i_2 + R_5 i_5 = E$

__Loop B:__   $R_4 i_4 - R_5 i_5 + R_6 i_6 = 0$

__Loop C:__   $R_2 i_2 - R_3 i_3 + R_4 i_4 = 0$

__Node a:__   $i_2 - i_4 - i_5 = 0$

__Node b:__   $i_1 - i_2 - i_3 = 0$

__Node c:__   $i_5 + i_6 - i_1 = 0$

We can express these equations as matrices:

\begin{equation*}
\underbrace{
\begin{bmatrix}
R_1 & R_2 & 0    & 0   &  R_5 & 0   \\
0   & 0   & 0    & R_4 & -R_5 & R_6 \\
0   & R_2 & -R_3 & R_4 &  0   & 0   \\
0   & 1   & 0    & -1  &  -1  & 0   \\
1   & -1  & -1   & 0   &  0   & 0   \\
-1  & 0   & 0    & 0   &  1   & 1   \\
\end{bmatrix}
}_{\mathbf{R}}
\cdot
\underbrace{
\begin{bmatrix}
i_1 \\
i_2 \\
i_3 \\
i_4 \\
i_5 \\
i_6 \\
\end{bmatrix}
}_{\mathbf{I}}
=\underbrace{
\begin{bmatrix}
E \\
0 \\
0 \\
0 \\
0 \\
0 \\
\end{bmatrix}
}_{\mathbf{V}}
\end{equation*}

```python
# Let's assign some values to our known variables:
E = 5

R1, R2, R3, R4, R5, R6 = 300, 500, 1000, 500, 1000, 500

R = np.array([[R1, R2,  0,   0,   R5, 0 ],
              [0,  0,   0,   R4, -R5, R6],
              [0,  R2, -R3,  R4,  0,  0 ],
              [0,  1,   0,  -1,  -1,  0 ],
              [1, -1,  -1,   0,   0,  0 ],
              [-1, 0,   0,   0,   1,  1]])

V = np.array([E, 0, 0, 0, 0, 0])
```

We can find $\mathbf{I}$ by matrix multiplication:
\begin{align*}
\mathbf{R}\cdot\mathbf{I}&=\mathbf{V}\\
\end{align*}

but we have to rearrange the equation.

\begin{align*}
\mathbf{I}  &=\frac{\mathbf{V}}{\mathbf{R}}\\
            &=\mathbf{R^{-1}}\cdot\mathbf{V}\\
\end{align*}

```python
R_inv = np.linalg.inv(R)

I = np.dot(R_inv, V)

print(I)
```

    [ 0.005  0.003  0.002  0.001  0.002  0.003]

Numpy conveniently solves linear matrix equations of the form `ax = b` so we don't need to rearrange.

```python
I = np.linalg.solve(R, V)
print(I)
```

    [ 0.005  0.003  0.002  0.001  0.002  0.003]

You can check your answer by confirming that the dot product, $\mathbf{R}\cdot\mathbf{I}$ is equal to $\mathbf{V}$.

```python
print(R @ I)
print(np.around(R @ I, 2))
```

    [  5.00000000e+00   0.00000000e+00   0.00000000e+00   0.00000000e+00
       3.46944695e-18  -8.67361738e-19]
    [ 5.  0.  0.  0.  0. -0.]

This process is a fast and simple way to solve simultaneous linear equations, for example when compared to eilimination methods.

### Systems of Equations Example: Rates of Change
A well known problem concerns an imaginary country with three cities, A, B and C.

At the end of each year, a fraction, $n$ of the people must leave each city.

Half of the people leaving a city move to one of the two options, and half to the other.

This gives us a system of simultaneous equations

\begin{align*}
\Delta A = \frac{Bn}{2} + \frac{Cn}{2} - An \\
\Delta B = \frac{An}{2} + \frac{Cn}{2} - Bn\\
\Delta C = \frac{An}{2} + \frac{Bn}{2} -Cn \\
\end{align*}

Matrices are a convenient way to represent this problem.

\begin{equation*}
\begin{bmatrix}
\Delta A \\
\Delta B \\
\Delta C \\
\end{bmatrix}
=
\underbrace{
\begin{bmatrix}
-n & \frac{n}{2} & \frac{n}{2} \\
\frac{n}{2} & -n & \frac{n}{2} \\
\frac{n}{2} & \frac{n}{2} & -n \\
\end{bmatrix}
}_{\mathbf{migration}}
\cdot
\underbrace{
\begin{bmatrix}
A \\
B \\
C \\
\end{bmatrix}
}_{\mathbf{population}}
\end{equation*}

Assigning values to $\mathbf{M}$ and $\mathbf{P}$.

\begin{equation*}
\begin{bmatrix}
\Delta a \\
\Delta b \\
\Delta c \\
\end{bmatrix}
=
\underbrace{
\begin{bmatrix}
-0.2 & 0.1 & 0.1 \\
0.1 & -0.2 & 0.1 \\
0.1 & 0.1 & -0.2 \\
\end{bmatrix}
}_{\mathbf{M}}
\cdot
\underbrace{
\begin{bmatrix}
A \\
B \\
C \\
\end{bmatrix}
}_{\mathbf{P}}
\end{equation*}

By solving the problem computationally, we can quickly solve multiple iterations.

In the example below there are two functions.

The first, `pop_change` computes the rate of change of the population.

The second, `pop_plot`:
- calls `pop_change` inside of it
- generates a line plot showing the population change over time

```python
# Initial population
P0 = np.array([190.0, 500.0, 30.0])

# Matrix migration or M
migration = np.full((3, 3), 0.01)
np.fill_diagonal(migration, -0.02)

def pop_change(P, M=migration):
    "Computes the change in population of three cities"
    return np.dot(P, M)

```

```python
def pop_plot(iterations, population):
    "Plots the population over the number of iterations specified"

    # make 2D array to append values to a column for each city
    record_pop = population.reshape((1,3))

    for i in range(iterations):

        # Increment population size
        population = population + pop_change(population)

        # Save the population data
        record_pop = (np.append(record_pop,
                                population.reshape((1,3)), # make 2D array
                                axis = 0))

    plt.plot(record_pop[:,0] , label="A")
    plt.plot(record_pop[:,1] , label="B")
    plt.plot(record_pop[:,2] , label="C")
    plt.xlabel("Years")
    plt.ylabel("Population")
    plt.legend();
    #plt.plot(record_pop)
    #print(record_pop)
```

```python
# Plot the change in the populations of the cities over time
pop_plot(iterations=120, population=P0)

print(f"total population = {np.sum(P0)}")
```

In the example above the array `migration` or `M` just describes the change in a quantity; the population of each city.

When a change in quantity is measured against a change in __time__ it describes the rate of change or flow of that quantity.

Rates are expressed as differential equations: the change in one parameter in relation to a change in another parameter. e.g.
 - velocity = change in position / change in time
 - acceleration = change in velocity / change in time

The Scipy library has easy to use tools for estimating the solutions to systems of (or single) differential equations like the example we just looked at.

The function `odeint` solves first order differential equations.

```python
from scipy.integrate import odeint
```

The function `odeint` takes a user-defined function as an argument.

This input function should simply define the derivative (e.g. rate of change) you want to solve.

For example, a function to calculate the derivative:

$\frac{dx}{dt} = t - x$

with initial value $x(0) = 1$

we can use:

```python
def dx_dt(x,t):
    return t - x
```

Two additional arguments are needed:
- the value(s) of t at which to evaluate x, starting with the initial value
- the initial value
- (if `dx_dt` takes any other argumemts they are entered as a tuple as the third argument)

```python
ts = np.linspace(0,5,100) # the value(s) of t at which to evaluate x
x0 = 1.0                  # the initial value

# odeint returns x at each value of t
xs = odeint(dx_dt, x0, ts)
```

The function `odeint` returns an nx1 array (2D column).

To plot the output we have to "flatten" it to a 1D array.

```python
# The function
#print(xs)
xs = np.array(xs).flatten()
#print(xs)

# plot the function
plt.xlabel("x")
plt.ylabel("t")
plt.plot(ts, xs);
```

We can use `odeint` to solve mulitple ordinary differential equations simultaneously (systems of ODEs).

In this case the value of each variable at the initial position should be input as a single data structure.

The function should output the rate of change of each variable as a single list.

For example we can use `odeint` to solve the population problem from earlier.

The function dP_dt is exactly the same as the function `change_pop` that we used  earlier. The only difference is that `dP_dt` returns the individual rates of chnage of each population as a list.

```python
migration = np.full((3, 3), 0.01)
np.fill_diagonal(migration, -0.02)

# Initial population
P0 = np.array([190.0, 500.0, 30.0])

# Time steps to evaluate
ts = np.arange(0, 150)

def dP_dt(P, t, M):

    dP_dt = np.dot(P, M)

    return [dP_dt[0],
            dP_dt[1],
            dP_dt[2]]

Ps = odeint(dP_dt, P0, ts, args=(migration,))

plt.plot(ts, Ps[:,0], label="A")
plt.plot(ts, Ps[:,1], label="B")
plt.plot(ts, Ps[:,2], label="C")
plt.xlabel("Time")
plt.ylabel("Population")
plt.legend();
```

For comparison, here is the original solution, generated using the user-defined `pop_plot`function.

```python
pop_init = np.array([190.0, 500.0, 30.0])
pop_plot(iterations=150, population=pop_init)
```

## Curve Fitting
When using functions to produce mathematical models of practical problems we often need to 'fit' the function to experimental data.

Python has several tools to find the curve which best describes the data set.

An example of this might be the data read from a sensor during an experiment.

In the experiment shown below, the bend in the arm is used to control the LED light.

The sensor used to detect the bend in the arm is made from conductive paint.

The resistance of the paint changes with the bend in the arm.

This change in resistance is detected as a voltage.

A microcontroller adjusts the light in repsonse to the voltage.

The circuit on which the sensor is based is very similar to the Wheatsone bridge circuit.

When calibrating a sensor like this, we need to find a relationship between the bend in the arm and the change in resisitance.

After that the voltage alone can be used to 'read' the bend in the arm.

This section will cover some useful Numpy features for curve fitting.

### Root mean square error (RMSE)
One way to quantify the fit between data and a model is to compute the RMSE.

__Error/residual:__ the difference between the observed value $y_i$ at $x$ and the modeled value $a(x_i)$

$$
\varepsilon_i = a(x_i) - y_i
$$

For $N$ data points, the *sum* of the squared errors is

$$S = \sum_{i=1}^{N}\varepsilon_i^2$$

The RMSE $E$ is computed as

$$E=\sqrt{\frac{1}{N}S}=\sqrt{\frac{1}{N}\sum{\varepsilon_i^2}}$$

The RMSE can be computed for any model and any data set  as an indicator of the "goodness of fit".

We will use it to assess the fit of the curves that we generate using Numpy.

The optimisation approach that involves minimising the sum of the squared errors (minimising the RMSE) is often referred to as a *least squares* approach.

```python
# Example data set
x = [0.000000000000000000e+00, 1.052631578947368363e+00, 2.105263157894736725e+00, 3.157894736842105310e+00,4.210526315789473450e+00, 5.263157894736841591e+00, 6.315789473684210620e+00,7.368421052631578760e+00,8.421052631578946901e+00,9.473684210526315042e+00,1.052631578947368318e+01,1.157894736842105132e+01,1.263157894736842124e+01,1.368421052631578938e+01,1.473684210526315752e+01, 1.578947368421052566e+01,1.684210526315789380e+01,1.789473684210526372e+01,1.894736842105263008e+01,2.000000000000000000e+01]
y = [7.445192947240600745e+01, 4.834835792411828947e+01, 6.873305436340778840e+01, 5.979576407972768948e+01,6.404530772390434379e+01,6.090548420541189500e+01, 7.157546008677115879e+01, 8.620253336570679892e+01, 1.138154622045899913e+02, 8.493639813028174501e+01, 9.783457330550828601e+01, 1.082064229481453594e+02, 1.063876210674365979e+02, 1.001971993955305038e+02, 1.061496321788094832e+02, 1.279575585921491836e+02, 1.556956405962417875e+02, 1.584164804859289859e+02, 1.753888794716459358e+02, 1.980941276403034124e+02]
x = np.array(x)
y = np.array(y)
```

### Fitting a Linear Function
A straight line can be fitted through a data set using the `linregress` function from `scipy.stats`.

```python
from scipy.stats import linregress
```

Example: `linregress(x, y)` returns the components of an equation of the form:
$$
y = mx + c
$$

It returns five values.

The first two are the gradient $m$ and y-intercept $c$ of the line.

The last three are statistical parameters that are outside of the scope of this class.

```python
m, c, r_value, p_value, std_err = linregress(x, y)

yfit = float(m) * x + c

# plot a scatter plot by setting 'o' as the marker
plt.plot(x, y, 'o', label='experiment data')

# plot the fitted linear function
plt.plot(x, yfit, label='fit');
```

The "goodness of fit" than then be assessed by calculating the RMSE.

Error, $\varepsilon_i = a(x_i) - y_i$

RMSE, $E=\sqrt{\frac{1}{N}\sum{\varepsilon_i^2}}$

```python
# error
e = (yfit - y)

# RMSE
rmse = np.sqrt(np.sum(e**2)/ len(y))

```

#### Adding a title
We can then a title.

```python
plt.plot(x, y, 'o', label='experiment data')
plt.plot(x, yfit, label='fit')

# print the RMSE as the title
plt.title('RMSE: '+str(rmse))
```

### Fitting a polynomial function
A polynomial curve can be fitted through a data set using the `polyfit` function from `numpy`.

The function arguments are:
- the two variables to analyse
- the order of the polynomial

The function returns:

the coefficients of each term of the polynomial.
 e.g. if we request a polynomial of order 2, we will have terms for, - x^2, x and x^0 and we should expect 3 return arguments.

```python
c, d, e = np.polyfit(x, y, 2)

f, g, h, i = np.polyfit(x, y, 3)

j, k, l, m, n = np.polyfit(x, y, 4)
```

The function `poly1D` can then be used to find the solution to y for all x, using the polynomial generated.

```python
yfit2 = np.poly1d([c, d, e])(x)
yfit3 = np.poly1d([f, g, h, i])(x)
yfit4 = np.poly1d([j, k, l, m, n])(x)

# Which produces the same output as
yfit2 = c*x**2 + d*x + e
yfit3 = f*x**3 + g*x**2 + h*x + i
yfit4 = j*x**4 + k*x**3 + l*x**2 + m*x + n
```

The RMSE can be found as before.

__NOTE__ It is neater to write a function to do this.

```python
# error
rmse2 = round(np.sqrt(np.sum((yfit2 - y)**2)/ len(y)), 2)
rmse3 = round(np.sqrt(np.sum((yfit3 - y)**2)/ len(y)), 2)
rmse4 = round(np.sqrt(np.sum((yfit4 - y)**2)/ len(y)), 2)
```

#### Adding a legend
We can add a legend to label each line on the graph.

Add a `label` when creating each plot.

Finally add the `legend` to the plot.

The argument `loc='best'` chooses the 'best' location for the legend to avoid obscuring your plotted data.

```python
# plot a scatter plot by setting 'o' as the marker
plt.plot(x, y, 'o', label='experiment data')

# plot the fitted 2nd order function
plt.plot(x, yfit2, label= f'2nd order RMSE{rmse2}')

# plot the fitted 3rd order function
plt.plot(x, yfit3, label= f'3rd order RMSE{rmse3}')

# plot the fitted 4th order function
plt.plot(x, yfit4, label= f'4th order RMSE{rmse4}')

# add a legend
plt.legend(loc='best')
```

### Fitting an arbitrary function
There are many Python functions available for curve fitting.

The function `curve_fit` can be used to fit a user-defined fuction.

We first define a function that we want to fit.

__Example 1__
```Python
def linear(x, a, b):
    return a * x + b
```
>`curve_fit(func, x, y)` will find values of `a` and `b` that fits the curve $y = ax + b$ to our data.

__Example 2__
```Python
def exponential(x, a, b):
    y = a * np.exp(b*x)
    return y
```
>`curve_fit(func, x, y)` will find values of `a` and `b` that fits the curve $y = a e^{bx}$ to our data.

`curve_fit` then returns two arguments:
 - an array with the optimised parameters
 - the covariance of the parameters (a statistical measure of accuracy)

```python
from scipy.optimize import curve_fit

def exponential(x, a, b):
    y = a * np.exp(b*x)
    return y

opt, cov = curve_fit(exponential, x, y)

# value of x for all x investigated
yfit = exponential(x, *opt)

plt.plot(x, y, 'o', label='experiment data')
plt.plot(x, yfit, 'r', label='fit')
plt.legend(loc='best')

rmse = np.sqrt(np.sum((yfit - y) ** 2) / len(y))
plt.title('RMSE: '+str(rmse));

print(f"y = {round(opt[0],2)} * e**({round(opt[1],2)}*x)")
```

## Vectorising Functions

Numpy functions applied to a single array, will be performed on each element in the array.

The function takes an array of values as an input argument.

```python
print(np.sqrt(a))
print(a ** (1/2))
```

    [[ 0.  0.  0.  0.]
     [ 0.  0.  0.  0.]]
    [[ 0.  0.  0.  0.]
     [ 0.  0.  0.  0.]]

For example, we can apply trigonometric functions, elementwise, to arrays, lists and tuples.

```python
x = np.array([0.0, np.pi/2, np.pi, 3*np.pi/2])
y = [0.0, np.pi/2, np.pi, 3*np.pi/2]
z = (0.0, np.pi/2, np.pi, 3*np.pi/2)

print(np.sin(x))
print(np.cos(y))
print(np.tan(z))
```

    [  0.00000000e+00   1.00000000e+00   1.22464680e-16  -1.00000000e+00]
    [  1.00000000e+00   6.12323400e-17  -1.00000000e+00  -1.83697020e-16]
    [  0.00000000e+00   1.63312394e+16  -1.22464680e-16   5.44374645e+15]

An array of values does not work as an input for all functions.

```python
def func(x):
    if x < 0:
        f = 2 * x
    else:
        f = 3 * x
    return f

x = np.array([2, -2])
# y = func(x) # Run this line after removing the # to see the error generated
```

This doesn't work because Python doesn't know what to do with the line

`if x < 0`

when `x` contains many values.

For some values of `x` the `if` statement may be `True`, for others it may be `False`.

A simple way around this problem is to vectorise the function.

We create a new function that is a *vectorized* form of the original function.

The new function and can be called with an array as an argument.

```python
funcvec = np.vectorize(func)

print(funcvec(x))
```

    [ 6 -4]

### Functions as function arguments
Functions that receive vectorised arguments are automatically applied to all elements of array inputs.

This is better exlpained with an example.

Recall the function `is_positive` from Seminar 5: Functions

```python
def is_positive(f, a):
    "Checks if the function value f(x) is positive"
    return f(a) > 0

# Apply is_positive to a non-vectorised function
print(is_positive(func, -3))
print()

# Apply is_positive to a vectorised function
print(x)
is_positive(funcvec, x)
```

    False

    [ 2 -2]

    array([ True, False], dtype=bool)

## Broadcasting

Another source of incompatibility that you are likely to encounter is in trying to use arrays with different shapes for arithmetic operations.

For example, you have one array that larger and another array that is smaller.

You may want to use the smaller array multiple times to perform an operation (such as a sum, multiplication, etc.) on the larger array.

This is achieved using the broadcasting mechanism.

The arrays can be broadcast together if all dimensions of the arrays are *compatible*

##### Dimensions are compatible when they are equal.

Consider the example below. `x` and `y` are the same shape, so we can addd them.

```python
x = np.ones((3,4))
print(x.shape)

y = np.full((3,4), 4)
print(y.shape)

# Add `x` and `y`
x + y
```

    (3, 4)
    (3, 4)

    array([[ 5.,  5.,  5.,  5.],
           [ 5.,  5.,  5.,  5.],
           [ 5.,  5.,  5.,  5.]])

##### Dimensions are compatible when the length of at least one of them is equal to 1.

```python
# 1 x 3 array
a = np.arange(1,4)

# integer
b = 2

# 1 x 3 array
result = a * b

print(a)
print()
print(b)
print()
print(result)
```

    [1 2 3]

    2

    [2 4 6]

In the dimension where `b` has size 1 and `a` has a size greater than 1 (i.e. 3), `b` behaves as if it were copied along that dimension.

```python
# 4 x 1 array
x = np.array([[0],
              [10],
              [20],
              [30]])

# 1 x 3 array
y = np.ones(3)

# 4 x 3 array
a = x * y

print(x)
print()
print(y)
print()
print(a)

```

    [[ 0]
     [10]
     [20]
     [30]]

    [ 1.  1.  1.]

    [[  0.   0.   0.]
     [ 10.  10.  10.]
     [ 20.  20.  20.]
     [ 30.  30.  30.]]

```python
# a: 4 x 3 array (see cell above)

# 1 x 3 array
b = np.arange(3)

# 4 x 3 array
result = a + b

print(a)
print()
print(b)
print()
print(result)
```

    [[  0.   0.   0.]
     [ 10.  10.  10.]
     [ 20.  20.  20.]
     [ 30.  30.  30.]]

    [0 1 2]

    [[  0.   1.   2.]
     [ 10.  11.  12.]
     [ 20.  21.  22.]
     [ 30.  31.  32.]]

The size of the output array is the maximum size along each dimension of the input arrays.

The 4x3 and 1x4 arrays shown in the cell below cannot be broadcast together.

The dimensions 3 and 4 are incompatible.

Note that if the array dimensions are incompatible, it will generate a ValueError.

Recall, the function `np.insert` that we used earlier.

An integer (length=1) can be broadcast into an array of any size.

```python
# Add items to an array
b = np.array([[1, 1],
              [2, 2],
              [3, 3]])

b = np.insert(b, 1, 4, axis=1)
print(b)
```

    [[1 4 1]
     [2 4 2]
     [3 4 3]]

Here are some examples of practical applications of broadcasting.

### Broadcasting Example: Calorie Calculator

Let's say we have a large data set; each datum is a list of parameters.

Example datum: a type of food and the the amount of fat, protein and carbohydrate in a serving of that food.

Our data set of food nutrients might look something like the table below:

|Food (100g) |Fat(g)|Protein (g)|Carbohydrate (g)|
|------------|----- |-----------|----------------|
|Almonds     |    49|         21|              22|
|Peas        |     0|          5|              14|
|Avocado     |    15|          2|               9|
|Kale        |     1|          3|              10|

By applying the following sclaing factors, we can calculate the number of calories in a food type due to fat, protein and carbohydrate:
 -  fat: 9 cal/g
 -  protein: 4 cal/g
 -  carbohydrate 4 cal/g

Using what we have studied so far, we could convert the table to calories using a loop:

```python

nutrients = np.array([[49, 21, 22],
                      [0,   5, 14],
                      [15,  2,  9],
                      [ 1,  3, 10]])

cal_convert = np.array([9, 4, 4])

calories = np.empty((4,3))

for index, value in enumerate(nutrients):
    calories[index] = value * cal_convert

```

However, it is faster and more concise to broadcast the two arrays together:

```python
nutrients = np.array([[49, 21, 22],
                      [0,   5, 14],
                      [15,  2,  9],
                      [ 1,  3, 10]])

cal_convert = np.array([9, 4, 4])

calories = nutrients * cal_convert

print(calories)
```

    [[441  84  88]
     [  0  20  56]
     [135   8  36]
     [  9  12  40]]

### Broadcasting Example: Vector Quantisation Algorithm
This is a simple algorithm used for catagorisation.

It determines which catagory a data point should belong to from its closest proximity to a set of values representing possible catagories.

Each value represents the mean of the corresponding catagory.

For example, colour quantisation is used in image processing reduces the number of distinct colors used in an image, while maintianing visual similarity to the original image.

CC BY-SA 3.0, https://commons.wikimedia.org/w/index.php?curid=1477223

In the plot below, each of the circles represents the mean height and weight of athletes grouped by type.

The square represents the height and weight of an athlete to be classified.

To find the closet point:
1. Use broadcasting to find the difference between the position of the __square__ and the position of each __circle__ in the x and y directions.

1. Find the distance, $d$ from the square, $s$ to each circle, $c$ using:
$d = \sqrt{(x_{c}-x_{s})^2 + (y_{c}-y_{s})^2}$

1. Choose the group corresponding to the minimum distance, $d_{min}$

```python
athlete = np.array([111.0,188.0])

categories = np.array([[102.0, 203.0],
                       [132.0, 193.0],
                       [45.0, 155.0],
                       [57.0, 173.0]])

# 1. broadcast
diff = categories - athlete

# 2. distance to each point (magnitude of values along axis 1 for each datum)
# dist = np.linalg.norm(diff,axis=1)
dist = np.sqrt(np.sum(diff**2,axis=1))

# 3. which group?
nearest = np.argmin(dist)
print(nearest)
```

    0

The nearest group is index 0 of the array `catagories`.

Based on mean height and weight, the athlete is most likely to be a basketball player.

## Resizing and Reshaping
We can change the size of an array in each dimension.

For example, you may want to edit the length of a dimension of an array to make it compatible with another array for broadcasting.

### Resizing
We can resize an array.

If the new array size is smaller, the original array will be trimmed to the new size.

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

b = np.resize(a,(2,1))
print(b)
print()

a.resize(2,1)
print(a)
```

    [[0]
     [1]]

    [[0]
     [1]]

If the new array size is larger, the extra space can either be filled with repeating copies of the original array.

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

b = np.resize(a,(4,4))
print(b)
```

    [[0 1 2 3]
     [0 1 2 3]
     [0 1 2 3]
     [0 1 2 3]]

or with zeros.

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

a.resize(4,4)
print(a)
```

    [[0 1 2 3]
     [0 0 0 0]
     [0 0 0 0]
     [0 0 0 0]]

### Reshaping
You can re-shape the array.
The new array must have the __same number of elements__ as the original array.

__Example:__ Using range to create a column vector:

```python
x = np.arange(0, 31, 10)
y = x.reshape((4,1))

# which can be written in one line as:
z = np.arange(0, 31, 10).reshape((4,1))

print(x)
print()
print(y)
print()
print(z)
```

    [ 0 10 20 30]

    [[ 0]
     [10]
     [20]
     [30]]

    [[ 0]
     [10]
     [20]
     [30]]

## Review Exercises

The folowing exercises are provided to practise what you have learnt in today's seminar.

The extension excercises cover topics that will not be included inthe exam but may be useful to you if you are familiar with using matrices and want to practise matrix manipulation using Python.

If you have not yet studied matrices, you can come back to this section when the mathematics used is more familiar to you.

### Review Exercise: 1D Arrays - Indexing and changing values.

In the cell below:

1. Create an array of zeros with length 25.
 Jump to Creating a Numpy Array of Zeros.

2. Change the first 10 values to 5.
 Jump to Changing items in an array.

3. Change the next 10 values to a sequence starting at 12 and increasig with steps of 2 to 30 - do this with one command.

4. Change the final 5 values to 30.

```python
# 1D array
```

```python
# Example Solution

# 1
a = np.zeros(25)

#2
a[: 10] = 5

#3
a[10: 20] = range(12, 31, 2)

#4
a[20:] = 30

print(a)
```

    [  5.   5.   5.   5.   5.   5.   5.   5.   5.   5.  12.  14.  16.  18.  20.
      22.  24.  26.  28.  30.  30.  30.  30.  30.  30.]

### Review Exercise: 1D Arrays - Correct the error

The code below, is supposed to:
 - change the last 5 values of the array x to the values [50, 52, 54, 56, 58]
 - print the result
There are some errors in the code.

Remove the comment markers and run the code to see the error message.

Then fix the code and run it again.

```python
#x = np.ones(10)
#x[5:] = np.arange(50, 62, 1)
#print(x)
```

```python
# Example Solution
import numpy as np
x = np.ones(10)
x[-5:] = range(50, 59, 2)
print(x)
```

    [  1.   1.   1.   1.   1.  50.  52.  54.  56.  58.]

### Review Exercise:  Creating 2D Arrays

Print:

__(A)__ An array with 3 rows and 2 columns (3 by 2), where all elements are 1.0.

Jump to Creating a Numpy Array.

__(B)__ An array of randomly generated integers, in the range 1 to 10, with 4 rows and one column (4 by 1).

Jump to Subpackages.

__(C)__ An array with 2 rows and 2 columns (2 by 2), where all elements are 4.

```python
# 2D arrays
```

```python
### Example Solution
# A
print(np.ones((3,2)))

# B
print(np.random.randint(1, 10, size=(4,1)))

# C
print(np.full((2,2), 4))
```

    [[ 1.  1.]
     [ 1.  1.]
     [ 1.  1.]]
    [[5]
     [4]
     [7]
     [1]]
    [[4 4]
     [4 4]]

### Review Exercise:  2D Arrays - Indexing and Basic Plotting

__(A)__ In the cell below, for the array `x`, write code to print:

* the first row of `x`
* the first column of `x`
* the third row of `x`
* the four values in the upper right hand corner of `x`
* the four values at the center of `x`
* the last two columns (third and fourth column) of `x`

Jump to Indexing into multi-dimensional arrays.

__(B)__ You have just *printed* the last two columns of `x`.

Now plot the third column of `x` against the fourth column of `x` as a line graph.

Jump to Introduction to Plotting.

```python
# 2D array
x = np.array([[4, 2, 1, 1],
              [2, 4, 2, 8],
              [2, 4, 3, 27],
              [4, 1, 4, 64]])
```

```python
# Example Solution

x = np.array([[4, 2, 1, 1],
              [2, 4, 2, 8],
              [2, 4, 3, 27],
              [4, 1, 4, 64]])

print(x[0])

print(x[:, 0])

print(x[2])

print(x[0:2, 2:4])

print(x[1:3, 1:3])

print(x[:, 2:4])

plt.plot(x[:, 2], x[:, 3])
```

### Review Exercise:  Boolean array indexing

Print the values of array `x` in the cell below that are...

__(A)__ ...greater than 3.

__(B)__ ...greater than 2 and less than 8.

__(C)__ ...less than 2 or greater than 8.

__(D)__ ...a multiple of `y`.

__(E)__ Change all odd numbers in x to 0

Jump to Boolean Array Indexing.

```python
x = np.arange(10)
y = 2
```

### Review Exercise: Curve Fitting

Using the example data set `a`,`b`:

__(A)__ Plot a scatter graph of `a` against `b`, with `a` on the horizontal axis and `b` on the vertical axis.

Jump to Scatter plot.

__(B)__ Fit __three different__ curves to the data. You can choose from:
 - a linear relationship. Jump to Fitting a linear relationship.
 - a polynomial curve. Jump to Fitting a polynomial curve.
 - an arbitrary function Jump to Fitting an arbitrary curve.

__(C)__ Plot each curve as a line of on the graph.

Jump to Line plot.

__(D)__ Find the route mean square error (RMSE) of each curve relative to the data set.

Jump to Route mean square error.

__(E)__ Display the RMSE of each curve as a figure legend.

Jump to Adding a legend.

__(F)__ Print the equation of the *optimal* curve i.e. the curve with the *smallest* RMSE.

```python
a = np.array([88438,45505,75127,115571,89911,87432,100083,85589,73104,86890,78580,70785,41050,57610,107537,59262,73038,87891,75368,111638,74911,71599,96774,79667,90725,93816,75859,64969,205688,71500,53098,71250,89615,94747,50400,63673,78257,72785,83015,150000,84699,67191,86298,117705,88935,89643,106678,97894,132164,59387,60684,96151,68794,74559,29430,88362,111792,57205,83651,87518,80129,86801,110761,63274,66143,110694,52590,59994,80460,103589,68298,59056,40294,161848,103100,86354,37428,43307,80792,77368,109159,71538,84783,86250,82900,74728,48597,75549,106942,102167,62708,60630,70273,84918,88693,74141,46627,119112,88260,97262,86095,110472,82734,84761,91715,103292,86339,147993,77560,100625,68094,78250,75426,86138,112344,115000,98846,90499,80029,61959,76779,68833,81026,66361,92737,76692,64974,103869,51951,108854,61038,75938,75346,40639,73156,80067,82322,52353,
              62832,207262,160106,77740,72011,167094,58458,41639,79528,66583,83993,138082,77366])
b = np.array([1.7,-0.4,0.5,2.6,1.4,1.5,1.5,1.7,-0.5,1.6,0.9,1.1,-1.7,0.3,1.8,0.5,1,1.9,0.1,2,1.7,1,1.2,1.5,1,1.1,1.2,0,2.6,1.4,-0.8,1.6,1.1,1.2,-1.4,-0.5,1.9,0,1.5,2.4,1.5,0.7,1.8,2,2.4,1.6,2,2.3,2,0.1,0.3,2.3,0,0,-1.7,1.9,2,0,0.9,1.3,0.4,1.6,2.3,-0.1,1.7,2.1,-0.9,0.1,1,1.9,0.4,-0.3,-2.4,2.7,1.3,2,-1.3,-1.5,0.7,1.1,2.3,1.1,0.7,0.9,1.1,0.1,-0.9,1.4,2.1,1.2,0.1,0.8,0.3,1.4,1.5,1,-0.5,2.4,0.9,1.5,1.6,1.2,1.3,1.8,0.8,1.8,1.9,2.6,1.5,1.8,1.8,0.6,0.7,1.2,1.5,2.5,1.1,1.6,1.6,1,0,0,1,0.5,1.7,0.6,0.1,1.7,0.2,2.1,0.1,0.9,0.8,-1.3,1.3,0.5,1.5,-0.6,1.2,2.4,2.6,1.1,0.8,2.5,-0.2,-2,0.1,0.1,1.6,2.6,1.2])

```

```python
### Example solution

a = np.array([88438,45505,75127,115571,89911,87432,100083,85589,73104,86890,78580,70785,41050,57610,107537,59262,73038,87891,75368,111638,74911,71599,96774,79667,90725,93816,75859,64969,205688,71500,53098,71250,89615,94747,50400,63673,78257,72785,83015,150000,84699,67191,86298,117705,88935,89643,106678,97894,132164,59387,60684,96151,68794,74559,29430,88362,111792,57205,83651,87518,80129,86801,110761,63274,66143,110694,52590,59994,80460,103589,68298,59056,40294,161848,103100,86354,37428,43307,80792,77368,109159,71538,84783,86250,82900,74728,48597,75549,106942,102167,62708,60630,70273,84918,88693,74141,46627,119112,88260,97262,86095,110472,82734,84761,91715,103292,86339,147993,77560,100625,68094,78250,75426,86138,112344,115000,98846,90499,80029,61959,76779,68833,81026,66361,92737,76692,64974,103869,51951,108854,61038,75938,75346,40639,73156,80067,82322,52353,
              62832,207262,160106,77740,72011,167094,58458,41639,79528,66583,83993,138082,77366])
b = np.array([1.7,-0.4,0.5,2.6,1.4,1.5,1.5,1.7,-0.5,1.6,0.9,1.1,-1.7,0.3,1.8,0.5,1,1.9,0.1,2,1.7,1,1.2,1.5,1,1.1,1.2,0,2.6,1.4,-0.8,1.6,1.1,1.2,-1.4,-0.5,1.9,0,1.5,2.4,1.5,0.7,1.8,2,2.4,1.6,2,2.3,2,0.1,0.3,2.3,0,0,-1.7,1.9,2,0,0.9,1.3,0.4,1.6,2.3,-0.1,1.7,2.1,-0.9,0.1,1,1.9,0.4,-0.3,-2.4,2.7,1.3,2,-1.3,-1.5,0.7,1.1,2.3,1.1,0.7,0.9,1.1,0.1,-0.9,1.4,2.1,1.2,0.1,0.8,0.3,1.4,1.5,1,-0.5,2.4,0.9,1.5,1.6,1.2,1.3,1.8,0.8,1.8,1.9,2.6,1.5,1.8,1.8,0.6,0.7,1.2,1.5,2.5,1.1,1.6,1.6,1,0,0,1,0.5,1.7,0.6,0.1,1.7,0.2,2.1,0.1,0.9,0.8,-1.3,1.3,0.5,1.5,-0.6,1.2,2.4,2.6,1.1,0.8,2.5,-0.2,-2,0.1,0.1,1.6,2.6,1.2])

# A
plt.plot(a, b,'o')

# B
c, d, e = np.polyfit(a, b, 2)

f, g, h, i = np.polyfit(a, b, 3)

j, k, l, m, n = np.polyfit(a, b, 4)

# C
yfit2 = np.poly1d([c, d, e])(a)
yfit3 = np.poly1d([f, g, h, i])(a)
yfit4 = np.poly1d([j, k, l, m, n])(a)

# yfit2 = c*a**2 + d*a + e
# yfit3 = f*a**3 + g*a**2 + h*a + i
# yfit4 = j*a**4 + k*a**3 + l*a**2 + m*a + n

plt.plot(a, yfit2,'o')
plt.plot(a, yfit3,'o')
plt.plot(a, yfit4,'o')

# D
# Define a function
def RMSE(ideal, real):
    return round(
        np.sqrt
        (np.sum
         ((ideal - real)**2)/ len(real))
        , 4)

# error
rmse2 = RMSE(yfit2, b)
rmse3 = RMSE(yfit3, b)
rmse4 = RMSE(yfit4, b)

# E
# You don't need to replot the data as shown here, you can just add the labels to your answer to part C
plt.plot(a, yfit2,'o', label = f"RMSE={rmse2}")
plt.plot(a, yfit3,'o', label = f"RMSE={rmse3}")
plt.plot(a, yfit4,'o', label = f"RMSE={rmse4}")
plt.legend(loc='best')

# F
print(np.min([rmse2, rmse3, rmse4]))
```

### Review Exercise:  Vectorising a function

Recall the function `is_positive` that we looked at in Seminar 5: Functions.

A similar function is shown in the cell below.

It takes an argument, `x` and returns:
- `True` if `x` is positive.
- `False` if `x` is negative.

Write a vectorised version of the function that will accept a Numpy array as an argument.

Jump to Vectorising a function.

Demonstrate your function works by calling it.

```python
def is_positive(x):
    "Checks if the function value f(x) is positive"
    return x > 0
```

```python
# Example Solution

def is_positive(x):
    "Checks if the function value f(x) is positive"
    return x > 0

is_positive_vec = np.vectorize(is_positive)

z = [1, 2, 3]

is_positive_vec(z)
```

    array([ True,  True,  True], dtype=bool)

### Review Exercise: Broadcasting
Use a single broadcasting operation to:
- add 2 to every element in the first column of `a`
- add 3 to every element in the second column of `a`
- subtract 1 from every element in the third column of `a`

```python
a = np.array([  [49, 21, 22],
                [0,   5, 14],
                [15,  2,  9],
                [ 1,  3, 10]])
```

```python
# Example Solution
print(a + np.array([2, 3, -1]))
```

    [[51 24 21]
     [ 2  8 13]
     [17  5  8]
     [ 3  6  9]]

### Extension Exercise:  Broadcasting and Mathematics with matrices - Simultaneous Equations, Polynomial Coeffcients

The equation of a parabola can be fully defined using three points that lie on the curve (provided the curve is not a  straight line).

The equation for a parabola is $y=ax^2+bx+c$, where $a$, $b$, and $c$ are constants.

__(A)__Given three points $(x_1,y_1)$, $(x_2,y_2)$, $(x_3,y_3)$,

three simultaneous equation may be written:
$$
\begin{split}
x_1^2a+x_1b+c&=y_1 \\
x_2^2a+x_2b+c&=y_2 \\
x_3^2a+x_3b+c&=y_3 \\
\end{split}
$$

Given that:

$(x_1,y_1)=(-2,2)$

$(x_2,y_2)=(1,-1)$

$(x_3,y_3)=(4,4)$

use matrix multiplication to find $a$,$b$ and $c$.

Jump to Simultaneous Equations.

Jump to Broadcasting.

Hint:

Build a 1D array, `y`,  containing values $y_1$, $y_2$ and $y_3$.

Build a 2D array, `x`, with which to multiply `y` to find a 1D array containing $a$, $b$ and $c$.

__(B)__ Show that your solution is correct by:
 - plotting a graph of $y=ax^2+bx+c$ against $x$ for range of x.
 - plotting the three points:
$(x_1,y_1)=(-2,2)$
$(x_2,y_2)=(1,-1)$
$(x_3,y_3)=(4,4)$
as points to check that they lie on the curve.

```python
# Solve the simultaneous equations to find a, b and c
```

```python
# Example solution

# 1D array, y
y = np.array([2, -1, 4])

# 2D array, x
x = np.array([[(-2**2), -2, 1],
              [(1**2),   1, 1],
              [(4**2),   4, 1]])

abc = np.linalg.solve(x, y)
print(abc)

# Another way to create matrix x without writing it out explicitly is to use:
# - broadcasting
# - the numpy function, power (raises first argument to power of second)
xs = np.array([-2, 1, 4])
xs = xs.reshape((3,1))
power = np.array([2, 1, 0])
x = np.power(xs, power)
```

    [ 0.8        -2.33333333  0.53333333]

### Extension Exercise:  Solving systems of first order differential equations - Predator & prey

The predator-prey equations represent a simplified model of two species interacting; for example:
 - rabbits $x$(prey)
 - population $y$: foxes (predators)

*The more foxes there are*, the more rabbits get eaten by foxes and the population of rabbits gets smaller...

...so there is less food for the foxes so the population of foxes gets smaller...

...so there is less danger to the rabbits, so the population fo rabbits gets bigger...

...so there is more food for the foxes, so the population of foxes gets bigger...

...*The more foxes there are*....

We can expect to see the two populations rise and fall cyclically with time.

The two related populations can be described by a pair of first-order non-linear ordinary differential equations.

\begin{align*}
\frac{dx}{dt}&=x(a-by)\\
\frac{dy}{dt}&=-y(c-dx)\\
\end{align*}

Where:

$a$ represents rabbit births

$b$ represents rabbit deaths due to foxes

$c$ represents fox deaths in the absence of rabbits

$d$ represents foxes kept alive by feeding on rabbits

($a,b,c,d$ are assumed to be positive)

__(A)__ Use `odeint` to estimate the solution to the couple differntial equations:
 - between time = 0 and time = 20 units (suggested time step = 0.1 unit)
 - where $a=b=c=d=1$
 - for initial conditions:
$x(0)=150$
$y(0)=100$

Jump to  Systems of Equations for Rates of Change.

__(B)__ Plot the two populations using a graph.
Include a legend to label the two populations.

__(C)__ Try changing the values of $a,b,c,d$ and the initial number of rabbits and foxes and observe the output in your graph.

```python
# Predator-prey equations
```

```python
# Example solution

# A
# Between time = 0 and time = 100
ts = np.linspace(0, 100, 100)

a,b,c,d = 1,1,1,1

# Define function that returns differential
def dP_dt(P, t):

    return [P[0]*(a - b*P[1]),
           -P[1]*(c - d*P[0])]

# Initial conditions
P0 = [150, 100]

# Solve using odeint
Ps = odeint(dP_dt, P0, ts)

# B
# plot the two populations using a graph
foxes = Ps[:,0]
rabbits = Ps[:,1]

plt.plot(ts, rabbits, label="rabbits")
plt.plot(ts, foxes,label="foxes")

t = ts[0:50:4]
r = rabbits[0:50:4]
f = foxes[0:50:4]

plt.xlabel("Time")
plt.ylabel("Population")
# Include a legend
plt.legend();

print(t)
print(r)
print(f)
```

# Summary

- Numpy arrays are a type of data structure optimised for numerical computations.
- Unlike data structure such as lists and tuples, an array can store data in multiple dimensions.
- Numpy array are a convenient way to represent vectors and matrices.
- Numpy funcions generally behave in an "elementwise" way.
- We can "vectorise" other functions so that they accept data structure as arguments.
- Broasdcasting is a useful tool for applying information in one numpy array to another without having to repeat or reshape the arrays to match one another.
- The matplotlib library can be used to quicky produce simple plots to visualise and check your solutions to mathematical problems.
- This includes fitting a curve or a relationship to a dataset.

###### Homework

1. __PULL__ the changes you made in-class today to your personal computer.
1. __COMPLETE__ any unfinished Review Exercises.
1. __PUSH__ the changes you make at home to your online repository.
1. Install __ffmpeg__...

#### Installing FFmpeg
FFmpeg is software for handling multimedia data.

You must have it installed to view the animatied figures we will produce in next week's seminar.

###### Installing FFmpeg on mac
Open a terminal.
Copy and paste the following command into the terminal to install homebrew (a linux-like package manager):
>`ruby -e "$(curl -fsSL https://raw.githubusercontent.com/Homebrew/install/master/install)"`

Copy and paste the following command into the terminal to install FFmpeg:
>`brew install ffmpeg`

To check ffmpeg has installed copy and paste the following command into the terminal:
>`ffmpeg -version`

If FFmpeg has installed a few lines of code will appear, starting with the version number which will be something like:
>`ffmpeg version 3.4 Copyright (c) 2000-2017 `

###### Installing FFmpeg on linux
Open a terminal.
Copy and paste the following commands into the terminal (one-by-one, pressing enter after one) to install ffmpeg:
>`sudo add-apt-repository ppa:kirillshkrogalev/ffmpeg-next
sudo apt-get update
sudo apt-get install ffmpeg`

To check ffmpeg has installed copy and paste the following command into the terminal:
>`ffmpeg -version`

If FFmpeg has installed a few lines of code will appear, starting with the version number which will be something like:
>`ffmpeg version 3.4 Copyright (c) 2000-2017 `

###### Installing FFmpeg on  windows

Open a terminal.

To find out what version of windows you have, copy and paste the following command into the terminal to see if your computer has a 32 or 64 bit CPU:
>`wmic os get osarchitecture`

Go to ffmpeg.zeranoe.com/builds/
 - Click the appropraite __Architecture__ for your computer.
 - Click Static __Linking__
 - Leave __Version__ as the default.
 - Click __Download Build__

Go to your computer's Downloads folder.
 - Right click on the newly downloaded .zip folder (it's name will start with `ffmpeg`).
 - Choose __Extract All__ from the drop-down menu.
 - If given the option to choose a location to extract the files to, choose your Program Files folder.
Otherwise, a non-.zip folder of the same name as the one you downloaded will appear in the Downloads folder.
Copy the folder to your Program Files folder.
 - Change rename of the folder you just copied into Program Files with the name:

`ffmpeg`

Go back to the terminal and copy and paste the following command into the terminal to add ffmpeg to the windows path:
>`PATH=C:\Program Files\ffmpeg\bin;%PATH%`

To check ffmpeg has installed copy and paste the following command into the terminal:
>`ffmpeg -version`

If FFmpeg has installed a few lines of code will appear, starting with the version number which will be something like:
>`ffmpeg version 3.4 Copyright (c) 2000-2017 `

Please do this before next week's class.

```python

```


text:
# EDOs de segundo orden: problemas de valores a la frontera

Todos los métodos vistos en las tres clases anteriores pretendían resolver **problemas de valor inicial**  en los que conocíamos el valor de la función $\mathbf{x}(t_0)$ y buscábamos conocer $\mathbf{x}(t)$ para $t \geq t_0$.

Sin embargo, muchas veces nos enfrentamos a problemas en los que no se conocen en todas las variables en el mismo tiempo, si no en en dos tiempos: los extremos del intervalo $[t_0,t_f]$. Nos interesa conocer $\mathbf{x}(t)$ en ese intervalo. En ese caso, decimos que estamos tratando con un problema de **valores a la frontera**.

## Preludio Matemático: Problema con condiciones en varios puntos

Existe todo un formalismo para manejar esta clase de problemas y para describirlos. Por simplicidad, nosotros nos concentraremos en problemas de sistemas de dos ecuaciones diferenciales de primer orden, obtenidos a partir de una **ecuación diferencial lineal de segundo orden**, que tiene la siguiente forma

$$
\begin{align}
a(t) \ddot{y}(t) + b(t) \dot{y}(t) &+  c(t) y(t) = f(t) \\
& \\
y(t_0) &= \alpha \\
y(t_f) &= \beta \\
\end{align}
$$

Con $a(t) \neq 0$, $b(t)$, $c
(t)$ y $f(t)$ funciones conocidas, $\alpha$ y $\beta$ constantes.

## Métodos de Shooting

### Planteamiento

Usando el truco para convertir la ecuación de segundo orden en un sistema de ecuaciones diferenciales ordinarias, sea $x_1(t) = y(t)$, $x_2(t) = \dot{y}(t)$. Así, el problema se convierte en resolver el sistema:

$$
\begin{align}
\dot{x_1} &= \dot{y} =  x_2 \\
\dot{x_2} &= \ddot{y} =  \frac{ f(t) -  b(t) x_2(t) -  c(t) x_1(t)  }{a(t)} \\
\end{align}
$$

Y las condiciones se vuelven

$$
\begin{align}
x_1(t_0) &= \alpha \\
x_1(t_f) &= \beta \\
\end{align}
$$

Para simplificar la notación, denotemos

$$
g(x_1,x_2,t) = \frac{ f(t) -  b(t) x_2(t) -  c(t) x_1(t)  }{a_2(t)} \\
$$

El sistema obtenido es muy parecido a los problemas de valores iniciales que ya sabemos resolver. La única diferencia es que no conocemos $x_2(t_0)$, si no $x_1(t_f)$. Sin embargo, **supogamos que si lo conozco**: sea $x_2(t_0)=z$ con $z$ una constante arbitraria. Entonces obtengo el sistema:

$$
\begin{align}
\dot{x_1} &= x_2 \\
\dot{x_2} &= g(x_1,x_2,t)  \\
& \\
x_1(t_0) &= \alpha \\
x_2(t_0) &= z \\
\end{align}
$$

El sistema ya cumple una de las condiciones de frontera requeridas, que $x_1(t_0) = \alpha$. Sin embargo, no sabemos si la otra condición original ($x_1(t_f)=\beta$) se cumple. Si ahora utilizamos uno de los métodos (no adaptativos, con paso de tiempo constante) para problemas de valores iniciales a dicho sistema, desde el tiempo $t_0$ hasta el tiempo $t_f$, usando un paso $h$ que **divide excatamente** al intervalo $[t_0,t_f]$, es decir, que cumple

$$
\frac{t_f - t_0}{h} = n \; \in \mathbb{N}
$$

los últimos valores obtenidos por el método corresponderán a  los valores $x_1(t_f)$ y $x_2(t_f)$.

Es claro que dichos valores dependerán de la condición incial $x_2(t_0)=z$. Así, denotemos al valor de $x_1(t_f)$ obtenido utilizando la condición incial $x_2(t_0)=z$ como $x_1(t_f ; z)$. Sabemos que dicho valor **segurmante no va a coincidir** con el valor real $\beta$. Así, podemos definir una función de diferencia $F(z)$, definida como

$$
F(z) = x_1(t_f ; z) - \beta
$$

Es claro que si existe un valor $z^*$ tal que $F(z^*) =0$ , entonces para ese $z^*$ la solución al problema de valor incial cumplirá la condición de frontera que queremos, por lo que para ese valor $z$ el problema de valor inicial se vuelve completamente equivalente al problema de valores a la frontera, por lo que **sus soluciones $x_{1}(t)$, $x_2(t)$ corresponden a las mismas funciones**

### ¿Cómo encontramos el valor $z^*$?

Nuevamente, notemos que el problema de encontrar $z^*$ tal que $F(z^*) =0$ corresponde a **resolver una ecuación trascendental** $F(z^*) =0$. Ya estudiamos previamente los métodos para resolver dicha ecuación, por lo que podemos utilizar cualquiera de ellos, pero utilizaremos en este caso el de Newton-Rhapson.

Así, los **métodos de Shooting** (*shooting* por que, al igual que con un cañon, vamos haciendo pruebas desde un punto inicial para ver si acertamos en el punto final) consisten en resolver la ecuación trascendental $F(z^*) =0$ y luego utilizar el valor $z^*$ para resolver el problema de valores iniciales que es equivalente al problema de valores a la frontera.

Noten que, por la definición de $F(z)$, para calcular el valor de la función en un punto **es necesario resolver numéricamente todo el sistema de ecuaciones diferenciales** con condición inicial $x_2(t_0)=z$. Esto va a volver  muy lento el procedimiento de resolver la ecuación trascendental.

### Ejercicio 1 y 2

Implementa una función `shootingRK4(alpha,beta,g,t_0,t_f,n,epsilon)` que tome como argumentos:

* `alpha`, `beta`: las condiciones a la frontera $\alpha$, $\beta$ del problema

* `g` la función **vectorial** $g(x_1,x_2,t)$

* `t_0`,`t_f` los tiempos inicial y final, `n` $>2$ un natural correspondiente el número de subintervalos que se van a tomar entre `t_0` y `t_f`. Eso **define** el paso de tiempo $h$ del sistema a resolver como

$$
\frac{t_f - t_0}{n} = h
$$

* `epsilon` $>0$, un  número real arbitrario que será nuestra tolerancia para el método de Newton-Rhapson.

La función va a implementar el método de Shooting resolviendo el sistema con la función `RK4` de la clase 16. La función debería proceder de la siguiente forma:

1. Primero definimos la función $F(z)$ que utilizando la función `RK4`, obtiene $x_1(t_f;z)$ y calcula $x_1(t_f;z) - \beta$

2. Ya definida $F(z)$, usar la función `miNewtonTol` de la clase 08 para aplicarle el método de Newton-Rhapson a la función $F$, utilizando una diferencia finita centrada y una tolerancia `epsilon`. Esto nos dará el valor $z^*$

3. Ya obtenido el valor de $z^*$, se debe de volver a utilizar `RK4` con la condición inicial $x_2(t_0) = z^*$ para resolver el sistema de ecuaciones.

Tu función debe de regresar el arreglo de tiempos `[t_0,...,t_n]` y las soluciones `[x_0,..,x_n]`.

### Ejercicio 3

Usa tu función `shootingRK4` para resolver los siguientes problemas de valores a la frontera

1. $\ddot{y} + 4y =0$, $y(0) = -2$, $y(\pi/4) = 10$

Respuesta: $y(t) = -2 \cos{(2t)} + 10 \sin{(2t)}$

2.  $t^2 \ddot{y} - 2 t \dot{y} + 2y =0$, $y(1) = -1$, $y(2)=1$

Respuesta: $y(t) = -\frac{5}{2}x + \frac{3}{2} x^2$

3. $\ddot{y} + 9y = \cos{x}$, $y(0) = 1$, $y(\pi/2)=-5/3$

Respuesta: $y(t) = \frac{1}{8}\cos{t} + \frac{7}{8}\cos{(3t)} + \frac{5}{3} \sin{(3t)}$

En una misma figura, compara tus soluciones con la solución analítica.

**Sugerencia** Ya que el método de Newton-Raphson necesita un valor inicial para $z$, puedes ayudarte de la solución analítica para proponer un valor de $z$ inicial cercano al real para que tu método de Newton no necesite demasiadas iteraciones.

## Métodos de diferencias finitas

Es claro que el método de shooting es **muy lento**, pues el tener que resolver toda la ecuación diferencial para evaluar $F$ puede ser muy tardado, además de que la solución obtenida siempre es aproximada, por lo que puede haber errores ahí que nos alejen de la condición verdadera y que nos den una mala solución.

Una alternativa para resolver problemas de valores a la frontera es utilizar diferencias finitas para aproximar las derivadas $\ddot{y}$ y $\dot{y}$. Utilizando una diferencia centrada para ambas derivadas, llegamos a la siguiente aproximación:

$$
\begin{align}
a(t) \frac{y(t+h) + y(t-h) - 2y(t)}{h^2} &+ b(t) \frac{y(t+h) - y(t-h)}{2h} +  c(t) y(t) = f(t) \\
& \\
y(t_0) &= \alpha \\
y(t_f) &= \beta \\
\end{align}
$$

Así, para un a partición regular de orden $n+1$ (i.e con $n$ subintervalos) del intervalo $[t_0,t_f]$, obtenemos los tiempos $t_0,\ldots,t_n$ que corresponden a un paso  $h = (t_f - t_0)/n$. Utilizando la notación de índices para el tiempo, el problema se convierte en:

$$
\begin{align}
a_i \frac{y_{i+1} + y_{i-1} - 2y_i}{h^2} &+ b_i \frac{y_{i+1} - y_{i-1}}{2h} +  c_i y_i = f_i \\
& \\
y_0 &= \alpha \\
y_n &= \beta \\
\end{align}
$$

Con $a_i = a(t_i)$, $b_i=b(t_i)$, $c_i = c(t_i)$, $f_i = f(t_i)$, $y_{i} = y(t_i)$

Podemos reordenar las ecuaciones para que quede de la siguiente forma:

$$
\begin{align}
\left(\frac{a_i}{h^2} - \frac{b_i}{2h}\right) y_{i-1} + \left(c_i - \frac{2a_i}{h^2} \right)& y_{i} +  \left(\frac{a_i}{h^2} + \frac{b_i}{2h}\right) y_{i+1} = f_i \\
& \\
y_0 &= \alpha \\
y_n &= \beta \\
\end{align}
$$

Si consideramos que la primera ecuación solo se cumple para $1 \leq i \leq n-1$, entonces el conjunto de ecuaciones representa un **sistema de $n+1$ ecuaciones lineales** de la forma $\mathbf{Ay} = \mathbf{b}$, para el que la solución $\mathbf{y} = (y_0,y_1,..,y_n)$ es la discretización de la función solución $y(t)$ a nuestra ecuación diferencial.

Dicho sistema puede resolverse utilizando los métodos (eliminación Gaussiana, Jacobi, Gauss-Siedel) que vimos en las clases de Álgebra Lineal.

### Ejercicio 4

Implementa una función `matriz(a,b,c,f,alpha,beta,t_0,t_f,n)`. Los argumentos de la función deben de ser

* `a`,`b`,`c`,`f` las funciones **escalares** $a(t)$,$b(t)$,$c(t)$ y $f(t)$

* `alpha`, `beta`: las condiciones a la frontera $\alpha$, $\beta$ del problema

* `t_0`,`t_f` los tiempos inicial y final, `n` $>2$ un natural correspondiente el número de subintervalos que se van a tomar. Eso define el paso de tiempo $h$ del sistema a resolver como

$$
\frac{t_f - t_0}{n} = h
$$

La función debe de regresar la matriz $A$, de $n+1 \times n+1$ y un vector $b$ de longitud $n+1$ que correspondan a la discretización del problema de valores a la frontera expuesto.

**Sugerencia:** la matriz $A$ debe de ser **tridiagonal**, es decir, los únicos elementos distintos de cero son elementos de la forma $A_{i,i}$, $A_{i+1,i}$ y $A_{i-1,i}$

### Ejercicio 5

Utiliza tu función `matriz` y alguno de los métodos de las clases de Álgebra lineal para resolver las ecuaciones del ejercicio 3. Compara tus soluciones con la solución analítica y con las soluciones obtenidas por el método de Shooting. ¿Es mejor este método o el método de Shooting?

**Sugerencia** si tienes dudas sobre si tu método de álgebra lineal converge para esa matriz, usa el operador `\` para apoyarte.

## Mecánica cuántica: la ecuación de Schrödinger independiente del tiempo

Un problema de valores a la frontera que surge comúnmente en física es la solución de la ecuación de Schrödinger dependiente del tiempo

La evolución en el tiempo de un sistema cuántico está definida por la ecuación de Schrödinger dependiente del tiempo:

$$
i \hbar \frac{\partial \psi}{\partial t} = H \psi
$$

Con $i$ la unidad imaginaria, $\hbar$ la constante reducida de Planck, $\psi$ la amplitud de probabildad del sistema y $H$ el hamiltoniano. La solución general a la ecuación es

$$
\psi = \sum_{k=1}^{\infty}A_k e^{-E_k t} \psi_k
$$

Con $\psi_k$ y $E_k$ las funciones propias y valores propios, respectivamente, del operador hamiltoniano, es decir, las que cumplen con la propiedad

$$
H \psi_k = E_k \psi_k
$$

A esa ecuación se le llama la **ecuación de Schrödinger independiente del tiempo**. Trabajando un sistema d una partícula en una dimensión, cuya energía potencial es un pozo infinito de la forma

$$
V(x) = \begin{cases}
0 & \text{si } 0 \leq x \leq 1 \\
\infty & \text{en otro caso}
\end{cases}
$$

Escribiendo el hamiltoniano de dicho sistema cuántico y haciendo una adimensionalización, la ecuación de Schrödinger independiente del tiempo se reduce al siguiente problema de valores iniciales

$$
\frac{1}{2} \frac{d^2 \psi_k}{dx^2}(x) + \epsilon_k \psi_k (x) = 0
 \\
\psi_{k} (0) = 0 \\
\psi_{k} (1) = 0 \\
$$

La solución analítica está dada por
$$
\psi_k(x) = \sqrt{2} \sin{(k\pi x)} \\
\epsilon_k = \frac{k^2 \pi^2}{2} \\
$$

### Ejercicio 6

Utiliza los valores de la solución analítica para $\epsilon_k$ y resuelve, utilizando el método de tu preferencia (Shooting o diferencias finitas), el problema de valor incial para obtener $\psi_k$. Compara tu solución con la solución analítica.

## La visión general

Resolver ecuaciones diferenciales ordinarias en la computadora es un problema fundamental de la física computacional por su relación íntima con mecánica clásica. No todos los métodos pueden resolver cualquier sistema, por lo que es importante buscar métodos con alta precisión para no tener errores como el del péndulo simple con el método de Euler

Existe un area dedicada a modelar los átomos de manera clásica. Así, resolver un sistema **termodinámico** (moleculas de aire en un cuarto, un fluido ideal) o **biológicos** (átomos en una proteína, moléculas en una membrana,etc) se convierte en resolver las ecuaciones de movimiento de todos los átomos del sistema. En 3 dimensiones, si tenemos $n$ átomos, habrá un sistema de $6n$ ecuaciones difernciales de primer orden que debemos de resolver para saber la evolución.

A esas simulaciones se les llama **dinámica molecular**. La dinámica molecular es una herramienta **revolucionaria** de la ciencia, pues permite resolver muchisísimos problemas. Hoy en día, se utlizan simulaciones de dinámica molecular para probar la eficiencia de fármacos y segurmante se están utilizando dichas simulaciones para probar y crear vacunas para el COVID19

```julia

```

ext: ipynb