euler = Plot[{x/3, -x/3}, {x, -1, 1}, Axes -> False, AspectRatio -> 1/3, PlotStyle -> {{Thick, Red}, {Thick, Red}}]; line = Graphics[{Thick, Line[{{-1, 0.2}, {-0.8, 0.2}, {-0.6, 0}, {-0.4, -0.1}, {-0.2, -0.0 .5}, {0, 0}, {0.2, 0.01}, {0.4, 0.1}, {0.6, 0.15}, {0.8, 0.1}, {1.0, 0.2}}]; point = Graphics[{PointSize[Large], Point[{{-1, 0.2}, {-0.8, 0.2}, {-0.6, 0}, {-0.4, -0.1}, {-0.2, -0.0 .5}, {0, 0}, {0.2, 0.01}, {0.4, 0.1}, {0.6, 0.15}, {0.8, 0.1}, {1.0, 0.2}}]}]; Show[line, euler, point]

Suppose that \( |f| \le M, \) where M is a positive constant. In this case the segments in the Euler polygon have slopes between -M and M and the polygon lies between two lines of slopes ±M through the point (x₀, y₀). It is natural to assume that the solution has the same property. This surmise can be proved without using the Euler polygon and without assuming that the solution is unique (Cauchy--Lipschitz theorem).

This algorithm can be accomplished either directly

Example: We illustrate 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 the 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:

\[ y' = x^2 y - 1.5 \,y , \qquad y(0) =1. \]

Since the given differential equation is linear, we can find its solution explicitly:

DSolve[{y'[x] == x*x*y [x] - 1.5*y[x], y[0] == 1}, y, x]

\[ y = e^{1.5\,x - x^3 /3} . \]

So we define its true solution in Mathematica:

true[x_] = E^(-1.5` x + 0.3333333333333333` x^3)

We substitute the slope function \( f(x,y) = x^2 y - 1.5 \,y \) into Euler's rule to obtain

\[ y_{n+1} = y_n + h\,f\left( x_n , y_n \right) = y_n + h \left( x_n^2 y_n - 1.5\, y_n \right) , \]

So application of Euler's rule for the above initial value problem, we will find an approximate solution at a discrete set of nodes 0 = x₀ < x₁ < ··· <x_N. The amount the independent variable changes each time is x_n - x_n-1 and is called the step size. Using the step length uniform, we get

\[ x_{n+1} = x_n + h = \left( n+1 \right) h \]

because we start with zero: 0 = x₀. The solution y(x) to the given initial value problem at the mesh point x_n is estimated by y_n. Factoring the slope function f(x_n,y_n), we get the linear difference equation of first order:

\[ y_{n+1} = y_n \left( 1 + h\,x_n^2 - 1.5\,h \right) = y_n \left( 1 + h^3\,n^2 - 1.5\,h \right) , \qquad n=0,1,2,\ldots ; \quad y_0 = 1. \]

We discussed these equations in section, and it has the explicit solution expressed via product of all intermediate terms:

\[ y_{n+1} = \prod_{k=0}^n \left( 1 + h^3\,k^2 - 1.5\,h \right) , \qquad n=0,1,2,\ldots . \]

Mathematica is able to solve the above difference equation, but it is represented through a special function.

sol = RSolveValue[{a[n + 1] == a[n]*(1 + n^2 *h^3 - 3/2*h), a[0] == 1}, a[n], n]

Let us check a few first terms:

\begin{align*} y_1 &= 1 - 1.5\,h , \\ y_2 &= \left( 1 - 1.5\,h \right) \left( 1 + h^3 - 1.5\,h \right) , \\ y_3 &= \left( 1 - 1.5\,h \right) \left( 1 + h^3 - 1.5\,h \right) \left( 1 + 4\,h^3 - 1.5\,h \right) , \\ y_4 &= \left( 1 - 1.5\,h \right) \left( 1 + h^3 - 1.5\,h \right) \left( 1 + 4\,h^3 - 1.5\,h \right) \left( 1 + 9\,h^3 - 1.5\,h \right), \\ \vdots & \qquad \vdots . \end{align*}

We plot these discrete values and compare with the true solution.

h = 0.1;
DiscretePlot[sol, {n, 0, 27}, PlotStyle -> {PointSize[0.02]}]
soldif = NDSolve[{y'[x] == x^2 *y[x] - 1.5*y[x], y[0] == 1}, y, {x, 0, 2.8}]
Plot[Evaluate[y[x]] /. soldif, {x, 0, 2.7}, PlotStyle -> {Black, Thick}, PlotRange -> {0, 3}]

Euler's discrete approximation	True solution

Now we are ready to use Mathematica. First, we do preprocessing by plugging in the numerical values:

f[x_, y_] := y*x^2 - 1.5*y; (* slope function *)
x0 := 0 (* starting point in x *)
y0 := 0 (* starting value for y *)
xf=2.0; (* Value of x at which y is desired *)
h = (xf - x0)/5.0 (* step size in x *)

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:

soln = DSolve[{f[x, y[x]] == y'[x], y[x0] == y0}, y, x]
plot = Plot[Evaluate[y[x] /. soln], {x, x0, xf},
PlotStyle -> {{Thickness[0.01], RGBColor[1, 0, 0]}}]

Now we implement Euler's algorithm:

X[0] = x0; Y[0] = y0;
X[1] = X[0] + h
Y[1] = Y[0] + f[X[0], Y[0]]*h
plot1 = ListPlot[{{X[0], Y[0]}, {X[1], Y[1]}}, Joined -> True,
PlotStyle -> {Thickness[0.005], RGBColor[0, 0, 1]}, DisplayFunction -> Identity]

The procedure is repeated using the next x value (x[2] = x[1] + h) and the y value approximation in the previous step.

X[2] = X[1] + h
Y[2] = Y[1] + f[X[1], Y[1]]*h
plot2 = ListPlot[{{X[0], Y[0]}, {X[1], Y[1]}, {X[2], Y[2]}}, Joined -> True,
PlotStyle -> {Thickness[0.005], RGBColor[0, 0, 1]}, DisplayFunction -> Identity]

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.

X[3] = X[2] + h
Y[3] = Y[2] + f[X[2], Y[2]]*h
plot3 = ListPlot[{{X[0], Y[0]}, {X[1], Y[1]}, {X[2], Y[2]}, {X[3], Y[3]}}, Joined -> True,
PlotStyle -> {Thickness[0.005], RGBColor[0, 0, 1]}, DisplayFunction -> Identity]

On the fourth step,we have

X[4] = X[3] + h
Y[4] = Y[3] + f[X[3], Y[3]]*h
plot4 = ListPlot[{{X[0], Y[0]}, {X[1], Y[1]}, {X[2], Y[2]}, {X[3], Y[3]}, {X[4], Y[4]}}, Joined -> True,
PlotStyle -> {Thickness[0.005], RGBColor[0, 0, 1]}, DisplayFunction -> Identity]

Finally, we plot Euler's approximations along with the explicit solution:

Show[plot, plot4, PlotRange -> Automatic]

In the above codes, a special option DisplayFunction was used. 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.

DisplayFunction -> $DisplayFunction default setting
DisplayFunction -> Identity generate no display
DisplayFunction -> f apply f to graphics objects to produce display

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.

The results of applying Euler's rule with h = 0.1

n	xn	yn	Exact solution	Absolute error
1	0.1	0.85	0.860995	0.0109949
2	0.2	0.72335	0.742796	0.0194464
3	0.3	0.617741	0.643393	0.0256518
4	0.4	0.530639	0.560645	0.030006
5	0.5	0.459534	0.492464	0.0329305
6	0.6	0.402092	0.436922	0.0348302
7	0.7	0.356254	0.392324	0.0360707
8	0.8	0.320272	0.357245	0.0369731
9	0.9	0.292729	0.330549	0.0378206
10	1.0	0.27253	0.311403	0.0388729

   ■

Example: Consider the initial value problem for the logistic equation

\[ \frac{{\text d}y}{{\text d}t} = y \left( 1-y \right) , \qquad y(0) = y_0 =0.5. \]

This problem has a unique solution

\[ y(t) = \frac{C\, e^t}{C\,e^t +1} , \qquad \mbox{where} \qquad C= \frac{y_0}{1-y_0} . \]

Using the following Mathematica code

y[0] = 0.5; h = 0.1
y[1] = y[0] + h*y[0]*(1 - y[0])
Do[y[n] = y[n - 2] + 2*h*y[n - 1]*(1 - y[n - 1]), {n, 2, 500}]
ListPlot[Table[{n*h, y[n]}, {n, 0, 500}]]

we observe instability

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

\[ y(x_{n+1}) = y(x_n ) + \int_{x_n}^{x_{n+1}} f(s, y(s))\,{\text d} s , \qquad n=0,1,2,\ldots . \]

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

true[x_] = Exp[x]/(1 + Exp[x])
y[0] = 0.5; h = 0.1; f[y_] = y^2 - y;
Do[y[n + 1] = y[n] + h f[y[n]], {n, 0, 500}]
y[1]

The results of applying Euler's rule with h = 0.1

n	xn	yn	Exact solution	Absolute error
1	0.1	0.475	0.524979	0.0499792
2	0.2	0.450063	0.549834	0.0997715
5	0.5	0.376852	0.622459	0.245607
10	1.0	0.266597	0.731059	0.464462
20	2.0	0.114573	0.880797	0.766224
30	3.0	0.0435246	0.952574	0.90905
50	5.0	0.00552563	0.993307	0.987782
100	10.0	0.0000286528	0.999955	0.999926
200	20.0	7.61083×10-10	1.	1.
400	40.0	5.3695×10-19	1.	1.

   ■

Example: We consider the following initial value problem for the Riccati equation

\[ y' = y^2 - x^2 , \qquad y(0) = 1/2 . \]

First, we find its explicit solution using mathematica standard command DSolve:

explicit = DSolve[{y'[x] == y[x]^2 - x^2 , y[0] == 1/2}, y, x] true[x_] = y[x] /. explicit Since the expression for the explicit solution involves special functions, we just evaluate its one particular value at x = 0.5 and plot the function true[0.5] {0.617769 + 2.77556*10^-17 I} Plot[true[x], {x, 0, 1}, PlotStyle -> Thick]

Now we apply Euler's rule using our subroutine euler:

euler[y^2 - x^2, {0, 1}, {0.5}, 10]

{{0, 0.5}, {0.1, 0.525}, {0.2, 0.551563}, {0.3, 0.577985}, {0.4, 0.602391}, {0.5, 0.622679}, {0.6, 0.636452}, {0.7, 0.640959}, {0.8, 0.633042}, {0.9, 0.609116}, {1., 0.565218}}

Therefore, we Euler's approximation to be 0.622679.

The results of applying Euler's rule with h = 0.1

n	xn	yn	Exact solution	Absolute error
1	0.1	0.525	0.525974	0.000973594
2	0.2	0.551563	0.552738	0.0011751
3	0.3	0.577985	0.578417	0.000432741
4	0.4	0.602391	0.600901	-0.00149051
5	0.5	0.62267	0.617769	-0.00490958
6	0.6	0.636452	0.626228	-0.0102239
7	0.7	0.640959	0.623052	-0.0179072
8	0.8	0.633042	0.604575	-0.028467
9	0.9	0.609116	0.566763	-0.0423522
10	1.0	0.565218	0.505431	-0.0597873

   ■

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

\[ \dot{y} = 2.3\, y -0.01\,y^2 -0.1\,y\, \int_0^t y(\tau )\,{\text d}\tau , \qquad y(0) =50. \]

Choosing the uniform grid \( t_k = k\,h , \quad k=0,1,2,\ldots , m; \) we integrate both sides to obtain the integral equation

\[ y(t ) = 50 + 2.3\,\int_0^t y(s)\,{\text d}s -0.01\,\int_0^t y^2 (s) \,{\text d}s -0.1\,\int_0^t y(s)\,{\text d}s \int_0^s y(\tau )\,{\text d}\tau \]

for each mesh point t = t_k. Since the double integral can be written as

\[ \int_0^t y(s)\,{\text d}s \int_0^s y(\tau )\,{\text d}\tau = \int_0^t \,{\text d}\tau \,\int_{\tau}^t {\text d}s \, y(s)\,y(\tau ) , \]

application of the left rectangular rule yields for the first mesh point t₁ = h:

\[ y(t_1 ) \approx 50 + 2.3\,h\,y(0) -0.01\,h\,y^2 (0) -0.1 \, h\, \int_0^{t_1} {\text d}s \, y(s)\,y(0 ) \approx 50 +2.3\,h\,50 -0.01\,h \,50^2 - 0.1\,h^2 \,50^2 . \]

So we choose the right hand-side as the approximate value y₁.

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

\[ y_{k+1} = y_k + 2.3\,h\,y_k -0.01\,h\,y^2_k -0.1 \, y_k \,h \int_0^{t_k} {\text d}s \, y(s)\qquad k=1,2,\ldots . \]

If the trapezoidal rule is used to approximate the integral, then this expression becomes

\[ y_{k+1} = y_k + 2.3\,h\,y_k -0.01\,h\,y^2_k -0.1 \, y_k \,h\, T_k (h) ,\qquad k=1,2,\ldots ; \]

where \( T_0 (h) = h^2 \,50^2 \) and

\[ T_{k} = T_{k-1} (h) + \frac{h}{2} \left( y_k + y_{k-1} \right) ,\qquad k=1,2,\ldots ; \]

Finally, we ask Mathematica to perform all calculations:

h := 0.001
T[0] := h^2 *50^2
y[0] := 50
y[1] := y[0] + 2.3*h*50 - 0.01*h*50^2 - 0.1*h^2*50^2
Do[ T[k] = T[k - 1] + h*(y[k] + y[k - 1])/2;
y[k + 1] = y[k] + 2.3*h*y[k] - 0.01*h*(y[k])^2 - 0.1*y[k]*h*T[k], {k, 1, 2000}]
ListPlot[Table[{k*h, y[k]}, {k, 1, 2000}]]

   ■

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:

Plot with some options:

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 actually y(0.1), then shift+return, and we have our nice answer.

First, we start with the 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. When thinking of parameters that 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 built-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:

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. \)

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 Mathematica command N. Of course you can readily change that or accomplish the corresponding value 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 illustrate 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:

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:

Return to Mathematica page
Return to the main page (APMA0330)
Return to the Part 1 (Plotting)
Return to the Part 2 (First Order ODEs)
Return to the Part 3 (Numerical Methods)
Return to the Part 4 (Second and Higher Order ODEs)
Return to the Part 5 (Series and Recurrences)
Return to the Part 6 (Laplace Transform)
Return to the Part 7 (Boundary Value Problems)