ODE questions

Euler’s method

  1. Use Euler’s method to solve the differential equation

\[\begin{equation} \frac{\mathrm{d}x}{\mathrm{d}t} = -x^3 + \mathrm{sin} t \end{equation}\]

with the initial condition \(x=0\) at \(t=0\).

  1. Plot \(x(t)\) up to \(t=100\), using 1000 steps.

Euler’s method II

  1. Use Euler’s method to solve the differential equation

\[\begin{equation} \frac{\mathrm{d}y}{\mathrm{d}x} = \mathrm{cos}(x)\mathrm{sin}(y) + \mathrm{cos}(y) \end{equation}\]

with the initial condition \(y=1\) at \(x=0\).

  1. Plot \(y(x)\) up to \(x=100\), using 1000 steps.

Second order Runge-Kutta method

  1. Use a second order Runge-Kutta method to solve the differential equation

\[\begin{equation} \frac{\mathrm{d}x}{\mathrm{d}t} = -x^3 + \mathrm{sin} t \end{equation}\]

with the initial condition \(x=0\) at \(t=0\).

  1. Plot \(x(t)\) up to \(t=100\), using 1000 steps. Investigate convergence with respect to the number of steps, and compare this to the performance of Euler’s method.

Modelling a non-linear pendulum

In the multivariable equations section we use Euler’s method to solve simultaneous first order ODEs. At the end of the tutorial you are also shown how to re-cast the second order ODE for a non-linear pendulum as two simultaneous first order ODEs. Using the given expressions on the tutorial page, and following the same method outlined as that outlined for the Strange Attractor, write a piece of code for modelling \(\theta\) as a function of time.

Tip: Combine the two variable \(\theta\) and \(\omega\) into a single vector \(\mathbf{r} = (\theta,\omega)\). The method will give us a value for \(\theta\) and \(\omega\), but we can ignore the \(\omega\) values generated if they are not needed.

The fourth order Runge-Kutta method

We can extend the approach outlined in [the Runge-Kutta tutorial] to higher orders. The most common method for the numerical solution of ODEs is the fourth-order Runge-Kutta method (RK4). This is accurate to terms of order \(h^4\) (equivalently, it carries an error or order \(h^5\)). It’s derivation is quite tedious, but the approach is the same as that for the second-order method: we Taylor expand around various points and take linear combinations that lead to the cancellation of terms in \(h^3\), \(h^4\) and so on.

The five resulting equations are:

\[\begin{equation} k_1 = hf(x,t) \end{equation}\]

\[\begin{equation} k_2 = hf(x+\frac{1}{2}k_1,t+\frac{1}{2}h) \end{equation}\]

\[\begin{equation} k_3 = hf(x+\frac{1}{2}k_2,t+\frac{1}{2}h) \end{equation}\]

\[\begin{equation} k_4 = hf(x+k_3,t+h) \end{equation}\]

\[\begin{equation} x(t+h) = x(t) + \frac{1}{6}(k_1+2k_2+2k_3+k_4) \end{equation}\]

  1. Use RK4 to solve the differential equation

\[\begin{equation} \frac{\mathrm{d}x}{\mathrm{d}t} = -x^3 + \mathrm{sin} t \end{equation}\]

with the initial condition \(x=0\) at \(t=0\).

  1. Plot \(x(t)\) for a number of steps \(N\) equal to 10, 20, 50 and 1000.

  2. How does convergence compare to that when using RK2? You may have already applied the RK2 method to this problem, as it was set as a question at the end of the tutorial.

Dynamical systems

The motion of two bodies of mass \(m_1\) and \(m_2\) attracted by gravity is the Kepler problem. In this case the moving bodies experience a force varying as the inverse square of the distance between them; for body one the distance to body two is

\[\begin{equation} \mathbf{r}_{12} = \mathbf{r}_1-\mathbf{r}_2, \end{equation}\]

whilst for body two the distance to body one is

\[\begin{equation} \mathbf{r}_{21} = \mathbf{r}_2-\mathbf{r}_1. \end{equation}\]

So that the corresponding forces are:

\[\begin{equation} \mathbf{F}_1 = -\frac{Gm_1m_2}{r_{12}}\hat{\mathbf{r_{12}}} \end{equation}\]

\[\begin{equation} \mathbf{F}_2 = -\frac{Gm_1m_2}{r_{21}}\hat{\mathbf{r_{21}}} \end{equation}\]

  1. Numerically compute and plot the orbits by integrating the equations of motion:

\[\begin{equation} m_1\mathbf{a}_1 = \mathbf{F}_1, \end{equation}\]

\[\begin{equation} m_2\mathbf{a}_2 = \mathbf{F}_2. \end{equation}\]

Use the Euler method for the integration:

\[\begin{equation} \mathrm{v}_1(t+\Delta t)= \mathrm{v}_1(t)+\Delta t \frac{\mathrm{F_1}}{m_1} \end{equation}\]

\[\begin{equation} \mathrm{v}_2(t+\Delta t)= \mathrm{v}_2(t)+\Delta t \frac{\mathrm{F_2}}{m_2} \end{equation}\]

\[\begin{equation} \mathrm{r}_1(t+\Delta t)= \mathrm{r}_1(t)+\Delta t \mathrm{v}_1(t) \end{equation}\]

\[\begin{equation} \mathrm{r}_2(t+\Delta t)= \mathrm{r}_2(t)+\Delta t \mathrm{v}_2(t) \end{equation}\]

Assume that the \(z\) axis is perpendicular to the plane of motion so that the motion is in two dimensions. Note that the force is not constant - it depends on \(r\), which couples together the equations for \(\mathbf{v}\) and \(\mathbf{r}\)

For the parameters assume \(G=1\), \(m_1=1\) and \(m_2=1\). For the initial conditions, assume that \(r_1 = [1,0]\), \(r_2 = [-1,0]\), \(v_1 = [0,0.5]\) and \(v_2 = [0,-0.5]\). Integrate up to the time \(t_\mathrm{max} = 20\) with a time step size of \(\delta t = 0.05\).

  1. Once you have working code, change the initial velocities and comment on what you observe.