Introduction to Initial Value Problems, Slides of Calculus

Download Introduction to Initial Value Problems and more Slides Calculus in PDF only on Docsity!

Chapter 2

Introduction to Initial Value

Problems

The purpose of this chapter is to study the simplest numerical methods for ap- proximating the solution to a first order initial value problem (IVP). Because the methods are simple, we can easily derive them plus give graphical interpretations to gain intuition about our approximations. Once we analyze the errors made in replacing the continuous differential equation by a difference equation, we will see that the methods only converge linearly which is quite slow. This is the motivation for looking at higher accurate methods in the next chapter. We will look at several numerical examples and verify the linear convergence of the methods and we will see that in certain situations one of the methods tends to oscillate and even“blow up” while the other always provides reliable results. This will motivate us to study the numerical stability of methods. We begin by looking at the prototype IVP that we consider in this chapter and the next. The differential equation for this IVP is first order and gives information on the rate of change of our unknown; in addition, an initial value for the unknown is specified. Later, in Chapter 4, we consider higher order IVPs and we will see that higher order IVPs can be written as a system of first order IVPs so that all the methods we study in this chapter and the next can be easily extended to systems. Before we investigate methods to approximate the solution of this prototype IVP we consider conditions which guarantee that the analytic solution exists, is unique and depends continuously on the data. In the sequel we will only be interested in approximating the solution to such problems. Once we have specified our prototype IVP we introduce the idea of approximating its solution using a difference equation. In general, we have to give up the notion of finding an analytic solution which gives an expression for the solution at any time and instead find a discrete solution which is an approximation to the exact solution at a set of finite times. The basic idea is that we discretize our domain, in this case a time interval, and then derive a difference equation which approximates

18 CHAPTER 2. INTRODUCTION TO INITIAL VALUE PROBLEMS

the differential equation in some sense. The difference equation is in terms of a discrete function and only involves differences in the function values; that is, it does not contain any derivatives. Our hope is that as the difference equation is imposed at more and more points (which much be chosen in a uniform manner) then its solution will approach the exact solution to the IVP. The simplest methods for approximating the solution to our prototype IVP are the forward and backward Euler methods which we derive by approximating the derivative in the differential equation at a point by the slope of a secant line. In § 2.2.2 we demonstrate the linear convergence of the method by introducing the concepts of local truncation error and global error. The important differences in explicit and implicit methods are illustrated by comparing these two Euler methods. In § 2.3 we present some models of growth/decay which fit into our prototype IVP and give results of numerical simulations for specific problems. In addition, we demonstrate that our numerical rate of convergence matches our theoretical rate. Lastly, we demonstrate numerically that the forward Euler method gives unre- liable results in certain situations whereas the backward Euler always appears to give reliable results. This leads us to introduce the concept of numerical instability and to see why this is required for a numerical simulation to converge to the exact solution.

2.1 Prototype initial value problem

A problem commonly encountered is an initial value problem where we seek a func- tion whose value is known at some initial time and whose derivative is specified for subsequent times. The following problems are examples of first order IVPs for y(t):

I :

