Ordinary differential equations

We can use calculus to study how quantities change in time. In the examples below we consider a second order linear ordinary differential equation (ODE) with constant coefficients.

Background

You might have previously encountered differentiation. Suppose that \(y\) is some function of \(x\).

Consider the first order ODE

\[ \frac{dy}{dx}=1. \]

Upon integration

\[ y(x)=x+C, \] where \(C\) is an integration constant.

What if \[ \frac{dy}{dx}=x \]

Can you integrate this ODE?

The above examples are first order ODEs because the highest order derivative is one. In a second order ODE the highest derivative is two. For example \[ \frac{d^2y}{dx^2}=-y \tag{1}\]

The solution to this ODE is

\[ y(x)=A\sin x + B \cos x\] where \(A\) and \(B\) are integration constants. Note that you can confirm this claim by twice differentiating the solution and substituting in Equation 1.

Newton’s Second Law

Second order ODEs arise in many different situations in applied mathematics. One of the most notable is in models of Newton’s Second Law \[ \mathbf{F}=m\mathbf{a}, \] where \(\mathbf{a}\) represents the acceleration of a particle, \(\mathbf{F}\) the resultant force acting on it and \(m\) the constant particle mass.

The acceleration can be written as the second order derivative of the particle position with respect to time, i.e. \[ \mathbf{a}=\frac{d^2\mathbf{r}}{dt^2}. \]

In physical situations the force term will be composed of terms that are relevant in a given situation.

Suppose a particle is subjected to a damping force and a spring-like force centred at the origin. We might write \[ \mathbf{F} = -\eta \frac{d \mathbf{r}}{dt}-k\mathbf{r} \] where \(\eta\) is a constant damping coefficient and \(k\) is a spring constant.

Hence we would obtain the second order ODE \[ m\frac{d^2\mathbf{r}}{dt^2}= -\eta \frac{d \mathbf{r}}{dt}-k\mathbf{r}. \]

Upon solving this equation we can describe the particle’s position.

Let’s consider a scalar second order ODE with constant coefficients. Let \(t\) represent time and \(y=y(t)\). Suppose that

\[ a\frac{d^2 y}{dt^2}+b\frac{dy}{dt}+cy=0, \tag{2}\] where \(a\), \(b\) and \(c\) are constants.

Consider the initial conditions

\[ y(t=0)=y_0 \quad \quad \frac{dy}{dt}\bigg|_{{t=0}}=d. \tag{3}\]

The app in Figure 1 allows you to explore the solution of the model (Equation 2 and Equation 3).

Can you identify the following three behaviours by adjusting model parameters:

periodic oscillations
damped oscillations
exponentially decaying but non-oscillatory solutions

Can you

describe model behaviour as \(a\rightarrow 0\).
explore model behaviour as \(c\rightarrow 0\).
describe how the oscillation period depends on model parameters.

Note that it is possible to solve Equation 2 and Equation 3 by seeking a solution of the form \[ y=e^{\lambda t}. \] This is left as an exercise.

#| standalone: true
#| components: [viewer]
#| viewerHeight: 800

from shiny import App, Inputs, Outputs, Session, render, ui
from shiny import reactive

import numpy as np
from pathlib import Path
import matplotlib.pyplot as plt
from scipy.integrate import odeint

app_ui = ui.page_fluid(
    ui.layout_sidebar(
        ui.panel_sidebar(
    ui.input_slider(id="a",label="a",min=0.01,max=3,value=0.1,step=0.001),
    ui.input_slider(id="b",label="b",min=0.0,max=15.0,value=10.0,step=0.1),             
    ui.input_slider(id="c",label="c",min=0.0,max=30.0,value=5.0,step=1.0),   
    ui.input_slider(id="y0",label="y(t=0)",min=0.0,max=20.0,value=5.0,step=1.0),
    ui.input_slider(id="y0p",label="dy/dt(t=0)",min=0.0,max=20.0,value=1,step=1.0),
    ui.input_slider(id="T",label="Simulation time",min=0.0,max=60.0,value=20.0,step=0.5),
              
          
            ),

        ui.panel_main(ui.output_plot("plot"),),
    ),
)

def server(input, output, session):
    
    @render.plot
    def plot():
        fig, ax = plt.subplots()
        #ax.set_ylim([-2, 2])
        # Filter fata
        
        
        a=float(input.a())
        b=float(input.b())
        c=float(input.c())
        y_0=float(input.y0())
        y_0_p=float(input.y0p())
        T=float(input.T())
        

        
        # Define rhs of LV ODEs
        def rhs_pop_model(y,t,a,b,c):
          rhs=np.zeros_like(y,dtype=float)

          z=y[1]

          #ay'' + by'+cy=0
          # y'=z
          # z'=y''=-(by'+cy)/a


          dy_dt=y[1]
          dz_dt=-(b*z+c*y[0])/a

          rhs[0]=dy_dt
          rhs[1]=dz_dt

          return rhs

        # Define discretised t domain
        t = np.linspace(0, T, 1000)

        # define initial conditions
        init_cond=[y_0,y_0_p]
        
        # Compute numerical solution of ODEs
        sol1 = odeint(rhs_pop_model, init_cond,t,args=(a,b,c))

        # Plot results
        y=sol1[:,0]
        yp=sol1[:,1]
        
        
        ax.plot(t, y)
        ax.set_xlabel('$t$')
        ax.set_ylabel('$y(t)$')

        plt.grid()
        #plt.show()
    
app = App(app_ui, server)

Figure 1

Note

At Dundee, core concepts from calculus (e.g. differential equations) are studied in the modules Maths 1A and Maths 1B and developed further in the modules Maths 2A and Maths 2B.

At Level 2 in the modules Computer algebra and dynamical systems and Introduction to Programming you would be introduced to techniques that are used to compute numerical solutions to differential equations.

At Level 3 in the module Differential Equations you would extend your knowledge of differential equations to include concepts such as Fourier Series and Partial Differential Equations. In the modules Mathematical Biology I and Mathematical Biology II you would also learn how to formulate and study mathematical models of biological systems.

You can find out more about these modules here.