ÐÓ°ÉÊÓƵ

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:

²â′Ìý=Ìýf(³æ,Ìý²â)

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(³æ,Ìý²â)
• The initialÌýx-±¹²¹±ô³Ü±ð,Ìýx_o
• The initialÌýy-±¹²¹±ô³Ü±ð,Ìý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(³æ,Ìý²â)
• x0, for the initialÌýx-±¹²¹±ô³Ü±ð,Ìýx_o
• y0, for the initialÌýy-±¹²¹±ô³Ü±ð,Ìý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}]

°Õ³ó±ðÌý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:

²â′Ìý=Ìý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(³æ,Ìý²â)
• x0, for the initialÌýx-±¹²¹±ô³Ü±ð,Ìýx_o
• y0, for the initialÌýy-±¹²¹±ô³Ü±ð,Ìý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:

²â′Ìý=Ìýf(³æ,Ìý²â)

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Ìý´Ú(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:

1. Read in the function and the initial values of all of the variables.

2. Initialize the solution list to contain the initial condition point.

3. 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.

4. 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Ìý´Ú(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."

5. 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,
ÌýÌýÌýÌý²â´Ç±ô»å=²â0,
ÌýÌýÌýÌý²õ´Ç±ô±ô¾±²õ³Ù=µ÷µ÷³æ0,²â0°¨°¨,
ÌýÌýÌýÌý³æ,²â,³ó
ÌýÌý°¨,

ÌýÌý³ó=±·°Ú(³æ²Ô-³æ0)/²õ³Ù±ð±è²õ±Õ;

ÌýÌýDo[ xnew=xold+h;
ÌýÌýÌýÌý²â²Ô±ð·É=²â´Ç±ô»å+³ó*(´Ú/.µ÷³æ-&²µ³Ù;³æ´Ç±ô»å,²â-&²µ³Ù;²â´Ç±ô»å°¨);
ÌýÌýÌýÌý²õ´Ç±ô±ô¾±²õ³Ù=´¡±è±è±ð²Ô»å°Ú²õ´Ç±ô±ô¾±²õ³Ù,µ÷³æ²Ô±ð·É,²â²Ô±ð·É°¨±Õ;
ÌýÌýÌýÌý³æ´Ç±ô»å=³æ²Ô±ð·É;
ÌýÌýÌýÌý²â´Ç±ô»å=²â²Ô±ð·É,
ÌýÌýÌýÌýµ÷²õ³Ù±ð±è²õ°¨
ÌýÌý±Õ;

ÌýÌý¸é±ð³Ù³Ü°ù²Ô°Ú²õ´Ç±ô±ô¾±²õ³Ù±Õ
]

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Ìý²Ô´Ç·É...