Math 431 • 9 March 1999

Euler's method for systems

⚠  If you did not tell Mathematica to evaluate automatically the initialization cells in this notebook when you opened it, be sure to evaluate the input cells in the first section below.

Procedure implementing Euler's Method for systems

Without essentially no modification whatsoever, the same function euler that implemented Euler's method for estimating solutions to first-order scalar IVPs works for systems of first-order IVPs!  (Lower-case y is changed to upper-case Y to suggest a vector-valued function.)

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

\(euler[F_, \ a_, \ Y0_, \ b_, \ n_]\  := 
	
      Module[{t, \ Y, \ Δt\  = \ \((b\  - \ a)\)\/n, \ i}, 

        		t[0]\  = \ a; \ Y[0]\  = \ Y0; 
		
        Do[
			t[i]\  = \ a\  + \ i\ Δt; 
			
          Y[i]\  = \ \ 
            Y[i\  - \ 1]\  + \ 
              Δt\ F[t[i\  - \ 1], \ Y[i\  - \ 1]], 

          			{i, \ 1, \ n}
		]; 
		
        Table[{t[i], \ Y[i]}, \ {i, \ 0, \ n}]
	]\)

Example: van der Pol equation

The differential equation

The ODE of this example is the "unforced" van der Pol equation

	\(TraditionalForm
      \`\(x'\)'\  - \ \((1\  - \ x\^2)\)\ x'\  + \ x\  = \ 0\)
	
with initial conditions

	x(0)  =  1,     x' (0)  =  1.

The equivalent system of first-order ODEs is

	\(TraditionalForm\`x'\  = \ \ \ \ \ \ \ y, \ \ \ \ \ \ \ \ 

      y'\  = \ \(-\ x\)\  + \ \((1\  - \ x\^2)\)\ y\)
	
with initial conditions:

	x(0)  =  1,  y(0)  =  1.

The vector field associated with the system of ODEs is:

F[t_,  {x_,  y_}]  := 
       {y,  -x  +  (1  -  (x^2)) y}

Applying Euler's method—crude approximation

Start with short t-interval with a small number of steps and a correspondingly modest increment Δt just to compare the results with paper-and-pencil calculations and to see how to manipulate the result of euler into a form for plotting.

estimate  =  euler[F,  0. ,  {1. ,  1. },  1,  4];

A good way to understand the result is to display it as a table consisting of t-values along with the corresponding estimates of (x,  y)-value pairs:

TableForm[estimate]

Here is a function to plot the the calculated, approximate x-values against t.  [More precisely, it "connects the dots" to plot a piecewise linear approximation to the graph of the exact solution x(t).]

\(txPlot[estimate_]\  := \ 
      Module[{xvalues, \ tvalues, \ txs}, \ 
        xvalues\  = \ Map[First, \ Last[Transpose[estimate]]]; \ 
        tvalues\  = \ First[Transpose[estimate]]; \ 
        txs\  = \ Transpose[{tvalues, \ xvalues}]; 
		
        ListPlot[txs, \ PlotJoined\  -> \ True, \ 
          PlotRange\  -> \ All, \ 
          AxesLabel\  -> \ {"t", \ "x"}]; ]\)

Use the function with the table of approximations from Euler's Method:

txPlot[estimate]

Explanation of txPlot

To understand how txPlot works, you may evaluate each expression in this subsection and examine its result.

Transpose[estimate]

xvalues  =  Map[First,  Last[Transpose[estimate]]]

tvalues  =  First[Transpose[estimate]]

txs  =  Transpose[{tvalues,  xvalues}]

Here is the plot of x vs. t:

ListPlot[txs,  PlotJoined  ->  True,  
      PlotRange  ->  All,  
      AxesLabel  ->  {"t",  "x"}];

The preceding steps were encapsulated into the function txPlot.

Here is a function to plot the calculated,  approximatey-values vs. t.

\(tyPlot[estimate_]\  := \ 
      Module[{yvalues, \ tvalues, \ tys}, \ 
        yvalues\  = \ Flatten[Map[Rest, \ Last[Transpose[estimate]]]]; \ 
        tvalues\  = \ First[Transpose[estimate]]; \ 
        tys\  = \ Transpose[{tvalues, \ yvalues}]; 
		
        ListPlot[tys, \ PlotJoined\  -> \ True, \ 
          PlotRange\  -> \ All, \ 
          AxesLabel\  -> \ {"t", \ "y"}]; ]\)

Use the function:

tyPlot[estimate]

Explanation of tyPlot

To understand how tyPlot works, you may evaluate each expression in the following cells and examine its result:

yvalues  =  Flatten[Map[Rest,  Last[Transpose[estimate]]]]

tys  =  Transpose[{tvalues,  yvalues}]

ListPlot[tys,  PlotJoined  ->  True,  
      PlotRange  ->  All,  
      AxesLabel  ->  {"t",  "y"}];

Here is a function for plotting the piecewise-linear estimate of the solution curve in the (x,  y)-phase plane:

phasePlot[xyestimate_]  :=  
      ListPlot[ Last[Transpose[xyestimate]],  PlotJoined  ->  True, 
         AspectRatio  ->  Automatic,  
        AxesLabel  ->  {"x",  "y"}];

Use the function:

phasePlot[estimate];

Explanation of phasePlot

To understand how phasePlot works, you may evaluate the following two cells and examine their results:

xypairs  =  Last[Transpose[estimate]]

ListPlot[xypairs,  PlotJoined  ->  True,  
      AspectRatio  ->  Automatic,  
      AxesLabel  ->  {"x",  "y"}];

Applying Euler's method—better approximation

Now use a longer interval [0,  2.5] with that same small number of steps:

estimate  =  euler[F,  0,  {1. ,  1. },  2.5,  10];

TableForm[estimate]

txPlot[estimate]

tyPlot[estimate]

phasePlot[estimate];

Applying Euler's method—long time interval

Next, a smaller step size and going all the way to t  =  25:

estimate  =  euler[F,  0,  {1. ,  1. },  25. ,  500];

phasePlot[estimate];

Exercise.  Investigate what happens to solution curves in the phase plane for various initial points.

Comparison with built-in numerical method NDSolve

Euler's Method is a rather crude way to numerically estimate solutions.  Mathematica's built-in function NDSolve uses sophisticated, state-of-the-art methods.   The next two cells use NDSolve to obtain an approximate solution to van der Pol's equation over a short t-interval.

approx  =  
      First[ NDSolve[{x'[t]  ==  y[t],  
            y'[t]  ==  -x[t]  +  (1  -  x[t]^2) y[t],  
            x[0]  ==  1,  y[0]  ==  1},  {x[t],  y[t]}, 
           {t,  0,  1}] ]

ParametricPlot[{x[t]  /.  approx,  y[t]  /.  approx}, 
       {t,  0,  1}];

Compare the rather smooth approximate solution curve just obtained to the jagged curve obtained from Euler's Method.

Exercise.  Use NDSolve for the t-interval [0,  25] and plot the result in the phase plane.  Compare the plot to the one that used Euler's Method.

Copyright © 1999 by Murray Eisenberg