Example Solution of Differential Equations

An ordinary differential equation (ODE), where the unknown is a function u(t), can be written in the generic form u'(t) = f (u(t),t). (7.3)

In addition, an initial condition, u(0) = u0, must be associated with this ODE to make the solution of (7.3) unique. The function f reflects an expression with u and/or t. Some important examples of ODEs and their corresponding forms of f are given below.

1. Exponential growth of money or populations:

where a is a given constant expressing the growth rate of u.

2. Logistic growth of a population under limited resources:

where a is the initial growth rate and R is the maximum possible value of u.

3. Radioactive decay of a substance:

where a is the rate of decay of u.

4. Body falling in a fluid:

where b models the fluid resistance, g is the acceleration of gravity, and u is the body's velocity (see Exercise 7.25 on page 405).

5. Newton's law of cooling:

where u is the temperature of a body, h is a heat transfer coefficient between the body and its surroundings, and s is the temperature of the surroundings.

Appendix B gives an introduction to ODEs and their numerical solution, and you should be familiar with that or similar material before reading on.

The purpose of the present section is to design and implement a class for the general ODE u' = f (u, t). You need to have digested the material about classes in Chapters 7.1.2, 7.3.1, and 7.3.2.

7.4.1 A Function for Solving ODEs

A very simple solution method for a general ODE on the form (7.3) is the Forward Euler method:

Here, Uk denotes the numerical approximation to the exact solution u at time tk, ^t is a time step, and if all time steps are equal, we have that tk = k^t, k = 0,..., n.

First we will implement the method (7.9) in a simple program, tailored to the specific ODE u' = u (i.e., /(u, t) = u in (7.3)). Our goal is to compute u(t) for t € [0,T]. How to perform this computation is explained in detail in Appendix B.2. Here, we follow the same approach, but using lists instead of arrays to store the computed u0,... ,un and to,..., tn values. The reason is that lists are more dynamical if we later introduce more sophisticated solution methods where ^t may change during the solution so that the value of n (i.e., length of arrays) is not known on beforehand. Let us start with sketching a "flat program":

# Integrate u'=u, u(0)=u0, in steps of dt until t=T u0 = 1 T = 3 dt = 0.1

u.append(u0) t.append(0) n = int(round(T/dt)) for k in range(n):

u.append(unew) tnew = t[-1] + dt t.append(tnew) from scitools.std import plot plot(t, u)

The new uk+1 and tk+1 values, stored in the unew and tnew variables, are appended to the u and t lists in each pass in the loop over k.

Unfortunately, the code above can only be applied to a specific ODE using a specific numerical method. An obvious improvement is to make a reusable function for solving a general ODE. An ODE is specified by its right-hand side function /(u, t). We also need the parameters u0, ^t, and T to perform the time stepping. The function should return the computed u0,..., un and t0,..., tn values as two separate arrays to the calling code. These arrays can then be used for plotting or data analysis. An appropriate function may look like def ForwardEuler(f, dt, u0, T):

Integrate u'=f(u,t), u(0)=u0, in steps of dt until t=T

u = []; t = [] # u[k] is the solution at time t[k] u.append(u0)

t.append(O) n = int(round(T/dt)) for k in range(n):

u.append(unew) tnes = t[-1] + dt t.append(tnes) return numpy.array(u), numpy.array(t)

Here, f(u, t) is a Python implementation of f (u, t) that the user must supply. For example, we may solve u' = u for t € (0, 3), with u(0) = 1, and At = 0.1 by the following code utilizing the shown ForwardEuler function:

# compare numerical solution and exact solution in a plot: from scitools.std import plot, exp u_exact = exp(t)

plot(t, u, 'r-', t, u_exact, 'b-', xlabel='t', ylabel='u', legend=('numerical', 'exact'), title="Solution of the ODE u'=u, u(0)=1")

Observe how easy it is to plot u versus t and also add the exact solution u = el for comparison.

7.4.2 A Class for Solving ODEs

Instead of having the numerical method for solving a general ODE implemented as a function, we now want a class for this purpose. Say the name of the class is ForwardEuler. To solve an ODE specified by a Python function f(u, t), from time t0 to some time T, with steps of size dt, and initial condition u0 at time t0, it seems convenient to write the following lines of code:

method = ForwardEuler(f, dt) method.set_initial_condition(u0, t0) u, t = method.solve(T)

The constructor of the class stores f and the time step dt. Then there are two basic steps: setting the initial condition, and calling solve to advance the solution to a specified time level. Observe that we do not force the initial condition to appear at t = 0, but at some arbitrary time. This makes the code more general, and in particular, we can call the solve again to advance the solution further in time, say to 2T:

method.set_initial_condition(u[-1] , t[-1]) u2, t2 = method.solve(2*T) plot(t, u, 'r-', t2, u2, 'r-')

