Preface

---

This is a tutorial made solely for the purpose of education and it was designed for students taking Applied Math 0330. It is primarily for students who have very little experience or have never used Mathematica before and would like to learn more of the basics for this computer algebra system. As a friendly reminder, don't forget to clear variables in use and/or the kernel.

Finally, the commands in this tutorial are all written in bold black font, while Mathematica output is in normal font. This means that you can copy and paste all commands into Mathematica, change the parameters and run them. You, as the user, are free to use the scripts for your needs to learn the Mathematica program, and have the right to distribute this tutorial and refer to this tutorial as long as this tutorial is accredited appropriately.

Euler's Methods

---

We start with the first numerical method for solving initial value problems that bears Euler's name. Leonhard Euler was born in 1707, Basel, Switzerland and passed away in 1783, Saint Petersburg, Russia. In 1738, he became almost blind in his right eye. Euler was one of the most eminent mathematicians of the 18th century, and is held to be one of the greatest in history. He is also widely considered to be the most prolific mathematician of all time. He spent most of his adult life in Saint Petersburg, Russia, except about 20 years in Berlin, then the capital of Prussia.

In 1768, Leonhard Euler (St. Petersburg, Russia) introduced a numerical method that is now called the Euler method or the tangent line method for solving numerically the initial value problem:

where f(x,y) is the given slope (rate) function, and \( (x_0 , y_0 ) \) is a prescribed point on the plane. Euler's method or rule is a very basic algorithm that could be used to generate a numerical solution to the initial value problem for first order differential equation. The solution that it produces will be returned to the user in the form of a list of points. The user can then do whatever one likes with this output, such as create a graph, or utilize the point estimates for other purposes. Euler's rule serves to illustrate the concepts involved in the advanced methods. It has limited usage because of the larger error that is accumulated as the process proceeds, so it requires a smaller step size. However, it is important to study because the error analysis is easier to understand.

To start, we need mesh or grid points, that is, the set of discrete points for independent variable at which we find approximate solutions. In other words, we will find approximate values of the unknow solution at these mesh points. For simplicity, we use uniform grid with fixed step length h; however, in practical applications, it is almost always not a constant. For convenience, we subdivide the interval of interest [a,b] (where usually \( a= x_0 \) is the starting point) into m equal subintervals and select the mesh/grid points

The considered initial value problem is assumed to have a unique solution \( y= \phi (x) , \) on the interval of interest {a,b], and its approximations at the grid points will be denoted by y_n, so we wish that \( y_n \approx \phi (x_n ) , \quad n=1,2, \ldots . \)

There are three main approaches (we do not discuss others in this section) to derive the Euler rule: either use finite difference approximation to the derivative or transfer the initial value problem to the Volterra integral equation and then truncate it, or apply Taylor series.

