### Solution


```python
# Special command for plotting in jupyter notebooks
%matplotlib inline
```

First import the required modules - `numpy` and `matplotlib`. 


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

Initially, it is good to have a look at our reference solution. Based on the given formula for the trajectory of the partical, we can work out the formula for its velocity as:

\begin{align}
v_x &= -\dfrac{1}{T_1}\exp \left( -\dfrac{t}{T_1} \right) 
\\
v_y &= -\dfrac{2}{T_2} \sin \left(\dfrac{t}{T_2}\right)
\end{align}

and acceleration as:

\begin{align}
a_x &= \dfrac{1}{T^2_1} \exp \left( -\dfrac{t}{T_1} \right) 
\\
a_y &= -\dfrac{2}{T^2_2} \cos\left(\dfrac{t}{T_2}\right)
\end{align}

We declare our constants as followed.


```python
# Declare the time constants
T1 = 1.0
T2 = 1.0/3.0
```

Then initialise the array to store our reference results. Here, our positions, velocities and accelerations of the particles in time are stored as 2D array of vectors - the first column corresponds to the $x$ component, and the second the $y$ component. 


```python
# Initialise array for position, velocity and acceleration
num_points = 500
r_ref = np.zeros((num_points, 2))
v_ref = np.zeros((num_points, 2))
a_ref = np.zeros((num_points, 2))

# Generate the reference time vector
t_ref = np.linspace(0, 2*np.pi, num_points)
```

The results can be calculated using the formula provided above.


```python
# Calculate the x and y coordinates of the particle
r_ref[:, 0] = np.exp(-t_ref/T1)
r_ref[:, 1] = 2.0*np.cos(t_ref/T2)

# Exact solution for velocity
v_ref[:, 0] = -(1.0/T1)*np.exp(-t_ref/T1)
v_ref[:, 1] = -(2.0/T2)*np.sin(t_ref/T2)

# Exact solution for acceleration
a_ref[:, 0] = (1.0/(T1**2))*np.exp(-t_ref/T1)
a_ref[:, 1] = -(2.0/(T2**2))*np.cos(t_ref/T2)
```

We can now plot the postion, velocity and acceleration fields.


```python
# Plot the solutions
# Initilise the figure
fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(24, 6))

# Plot the position
ax1.plot(r_ref[:, 0], r_ref[:, 1], linewidth=2.0)
ax1.set_title('Position')
ax1.set_xlabel('x')
ax1.set_ylabel('y')

# Plot the velocity
ax2.plot(v_ref[:, 0], v_ref[:, 1], linewidth=2.0)
ax2.set_title('Velocity')
ax2.set_xlabel('x')
ax2.set_ylabel('y')

# Plot the acceleration
ax3.plot(a_ref[:, 0], a_ref[:, 1], linewidth=2.0)
ax3.set_title('Acceleration')
ax3.set_xlabel('x')
ax3.set_ylabel('y');
```

### Numerical differentiation

#### Forward difference method

Now we assume that we only have a finite number of points for the location of the particle, `r`, in time. For this solution, we assume that we only have 20 data points, i.e. `num_points = 20`, for the particle's trajectory in time interval $t = 0 - 2 \pi$s.


```python
# Generate our finite time vector, with 20 points
num_points = 20
t = np.linspace(0, 2*np.pi, num_points)

# Calculate the position used for the numerical term
r_num = np.zeros((num_points, 2))
r_num[:, 0] = np.exp(-t/T1)
r_num[:, 1] = 2.0*np.cos(t/T2)
```

To calculate the particles velocity and acceleration we need to calculate the first and second derivatives of its position, `r`, with respect to time. Reminds ourselve that here we only have discrete number of points, rather than a continuous function. Hence, a numerical scheme is required.

Python offers the `np.diff()` comand to calculate the difference between consecutive elements in a `numpy` array. The resulting vector is one element shorter than the original, e.g. if `x=[1, 4, 3, 5]`, then `np.diff(x)` returns `[3, -1, 2]`. 

