An Intuitive Elucidation of the Forward Euler Numerical Method

Jason D. Morken1,2,*

1. School of Mathematical and Statistical Sciences, Arizona State University, Tempe 2. Department of Chemistry and Biochemistry, Arizona State University, Tempe

October, 2015

Recall that we can compute the slope m between any two points (t0, y0) and (t1, y1) using the formula

m = y1 −y0 t1 − t0

. (1)

Now, recall that if we choose any two real numbers, say 0.1 and 0.0001 for example, there are infinitely many real numbers between those two numbers. The differential equations we are dealing with are defined continuously over subsets of R (the set of real numbers). A computer cannot store an infinite number of points and therefore we must discretize the interval on which we wish to compute our numerics. By discretize, we mean represent the continuous interval as an array (vector) containing a discrete number of elements (as opposed to an infinite number of elements).

Consider a one-dimensional ODE initial value problem (IVP) dy/dt = f(t, y) with y(t0) = y0, and where t, y ∈ R (i.e., where variables t and y are defined over the set of real numbers). If we wish to solve the IVP numerically over some subset of the real numbers, we must discretize the time interval on which we wish to solve the differential equation numerically. Then the numerical method will compute approximate solution values (to the IVP) at the corresponding time values in our discretization.

Suppose we wish to take N timesteps. This means our discretization will contain N + 1 time points. Then ∆t = (tN − t0)/N. NOTE: it is common practice to denote ∆t by h, and this is the notation used in the protocol. Notice the inverse proportionality between N and ∆t (or N and h). Hence, if ∆t is decreased by a factor of 10, then N consequently increases by a factor of 10. With these definitions, we have that tn = t0 + n∆t for each n = 0, 1, …, N. Therefore, our time vector t will take the form t = [t0, t1, …, tN ]. Then the numerical method computes solution approximations at each time value in the vector t and thus produces a vector y of the form y = [y0, y1, …, yN ] which inherits the dimensions of the vector t.

For the purposes of this class, adjacent elements in the vector t will be of uniform distance apart ∆t. Then t1 = t0 + ∆t. Therefore, t1 − t0 = ∆t and thus equation (2) becomes

m = y1 −y0

∆t . (2)

Now, notice that y1 approximates the actual solution value y(t1) = y(t0 + ∆t). However, since we are given the initial condition y(t0) = y0, then y0 in our solution vector is known and is exact. Therefore, we have that

m ≈ y(t0 + ∆t) −y(t0)

∆t . (3)

1

This is starting to look similar to the limit definition of the derivative

dy(t0)

dt = lim

∆t→0

y(t0 + ∆t) −y(t0) ∆t

. (4)

Therefore, if we choose ∆t sufficiently small, we should have that

dy(t0)

dt ≈

y(t0 + ∆t) −y(t0) ∆t

. (5)

Writing this in terms of elements in our solution vector, we get

dy(t0)

dt = f(t0, y0) ≈

y1 −y0 ∆t

. (6)

Now, we will always know the expression for the function f which simply represents the right-hand side of the ODE. Therefore, f(t0, y0) is a known value. Recall that y0 is also a known value. Then the only unknown here is y1 which can be computed by rearranging equation (6) to get

y1 = y0 + ∆tf(t0, y0). (7)

If ∆t is sufficiently small, equation (7) should give us a fairly accurate approximation of the actual solution value y(t1). After we compute y1, we can compute y2 using

y2 = y1 + ∆tf(t1, y1).

Continuing on in this fashion we arrive at the forward Euler method

yn+1 = yn + ∆tf(tn, yn), n = 0, 1, 2, …, N. (8)

Therefore, we have a vector containing a discrete number of exact time values and a vector of corresponding solution approximations at each time value where yn represents the numerical approximation of the exact (analytical) solution value y(tn) for each n = 1, 2, …, N, where n = 0 is excluded since y0 = y(t0) is known.

Final note: although, analytically, we represent elements in the vectors t and y using indices n = 0, 1, …, N, MATLAB cannot index vectors starting at n = 0. Therefore, the vector t = [t0, t1, …, tN ] must be represented in MATLAB as t = [t(1),t(2),…,t(N+1)]. Further, the output of euler.m will be two column vectors t and y as opposed to row vectors which we often use in our analytical discussion.

This document was created in LATEX.

## Needs help with similar assignment?

We are available 24×7 to deliver the best services and assignment ready within 3-4 hours? Order a custom-written, plagiarism-free paper

Get Answer Over WhatsApp

Order Paper Now