{ (^) y′(t) = sin πt 0 < t ≤ 4 y(0) = 0.

II :

{ (^) y′(t) + y (^2) (t) = t 2 < t ≤ 10 y(2) = 1. (2.1) Clearly, these examples are special cases of the following general IVP.

General IVP: find y(t) satisfying

dy dt

= f (t, y) t 0 < t ≤ T (2.2a)

y(t 0 ) = y 0. (2.2b)

Here f (t, y) is the given derivative of y(t) and y 0 is the known value at the initial time t 0. For example, for the IVP II in (2.1) we have f (t, y) = t − y^2 , t 0 = 2, T = 10 and y 0 = 1. Note that both linear and nonlinear differential equations are included in the general equation (2.2a). For certain choices of f (t, y) we can find an analytic solution to (2.2). In the simple case when f = f (t), i.e., f is a function of t and not both t and y, we can

20 CHAPTER 2. INTRODUCTION TO INITIAL VALUE PROBLEMS

We can always verify that we haven’t made an error in determining the solution by demon- strating that it satisfies the differential equation. Here we have

y(t) = Ce−^

t 22 ⇒ y′(t) = C − 2 t 2 e−^

t 22 = −t

( Ce−^

t 22 ) = −ty(t)

so the equation is satisfied. To determine a unique solution we impose the value of y(t) at some point; here we set y(0) = 2 to get the particular solution y(t) = 2e−t (^2) / 2 because

y(0) = 2, y(t) = Ce−^

t 22 ⇒ 2 = Ce^0 ⇒ C = 2.

Even if we are unable to determine the analytic solution to (2.2), we can still gain some qualitative understanding of the behavior of the solution. This is done by the visualization technique of plotting the tangent line to the solution at numerous points (t, y); recall that the slope of the tangent line to the solution curve is given and is just f (t, y). Mathematical software with graphical capabilities often provide commands for automatically drawing a direction field with arrows which are scaled to indicate the magnitude of the slope; typically they also offer the option of drawing some solutions or streamlines. Using direction fields to determine the behavior of the solution is illustrated in the following example.

Example 2.2. Direction fields

Draw the direction fields for the ODE

y′(t) = t^2 + y(t) 0 < t < 4

and indicate the specific solution which satisfies y(0) = 1. At each point (t, y) we draw the line with slope t^2 + y; this is illustrated in the figure below where numerous streamlines have been sketched. To thread a solution through the direction field start at a point and follow the solution, remembering that solutions don’t cross and that nearby tangent lines should be nearly the same. To see which streamline corresponds to the solution with y(0) = 1 we locate the point (0, 1) and follow the tangents; this solution is indicated by a thick line in the direction field plot below. If a different initial condition is imposed, then we get a different streamline.

• 2 - 1 1 2
  • 4
  • 2

2

4

2.1. PROTOTYPE INITIAL VALUE PROBLEM 21

Before we discuss methods for approximating the solution of the IVP (2.2) we first need to ask ourselves if our prototype IVP actually has an analytic solution, even if we are unable to find it. We are only interested in approximating the solution to IVPs which have a unique solution. However, even if we know that a unique solution exists, we may still have unreliable numerical results if the solution of the IVP does not depend continuously on the data. If this is the case, then small changes in the data can cause large changes in the solution and thus roundoff errors in our calculations can produce meaningless results. In this situation we say the IVP is ill- posed or ill-conditioned, a situation we would like to avoid. Luckily, most differential equations that arise from modeling real-world phenomena are well-posed. The conditions that guarantee well-posedness of a solution to (2.2) are well known and are presented in Theorem 2.1. Basically the theorem requires that the derivative of y(t) (given by f (t, y)) be continuous and, moreover, this derivative is not allowed to change too quickly as y changes. A basic problem in calculus is to determine how much a continuous function changes as the independent variables change; clearly we would like a function to change a small amount as an independent variable changes but this is not always the case. The concept of Lipschitz continuity 1 gives a precise measure of this “degree of continuity”. To understand this concept first think of a linear function g(x) = ax + b and consider the effect changing x has on the dependent variable g(x). We have

|g(x 1 ) − g(x 2 )| = |ax 1 + b − (ax 2 + b)| = |a| |x 1 − x 2 |.

This says that as the independent variable x varies from x 1 to x 2 the change in the dependent variable g is governed by the slope of the line, i.e., a = g′(x). For a general function g(x) Lipschitz continuity on an interval I requires that the magnitude of the slope of the line joining any two points x 1 and x 2 in I must be bounded by a real number. Formally, a function g(x) defined on a domain D ⊂ R^1 is Lipschitz continuous on D if for any x 1 6 = x 2 ∈ D there is a constant L such that |g(x 1 ) − g(x 2 )| ≤ L|x 1 − x 2 | ,

or equivalently |g(x 1 ) − g(x 2 )| |x 1 − x 2 |

≤ L.

Here L is called the Lipschitz constant. This condition says that we must find one constant L which works for all points in the domain. Clearly the Lipschitz constant is not unique; for example, if L = 5, then L = 5. 1 , 6 , 10 , 100 , etc. also satisfy the condition. If g(x) is differentiable then an easy way to determine the Lipschitz constant is to find a constant such that |g′(x)| ≤ L for all x ∈ D. The linear function g(x) = ax + b is Lipschitz continuous with L = |a| = |g′(x)|. Lipschitz continuity is a stronger condition than merely saying the function is continuous so a Lipschitz continuous function is always continuous but the converse is not true. For example, the function g(x) =

x is continuous on D = [0, 1] but is not Lipschitz continuous on D because g′(x) = 1/(

x) is not bounded on [0, 1]. (^1) Named after the German mathematician Rudolf Lipschitz (1832-1903).

2.2. DISCRETIZATION 23

the solution where we give up finding a formula for the solution at all times and instead find an approximation at a set of distinct times. Probably the most obvious approach to discretizing a differential equation is to approximate the derivatives in the equation by difference quotients to obtain a difference equation which involves only differences in function values. The solution to the difference equation will not be a continuous function but rather a discrete function which is defined over a finite set of points. When plotting the discrete solution one often draws a line through the points to get a continuous curve but remember that interpolation must be used to determine the solution at points other than grid points. In this chapter we concentrate on two of the simplest difference equations for (2.2a). Because the difference equation is defined at a finite set of points we first discretize the time domain [t 0 , T ]; alternately, if our solution depended on the spatial domain x instead of t we would discretize the given spatial interval. For now we use N + 1 evenly spaced points tn, n = 0, 1 , 2 ,... , N

t 1 = t 0 + ∆t, t 2 = t 0 + 2∆t, · · · , tN = t 0 + N ∆t = T ,

where ∆t = (T − t 0 )/N is called the step size or time step. Our task is to find a means for approximating the solution at each of these discrete values and our hope is that as we perform more calculations with N getting large, or equivalently ∆t → 0 , our approximate solution will approach the exact solution in some sense. In the left plot in Figure 2.1 we plot an exact solution (the continuous curve) to a specific IVP and a discrete approximation for ∆t =

1. 5. The approximate solution is plotted only at the points where the difference equation is enforced. From this plot we are unable to say if our discrete solution appears to approach the continuous solution; it is important to realize that when we test our code results at a single time step do not confirm that the results are correct. However, in the figure on the right we plot this discrete solution plus three additional discrete solutions where we decrease ∆t by half for each approximation. By observing the plot and using the “eyeball norm” we can convince ourselves that as ∆t → 0 our discrete solution approaches the analytic solution. One of our goals is to make this statement precise and to determine the rate at which our approximate solution converges to the exact solution.

2.2.1 The Euler methods

We will see that there are many approaches to deriving discrete methods for our IVP (2.2) but the two simplest methods use the slope of a secant line to approximate the derivative in (2.2a). These methods are called the forward and backward Euler methods named after Leonhard Euler.^2 The methods can be derived from several different viewpoints but here we use the secant line approximation

y′(t) ≈

y(tn+1) − y(tn) ∆t

for t ∈ [tn, tn+1]. (2.4)

(^2) Euler (1707-1783) was a Swiss mathematician and physicist.

24 CHAPTER 2. INTRODUCTION TO INITIAL VALUE PROBLEMS

1 2 3 4

Ê

Ê Ê (^) Ê (^) Ê Ê^ Ê

Ê

Ê

‡ ‡ ‡ ‡ ‡ ‡ (^) ‡ (^) ‡ ‡ (^) ‡ ‡ (^) ‡ (^) ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡^ ‡^ ‡^

‡ ‡^

‡^ ‡

‡^ ‡

‡ ‡

‡ ‡

‡

Ï Ï Ï Ï Ï Ï ÏÏÏÏ ÏÏÏÏÏÏ ÏÏÏÏÏÏÏÏÏÏÏÏÏÏÏÏÏÏÏÏÏÏÏÏÏÏ

ÏÏÏÏ ÏÏÏ Ï^ Ï^ Ï

Ï^ Ï Ï^ Ï Ï^ Ï

Ï^ Ï

Ï^ Ï

Ï^ Ï

Ú Ú Ú Ú Ú (^) Ú (^) Ú Ú Ú Ú Ú^ Ú

Ú Ú

Ú

Ú

Ú

1 2 3 4

Figure 2.1: The exact solution to an IVP is shown as a solid curve. In the figure on the left a discrete solution using ∆t = 0. 5 is plotted. From this plot, it is not possible to say that the discrete solution is approaching the exact solution. However, in the figure on the right the discrete solutions for ∆t = 0. 5 , 0. 25 , 0. 125 , and 0. are plotted. From this figure, the discrete approximations appear to be approaching the exact solution as ∆t decreases.

Using this difference quotient to approximate y′(tn) gives one method and using it to approximate y′(tn+1) gives the other method. When we write a difference equation we need to use different notation for the exact solution y(t) and its discrete approximation; to this end, we let Y n^ denote the approximation to y(tn). Clearly Y 0 = y 0 which is the given initial condition (2.2b). We now want to write a difference equation which will allow us to calculate Y 1 ≈ y(t 1 ). We use the differential equation evaluated at t 0 , i.e., y′(t 0 ) = f (t 0 , y 0 ), and the approximation for y′(t 0 ) from (2.4) with n = 0 and t = t 0 , i.e., y′(t 0 ) ≈ (Y 1 − Y 0 )/∆t, to get the difference equation

Y 1 − Y 0 ∆t

= f (t 0 , Y 0 ).

We have a starting point Y 0 = y 0 from our initial condition and thus we can solve for our approximation to y(t 1 ) from

Y 1 = Y 0 + ∆tf (t 0 , Y 0 ). (2.5)

Once Y 1 is obtained we can repeat the process to obtain a difference equation for Y 2. This procedure is known as the forward Euler method and is defined by the following formula.

Forward Euler: Y n+1^ = Y n^ + ∆tf (tn, Y n) , n = 0, 1 , 2 ,... , N − 1 (2.6)

The term “forward” is used in the name because we write the equation at the point tn and difference forward in time to tn+1; note that this implies that the given slope f is evaluated at the known point (tn, Y n).

26 CHAPTER 2. INTRODUCTION TO INITIAL VALUE PROBLEMS

leads to the difference equation

Y 1 − Y 0 ∆t

= f (t 1 , Y 1 ) ⇒ Y 1 = Y 0 + ∆tf (t 1 , Y 1 ).

It is important to realize that this equation is inherently different from (2.5) because we must evaluate f (t, y) at the unknown point (t 1 , Y 1 ). In general, this leads to a nonlinear equation to solve for each Y n^ which can be computationally expensive. For example, if we have f (t, y) = ety^ then the equation for Y 1 is

Y 1 = Y 0 + ∆tet^1 Y^

1

which is a nonlinear equation for the unknown Y 1. This method is called the back- ward Euler method because we are writing the equation at tn+1 and differencing backwards in time to tn.

Backward Euler: Y n+1^ = Y n^ + ∆tf (tn+1, Y n+1) , n = 0, 1 , 2 ,... , N − 1 (2.7)

The difference between the forward and backward Euler schemes is so impor- tant that we use this characteristic to broadly classify methods. The forward Euler scheme given in (2.6) is called an explicit scheme because we can write the un- known explicitly in terms of known values. The backward Euler method given in (2.7) is called an implicit scheme because the unknown is written implicitly in terms of known values and itself. The terms explicit/implicit are used in the same manner as explicit/implicit differentiation. In explicit differentiation a function to be differentiated is given explicitly in terms of the independent variable such as y(t) = t^3 + sec t; in implicit differentiation the function is given implicitly such as y^2 + sin y − t^2 = 4 and we want to compute y′(t). In the exercises you will get practice in identifying schemes as explicit or implicit. In the case when f (t, y) is linear in y we can get a general formula for Y n^ in terms of Y 0 and ∆t for both the forward and backward Euler methods by applying the formulas recursively. This means that we can compute an approximation to y(tn) without computing approximations at y(tn− 1 ),... , y(t 1 ) which we normally have to do. The following example illustrates this for the forward Euler method and in the exercises you are asked to find the analogous formula for the backward Euler method. In the next example we fix the time step and compare the relative error for a range of final times; in Exercise 2.5 we fix the final time and reduce the time step.

Example 2.3. General solution to forward euler difference equation for lin- ear equations

Consider the IVP y′(t) = −λy 0 < t ≤ T, y(0) = y 0

2.2. DISCRETIZATION 27

whose exact solution is y(t) = y 0 e−λt. Find the general solution for Y n^ in terms of Y 0 and ∆t for the forward Euler method.

For the forward Euler method we have

Y 1 = Y 0 + ∆t(−λY 0 ) = (1 − λ∆t)Y 0.

Similarly Y 2 = (1 − λ∆t)Y 1 = (1 − λ∆t)^2 Y 0

Continuing in this manner gives

Y n^ = (1 − λ∆t)nY 0.

Example 2.4. Error in the forward euler method for fixed time

In this example, we fix the time step ∆t and compare the relative error for a range of times using the formula in Example 2.3. Set λ = 5, y 0 = 2 and ∆t = 1/ 20 in Example 2.3 and compute the relative error at t = 0. 2 , 0. 4 ,... , 1. 0.

For this choice of ∆t and λ the forward Euler formula we derived in Example 2.3 reduces to Y n^ = 2(0. 95 n). We give the relative error as a percent; it is computed by taking the actual error, normalizing by the exact solution and converting to a percent. The table below gives the approximations.

t Percentage Error

0.2 2.55% 0.4 5.04% 0.6 7.47% 0.8 9.83% 1.0 12.1%

As can be seen from the table the relative error increases as the time increases which we expect because the errors are being accumulated; this will be discussed in detail in the next section. It is important to compute a relative error because if the exact solution is near zero then the absolute error may be small while the relative error is large.

Example 2.5. Error in the forward euler method as time step decreases

In this example we fix the time at which the error is calculated and vary ∆t using the formula derived in Example 2.3 for the forward Euler method with λ = 5 and y 0 = 2. We compute the absolute and relative errors at t = 1 for ∆t = 1/ 10 ,... , 1 / 320 and discuss the results. Then we determine how many times do have to halve the time step if we want the relative error to be less than 1%.

We use the expressions for Y n^ from Example 2.3 to determine the approximations to y(1) = 2e−^5 ≈ 1 .3476 10−^2.

2.2. DISCRETIZATION 29

t t 0 t 1 t 2 t 3

y(t)

s Y 0 = y(t 0 )

Y 1 s

Y (^) s 2

s

Y 3

c y(t 1 )

c y (^) (t 2 )

y(t 3 ) c

Figure 2.3: The exact solution and the discrete solution agree at t 0. At t 1 the error |Y 1 − y(t 1 )| is due to approximating the derivative in the ODE by a difference quotient. At t 2 the error |Y 2 − y(t 2 )| is due to approximating the derivative in the ODE and the fact that the starting value, Y 1 , does not lie on the solution curve as Y 0 did.

and is the difference in y(t 2 ) and Y˜ 2 = y(t 1 ) + ∆tf

t 1 , y(t 1 )

, i.e., the remainder when the exact solution is substituted into the difference equation at t 1. The global error at t 2 is illustrated in the plot on the left of the figure and is just the difference in y(t 2 ) and Y 2.

We now want to determine the local truncation error for the forward Euler method so we substitute the exact solution to (2.2a) into the difference equation (2.6) and calculated the remainder. If τn+1 denotes the local truncation error at the (n + 1)st time step then

τn+1 = y(tn+1) −

[

y(tn) + ∆tf

tn, y(tn)

)]

Our strategy is to first quantify the local truncation error in terms of ∆t and then use this result to determine the global error. In order to to combine terms in (2.8) we need all terms to be evaluated at the same point

tn, y(tn)

. The only term not at this point is the exact solution y(tn+1) so we use a Taylor series with remainder (see (??) in Appendix) because the expansion is in terms of y(tn) and its derivatives at tn. We have

y(tn+1) = y(tn + ∆t) = y(tn) + ∆ty′(tn) +

(∆t)^2 2!

y′′(ξi) ξi ∈ (tn, tn+1).

Substituting this into the expression (2.8) for the truncation error yields

τn+1 =

[

y(tn) + ∆tf

tn, y(tn)

(∆t)^2 2!

y′′(ξi)

]

[

y(tn) + ∆tf

tn, y(tn)

)]

