Data Science and Computing with Python for Pilots and Flight Test Engineers

Let \(y(x)\) be a function of \(x\) and let \(y'(x) := \frac{dy}{dx}\) denote its first derivative with respect to \(x\).

We wish to solve a first-order differential equation of the form:

$$ y'(x) = f(y(x)) $$

where \(f(y)\) is an arbitrary function of \(y\). (If your differential equation is not of this form yet, you must use algebraic manipulations to get it into this form, i.e. to isolate \(y’\) on the left-hand side. If the differential equation is of a higher order, you have to turn it into a system of first-order equations first, and then apply this technique (see next lesson).

Let \(y_0\) be the value of \(y(x)\) at \(x_0\). This is called the initial condition.

Numerical integration proceeds with small steps in \(x\) (you can think of parameter \(x\) as time). Assume that at step \(x_i\) the value of \(y\) is \(y(x_i)\). We can then compute

$$ y'(x_i) = f(y(x_i)) $$

A small time step \(\Delta x\) later, at \(x_{i+1}=x_i+\Delta x\), \(y(x)\) will have the value

$$ y(x_{i+1}) = y(x_i) + y'(x_i) \cdot \Delta x $$

All we need to get this started is to specify an initial condition, i.e. to specify at some \(x_0\) what the value of \(y(x_0)\) was, and the rest follows from the above equations.

We shall implement this now and write our own crude numerical integrator for differential equations:

import numpy as np
import matplotlib.pyplot as plt

def differential_equation_integrator(function, x, y, dx, number_of_steps):
    """ A simple numerical integrator. """
    x_list = [x]
    y_list = [y]
    
    
    for i in range(0, number_of_steps):
        y_prime = function(y) # 
        y = y + y_prime * dx # perform a step in x
        x = x + dx
        
        x_list.append(x)
        y_list.append(y)
        
    return np.array(x_list), np.array(y_list)
    

Now run the integrator on some sample problem. First we need to define the function that represents the differential equation we wish to solve:

def diff_equation(y):
    """ Differential equation for function y. """
    y_prime = (1.0/y**2) # change this depending on the differential equation you wish to solve.
    return y_prime

Then we need to run the actual numerical integration of this equation and plot the result:

# Specify initial conditions, i.e. value of y(x_0) at initial point x_0:
x = 1.0
y = 0.5

# Specify step size in integration and number of steps to be taken
dx = 0.0001 # integration step - make it smaller to make the integration more accurate (but it will also take longer to compute)
number_of_steps = 1000 # together with step size this will determine which value of $x$ will be reached during integration.

x, y = differential_equation_integrator(diff_equation, x, y, dx, number_of_steps)
    
    
# Plot the solution, i.e. y(x) as a function of x:
plt.plot(x, y, "r-")    
plt.xlabel('x')
plt.ylabel('y')

plt.show()

As always, copy and paste the above code cells into your Jupyter notebook and run. For simplicity, the above solution is for first-order ordinary differential equations with constant coefficients. It is left as an exercise to the reader to generalize the code for cases, where the coefficients in the equation depend on \(x\), i.e. when diff_equation(y, x), and not just diff_equation(y).