# 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
t, tau, v0 = symbols("t tau v0")
x = Function("x")
```

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$.

4c9f19115750c2ebd7c52eb01ce2336d846341ce
ipynb
Jupyter Notebook
Lecture01.ipynb
garth-wells/IA-maths-Ipython
7497ee12f2a36d6ec6025ccdb3d2ce6387dc7c9a
MIT
2017-05-21T16:30:35.000Z
2021-12-14T00:23:51.000Z
Lecture01.ipynb
garth-wells/IA-maths-Ipython
7497ee12f2a36d6ec6025ccdb3d2ce6387dc7c9a
MIT
2016-10-10T09:06:41.000Z
2016-11-04T11:05:32.000Z
Lecture01.ipynb
garth-wells/IA-maths-Jupyter
7497ee12f2a36d6ec6025ccdb3d2ce6387dc7c9a
MIT
2017-05-15T08:22:51.000Z
2021-11-02T05:27:59.000Z
Qwen/Qwen-72B
1. YES
2. YES