30 CHAPTER 2. INTRODUCTION TO INITIAL VALUE PROBLEMS

t t 0 t 1 t 2

y(t) y(t)

Y 0 r

Y 1 r

Y (^2) r

r

Y 2 = Y 1 + ∆tf (t 1 , Y 1 )

 global error

b y(t 1 )

b y(t 2 )

t t 0 t 1 t 2

Y 0 r

b y(t 1 )

b y(t 2 )



Y˜ 2 r  local error

Y^ ˜ 2 = y(t 1 ) + ∆tf

t 1 , y(t 1 )

Figure 2.4: A comparison of the global error (left figure) and the local truncation error (right figure) at t 2 for the forward Euler method. The global error is the total error made whereas the local truncation error is the error due to the discretization of the differential equation.

(∆t)^2 2!

y′′(ξi) ,

where we have used the differential equation at tn, i.e., y′(tn) = f

tn, y(tn)

. If y′′(t) is bounded on [0, T ], say |y′′(t)| ≤ M and T = t 0 + N ∆t, then we have

τ = max 1 ≤n≤N |τn| ≤

M

(∆t)^2. (2.9)

We say that the local truncation error for Euler’s method is order (∆t)^2 which we write as O

∆t^2

. This says that the local error is proportional to the square of the step size; i.e., it is a constant times the square of the step size. This means that if we compute the local error for ∆t then the local error using ∆t/ 2 will be reduced by approximately (1/2)^2 = 1/ 4. Remember, however, that this is not the global error but rather the error made because we have used a finite difference quotient to approximate y′(t). We now turn to estimating the global error in the forward Euler method. We should expect to only be able to find an upper bound for the error because if we can find a formula for the exact error, then we can calculate this and add it to the approximation to get the exact solution. The proof for the global error for the forward Euler method is a bit technical but it is the only global error estimate that we will derive because the methods we consider will follow the same relationship between the local and global error as the Euler method.

