Euler Method - Mathematica Implementation Part 3

Numerical Methods for Solving Differential Equations

Euler's Method

Using the Method with Mathematica

(continued from last page...)

Let's build at a very basic program that could be used to generate a numerical solution to a first order initial value problem of the form:

y′ = f(x, y)

y(x_o) = y_o

using Euler's Method. The solution that it produces will be returned to the user in the form of a list of points.

Designing the Program's Output

The form of the output of our program is, perhaps, the most important part of our program design. After all, it's for this that the user of our program is doing everything. (Even if that user happens to be us.)

You'll notice that in the program description given above, we are told that the "solution that it produces will be returned to the user in the form of a list of points." Well, that's easy enough to handle!

Based on what we already know about Mathematica's handling of lists and points, we would expect an typical solution to look something like this:

{{0,0},{0.25,0},{0.5,0.0625},{0.75,0.21875},{1,0.515625}}

The user can then do whatever he likes with this output, such as create a graph, or utilize the point estimates for other purposes.

Designing the Program's Input

So now that we have a goal, in the form of the program's output form, 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 function f(x, y)
The initial x-value, x_o
The initial y-value, y_o
The final x-value
The number of subdivisions to be used when chopping the solution interval into individual jumps

I think that's about it! 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 function f(x, y)
x0, for the initial x-value, x_o
y0, for the initial y-value, y_o
xn, for the final x-value
steps, for the number of subdivisions

Now let's think about the form of command we'd like to be able to give Mathematica whenever we want the Euler Method numerical solution to a problem. It would be nice if the new program that we are building behaved in a similar fashion to the way Mathematica's built in functions behave. Think about the way the Plot command that we've used so many times is structured. An example of it's use is:

Plot[x^2+3x, {x,-2,5}]

The parameters that we pass to the Plot command here are:

x^2+3x, the function which is to be plotted
{x,-2,5}, a list consisting of three parts:
- x, the variable to be used for the plot (any other variable names appearing in the function should previously have been given values)
- -2, the starting value of the variable's interval for plotting
- 5, the ending value of the variable's interval for plotting

Maybe we could design our program, let's call it euler, to receive its parameters in a similar way. Say we were told to solve the preliminary example, back in the introduction:

y′ = x + 2y

y(0) = 0

numerically, finding a value for the solution at x = 1, and using 4 steps. Doesn't the following command look like the kind of structure Mathematica would use:

euler[x+2y,{x,0,1},{y,0},4]

if it had a built-in routine for using Euler's Method?

Note that we don't capitalize our function's name, euler, since capitalized names are reserved, by convention, for Mathematica's built in commands. Let's look over the parameters:

x+2y is the function f(x, y) that comes from the initial value problem's differential equation
{x,0,1} is a list consisting of three parts:
- x is the independent variable's name
- 0 is the independent variable's initial value, x_o, which comes from the initial condition of the initial value problem
- 1 is the independent variable's final value, which comes from the description of what is required
{y,0} is a list consisting of two parts:
- y is the dependent variable's name
- 0 is the dependent variable's initial value, y_o, which comes from the initial condition of the initial value problem
4 is the number of steps to be used when chopping up the solution interval into small jumps

Is this making sense? I hope so, since I'm trying to take it very slowly for you beginners out there. (I know, you expert programmers are getting bored. Patience!)

Using what we just designed as an input form for solving this particular problem, let's now generalize. Earlier we gave names to the parameters our program would need. I'll repeat them here:

f, for the function f(x, y)
x0, for the initial x-value, x_o
y0, for the initial y-value, y_o
xn, for the final x-value
steps, for the number of subdivisions

Remember? If we put these into a kind of "input template", based on the last example, the general command for executing our program might look something like this:

euler[f,{x,x0,xn},{y,y0},steps]

Do you agree? Of course, we still have no code written to make the program do anything, but at least we know what the user of our program is going to type in as input.

Getting from Input to Output

So now we come to the hard part! How do we take the program user's input, and produce the expected output? What should be the "guts" of our program?

Perhaps we could describe these "guts" in English before we try doing any actual coding. Recall that we summarized Euler's Method back in the Introduction as follows:

Summary of Euler's Method

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

y′ = f(x, y)

y(x_o) = y_o

we decide upon what interval, starting at the initial condition, we desire to find the solution. We chop this interval into small subdivisions of length h. Then, using the initial condition as our starting point, we generate the rest of the solution by using the iterative formulas:

x_n+1 = x_n + h

y_n+1 = y_n + h f(x_n, y_n)

to find the coordinates of the points in our numerical solution. We terminate this process when we have reached the right end of the desired interval.

Breaking this summary down into small steps we might describe it like this:

Read in the function and the initial values of all of the variables.
Initialize the solution list to contain the initial condition point.
Calculate the step-size, h. This can be found by taking the length of the solution interval, and dividing it by the number of steps to be used.
Find new points and append them to the solution list by using the Euler iteration formulas repeatedly, as follows:
1. Find the "new" x-value from the "old" (previous) x-value, using the formula
  x_n+1 = x_n + h.
2. Find the "new" y-value from the "old" (previous) y-value, using the formula
  y_n+1 = y_n + h f(x_n, y_n).
3. Append the ordered pair consisting of the "new" x-value and y-value to the list of points making up the solution.
4. Reset the the "new" x-value and y-value to be the "old" x-value and y-value, in readiness for the next repetition of the "loop."
Return the completed list of points to the user.

The Program Itself

Now we have the job of translating this list of steps from English into Mathematica commands. The result of the translation looks something like this:

euler[f_,{x_,x0_,xn_},{y_,y0_},steps_]:=

Block[{ xold=x0,
    yold=y0,
    sollist={{x0,y0}},
    x,y,h
  },

  h=N[(xn-x0)/steps];

  Do[ xnew=xold+h;
    ynew=yold+h*(f/.{x->xold,y->yold});
    sollist=Append[sollist,{xnew,ynew}];
    xold=xnew;
    yold=ynew,
    {steps}
  ];

  Return[sollist]
]

Ooh! What a lot of new commands! It looks like we need to carefully sift through them all in order to fully understand how the program works.

Let's go and analyze the program now...

$Compass$

If you're lost, impatient, want an overview of this laboratory assignment, or maybe even all three, you can click on the compass button on the left to go to the table of contents for this laboratory assignment.