With `dx = np.diff(x)` and `dt = np.diff(t)` the division `dx/dt` gives an approximation of the particle’s velocity in the x direction. The second derivative can be calculated using: 

    np.diff(x,2)/(np.power(dt[0:num_points - 2], 2)

where `np.diff(x,2)` is the second order difference of `x` – using the above example `np.diff(x,2)` returns `[-4, 3]`. 

For our case, we need to apply the differentiation for the $x$ and $y$ direction of the position, which are stored in the first and the second column of `r`, respectively.


```python
# Initialise array for velocity and acceleration
# calculated using forward difference method
v_f = np.zeros((num_points-1, 2))
a_f = np.zeros((num_points-2, 2))

# Forward difference method to calculate velocity
v_f[:, 0] = np.diff(r_num[:, 0]) / np.diff(t)
v_f[:, 1] = np.diff(r_num[:, 1]) / np.diff(t)
# And then the acceleration
# Calculate acceleration
a_f[:, 0] = np.diff(r_num[:, 0], 2) / (np.power(np.diff(t)[0:num_points-2], 2))
a_f[:, 1] = np.diff(r_num[:, 1], 2) / (np.power(np.diff(t)[0:num_points-2], 2))
```

We can look at our forward difference scheme, with 20 data points, compared to the reference - 'smooth' - result. 


```python
# Generate the subplots
fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(24, 6))

# Plot the position
ax1.plot(r_num[:, 0], r_num[:, 1], linewidth=2.0, label='Forward')
ax1.plot(r_ref[:, 0], r_ref[:, 1], linewidth=2.0, label='Reference')
ax1.set_title('Position')
ax1.legend()

# Plot the velocity
ax2.plot(v_f[:, 0], v_f[:, 1], linewidth=2.0, label='Forward')
ax2.plot(v_ref[:, 0], v_ref[:, 1], linewidth=2.0, label='Reference')
ax2.set_title('Velocity')
ax2.legend(loc=2)

# Plot the acceleration
ax3.plot(a_f[:, 0], a_f[:, 1], linewidth=2.0, label='Forward')
ax3.plot(a_ref[:, 0], a_ref[:, 1], linewidth=2.0, label='Reference')
ax3.set_title('Acceleration')
ax3.legend();
```

#### Central difference method

In the forward difference method, the current velocity is calculated using the current and the following position. The current acceleration is calculated using the previous and next velocity. This leads to a *mismatch* between velocity on one hand and position and acceleration on the other by one timestep. The centre difference method allows to calculate velocities and accelerations at the same positions over the interval `[1, N-1]`. This method is implemented below, and full understanding is not required. 


```python
# Initialise array for velocity and acceleration
# calculated using central difference method
v_c = np.zeros((num_points-1, 2))
a_c = np.zeros((num_points-2, 2))

# Positions and time difference between adjacent
# elements of position and time. The argument axis=0
# indicates that we carry out np.diff on column-wise 
# of r
dr = np.diff(r_num, axis=0)
dt = np.diff(t)

# Centre difference method
for i in range(1, num_points-1):
    v_c[i-1, :] = (dr[i, :] + dr[i-1, :])/(dt[i] + dt[i-1])
    a_c[i-1, :] = (dr[i, :]/dt[i] - dr[i-1, :]/dt[i-1])/(0.5*(dt[i] + dt[i-1]))
```

We can quantitatively examine the difference between the forward and central difference scheme to calculate the velocity and the acceleration. Again, here we stick to only 20 data points as before.


```python
# Generate the subplots
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(16, 6))

# Plot the velocity
ax1.plot(v_ref[:, 0], v_ref[:, 1], label='Reference', lw=3.0, ls='--')
ax1.plot(v_f[:, 0], v_f[:, 1], label='Forward')
ax1.plot(v_c[:, 0], v_c[:, 1], label='Central')
ax1.set_title('Velocity')
ax1.legend(loc=2)

# Plot the acceleration
ax2.plot(a_ref[:, 0], a_ref[:, 1], label='Reference', lw=3.0, ls='--')
ax2.plot(a_f[:, 0], a_f[:, 1], label='Forward')
ax2.plot(a_c[:, 0], a_c[:, 1], label='Central')
ax2.set_title('Acceleration')
ax2.legend();
```

#### A programmer approach - with interactive!

We can define functions that do our calculations for the forward and central difference method as followed. 


```python
def forward_method(r, t):
    n = len(t)
    
    # Initialise array for velocity and acceleration
    # calculated using forward difference method
    v_for = np.zeros((n-1, 2))
    a_for = np.zeros((n-2, 2))

    # Forward difference method to calculate velocity
    v_for[:, 0] = np.diff(r[:, 0]) / np.diff(t)
    v_for[:, 1] = np.diff(r[:, 1]) / np.diff(t)
    
    # And then the acceleration
    # Calculate acceleration
    a_for[:, 0] = np.diff(r[:, 0], 2) / (np.power(np.diff(t)[0:n-2], 2))
    a_for[:, 1] = np.diff(r[:, 1], 2) / (np.power(np.diff(t)[0:n-2], 2))

    return v_for, a_for
```


```python
def central_method(r, t):
    # Initialise array for velocity and acceleration
    # calculated using central difference method
    n = len(t)
    v_cen = np.zeros((n, 2))
    a_cen = np.zeros((n, 2))

    # Positions and time difference between adjacent
    # elements of position and time
    dr = np.diff(r, axis=0)
    dt = np.diff(t)

    # Centre difference method
    for i in range(1, n-2):
        v_cen[i-1, :] = (dr[i, :] + dr[i-1, :])/(dt[i] + dt[i-1])
        a_cen[i-1, :] = (dr[i, :]/dt[i] - dr[i-1, :]/dt[i-1]) / (0.5*(dt[i] + dt[i-1]))

    return v_cen, a_cen
```

By doing so, we can easily calculate the results achieved from either scheme. We can call these functions for the interactive widget, which illustrates the dependence of these two methods on the number of available data points. Full understanding of the code is not required, yet curious students are encouraged to explore. 


```python
# Interactive notebook to compare the two methods
from ipywidgets import *
```


```python
def compare_difference(num_points):
    t = np.linspace(0, 2*np.pi, num_points)

    # Calculate the position used for the numerical term
    r_num = np.zeros((num_points, 2))
    r_num[:, 0] = np.exp(-t/T1)
    r_num[:, 1] = 2.0*np.cos(t/T2)

    # Forward difference method
    v_for, a_for = forward_method(r_num, t)

    # Central difference metho
    v_cen, a_cen = central_method(r_num, t)

    # Generate the subplots
    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(16, 6))

    # Plot the velocity
    ax1.plot(v_ref[:, 0], v_ref[:, 1], label='Reference', lw=3.0, ls='--')
    ax1.plot(v_for[:, 0], v_for[:, 1], label='Forward')
    ax1.plot(v_cen[:, 0], v_cen[:, 1], label='Central')
    ax1.set_title('Velocity')
    ax1.legend(loc=2)

    # Plot the acceleration
    ax2.plot(a_ref[:, 0], a_ref[:, 1], label='Reference', lw=3.0, ls='--')
    ax2.plot(a_for[:, 0], a_for[:, 1], label='Forward')
    ax2.plot(a_cen[:, 0], a_cen[:, 1], label='Central')
    ax2.set_title('Acceleration')
    ax2.legend()
```


```python
interact(compare_difference, num_points=IntSlider(
        min=10, max=200, step=10, value=20, description='No. points'));
```

# Lecture 1: solving ordinary differential equations

This lecture introduces ordinary differential equations, and some techniques for solving first order equations. This notebook uses computer algebra via [Sympy](<http://sympy.org/>) to solve some ODE examples from the lecture notes.

# Importing SymPy

To use Sympy, we first need to import it and call `init_printing()` to get nicely typeset equations:


```python
import sympy
from sympy import symbols, Eq, Derivative, init_printing, Function, dsolve, exp, classify_ode, checkodesol

# This initialises pretty printing
init_printing()
from IPython.display import display

# Support for interactive plots
from ipywidgets import interact

# This command makes plots appear inside the browser window
%matplotlib inline
```

# Example: car breaking

During braking a car’s velocity is given by $v = v_{0} e^{−t/\tau}$. Calculate the distance travelled.

We first define the symbols in the equation ($t$, $\tau$ and $v_{0}$), and the function ($x$, for the displacement):


```python

```

Next, we define the differential equation, and print it to the screen for checking:


```python
eqn = Eq(Derivative(x(t), t), v0*exp(-t/(tau)))
display(eqn)
```

The `dsolve` function solves the differential equation symbolically:


```python
x = dsolve(eqn, x(t))
display(x)
```

where $C_{1}$ is a constant. As expected for a first-order equation, there is one constant.

SymPy is not yet very good at eliminating constants from initial conditions, so we will do this manually assuming that $x = 0$ and $t = 0$:


```python
x = x.subs('C1', v0*tau)
display(x)
```

Specifying a value for $v_{0}$, we create an interactive plot of $x$ as a function of the parameter $\tau$:


```python
x = x.subs(v0, 100)

def plot(τ=1.0):
    x1 = x.subs(tau, τ)

    # Plot position vs time
    sympy.plot(x1.args[1], (t, 0.0, 10.0), xlabel="time", ylabel="position");

interact(plot, τ=(0.0, 10, 0.2));
```


<p>Failed to display Jupyter Widget of type <code>interactive</code>.</p>
<p>
  If you're reading this message in Jupyter Notebook or JupyterLab, it may mean
  that the widgets JavaScript is still loading. If this message persists, it
  likely means that the widgets JavaScript library is either not installed or
  not enabled. See the <a href="https://ipywidgets.readthedocs.io/en/stable/user_install.html">Jupyter
  Widgets Documentation</a> for setup instructions.
</p>
<p>
  If you're reading this message in another notebook frontend (for example, a static
  rendering on GitHub or <a href="https://nbviewer.jupyter.org/">NBViewer</a>),
  it may mean that your frontend doesn't currently support widgets.
</p>



# Classification

We can ask SymPy to classify our ODE, e.g. show that it is first order):