32 CHAPTER 2. INTRODUCTION TO INITIAL VALUE PROBLEMS

Step III. Recognize that the sum in (2.11) is a geometric series whose sum is known to get En ≤

τ ∆tL

[(1 + ∆tL)n^ − 1]. (2.12)

Step IV. Use the Taylor series expansion of e∆tL^ near zero to bound (1 + ∆tL)n by en∆tL^ which in turn is less than eT L.

Step V. Use the bound (2.9) for τ to get the final result

En ≤

M ∆t^2 2∆tL (eT L^ − 1) = C∆t where C =

M

2 L

(eT L^ − 1). (2.13)

We now give the details for each step. For the first step we use the fact that the local truncation error is the remainder when we substitute the exact solution into the difference equation; i.e.,

τn = y(tn) − y(tn− 1 ) − ∆tf

tn− 1 , y(tn− 1 )

To get the desired expression for En we solve for y(tn) in the above expression, substitute into the definition for En and use the triangle inequality; then we use the forward Euler scheme for Y n^ and (2.3). We have

En =

∣[τn + y(tn− 1 ) + ∆tf (tn− 1 , y(tn− 1 ))]^ − Y n

∣τn + y(tn− 1 ) + ∆tf (tn− 1 , y(tn− 1 ))^ − [Y n−^1 + ∆tf (tn− 1 , Y n−^1 )]∣∣

