Solving simple differential equations numerically

Differential equations arise in situations where a quantity evolves, usually over time, according to a given relationship. They are extremely common in engineering and physics, and appear quite naturally. One of the classic examples of a (very simple) differential equation is the law of cooling devised by Newton. The temperature of a body cools at a rate proportional to the current temperature. Mathematically, this means that we can write the derivative of the temperature T of the body at time t > 0 using the differential equation

where k is a positive constant that determines the rate of cooling. This differential equation can be solved analytically by first "separating the variables" and then integrating and rearranging. After performing this procedure, we obtain the general solution

where T₀ is the initial temperature.

In this recipe, we will solve a simple ordinary differential equation numerically using the solve_ivp routine from SciPy.

Getting ready

We will demonstrate the technique for solving a differential equation numerically in Python using the cooling equation described previously since we can compute the true solution in this case. We take the initial temperature to be T₀= 50 and k = 0.2. Let's also find the solution for t values between 0 and 5.

A general (first order) differential equation has the form

where f is some function of t (the independent variable) and y (the dependent variable). In this formula, T is the dependent variable and f(t, T) = -kt. The routines for solving differential equations in the SciPy package require the function f and an initial value y₀and the range of t values where we need to compute the solution. To get started, we need to define our function f in Python and create the variables y₀and t range ready to be supplied to the SciPy routine:

def f(t, y):
    return -0.2*y

t_range = (0, 5)

Next, we need to define the initial condition from which the solution should be found. For technical reasons, the initial y values must be specified as a one-dimensional NumPy array:

T0 = np.array([50.])

Since, in this case, we already know the true solution, we can also define this in Python ready to compare to the numerical solution that we will compute:

def true_solution(t):
    return 50.*np.exp(-0.2*t)

How to do it...

Follow these steps to solve a differential equation numerically and plot the solution along with the error:

1. We use thesolve_ivp routine from the integrate module in SciPy to solve the differential equation numerically. We add a parameter for the maximum step size, with a value of 0.1, so the solution is computed at a reasonable number of points:

sol = integrate.solve_ivp(f, t_range, T0, max_step=0.1)

1. Next, we extract the solution values from the solobject returned from the solve_ivpmethod:

t_vals = sol.t
T_vals = sol.y[0, :]

1. Next, we plot the solution on a set of axes as follows. Since we are also going to plot the approximation error on the same figure, we create two subplots using the subplots routine:

fig, (ax1, ax2) = plt.subplots(1, 2, tight_layout=True)
ax1.plot(t_vals, T_vals)
ax1.set_xlabel("$t$")
ax1.set_ylabel("$T$")
ax1.set_title("Solution of the cooling equation")

This plots the solution on a set of axes displayed in the left-hand side of Figure 3.1.

1. To do this, we need to compute the true solution at the points that we obtained from the solve_ivp routine, and then calculate the absolute value of the difference between the true and approximated solutions:

err = np.abs(T_vals - true_solution(t_vals))

1. Finally, on the right-hand side of Figure 3.1, we plot the error in the approximation with a logarithmic scale on the y axis. We can then plot this on the right-hand side with a logarithmic scale y axis using the semilogy plot command as we saw in Chapter 2, Mathematical Plotting with Matplotlib:

ax2.semilogy(t_vals, err)
ax2.set_xlabel("$t$")
ax2.set_ylabel("Error")
ax2.set_title("Error in approximation")

The left-hand plot in Figure 3.1 shows decreasing temperature over time, while the right-hand plot shows that the error increases as we move away from the known value given by the initial condition:

How it works...

Most methods for solving differential equations are "time-stepping" methods. The pairs (t_i, y_i) are generated by taking small t steps and approximating the value of the function y. This is perhaps best illustrated by Euler's method, which is the most basic time-stepping method. Fixing a small step size h > 0, we form the approximation at the ith step using the formula

starting from the known initial value y₀. We can easily write a Python routine that performs Euler's method as follows (there are, of course, many different ways to implement Euler's method; this is a very simple example):

1. First, we set up the method by creating lists that will store the t values and y values that we will return:

def euler(func, t_range, y0, step_size):
    """Solve a differential equation using Euler's method"""
    t = [t_range[0]]
    y = [y0]
    i = 0

1. Euler's method continues until we hit the end of the t range. Here, we use a while loop to accomplish this. The body of the loop is very simple; we first increment a counter i, and then append the new t and y values to their respective lists:

    while t[i] < t_range[1]:
        i += 1
        t.append(t[i-1] + step_size)  # step t
        y.append(y[i-1] + step_size*func(t[i-1], y[i-1]))  # step y
    return t, y

The method used by the solve_ivp routine, by default, is the Runge-Kutta-Fehlberg method (RK45), which has the ability to adapt the step size to ensure that the error in the approximation stays within a given tolerance. This routine expects three positional arguments: the function f, the t range on which the solution should be found, and the initial y value (T₀in our example). Optional arguments can be provided to change the solver, the number of points to compute, and several other settings.

The function passed to the solve_ivp routine must have two arguments as in the general differential equation described in Getting ready section. The function can have additional arguments, which can be provided using the args keyword for the solve_ivp routine, but these must be positioned after the two necessary arguments. Comparing the euler routine we defined earlier to the solve_ivp routine, both with a step size of 0.1, we find that the maximum true error between the solve_ivp solution is in the order of 10^-6, whereas the euler solution only manages an error of 31. The euler routine is working, but the step size is much too large to overcome the accumulating error.

The solve_ivp routine returns a solution object that stores information about the solution that has been computed. Most important here are the t and y attributes, which contain the t values on which the solution y is computed and the solution y itself. We used these values to plot the solution we computed. The yvalues are stored in a NumPy array of shape (n, N), where nis the number of components of the equation (here, 1), and Nis the number of points computed. The yvalues held in solare stored in a two-dimensional array, which in this case has 1 row and many columns. We use the slice y[0, :]to extract this first row as a one-dimensional array that can be used to plot the solution in step 4.

We use a logarithmically scaled y axis to plot the error because what is interesting there is the order of magnitude. Plotting it on a non-scaled y axis would give a line that is very close to the x axis, which doesn't show the increase in the error as we move through the t values. The logarithmically scaled y axis shows this increase clearly.

There's more...

The solve_ivp routine is a convenient interface for a number of solvers for differential equations, the default being the Runge-Kutta-Fehlberg (RK45) method. The different solvers have different strengths, but the RK45 method is a good general-purpose solver.

See also

For more detailed instructions on how to add subplots to a figure in Matplotlib, see theAdding subplots recipe from Chapter 2, Mathematical Plotting with Matplotlib.