The initial condition of this second simulation is the final u and t values of the first simulation. To plot the complete solution, we just plot the individual simulations.

The task now is to write a class ForwardEuler that allow this type of user code. Much of the code from the ForwardEuler function above can be reused, but it is reorganized into smaller pieces in a class. Such reorganization is known as refactoring (see also Chapter 3.6.2). An attempt to write the class appears below.

class ForwardEuler: tttttt

by the ForwardEuler method.

Class attributes: t: array of time values u: array of solution values (at time points t) k: step number of the most recently computed solution f: callable object implementing f(u, t) dt: time step (assumed constant) ttt def __init__(self, f, dt):

self.f, self.dt = f, dt def set_initial_condition(self, u0, t0=0):

self.u = [] # u[k] is solution at time t[k] self.t = [] # time levels in the solution process self.u.append(float(u0)) self.t.append(float(t0)) self.k =0 # time level counter def solve(self, T):

"""Advance solution in time until t <= T."""

unew = self.advance()

self.u.append(unew)

return numpy.array(self.u), numpy.array(self.t)

def advance(self):

"""Advance the solution one time step."""

# load attributes into local variables to

# obtain a formula that is as close as possible

We see that we initialize two lists for the uk and tk values at the time we set the initial condition. The solve method implements the time loop as in the ForwardEuler function. However, inside the time loop, a new uk+1 value (unew) is computed by calling another class method, self. advance, which here implements the numerical method (7.9).

Changing the numerical method is just a matter of changing the advance function only. We could, of course, put the numerical updating formula explicitly inside the solve method, but changing the numerical method would then imply editing internals of a method rather than replacing a complete method. The latter task is considered less error-prone and therefore a better programming strategy.

7.4.3 Verifying the Implementation

We need a problem where the exact solution is known to check the correctness of class ForwardEuler. Preferably, we should have a problem where the numerical solution is exact such that we avoid dealing with approximation errors in the Forward Euler method. It turns out that if the solution u(t) is linear in t, then the Forward Euler method will reproduce this solution exactly. Therefore, we choose u(t) = at + uo, with (e.g.) a = 0.2 and u0 = 3. The corresponding f is the derivative of u, i.e., f (u, t) = a. This is obviously a very simple right-hand side without any u or t. We can make f more complicated by adding something that is zero, e.g., some expression with u - at - u0, say (u - at - u0)4, so that f (u, t) = a + (u - at - u0)4.

We implement our special f and the exact solution in two functions _f1 and _u_solution_f1:

Testing the code is now a matter of performing the steps u0 = 3 dt = 0.4 T = 3

method = ForwardEuler(_f1, dt) method.set_initial_condition(u0, 0) u, t = method.solve(T) u_exact = _f1_solution(t) print 'Numerical:\n', u print 'Exact:', '\n', u_exact

The output becomes

Numerical:

[ 3. 3.08 3.16 3.24 3.32 3.4 3.48 3.56 3.64] Exact:

showing that the code works as it should in this example.

7.4.4 Example: Logistic Growth

A more exciting application is to solve the logistic equation (B.23), with the f (u, t) function specified in (7.5).

First we may create a class for holding information about this problem: a, R, the initial condition, and the right-hand side. We may also add a method for printing the equation and initial condition. This problem class can then be expressed as class Logistic:

Problem class for a logistic ODE

Return f(u,t) for the logistic ODE

Return ODE and initial condition

return "u'(t) = %g*u*(1 - u/°/.g)\nu(0)=°/.g" % \ (self.alpha, self.R, self.u0)

Running a case with a = 0.2, R = 1, u(0) = 0.1, ^t = 0.1, and simulating up to time T = 40, can be performed in the following function:

def logistic():

method = ForwardEuler(problem, dt) method.set_initial_condition(problem.u0, 0) u, t = method.solve(T)

from scitools.std import plot, hardcopy, xlabel, ylabel, title plot(t, u)

title('Logistic growth: alpha=0.2, dt=°/g, /d steps' \

The resulting plot is shown in Figure 7.3. Note one aspect of this function: the "star import", as in from scitools.std import *, is not allowed inside a function (or class method for that sake), so we need to explicitly list all the functions we need to import. (We could, as in the previous example, just import plot and rely on keyword arguments to set the labels, title, and output file.)

The ForwardEuler class is further developed in Chapter 9.4, where it is shown how we can easily modify the class to implement other numerical methods. In that chapter we extend the implementation to systems of ODEs as well.

Logistic growth: alpha=0.2, dt=0.1, 400 steps

Logistic growth: alpha=0.2, dt=0.1, 400 steps

Fig. 7.3 Plot of the solution of the ODE problem u' = 0.2u(1 - u), «(0) = 0.1.
0 0

Post a comment