≤

∣τn| + |y(tn− 1 ) − Y n−^1 | + ∆t|f (tn− 1 , y(tn− 1 ))^ − f (tn− 1 , Y n−^1 )

≤ |τn| + En− 1 + ∆tL|y(tn− 1 ) − Y n−^1 | = |τn| + (1 + ∆tL)En− 1.

In the final step we have used the Lipschitz condition (2.3) which is a hypothesis of the theorem. Since |τn| ≤ τ , we have the desired result. For the second step we apply (2.10) recursively

En ≤ KEn− 1 + τ ≤ K[KEn− 2 + τ ] + τ = K^2 En− 2 + (K + 1)τ ≤ K^3 En− 3 + (K^2 + K + 1)τ ≤ · · ·

≤ KnE 0 + τ

n∑− 1

i=

Ki^.

Because we assume for analysis that there are no roundoff errors, E 0 = |y 0 −Y 0 | = 0 we are left with τ

∑n− 1 i=0 K

i. For the third step we simplify the sum by noting that

it is a geometric series of the form

∑n− 1 k=0 ar k (^) with a = τ and r = K. From calculus

we know that the sum is given by a(1 − rn)/(1 − r) so that if we use the fact that K = 1 + ∆tL we arrive at the result (2.12)

En ≤ τ

