<jupyter_start><jupyter_text>Lab 8Last week we used the improved Euler method to numerically solve a DE of the form\begin{align*}\frac{\mathrm{d}\mathbf{y}}{\mathrm{d}t} = f(\mathbf{y}, t),\end{align*}where $\mathbf{y}$ is a vector of dependent variables and $f$ is a linear function of $\mathbf{y}$.This week we will use the SciPy `odeint` function in place of our implmentation. There are several advantages to using SciPy code, including that it is well-tested, provides a greater range of features, and uses more advanced techniques. One disadvantage is that it makes our code dependendent on SciPy, but it's worth it.We will also show how you can approximate solutions to systems of DEs that you may not be able to solve analytically. Setup<jupyter_code>from numpy import arange, empty, exp, array, linspace, isclose, allclose, sin, pi
from scipy.integrate import odeint
from scipy.optimize import minimize
from plotly.offline import init_notebook_mode
from plotly import graph_objs as go<jupyter_output><empty_output><jupyter_text>A Nonlinear DEImagine that the number of rabbits in an area as a function of time is modelled by $r(t)$ and the number of foxes in the same area at the same time is given by $f(t)$. We can modelled their linked fates using the coupled DEs\begin{align*}\frac{\mathrm{d}r}{\mathrm{d}t} &= 0.08 r - 0.0004 rf, \\\frac{\mathrm{d}f}{\mathrm{d}t} &= -0.09 f + 0.0005 rf,\end{align*}with initial conditions $r_0 = 300$ and $f_0 = 100$. See Section 6.3 (and the rest of Chapter 6) if you are curious about where this equation comes from.For us, it is important to note that- this system of DEs cannot be represented as matrix multiplied by a $(r\, f)^\intercal$ vector and- we have not studied how to solve this equation.We can still define it as a Python function, however. Note that in the following function `rf` is a variable that contains $r$ and $f$ as a vector (actually a NumPy `array`, but it is one-dimensional so we call it a vector).<jupyter_code>def drfdt(rf, t):
    r, f = rf
    drdt = 0.08*r - 0.0004*r*f
    dfdt = -0.09*f + 0.0005*r*f
    return array([drdt, dfdt])<jupyter_output><empty_output><jupyter_text>`odeint` works almost exactly the same as `euler_improved` did, but it takes the values of $t$ for which you would like $r$ and $f$ as inputs. Note that the initial values are input as `[300, 100]`, and that the first element of $t$ must correspond to the time of those initial values (in this case `t[0] == 0`).<jupyter_code>t = linspace(0, 100, 100)
rf = odeint(drfdt, [300, 100], t)
fig = go.Figure()
fig.add_trace(go.Scatter(x=t, y=rf[:,0], name='rabbits'))
fig.add_trace(go.Scatter(x=t, y=rf[:,1], name='foxes'))
fig.show('png')<jupyter_output><empty_output><jupyter_text>You do not have to specify step sizes for `odeint` (like you did for `euler_improved`). It figures that out for itself. In fact, if you want to know how many rabbits and foxes that you will have after 40 years, you can call it with just that input:<jupyter_code>rf = odeint(drfdt, [300, 100], [0, 40])
print(rf[-1, 0], 'rabbits')
print(rf[-1, 1], 'foxes')<jupyter_output>70.757766160244 rabbits
208.23878627992096 foxes<jupyter_text>If you wanted to figure out exactly when the number of rabbits reaches its lowest ebb, you can minimize the result numerically, again using SciPy.<jupyter_code>def rabbits(t):
    return odeint(drfdt, [300, 100], [0, t])[-1, 0]