```python
classify_ode(eqn)
```




    ('separable',
     '1st_exact',
     '1st_linear',
     'Bernoulli',
     '1st_power_series',
     'lie_group',
     'nth_linear_constant_coeff_undetermined_coefficients',
     'nth_linear_constant_coeff_variation_of_parameters',
     'separable_Integral',
     '1st_exact_Integral',
     '1st_linear_Integral',
     'Bernoulli_Integral',
     'nth_linear_constant_coeff_variation_of_parameters_Integral')



# Parachutist

Find the variation of speed with time of a parachutist subject to a drag force of $kv^{2}$.

The equations to solve is

$$
\frac{m}{k} \frac{dv}{dt} = \alpha^{2} - v^{2}
$$

where $m$ is mass, $k$ is a prescribed constant, $v$ is the velocity, $t$ is time and $\alpha^{2} = mg/k$ ($g$ is acceleration due to gravity).

We specify the symbols, unknown function $v$ and the differential equation


```python
t, m, k, alpha = symbols("t m k alpha")
v = Function("v")
eqn = Eq((m/k)*Derivative(v(t), t), alpha*alpha - v(t)*v(t))
display(eqn)
```

First, let's classify the ODE:


```python
classify_ode(eqn)
```




    ('separable', '1st_power_series', 'lie_group', 'separable_Integral')



We see that it is not linear, but it is separable. Using `dsolve` again,


```python
v = dsolve(eqn, v(t))
display(v)
```

SymPy can verify that an expression is a solution to an ODE:


```python
print("Is v a solution to the ODE: {}".format(checkodesol(eqn, v)))
```

    Is v a solution to the ODE: (True, 0)


Try adding the code to plot velocity $v$ against time $t$.


```python

```


```python

```

4ca4868b6192c277863ac587523769fda0ef7cc0
ipynb
Jupyter Notebook
.ipynb_checkpoints 2/Assignment solution-checkpoint.ipynb
hphilamore/ILAS_PyEv2019
b3c8ebe00d6795f67879a50ce6ef517353b069c1
MIT
.ipynb_checkpoints 2/Assignment solution-checkpoint.ipynb
hphilamore/ILAS_PyEv2019
b3c8ebe00d6795f67879a50ce6ef517353b069c1
MIT
.ipynb_checkpoints 2/Assignment solution-checkpoint.ipynb
hphilamore/ILAS_PyEv2019
b3c8ebe00d6795f67879a50ce6ef517353b069c1
MIT
Qwen/Qwen-72B
1. YES
2. YES