1. We start with finite difference method. If we approximate the derivative in the left-hand side of the differential equation y' = f(x,y) by the finite difference
   \[ y' (x_n ) \approx \frac{y_{n+1} - y_n}{h} \]
   on the small subinterval \( [x_{n+1} , x_n ] , \) we arrive at the Euler's rule when right-hand side is evaluated at x = x_n.
   \begin{equation*} y_{n+1} = y_n + (x_{n+1} - x_n ) f( x_n , y_n ) \qquad \mbox{or} \qquad y_{n+1} = y_n + h f_n , \end{equation*}
   where the following notations are used: \( h=x_{n+1} - x_n \) is the step length (which is assumed to be constant for simplicity), \( f_n = f( x_n , y_n ) \) is the value of slope function at mesh point, and \( y_n \) denotes the approximate value of the actual solution \( y=\phi (x) \) at the point \( x_n \) (\( n=1,2,\ldots \) ).
2. Integrating the identity y'(x) = f(x,y(x)) between x_n and x_n+1, we obtain
   \[ y(x_{n+1}) - y(x_n ) = \int_{x_n}^{x_{n+1}} f(t, y(t))\,{\text d} t . \]
   Approximating the integral by the left rectangular (very crude) rule for numerical integration (length of interval times value of integrand at left end point) and identifying \( y(x_n ) \) with y_n, we again obtain the Euler formula.


3. We may assume the possibility of expanding the solution in a Taylor series around the point x_n:
   \[ y(x_{n+1}) = y(x_n +h) = y(x_{n}) + h\, f(x_n, y(x_n )) + \frac{1}{2}\, h^2 \,y'' (x_n ) + \cdots . \]
   The Euler formula is the result of truncating this series after the linear term in h.

Each of these interpretations points the way to a class of generalizations of Euler's method that are discussed later. It is interesting to note that the generalization indicated by (i) which seems to be the most straightforward, has proved to be the least fruitful of the three.

When the slope function is evaluated at the right end point, we come to the so called backward Euler rule:

However, this is not the only one approximation possible, one may consider the central difference: which usually gives more accurate approximation. Using central difference, we get two-step approximation: which requires two starting points to solve this recurrence of second order. So one can use standard Euler rule to determine y₁: \( y_1 = y_0 + h\, f\left( x_0 , y_0 \right) . \) However, this numerical algorithm may produce instability and we may observe numerical solution that is different from the true solution. Such solution is called "ghost solution." This phenomenon is caused by roundoff errors.

To avoid the instability and to prevent the appearance of ghost solutions, many refine numerical procedures to integrate the ordinary differential equations have been proposed. In practical calculations, the central difference scheme is used when evaluations of slope function in intermediate points (as, for instance, required by Runge--Kutta methods) are not permitted. We mention two remedies to remove ghost solutions.

Consider the following mixed difference scheme parameterized by μ, 0 <μ<1, for the first order differential equation: If we set μ = 1, then we have Euler's rule. When μ =0, the scheme is nothing but the central difference scheme. This algorithm does not produce ghost solutions when μ > h/2.

Let φ be a real-valued function on \( \mathbb{R} \) that satisfies the property:

There exists a variety of functions φ that satisfy the above condition, e.g., \( \phi (h) = 1 - e^{-h} \quad \Longrightarrow \quad \phi (hq)/q = \left( 1- e^{-hq} \right) /q . \) Then the following nonstandard schemes are stable:

• nonstandard explicit Euler method given by
  \[ \frac{y_{k+1} - y_k}{\phi (hq)/q} = f(y_k ), \qquad y_0 = y(0), \quad k=0,1,2,\ldots . \]
• nonstandard implicit Euler method given by
  \[ \frac{y_{k+1} - y_k}{\phi (hq)/q} = f(y_{k+1} ), \qquad y_0 = y(0), \quad k=0,1,2,\ldots . \]

Example: Consider the initial value problem for the logistic equation

This problem has a unique solution Using the following Mathematica code we observe instability

Another approach is based on transferring the given IVP \( y' = f(x,y) , \quad y(x_0 ) = y_0 \) to equivalent integral equation

Application of left Riemann approximation, we come to the Euler's rule.

Example: We demonstrate how the Euler method works on an example of a linear equation when all calculations become transparent. Note that Mathematica attempts to solve the ODE using its function. A solution may or may not be possible using this function. In case when a solution may not be known or determined, Euler's method provides an alternative to approximate the solution numerically. Consider the initial value problem:

First, we do preprocessing by plugging in the numerical values: We we attempt to find an approximation for x=2 with a step size of 0.4. Then we find the explicit expression for the formulated above IVP using standard Mathematica command and plot it: Now we implement Euler's algorithm: The procedure is repeated using the next x value (x[2] = x[1] + h) and the y value approximation in the previous step. The procedure on the third step is repeated using the next x value (x[3] = x[2] + h) and the y value approximation in the previous step. On the fourth step,we have Finally, we plot Euler's approximations along with the explicit solution:

In the above codes, a special option DisplayFunction was used. Actually, you can either remove this option or replace with a standard one: DisplayFunction -> $DisplayFunction. All Mathematica graphics functions such as Show and Plot have an option DisplayFunction,which specifies how the Mathematica graphics and sound objects they produce should actually be displayed. The way this works is that the setting you give for DisplayFunction is automatically applied to each graphics object that is produced.

Within the Mathematica kernel, graphics are always represented by graphics objects involving graphics primitives. When you actually render graphics, however, they must be converted to a lower-level form which can be processed by a Mathematica front end, such as a notebook interface, or by other external programs. The standard low-level form that Mathematica uses for graphics is PostScript. The Mathematica function Display takes any Mathematica graphics object, and converts it into a block of PostScript code. It can then send this code to a file, an external program, or in general any output stream.

We combine all steps into subroutine: In the above code,
n is number of steps
x0 = initial condition for x
y0 = initial condition for y
xf = final value of x at which y is desired.
f = slope function for the differential equation dy/dx = f(x,y) This loop calculates the following
AV = approximate value of the ODE solution using Euler's Method by calling the module EulerMethod
Et = true error
et = absolute relative true error percentage
Ea = approximate error
ea = absolute relative approximate error percentage
sig = least number of significant digits correct in approximation

The following code uses Euler's Method to calculate intermediate step values for the purpose of displaying the solution over the range specified. The number of steps used is the maximum value, n (in this example, n=128).

Finally, we plot the exact solution along with Euler's approximation:

Now we compare the value at the final point (denoted by xf) with a step size

We plot dependences of errors with the step size.

Example: To demonstrate the latter approach, we consider the initial value problem for the integro-differential equation

Choosing the uniform grid \( t_k = k\,h , \quad k=0,1,2,\ldots , m; \) we integrate both sides to obtain the integral equation for each mesh point t = t_k. Since the double integral can be written as application of left rectangular rule yields for the first mesh point t₁ = h: So we choose the right hand-side as the approximate value y₁.

For the general step in Euler's rule, we have

If the trapezoidal rule is used to approximate the integral, then this expression becomes where \( T_0 (h) = h^2 \,50^2 \) and Finally, we ask Mathematica to perform all calculations:

There are several ways to implement the Euler numerical method for solving initial value problems. We demonstrate its implementations in a series of codes. To start, define the initial point and then the slope function:

Next, define the step size:

Now we define the Euler method itself:
Create the table of approximations using Euler's rule:

Plot with some options:

or
or

Another way is to make a loop.

First of all, it is always important to clear all previous assignment to all the variables that we are going to use, so we have to type: Clear[y]

The basic structure of the loop is:
Do[some expression with n, {n, starting number, end number}]
or with option increment, denoted by k:
Do[some expression with n, {n, starting number, end number, k}]

The function of this Do loop is to repeat the expression, with n taking values from “starting number” to “ending number,” and therefore repeat the expression (1+ (ending number) - (starting number))  many times. For our example, we want to iterate 99 steps, so n will go from 0, 1, 2, ..., until 99, which is 100 steps in total. There is just one technical issue: we must have n as an integer. Therefore we let y[n] denote the actual value of y(x) when x = n*h. This way everything is an integer, and we have our nice do loop.
Finally put y[10] , which is actualy y(0.1), then shift+return, and we have our nice answer.

First, we start with output of our program, which is, perhaps, the most important part of our program design. You'll notice that in the program description, we are told that the "solution that it produces will be returned to the user in the form of a list of points. The user can then do whatever one likes with this output, such as create a graph, or utilize the point estimates for other purposes.

Next, we must decide upon the information that must be passed to the program for it to be able to achieve this. In programming jargon, these values are often referred to as parameters. So what parameters would an Euler program need to know in order to start solving the problem numerically? Quite a few should come to mind:

The slope function f(x, y)
The initial x-value, x₀
The initial y-value, y₀
The final x-value,
The number of subdivisions to be used when chopping the solution interval into individual jumps or the step size.

To code these parameters into the program, we need to decide on actual variable names for each of them. Let's choose variable names as follows:

f, for the slope function f(x, y)
x₀,      for the initial x-value, x₀
y₀,      for the initial y-value, y₀
x_n,      for the final x-value
Steps, for the number of subdivisions
h,        for step size

There are many ways to achieve this goal. However, it is natural to have the new subroutine that looks similar to build-in Mathematica's commands. So our program might look something like this:

In order to use Euler's method to generate a numerical solution to an initial value problem of the form:

We have to decide upon what interval, starting at the initial point x₀, we desire to find the solution. We chop this interval into small subdivisions of length h, called step size. Then, using the initial condition \( (x_0,y_0) \) as our starting point, we generate the rest of the solution by using the iterative formulas:

to find the coordinates of the points in our numerical solution. The algorithm is terminated when the right end of the desired interval has been reached.

If you would like to use a built-in Euler’s method, unluckily there is none. However, we can define it ourselves. Simply copy the following command line by line:

Now we have our euler function:   euler[f(x,y), {x,x0,x1},{y,y0},steps]
Then this script will solve the differential equation y’=f(x,y), subject to the initial condition y(x0)=y0, and generate all values between x₀ and x1. The number of steps for the Euler’s method is specified with steps.
To solve the same problem as above, we simply need to input:

Example: Consider the initial value problem \( y' = x+y, \quad y(0) =1. \)

As you can see, the output is really a big table of values, and you can really just read the last one off to get y(1). If the only one final value is needed, then Return command should be replaced with

We can plot the output as follows

Or we can plot it as

The following code, which uses a slightly different programming paradigm, implements Euler's method for a system of differential equations:

And the usage message is:

euler::usage = "euler[F, t0, Y0, b, n] gives the numerical solution to {Y' == F[t, Y], Y[t0] == Y0} over the interval\n
[t0, b] by the n-step Euler's method. The result is in the form of a table of {t, Y} pairs."

Note that this function uses an exact increment h rather than converting it explicitly to numeric form using N. Of course you can readily change that or accomplish the corresponding thing by giving an approximate number for one of the endpoints.

Next we plot the points by joining them with a curve:

Another way without writing a subroutine:

The results can also be visualized by connecting the points:

Next, we demonstrate application of Function command:

Here is another steamline of the Euler method (based on Mathematica function), which we demonstrate in the following

Example: Consider the initial value problem \( y' = y^2 - 3\, x^2 , \quad y(0) =-1. \)

Now we can plot our results, making sure we only refer to values of x and y that we have defined:

Finally, we plot Euler approximations along with the actual solution:

Now we are going to repeat the problem, but using Mathematica lists format instead. The idea is still the same---we set a new x value, compute the corresponding y value, and save the two quantities in two lists. In particular, in a list we are simply storing values; In previously presented code, we were actually storing rules, which are more complicated to store and evaluate. We start out with just one value in each list (the initial conditions). Then we use the Append command to add our next pair of values to the list.

Now we can plot our results:

Euler's method is first-order accurate, so that errors occur one order higher starting at powers of \( h^2 . \)
Euler's method is implemented in NDSolve as ExplicitEuler:

 

 

II. Backward Euler Method

---

 

Backward Euler formula: where h is the step size (which is assumed to be fixed, for simplicity) and \( f_{n+1} = f \left( x_{n+1}, y_{n+1} \right) . \)

Example: Consider the following initial value problem:

Here is the Mathematica code that solve this problem:
Then we plot solution:

NSolve has to spend time to compute all roots to the equation (which can be computationally expensive). FindRoot does a pretty fast search looking for only a single root, so it is quick for complex equations.
If you care about all possible roots, or if you have no clue where the roots of the equation may be, FindRoot is a terrible choice. If you only care about a single root and have a rough idea of where it might be, though, FindRoot will find it quickly.