result = minimize(rabbits, 40)
print(result)<jupyter_output>fun: 70.67226332711697
 hess_inv: array([[3.23550154]])
      jac: array([-3.81469727e-06])
  message: 'Optimization terminated successfully.'
     nfev: 10
      nit: 4
     njev: 5
   status: 0
  success: True
        x: array([40.73864735])<jupyter_text>You get a heap of diagnostic information from `minimize`, but the most important thing is that "Optimization terminated successfully" and that the value is stored in `x`:<jupyter_code>print('rabbits rebounded after', result.x[0], 'years')<jupyter_output>rabbits rebounded after 40.738647346149506 years<jupyter_text>ExercisesWe will now numerically solve a differential equation that we cannot solve analytically. Pendulum motion can be defined as\begin{align*}\frac{\mathrm{d}^2\theta}{\mathrm{d}t^2} = -\frac{g}{L}\sin\theta,\end{align*}where $g$ is gravitational acceleration, $L$ is the length of the pendulum, and $\theta$ is theangle the pendulum makes with the vertical as shown in the figure.Figure 1: A simple pendulumThe above equation is a second-order nonlinear differential equation and we don’t have away to solve this equation algebraically. That is, we can’t use the characteristic equationmethod or method of undetermined coefficients to solve this equation.We can, however, convert it into a system of first-order DEs and find an approximate solution using `odeint`. By setting $v=\mathrm{d}\theta/\mathrm{d}t$, we obtain the equivalent formulation\begin{align*}\frac{\mathrm{d}\theta}{\mathrm{d}t} &= v \\\frac{\mathrm{d}v}{\mathrm{d}t} &= -\frac{g}{L}\sin\theta.\end{align*}Suppose that $g/L = 10$. Write a function to calculate the vector$\left(\begin{array}{1}\frac{\mathrm{d}\theta}{\mathrm{d}t}\\\frac{\mathrm{d}v}{\mathrm{d}t}\end{array}\right)$ as a function of $\left(\begin{array}{1}\theta\\v\end{array}\right)$ and $t$. (Note that NumPy does not care whether you return a column or row vector. Usually it figures out what you mean from the context.)<jupyter_code>def dthetavdt(thetav, t):
    z, v = thetav
    dzdt = v
    dvdt = -10*sin(z)
    return array([dzdt, dvdt])<jupyter_output><empty_output><jupyter_text>The following cell should run without errors.<jupyter_code>assert allclose(dthetavdt([pi/2, 0], 0), [0, -10.0])
assert allclose(dthetavdt([pi/12, 0], 0), [0, -2.5881904510252074])
assert allclose(dthetavdt([0, 1], 0), [1, 0])<jupyter_output><empty_output><jupyter_text>Use `odeint` to plot $\theta$ on the interval $0\leq t \leq 2.5$ when $\theta_0=\frac{\pi}{12}$ and $v_0=0$.<jupyter_code>t = linspace(0, 2.5, 100)
thetav = odeint(dthetavdt, [pi/12, 0], t)
fig = go.Figure()
fig.add_trace(go.Scatter(x=t, y=thetav[:,0], name='Theta'))
fig.add_trace(go.Scatter(x=t, y=thetav[:,1], name='Velocity'))
fig.show('png')<jupyter_output><empty_output><jupyter_text>In the following cell, calculate and print the values for $\theta(2.5)$ and $v(2.5)$ when $\theta_0=\frac{\pi}{12}$ and $v_0=0$.<jupyter_code>thetav = odeint(dthetavdt, [pi/12, 0], [0, 2.5])
print(thetav[-1, 0], 'Theta')
print(thetav[-1, 1], 'Velocity')<jupyter_output>-0.004681934744160246 Theta
-0.8253870220740621 Velocity<jupyter_text>Plot $\theta$ on the interval $0\leq t \leq 2.5$ when $\theta_0=\frac{\pi}{2}$ and $v_0=0$.<jupyter_code>t = linspace(0, 2.5, 100)
thetav = odeint(dthetavdt, [pi/2, 0], t)
fig = go.Figure()
fig.add_trace(go.Scatter(x=t, y=thetav[:,0], name='Theta'))
fig.add_trace(go.Scatter(x=t, y=thetav[:,1], name='Velocity'))
fig.show('png')<jupyter_output><empty_output><jupyter_text>In the following cell, calculate and print the values for $\theta(2.5)$ and $v(2.5)$ when $\theta_0=\frac{\pi}{2}$ and $v_0=0$.<jupyter_code>thetav = odeint(dthetavdt, [pi/2, 0], [0, 2.5])
print(thetav[-1, 0], 'Theta')
print(thetav[-1, 1], 'Velocity')<jupyter_output>1.4510997040393572 Theta
-1.5453873087645622 Velocity