1 − Kn 1 − K

= τ

Kn^ − 1 K − 1

τ ∆tL

[

(1 + ∆tL)n^ − 1

]

2.3. NUMERICAL COMPUTATIONS 33

To justify the fourth step we know that for real z the Taylor series expansion ez^ = 1 + z + z^2 /2! + · · · near zero implies that 1 + z ≤ ez^ so that (1 + z)n^ ≤ enz^. If we set z = ∆tL we have (1 + ∆tL)n^ ≤ en∆tL^ so that

En ≤

τ ∆tL

(en∆tL^ − 1).

For the final step we know from the hypothesis of the theorem that |y′′(t)| ≤ M so τ ≤ M ∆t^2 / 2. Also n in En is the number of steps taken from t 0 so n∆t = tn ≤ T where T is the final time and so en∆tL^ ≤ eT L. Combining these results gives the desired result (2.13). 

In general, the calculation of the local truncation error is straightforward (but sometimes tedious) whereas the proof for the global error estimate is much more involved. However, for the methods we consider, if the local truncation error is O(∆tr^ ) then we expect the global error to be one power of ∆t less, i.e., O(∆tr−^1 ). One can show that the backward Euler method also has a local truncation error of O(∆t^2 ) so we expect a global error of O(∆t). You are asked to demonstrate the local truncation rigorously in the exercises. The calculation is slightly more complicated than for the forward Euler method because the term f

