Math 331.1 • 16 September 2005

Euler's Method

Copyright © 2005 Murray Eisenberg

Preparation (technical stuff)

<< Utilities`Notation`

Off[General::spell1]
Symbolize[(Subscript[t, 0])]
Symbolize[(Subscript[y, 0])]
On[
      General::spell1]

<< Graphics`

thicker  =  Thickness[0.0075];
textFormat  =  {FontWeight  ->  "Bold",  
          FontSize  ->  16};

Euler's Method procedure

Documentation

\(\(\(euler::usage\  = \ \*""euler[f, (Subscript[t, 0]), (Subscript[y, 0]), b, n] gives the numerical solution to {y' == f[t, y], y[(Subscript[t, 0])] == (Subscript[y, 0])} over the interval [(Subscript[t, 0]), b] by the n-step Euler's method.  The result is in the form of a table of {t, y} pairs."";\)\(	\)\)\)

Definition:

\(euler[f_, \ \(t\_0\) : _, \ \(y\_0\) : _, \ b_, \ n_]\  := 
	
      Module[{t, \ y, \ h, \ 
          i}, \ 
		h\  = \ \((b\  - \ t\_0)\)\/n; 
		t[0]\  = \ t\_0; \ 
        y[0]\  = \ 
          y\_0; 
		Do[
			t[i]\  = \ 
            t\_0\  + \ i\ h; 
			y[i]\  = \ 
            f[t[i\  - \ 1], \ y[i\  - \ 1]]\ h\  + \ 
              y[i\  - \ 1], 
			{i, \ 1, \ 
            n\ }
		]; 
		Table[{t[i], \ y[i]}, \ {i, \ 0, \ 
            n}]
	]\)

Example

We shall apply the procedure euler defined above in order to solve

	y'  =  
        (t^2)  -  y,            y(0)  =  1

for 0   ≤  t  ≤  2 numerically by Euler's method using 10 steps.  Then we shall compare that solution to both the exact solution obtained from Mathematica's DSolve function and to the numerical solution provided by Mathematica's NDSolve function

The right-hand-side f(t,  y)

Clear[f]
f[t_,  y_]   :=  (t^2)  -  y

Applying Euler's Method

{TraditionalForm[(Subscript[t, i])],  TraditionalForm[(Subscript[y, i])]}

eulerTable  =  euler[f,  0. ,  1. ,  2. ,  10] ;
TableForm[eulerTable,  
      TableHeadings  ->  {None,  
          TraditionalForm  /@  {(Subscript[t, i]),  (Subscript[y, i])}}]

"Connect the dots" to obtain a piecewise linear curve:

DisplayTogether[
        eulerDots  =  
          ListPlot[eulerTable,  PlotStyle  ->  PointSize[0.025],  
            PlotRange  ->  {0,  2},  
            AspectRatio  ->  Automatic,  
            TextStyle  ->  textFormat],  
        eulerCurve  =  
          ListPlot[eulerTable,  PlotJoined  ->  True,  
            PlotStyle  ->  thicker,  
            AspectRatio  ->  Automatic]];

Comparing with exact solution from built-in DSolve

soln  =  
      DSolve[{y'[t]  ==  f[t,  y[t]],  y[0]  ==  1},  y[t],  
        t]

actualCurve  =  
        Plot[y[t]  /.  soln,  {t,  0,  2},  
          PlotStyle  ->  {Red,  thicker},  
          PlotRange  ->  {0,  2},  
          AspectRatio  ->  Automatic,  
          TextStyle  ->  textFormat];

DisplayTogether[eulerDots,  eulerCurve,  actualCurve];

Comparing with built-in NDSolve

ndSoln  =  
      NDSolve[{y'[t]  ==  f[t,  y[t]],  y[0]  ==  1},  
        y,  {t,  0,  2}]

ndCurve  =  
        Plot[y[t]  /.  ndSoln,  {t,  0,  2},  
          PlotStyle  ->   {Magenta,  thicker},  
          PlotRange  ->  {0,  2},  
          AspectRatio  ->  Automatic,  
          TextStyle  ->  textFormat];

DisplayTogether[eulerDots,  eulerCurve,  ndCurve];