A Comparison of Different Numerical Schemes to Solve Nonlinear First Order ODE ()
1. Introduction
Consider the nonlinear first order ordinary differential equation
(1)
Ordinary differential equations occur in many scientific disciplines, including physics, chemistry, biology, and economics. In mathematics, an ordinary differential equation (ODE) is a differential equation containing one or more functions of one independent variable and the derivatives of those functions. Ordinary differential equations (ODEs) arise in many contexts of mathematics and social and natural sciences. Mathematical descriptions of change use differentials and derivatives. Various differentials, derivatives, and functions become related via equations, such that a differential equation is a result that describes dynamically changing phenomena, evolution, and variation. Some ODEs can be solved explicitly in terms of known functions and integrals. When that is not possible, the equation for computing the Taylor series of the solutions may be useful. For applied problems, numerical methods for ordinary differential equations can supply an approximation of the solution. In addition, some methods in numerical partial differential equations convert the partial differential equation into an ordinary differential equation, which must then be solved. Euler method (also called forward Euler method) is a first-order numerical procedure for solving ordinary differential equations (ODEs) with a given initial value. It is the most basic explicit method for numerical integration of ordinary differential equations and is the simplest Runge-Kutta method [1] [2] [3]. The backward Euler method (or implicit Euler method) is one of the most basic numerical methods for the solution of ordinary differential equations, the backward Euler method has order one. This means that the local truncation error (defined as the error made in one step) is
. The error at a specific time t is
. It is similar to the (standard) Euler method, but differs in that it is an implicit method. This differs from the (forward) Euler method in that the latter uses
in place of
. The backward Euler method is an implicit method: the new approximation
appears on both sides of the equation, and thus the method needs to solve an algebraic equation for the unknown
. For non-stiff problems, this can be done with fixed-point iteration [4] [5] [6] [7]. This paper is organized as follows: In Chapter 2 we presented three simplest numerical stencils. In Chapter 3 we introduce a numerical formulation of ODE and apply examples.
2. Three Simplest Numerical Stencils
2.1. Forward Euler Stencil
The purpose of this paper is to derive the three simplest numerical stencils to solve the first order equation
(2)
In this expression, f is assumed to be a known function of the independent variable x and the function that we are trying to solve for
. The simplest numerical stencils to solve this equation will give us an approximation to y at some point
given some knowledge of y at
. All of these stencils are based on the Taylor series approximation for
about
to linear order:
(3)
Let us define h as the difference between x and X, and then get rid of x in the above expression:
(4)
Now, we can remove the first derivative of y by making use of the differential equation:
(5)
The forward Euler algorithm involves discarding the terms of order h2 and higher in this expression and hence obtaining an approximation to
in terms of
. We can accomplish this by converting the above series into a polynomial using the convert/polynom command. Also, we can take
as the ith point on an evenly spaced lattice
, where h is the lattice spacing and i is an integer; it then follows that “
”. Furthermore, we label the numeric approximation of
at
a “
and approx;
”. Implementing these steps and notational changes:
(6)
The last expression is the forward Euler stencil for solving ode. It is called an explicit algorithm because the quantity we wish to calculate
is given explicitly in terms of
.
2.2. Backward Euler Stencil
The backward Euler stencil is obtained in a similar fashion, except we identify h,
, and
differently:
(7)
This stencil is explicit because it will not in general be possible to algebraically isolate what we want to calculate, i.e.
, except for very specific forms of f. We can rearrange the stencils to better demonstrate their geometric meaning:
(8)
These forms illustrate that the forward Euler stencil (first equation) is a simple approximation of the first derivative of
in the interval
using ode evaluated at the left hand side of the interval. Conversely, the backward Euler stencil uses the differential equation to evaluate the derivative at the righthand side of the interval [8] [9].
2.3. Trapezoidal Method
There is another stencil we can derive that is the average of these two approaches:
(9)
This is called the trapezoidal method, and it uses ode evaluated at both ends of the interval to approximate
. Like the backward Euler stencil, it represents an implicit scheme.
3. Numerical Formulation
As an example, we can look at what these stencils look like for the special case of
, where
is a constant:
For this special case, we see that it is possible to isolate
for each stencil. Let’s do this using a loop:
Now, lets go back to the general case by undefined f:
We now want to determine what the error is in each stencil. To do this, we need to rewrite each of them in the standard form
therefore,
Note that these relations define our numeric stencils in these sense that they give exact relations between the approximations
and
. But if we replace the approximations with the exact values of
they are supposed to represent, the above will represent only approximate equalities. That is, if we make the changes
,
,
,
, the left hand sides of the above equations will only be approximately equal to zero. We call the magnitude of the discrepancy the “one step error” or “local error” in the stencil. Hence, the one step error in the various stencils will be:
To estimate the magnitudes of these errors, we expand each of the above expression as a power series about
(since we are implicitly assuming h is a small quantity):
where,
where,
Notice that first, second and third derivatives of y appear explicitly in these expressions. The first derivative terms can be removed by making use of ode. The second derivative can also be removed by examining the derivative of ode:
and
Similarly, the third derivative can be removed
We now substitute our formulae for the derivatives of y into the error expressions:
(10)
where,
(11)
where,
(12)
where,
We see that the local error for the forward and backward schemes is
. Furthermore, the leading order terms for those stencils are the negative of the other one. Since the trapezoidal scheme is the average of the forward and backward algorithms, this explains why the leading order error for the trapezoidal scheme is
. In other words, the trapezoidal scheme is more accurate than the other two. Note that if all we are after is the magnitude of the leading order terms in the errors, we can just expand the above in a low order series:
(13)
(14)
(15)
We now turn our attention to actual numerical algorithms employing these stencils. Since the forward Euler scheme is explicit, it is particularly easy to develop some code for it that works for arbitrary functions f. It is useful to rewrite the stencil with
isolated on the LHS:
Now we can transcribe this as a mapping that takes the step size and the old values of
and
and returns the new value
We could have equivalently accomplished this by using the unapply command, which converts expression into mappings without having to re-write them manually:
Here is a procedure that calculates the numerical solution of ode in the interval
once the form of
has been fixed. Its arguments are initial data in the form
and the number of steps N. The output is a list of points, which we compare in the plot to the output generated by numeric for the same problem.
Example: Solve the equation
Solution:
Of course Figure 1 it is really comparing the results of two numeric approximations to the true solution of the problem. It is also useful to compare numeric answers to analytic results (if available). For the above choice of
there is an analytic solution in terms of Airy functions:
Figure 1. Comparison between forward euler and dsolve numeric.
(16)
where
Here is some code that generates a movie comparing how the numerical solution converges to the analytic solution as the number of cells is increased:
We would also like to examine the performance of the backward and trapezoidal stencils see Figure 2, so let’s make a simple choice of
that admits an analytic solution allows each stencil to be solve for
explicitly:
Figure 2. Comparison between backward and trapezoidal and dsolve numeric.
As for the above code, it will be useful to realize each stencil as a mapping with arguments
and
. We can do this using a loop and the unapply command:
Conversely, we could have accomplished the same thing in one line by using the map command, which applies the same mapping onto each element of a list (in this case we are applying operations onto each element of Stencils to generate a new list stencil):
It is interesting to note that even though each of the stencils give different expressions for
, they actually agree with one another to order (
). To see this, let’s perform some Taylor series expansions:
Now, here is a procedure Euler solution that calculates a numeric solution for
for
using N cells and assuming
and a given value of
. The value of choice dictates which stencil is used: 1 for forward Euler, 2 for backward Euler, and 3 for trapezoidal. Here is a plot of the numeric output for each stencil compared to the actually solution (because we are going to plot 4 curves on our graph). It seems as if the trapezoidal method performs much better than the other two.
,
,
.
We can quantify exactly how well the various stencils are doing by comparing the numeric and actual values for y at some fixed value of x, say
. We call this the global error in the numeric solution at
, which we denote by
In this procedure, note that the
suffix after Euler solution
has the effect of picking out the last element of the list (which is itself a list of 2 quantities), and then picking out the second element of that sub-list (which is the numerical answer for
. We now generate a log-log plot of the global errors for each stencil as a function of N.
For
, the error curves look linear on this log-log plot. This implies that above some threshold number of steps, the errors obey an approximate power law
where a and b are constants. We can determine the values of these parameters by fitting a power law to our error data using statistics. Figure 3 and Figure 4 illustrate four solutions Forward Euler, Backward Euler, Trapezoidal and analytic solution and their errors.
(PowerFit). We first need to regenerate our data for
, the range for which we believe a power law ought to be valid. Then, we use the Statistics (PowerFit) command to perform the fit. Note that this function requires the input to be a matrix, which is why we use the convert Matrix structure. The raw output
Figure 3. Four solutions forward euler, backward euler, trapezoidal and analytic solution and dsolve numeric.
Figure 4. For
, the error curve and dsolve numeric.
of the fitting algorithms is shown for each stencil, as well as the value of b in the power-law
.
.
.
.
We see that for asymptotically large numbers of steps
, the global errors are
for the forward and backward Euler stencils, and
for the trapezoidal method. This confirms the general expectation that for a stencil with one step error
, we expect the global error to obey
Recall that we showed in (13, 14, 15) above that
for forward and backward Euler, and
for the trapezoidal algorithm.
4. Conclusion
In this paper, we have investigated the solution of ODEs using three approaches. The first one, is the Forward Euler method, the second is Backward Euler method and the third one is a numerical method based on Trapezoidal method. In general, we find the numerical method of Trapezoidal method, is a good technique that helped us to find an approximation of the exact solution with small error. In the last, we note the Trapezoidal method.
Acknowledgements
The author would like to thank the reviewers for their valuable comments which have improved the paper.