tn+1, y(tn+1)

must be expanded in a Taylor series with two independent variables whereas for the forward Euler method we did not have to expand f

tn, y(tn)

because it was already at the desired point

tn, y(tn)

In the following section we look at some specific examples of the IVP (2.2) and use both forward and backward Euler methods; we will demonstrate that our numerical rate of convergence agrees with the theoretical rate. However, we should keep in mind that O(∆t) is a very slow rate of convergence and ultimately we need to derive methods which converge more quickly to the solution.

2.3 Numerical computations

In this section we provide some numerical simulations for IVPs of the form (2.2) using both the forward and backward Euler methods. We begin by looking at some examples of IVPs which include exponential, logistic and logarithmic growth/decay. Before providing the simulation results we discuss the computer implementation of both the forward and backward Euler methods. For the simulations presented we choose problems with known analytic solutions so that we can compute numerical rates of convergence and compare with the theoretical result given in Theorem 2.2. To do this we compute approximate solutions for a sequence of step sizes where ∆t → 0 and then compute a numerical rate of convergence using § 1.6. As expected, the numerical rate of convergence for both methods is linear; however, we will see that the forward Euler method does not provide reliable results for all choices of ∆t for some problems. Reasons for this failure will be discussed in § ??.

2.3. NUMERICAL COMPUTATIONS 35

