# Python and Jupyter notebooks: part 01

Last revised: 01-Jan-2020 by Dick Furnstahl [furnstahl.1@osu.edu]

This is a whirlwind tour of just the minimum we need to know about Python and Jupyter notebooks to get started doing mechanics problems.  We'll add more features and details as we proceed.

A Jupyter notebook is displayed on a web browser on a computer, tablet (e.g., IPad), or even your smartphone.  The notebook is divided into *cells*, of which two types are relevant for us:
* Markdown cells: These have headings, text, and mathematical formulas in $\LaTeX$ using a simple form of HTML (language of webpages) called markdown.
* Code cells: These have Python code (or other languages, but we'll stick to Python).

Either type of cell can be selected with your cursor and will be highlighted in color when active.  You evaluate an active cell with shift-return (as with Mathematica) or by pressing `Run` on the toolbar.  Some notes:
* When a new cell is inserted, by default it is a Code cell and will have `In []:` in front.  You can type Python expressions or entire programs in a cell.  How you break up code between cells is your choice and you can always put Markdown cells in between.  When you evaluate a cell it gets the next number, e.g., `In [5]:`.
* On the menu bar is a pulldown menu that lets you change back and forth between Code and Markdown cells.  Once you evaluate a Markdown cell, it gets formatted (and has a blue border).  To edit the Markdown cell, double click in it. 

**Try double-clicking on this cell and then shift-return.**  You will see that a bullet list is created just with an asterisk and a space at the beginning of lines (without the space you get *italics* and with two asterisks you get **bold**).  **Double click on the title header above and you'll see it starts with a single #.**  Headings of subsections are made using ## or ###.  See this [Markdown cheatsheet](https://github.com/adam-p/markdown-here/wiki/Markdown-Cheatsheet) for a quick tour of the Markdown language (including how to add links!).

**Now try turning the next (empty) cell to a Markdown cell and type:** `Einstein says $E=mc^2$` **and then evaluate it.**  This is $\LaTeX$! (If you forget to convert to Markdown and get `SyntaxError: invalid syntax`, just select the cell and convert to Markdown with the menu.)


```python

```

The menus enable you to rename your notebook file (always ending in `.ipynb`) or `Save and Checkpoint` to save the changes to your notebook.  You can insert and delete cells (use the up and down arrows in the toolbar to easily move cells).  You will often use the `Kernel` menu to `Restart` the notebook (and possibly clear output).

## Ok, time to try out Python expressions and numpy

We can use the Jupyter notebook as a super calculator much like Mathematica and Matlab.  **Try some basic operations, modifying and evaluating the following cells, noting that exponentiation is with** `**` **and not** `^`.


```python
 1 + 1  # Everything after a number sign / pound sign / hashtag) 
        #  is a comment
```




    2




```python
3.2 * 4.713
```




    15.081600000000002



If we want a floating point number (which will be the same as a `double` in C++), we *always* include a decimal point (even when we don't have to) while a number without a decimal point is an integer. (Note: in the Python we use there is no problem with 1/2 being evaluated as zero, as in C++, but we will still keep the convention of including decimal points.)


```python
3.**2
```




    9.0



We can define integer, floating point, and string variables, perform operations on them, and print them.  Note that we don't have to predefine the type of a variable (unlike C++) and we can use underscores in the names (unlike Mathematica).  **Evaluate the following cells and then try your own versions.** 


```python
x = 5.
print(x)
x   # If the last line of a cell returns a value, it is returned as output.
```

    5.0





    5.0




```python
y = 3.*x**2 - 2.*x + 7.
print('y = ', y)           # Strings delimited by ' 's (or " "s)
print(f'y = {y:.0f}')      # Just a preview: more on formatting later 
print(f'y = {y:.2f}')      #  (note: this uses an "fstring")
```

    y =  72.0
    y = 72
    y = 72.00



```python
first_name = 'Dick'     # Strings delimited by ' 's
last_name = 'Furnstahl'
full_name = first_name + ' ' + last_name  # you can concatenate strings 
print(full_name)
```

    Dick Furnstahl


Ok, how about square roots and trigonometric functions and ...


```python
ans = np.sqrt(2)
print(f'The answer is {ans: .4f}')
```

    The answer is  1.4142



```python
np.sin(np.pi)
```




    1.2246467991473532e-16



We need to `import` these functions through the numpy library. There are other choices, but numpy works with the arrays we will use.  Note: *Never* use `from numpy import *` instead of `import numpy as np`.  Here `np` is just a abbreviation for numpy (which we can choose to be anything, but `np` is conventional).


```python
import numpy as np
```


```python
print(np.cos(0.))
```

    1.0


Now functions and constants like `np.sqrt` and `np.pi` will work.  **Go back and fix the square root and sine.**

### Debugging aside . . .

Suppose you try to import and it fails (**go ahead and evaluate the cell**):


```python
import numpy
```

When you get a `ModuleNotFoundError`, the first thing to check is whether you have misspelled the name. Try using Google, e.g., search for "python numpie". In this case (and in most others), Google will suggest the correct name (here it is numpy).  If the name does exist, check whether it sounds like the package you wanted.

If you have the correct spelling, check whether you have installed the relevant package.  If you installed Python with Anaconda (which we will assume you did -- if not, do it!), then use `conda list`, e.g., `conda list numpy` in a Terminal window (on a Mac or Linux box) or in an Anaconda Prompt window (on a Windows PC).

### numpy arrays

At first, we will most often use numpy arrays so we'll start with those.  They are *like* lists delimited by square brackets, i.e., `[]`s, and we will construct them with `np.arange(min, max, step)` to get an array from `min` to `max` in steps of `step`. Examples:


```python
t_pts = np.arange(0., 10., .1)
t_pts
```




    array([0. , 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1. , 1.1, 1.2,
           1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2. , 2.1, 2.2, 2.3, 2.4, 2.5,
           2.6, 2.7, 2.8, 2.9, 3. , 3.1, 3.2, 3.3, 3.4, 3.5, 3.6, 3.7, 3.8,
           3.9, 4. , 4.1, 4.2, 4.3, 4.4, 4.5, 4.6, 4.7, 4.8, 4.9, 5. , 5.1,
           5.2, 5.3, 5.4, 5.5, 5.6, 5.7, 5.8, 5.9, 6. , 6.1, 6.2, 6.3, 6.4,
           6.5, 6.6, 6.7, 6.8, 6.9, 7. , 7.1, 7.2, 7.3, 7.4, 7.5, 7.6, 7.7,
           7.8, 7.9, 8. , 8.1, 8.2, 8.3, 8.4, 8.5, 8.6, 8.7, 8.8, 8.9, 9. ,
           9.1, 9.2, 9.3, 9.4, 9.5, 9.6, 9.7, 9.8, 9.9])



If we give a numpy array to a function, each term in the list is evaluated with that function:


```python
x = np.arange(1., 5., 1.)
print(x)
print(x**2)
print(np.sqrt(x))
```

    [1. 2. 3. 4.]
    [ 1.  4.  9. 16.]
    [1.         1.41421356 1.73205081 2.        ]


We can pick out elements of the list.  **Why does the last one fail?** 


```python
# List index is out range.
# List elements begin with indices at 0.
print(x[0])
print(x[3])
print(x[2])
```

    1.0
    4.0
    3.0


An alternative to `np.arange` is `np.linspace(min, max, number)` to get an array from `min` to `max` with `number` elements.  Example: 


```python
u_pts = np.linspace(0., 10., 10)
u_pts
```




    array([ 0.        ,  1.11111111,  2.22222222,  3.33333333,  4.44444444,
            5.55555556,  6.66666667,  7.77777778,  8.88888889, 10.        ])



**Change the last statement for `u_pts` to get an array from 0 to 10 spaced by 1.**

## Getting help

You will often need help identifying the appropriate Python (or numpy or scipy or ...) command or you will need an example of how to do something or you may get an error message you can't figure out.  In all of these cases, Google (or equivalent) is your friend. Always include "python" in the search string (or "numpy" or "matplotlib" or ...) to avoid getting results for a different language. You will usually get an online manual as one of the first responses if you ask about a function; these usually have examples if you scroll down. Otherwise, answers from *Stack Overflow* queries are your best bet to find a useful answer.

## Functions

There are many Python language features that we will use eventually, but in the short term what we need first are functions.  Here we first see the role of *indentation* in Python in place of {}s or ()s in other languages.  We'll always indent four spaces (never tabs!).  We know a function definition is complete when the indentation stops. 

To find out about a Python function or one you define, put your cursor on the function name and hit shift+Tab+Tab. **Go back and try it on `np.arange`.**  


```python
# Use "def" to create new functions.  
#  Note the colon and indentation (4 spaces).
def my_function(x):
    """This function squares the input.  Always include a brief description
        at the top between three starting and three ending quotes.  We will
        talk more about proper documentation later.
        Try shift+Tab+Tab after you have evaluated this function.
    """
    return x**2

print(my_function(5.))

# We can pass an array to the function and it is evaluated term-by-term.
x_pts = np.arange(1., 10., 1.)
print(my_function(x_pts))
```

    25.0
    [ 1.  4.  9. 16. 25. 36. 49. 64. 81.]



```python
# Two variables, with a default for the second
def add(x, y=4.):
    """Add two numbers."""
    print("x is {} and y is {}".format(x, y))
    return x + y  # Return values with a return statement

# Calling functions with parameters
print('The sum is ', add(5, 6))  # => prints out "x is 5 and y is 6" and returns 11

# Another way to call functions is with keyword arguments
add(y=6, x=5)  # Keyword arguments can arrive in any order.

```

    x is 5 and y is 6
    The sum is  11
    x is 5 and y is 6





    11



**How do you explain the following result?**


```python
add(2)
```

    x is 2 and y is 4.0





    6.0



### Debugging aside . . .

There are two bugs in the following function.  **Note the line where an error is first reported and fix the bugs sequentially (so you see the different error messages).**


```python
def hello_function():
    msg = "hello, world!"
    print(msg)
    return msg
```

## Plotting with Matplotlib

Matplotlib is the plotting library we'll use, at least at first.  We'll follow convention and abbreviate the module we need as `plt`.  The `%matplotlib inline` line tells the Jupyter notebook to make inline plots (we'll see other possibilities later).


```python
%matplotlib inline

import numpy as np
import matplotlib.pyplot as plt
```

Procedure we'll use to make the skeleton plot:
0. Generate some data to plot in the form of arrays.
1. Create a figure;
2. add one or more subplots;
3. make a plot and display it.


```python
t_pts = np.arange(0., 10., .1)     # step 0.
x_pts = np.sin(t_pts)  # More often this would be from a function 
                       #  *we* write.

my_fig = plt.figure()              # step 1.
my_ax = my_fig.add_subplot(1,1,1)  # step 2: rows=1, cols=1, 1st subplot
my_ax.plot(t_pts, x_pts)           # step 3: plot x vs. t
```

NOTE: When making just a single plot, you will more usually see steps 1 to 3 compressed into `plt.plot(t_pts, np.sin(t_pts))`.  *Don't do this.*  It saves a couple of lines but restricts your ability to easily extend the plot, which is what we want to make easy. 

We can always go back and dress up the plot:


```python
my_fig = plt.figure()
my_ax = my_fig.add_subplot(1,1,1)  # nrows=1, ncols=1, first plot
my_ax.plot(t_pts, x_pts, color='blue', linestyle='--', label='sine')

my_ax.set_xlabel('t')
my_ax.set_ylabel(r'$\sin(t)$')  # use $...$ to get LaTeX and r to render it
my_ax.set_title('Sine wave')

# here we'll put the function in the call to plot!
my_ax.plot(t_pts, np.cos(t_pts), label='cosine')  # just label the plot

my_ax.legend();  # turn on legend (try removing the semicolon)

```

Now make two subplots:


```python
y_pts = np.exp(t_pts)         # another function for a separate plot

fig = plt.figure(figsize=(10,5))  # allow more room for two subplots

# call the first axis ax1
ax1 = fig.add_subplot(1,2,1)  # one row, two columns, first plot
ax1.plot(t_pts, x_pts, color='blue', linestyle='--', label='sine')
ax1.plot(t_pts, np.cos(t_pts), label='cosine')  # just label the plot
ax1.legend()

ax2 = fig.add_subplot(1,2,2)  # one row, two columns, second plot
ax2.plot(t_pts, np.exp(t_pts), label='exponential')  
ax2.legend();

```

## Solving ordinary differential equations (ODEs)

Newton's 2nd Law for one particle:

$\begin{align}
  \mathbf{F} = m\mathbf{a} = m\frac{d\mathbf{v}}{dt} = m\frac{d^2\mathbf{x}}{dt^2}
\end{align}$ 

is a first-order differential equation in the velocity vector and a second-order differential equation in the position vector. Here we assume that $\mathbf{F} = \mathbf{F}(\mathbf{x}, \mathbf{v}, t)$ does not have higher derivatives. So we need to know how to solve such differential equations, usually with initial conditions (as opposed to boundary conditions). 

Here is how to do it with the Scipy function `odeint`.  (Note that `solve_ivp` is now preferred to `odeint`; we'll switch to that later but it is easier to get started with `odeint`.)

### First-order ODE


```python
# Import the required modules
import numpy as np
import matplotlib.pyplot as plt

from scipy.integrate import odeint  # Get only odeint from scipy.integrate
# for the future:
from scipy.integrate import solve_ivp   # Now preferred to odeint
```

Let's try a one-dimensional first-order ODE, say:

$\begin{align}
\quad 
\frac{dv}{dt} = -g, \quad \mbox{with} \quad v(0) = 10
\end{align}$

in some appropriate units (we'll use MKS units by default).  This ODE can be separated and directly integrated:

$\begin{align}
  \int_{v_0=10}^{v} dv' = - g \int_{0}^{t} dt'
  \quad\Longrightarrow\quad
    v - v_0 = - g (t - 0)
  \quad\Longrightarrow\quad
   v(t) = 10 - gt
\end{align}$




```python
# Define a function which calculates the derivative
def dv_dt(v, t, g=9.8):
    """Returns the right side of a simple first-order ODE with default g."""
    return -g   

t_pts = np.linspace(0., 10., 101)     # 101 points from t=0 and t=10.
v_0 = 10.0  # the initial condition
v_pts = odeint(dv_dt, v_0, t_pts)  # odeint( function for rhs, 
                                   #         initial value of v(t),
                                   #         array of t values )

```


```python

```

Let's check the output $v(t)$.  Does it make sense?


```python
print(v_pts[0:10])
v_pts.shape
```

    [[10.  ]
     [ 9.02]
     [ 8.04]
     [ 7.06]
     [ 6.08]
     [ 5.1 ]
     [ 4.12]
     [ 3.14]
     [ 2.16]
     [ 1.18]]





    (101, 1)



Now plot $v(t)$ calculated numerically and the analytic solution:


```python
g = 9.8
v_pts_exact = v_0 - g*t_pts

fig = plt.figure(figsize=(8,4))

ax = fig.add_subplot(1,2,1)
ax.plot(t_pts, v_pts, label='numerical', color='blue', lw=4)
ax.plot(t_pts, v_pts_exact, label='exact', color='red')
ax.set_xlabel('t [sec]')
ax.set_ylabel('v(t) [meters/sec]')
ax.legend();

diff = v_pts_exact - v_pts.flatten()
ax2 = fig.add_subplot(1,2,2)
ax2.semilogy(t_pts, diff, label='difference', color='red')
ax2.set_xlabel('t [sec]')
ax2.set_ylabel('v(t) [meters/sec]')
ax2.legend();

fig.tight_layout()
```

We adjusted the colors and linewidths so you could see that the lines are on top of each other. **Why are the differences around $10^{-15}$?**

Try solving instead:

$\begin{align}
\quad 
\frac{dx}{dt} = -x \quad \mbox{with} \quad x(0) = 1
\end{align}$

both analytically and numerically.


```python
def dx_dt(x, t):
    """Function returns the right side of the above ODE."""
    return -x

t_pts = np.linspace(0, 10, 101)
x_0 = 1
x_pts = odeint(dx_dt, x_0, t_pts)

print(x_pts[0:10])
print(x_pts.shape)
```

    [[1.        ]
     [0.90483745]
     [0.81873077]
     [0.7408182 ]
     [0.67032006]
     [0.60653067]
     [0.54881165]
     [0.49658532]
     [0.44932898]
     [0.40656968]]
    (101, 1)


### Second-order ODE

Suppose we have a second-order ODE such as:

$$
\quad y'' + 2 y' + 2 y = \cos(2x), \quad \quad y(0) = 0, \; y'(0) = 0
$$

We can turn this into two first-order equations by defining a new dependent variable equal to $y'$ (we'll call it $z$):

$$
\quad z \equiv y' \quad \Rightarrow \quad z' + 2 z + 2y = \cos(2x) 
  \quad\mbox{or}\quad z' = -2z -2y + \cos(2x), \quad\mbox{with}
\quad z(0)=y(0) = 0.
$$

Now introduce the vector 

$$
  \mathbf{U}(x) = \left(\begin{array}{c}
                         y(x) \\
                         z(x)
                        \end{array}
                  \right)
        \quad\Longrightarrow\quad
    \frac{d\mathbf{U}}{dx} = 
                             \left(\begin{array}{c}
                                    y' \\
                                    z'
                                   \end{array}
                             \right) 
                           =  
                             \left(\begin{array}{c}
                                    z \\
                                    -2 y' - 2 y + \cos(2x)
                                   \end{array}
                             \right)
              \quad\mbox{with}\quad
    \mathbf{U}(0) = \left(\begin{array}{c}
                         0 \\
                         0
                        \end{array}
                  \right)
$$

We can solve this system of ODEs using `odeint` with lists, as follows:


```python
# Define a function for the right side
def dU_dx(U, x):
    """Right side of the differential equation to be solved.
    U is a two-component vector with y=U[0] and z=U[1]. 
    Thus this function should return [y', z']
    """
    return [U[1], -2*U[1] - 2*U[0] + np.cos(2*x)]

# initial condition U_0 = [y(0)=0, z(0)=y'(0)=0]
U_0 = [0., 0.]

x_pts = np.linspace(0, 15, 200)  # Set up the mesh of x points
U_pts = odeint(dU_dx, U_0, x_pts)  # U_pts is a 2-dimensional array
y_pts = U_pts[:, 0]   # Ok, this is tricky.  For each x, U_pts has two 
                      #  components.  We want the upper component for all
                      #  x, which is y(x).  The : means all of the first 
                      #  index, which is x, and the 0 means the first
                      #  component in the other dimension.
                      # What is U_pts[:, 1]? (Try plotting it.)
```


```python
fig = plt.figure()
ax = fig.add_subplot(1,1,1)
ax.plot(x_pts, y_pts)
ax.set_xlabel('x')
ax.set_ylabel('y')
```

## Widgets!

A widget is an object such as a slider or a check box or a pulldown menu.  We can use them to make it easy to explore different parameter values in a problem we're solving, which is invaluable for building intuition.  They act on the argument of a function.  We'll look at a simple case here but plan to explore this much more as we proceed.

The set of widgets we'll use here (there are others!) is from `ipywidgets`; we'll conventionally import the module as `import ipywidgets as widgets` and we'll also often use `display` from `Ipython.display`.


```python
import numpy as np
import matplotlib.pyplot as plt

import ipywidgets as widgets
from IPython.display import display

%matplotlib inline
```

The simplest form is to use `interact`, which we pass a function name and the variables with ranges.  By default this makes a *slider*, which takes on integer or floating point values depending on whether you put decimal points in the range. **Try it! Then modify the function and try again.**


```python
# We can do this to any function
def test_f(x=5.):
    """Test function that prints the passed value and its square.
       Note that there is no return value in this case."""
    print ('x = ', x, ' and  x^2 = ', x**2)
    
widgets.interact(test_f, x=(0.,10.));
```


    interactive(children=(FloatSlider(value=5.0, description='x', max=10.0), Output()), _dom_classes=('widget-inte…



```python
# Explicit declaration of the widget (here FloatSlider) and details
def test_f(x=5.):
    """Test function that prints the passed value and its square.
       Note that there is no return value in this case."""
    print ('x = ', x, ' and  x^2 = ', x**2)
    
widgets.interact(test_f, 
                 x = widgets.FloatSlider(min=-10,max=30,step=1,value=10));
```


    interactive(children=(FloatSlider(value=10.0, description='x', max=30.0, min=-10.0, step=1.0), Output()), _dom…


Here's an example with some bells and whistles for a plot.  **Try making changes!**


```python
def plot_it(freq=1., color='blue', lw=2, grid=True, xlabel='x', 
            function='sin'):
    """ Make a simple plot of a trig function but allow the plot style
        to be changed as well as the function and frequency."""
    t = np.linspace(-1., +1., 1000)  # linspace(min, max, total #)

    fig = plt.figure(figsize=(8,6))
    ax = fig.add_subplot(1,1,1)

    if function=='sin':
        ax.plot(t, np.sin(2*np.pi*freq*t), lw=lw, color=color)
    elif function=='cos':
        ax.plot(t, np.cos(2*np.pi*freq*t), lw=lw, color=color)
    elif function=='tanh':
        ax.plot(t, np.tanh(2*np.pi*freq*t), lw=lw, color=color)

    ax.grid(grid)
    ax.set_xlabel(xlabel)
    
widgets.interact(plot_it, 
                 freq=(0.1, 60.), color=['blue', 'red', 'green'], 
                 lw=(1, 10), xlabel=['x', 't', 'dog'],
                 function=['sin', 'cos', 'tanh'])
    
```


    interactive(children=(FloatSlider(value=1.0, description='freq', max=60.0, min=0.1), Dropdown(description='col…





    <function __main__.plot_it(freq=1.0, color='blue', lw=2, grid=True, xlabel='x', function='sin')>




```python

```

dd7351cf6a6385f63bb379e80e36f1157769ef5e
ipynb
Jupyter Notebook
2020_week_1/5300_Jupyter_Python_intro_01.ipynb
CLima86/Physics_5300_CDL
d9e8ee0861d408a85b4be3adfc97e98afb4a1149
MIT
2020_week_1/5300_Jupyter_Python_intro_01.ipynb
CLima86/Physics_5300_CDL
d9e8ee0861d408a85b4be3adfc97e98afb4a1149
MIT
2020_week_1/5300_Jupyter_Python_intro_01.ipynb
CLima86/Physics_5300_CDL
d9e8ee0861d408a85b4be3adfc97e98afb4a1149
MIT
Qwen/Qwen-72B
1. YES
2. YES