Demo #4
Introduction to Numerical Methods:
Euler's Method

In this demo we show how the simple Euler algorithm can be implemented by using Mathematica as a programming language. While this is convenient for this course, you should note that a comparable computer program written in a programming language such as Fortran or C would run much faster, because of particular details and structure of the Mathematica environment. Thus, a complex numerical algorithm that demands a large amount of computer time to run would be best written in one of these other languages. For the present course, however, much of what we would like to do can be accomplished simply and efficiently in the Mathematica environment.

Our objective in the present demo is to solve the following nonlinear ODE:

(dy/dx)=(x^2) - (y^2)

subject to the initial condition y[0] = 1.

Our first step is to define the derivative function for the ODE that we wish to solve:

f[x_,y_] := x^2 - y^2

To approximate the solution with IV y(0)=1 on the interval [0,2], we set the following initial values (the semicolons suppress the output):

x0 = 0.;
y0 = 1.;
xf = 2.;

Next, we set the number of steps and the stepsize:

n = 10
h = ( xf - x0 )/n

The algorithm is started by initializing x and y to the values set by the initial condition:

x = x0
y = y0

Euler's method is implemented by stepping away from the initial condition with the following Do loop:

Do[ k = f[x,y];			(* the left-hand slope *)
	y = y + h*k;		(* Euler step to update y *)
	x = x + h;			(* update x	*)
	Print[ x,"    ",y], (* display current results *)
	{i,1,n} ]

We save these in an array for future comparison to "exact" results:

euler = {{0.2,0.8},{0.4,0.68},{0.6,.61952},
          {0.8,.614759},{1.,.667173},{1.2,.778149},
          {1.4,.945046},{1.6,1.15842},{1.8,1.40203},
          {2.,1.65689}}

Now, to see how accurate our simple Euler routine is, we generate an "exact" numerical solution (to higher accuracy than the six digits shown above):

Clear[x,y]
NDSolve[ {y'[x] == f[x,y[x]],y[0]==1},y,{x,0,2}]

plt1 = Plot[Evaluate[ y[x] /. %],{x,0,2}]

Now, let's see how our accurate our Euler method with the (coarse) step size of h=0.2 is:

plt2 = ListPlot[euler]

Show[ plt1,  plt2]

Thus, we see that the Euler method works surprisingly well even with a coarse mesh of points. By doubling the number of points to n = 20 (and thereby halving the step size h), the Euler and NDSolve results would be indistinguishable.

As a general rule of thumb, one halves the value of h repeatedly until convergence of the y values across the interval is achieved to some prescribed accuracy.