if it is positive; otherwise it is the rate of decay. In the exercises you are asked to explore an exponential growth model for bread mold.

Logistic growth and decay

The previous model of population growth assumes there is an endless supply of resources and no predators. Logistic growth of a population attempts to incorporate resource availability by making the assumption that the rate of population growth (i.e., the proportionality constant) is dependent on the population density. Figure ?? compares exponential growth and logistic growth; clearly exponential growth allows the population to grow in an unbounded manner whereas logistic growth requires the population to stay below a fixed amount K which is called the carrying capacity of the population. When the population is considerably below this threshold the two models produce similar results. The logistic model we consider restricts the growth rate in the following way

r = r 0

p K

where K is the maximum allowable population and r 0 is a given growth rate for small values of the population. As the population p increases to near the threshold value K then (^) Kp becomes closer to one (but less than one) and so the term (1 − (^) Kp ) gets closer to zero and the growth rate decreases because of fewer resources; the limiting value is when p = K and the growth rate is zero. However when p is small compared with K, the term (1 − (^) Kp ) is near one and it behaves like exponential growth with a rate of r 0. Assuming the population at any time is proportional to the current population using the proportionality constant (2.15), our differential equation becomes

p′(t) = r 0

p(t) K

p(t) = r 0 p(t) − r 0 K

p^2 (t)

along with p(t 0 ) = p 0. This equation is nonlinear in the unknown p(t) due to the p^2 (t) term and is more difficult to solve than the exponential growth equation. However, it can be shown that the solution is

p(t) = Kp 0 (K − p 0 )e−r^0 t^ + p 0

which can be verified by substitution into the differential equation. We expect that as we take the limt→∞ p(t) we should get the threshold value K. Clearly this is true because limt→∞ e−r^0 t^ = 0.

Logarithmic growth

Logarithmic growth has a period of rapid increase followed by a period where the growth rate slows. This is in contrast to exponential growth which has a slow initial growth rate and then grows more rapidly. Unlike logistic growth, there is no

36 CHAPTER 2. INTRODUCTION TO INITIAL VALUE PROBLEMS

upper bound for the growth. The IVP whose solution obeys the logarithmic growth y(t) = a + b ln t is

y′(t) =

b t

1 < t ≤ T, y(1) = a. (2.17)

An application of logarithmic for modeling strength of earthquakes is explored in the exercises.

2.3.2 Computer implementation

The computer implementation of the forward Euler method requires a single time loop. The information which changes for each IVP is the interval [t 0 , T ], the initial condition, the given slope f (t, y) and the exact solution for the error calculation. We incorporate separate functions for f (t, y) and the exact solution and input the other variables as well as the step size ∆t. The solution should not be stored for all time; if it is needed for plotting, etc., then the time and solution should be written to a file to be used later. For a single equation it does not take much storage to keep the solution at each time step but when we encounter systems and problems in higher dimensions the storage to save the entire solution can become prohibitively large so it is not good to get in the practice of storing the solution at each time. The following pseudocode gives an outline of one approach to implement the forward Euler method using a uniform time step. Note that in this algorithm we have specified N , the number of time steps and computed ∆t by (T − t 0 )/N. The only problem with this is when the time interval is, for example, π, then N ∆t may not be exactly equal to π due to roundoff so for the error you must compare the approximate solution at computed time t with the exact solution there rather than at T. Alternately, if we specify the initial and final times along with ∆t, then N can be computed from N = (T − t 0 )/∆t. Even then N ∆t may not be exactly T so care should be taken.

Algorithm 1. forward euler method

Define: external function for the given slope f (t, y) and exact solution for error calculation

Input: the initial time, t 0 , the final time T , the initial condition y = y 0 , and the number of uniform time steps N

∆t = (T − t 0 )/N

t = t 0

for i=1:N

m = f (t, y) y = y + ∆t m t = t + dt output t, y