Errors in Numerical Procedures

Since numercial methods for solving initial value initial value problems only yield approximate solutions, it is very important in 
practice to understand how well a numerical solution 
approximates the actual solution of a differential equation.
That is, one needs useful upper bounds for `truncation error.'

In \(In\ Euler' s\ method\ 
          \((and\ its\ many\ refinements, \ 
            such\ as\ the\ 
Runge - Kutta\ method)\), \ 
        there\ are\ two\ main\ observations\ to\)\(make . 
	
\((1)\)\ The\ quality\ of\ the\ approximate\ solution\ 
          depends\ on\ the\ step\)size h,  
            a smaller step size usually yielding a better approximation .   Therefore,  
        h will figure prominently as aparameter in any reasonable upper bound for 
              error .   In Euler' s method,  
        one generally expects a (global) error that is no\(larger\ than\ Ch\ where\ C\ is\ a\ 
              constant . 
		
\((2)\)\ Local\ truncation\ error\ 
              \((that\ is, \ 
                the\ error\ incurred\ in
the\ first\ step\ of\ a\ method)\)\ 
              is\ easier\ to\ understand\ than\) \(global\ truncation\ error\ 
          \((the\ error\ that\ accumulates\ in
all\ the\ steps, \ 
            from\ first\ to\ last, \ 
            of\ the\ approximation . 
To\ appreciate\ the\ difference, \ 
            note\ that\ local\ truncation
error\ in\ Euler' s\ method\ is\ 
              bounded\ by\ Ch\^2\ where\ C\ is\ roughly
half\ the\ second\ 
              derivative\ of\ the\ solution\ function; \ 
            global
truncation\ error\ is\ also\ bounded\ by\ Ch, \ 
            but\ a\ useful\ formula
for\ C\ is\ much\ harder\ to\ 
              \(find . \)\)\)

Local Truncation Error

In what follows we will consider the initial value problem 
                              y'(x) = y(x),     y(0) = 1,
which has the solution y =(e^x)

\(f[x_, y_]\  := \ y
\)\(EulerMethod[\ f_, {x0_, \ y0_}, \ h_, \ n_\ ]\  := \ 

		
      Module[
		\ \ \ {EMStep}, 
		\ \ \ 
        Interpolation[
\ \ \ 	\ \ \ \ \ \ \ \ \ 
          EMStep[{x_, y_}]\  := \ N[{x + h, y + h*f[x, y]}]; 

          		\ \ \ \ \ \ \ \ \ \ NestList[EMStep, {x0, y0}, n], 

          			\ \ \ \ \ \ \ InterpolationOrder\  -> \ 1
		\ \ \ ]

        	\ \ ]\ \ \ \ \ \ \ \ \ \ \ \ \ \)

h=1:

ya =  EulerMethod[f, {0, 1}, 1, 1]err  =  Abs[ya[1]  -  Exp[1]]

which is significantly larger than 1/2.  The problem is that we
really need to take account of the second derivative of y for all
values of x between 0 and 0+h. The maximum value of y''(x) 
on the interval [0,1] is  e, and e/2 is about 1.36, so we're 
actually doing OK.  Let's try h=.1:

ya =  EulerMethod[f, {0, 1},  .1, 1]err  =  Abs[ya[ .1]  -  Exp[ .1]]

Better.  Let's finish with something *really* small.

ya  =  EulerMethod[f, {0, 1},  .000001, 1]err  =  Abs[ya[ .000001]  -  Exp[ .000001]]

Global Truncation Error

Now we will retry the same values of h but always compare the 
approximate and actual solutions at x = 1.  For h = 1, this
is already done.  For h = .1, we obtain

ya =  EulerMethod[f, {0, 1},  .1, 11]err  =  Abs[ya[1]  -  Exp[1]]

The global error seems roughly the same size as the step size now.
With h=.001, we get

ya =  EulerMethod[f, {0, 1},  .001, 1000]err  =  Abs[ya[1]  -  Exp[1]]

which is again about the same size as h.  Let's try h=.00001:

ya =  EulerMethod[f, {0, 1},  .00001, 100000]err  =  Abs[ya[1]  -  Exp[1]]

Not bad,  
    but I wouldn' t want to wait around for many more decimal\(TextForm
    \`places\ of\ \(accuracy!\)\ \ Euler' s\ method\ turns\ out\ to\ be\ 
      rather\)expensive computationally .   Since global error is 
    proportional to step size,  
    it would take about (10^20)\(TextForm
    \`operations\ to\ grind\ out\ 20\ digits\ of\ \(accuracy!\)\ \ More\ 
      refined\)methods,  
    such as the Runge - Kutta method are much more practical for computational purposes .

Unstable Equilibria and Global Truncation Error

Mathematica provides a routine --called `NDSolve'-- 
for solving ODE's numerically,and it is much better than the 
home-made Euler's method routine we've been using.  It uses the
same syntax as DSolve, except that it requires a range for the 
independent variable.  That is, to solve the above ODE, we type

ya[x_]  =  
      First[y[x] /. 
          NDSolve[ {y'[x]  ==  y[x],  y[0]  ==  1}, 
            y[x], {x, 0, 1}]]

ya[1]

which is a pretty good approximation for e.  Of course, NDSolve
has options such as AccuracyGoal that you can set to improve the number of accurate digits in the approximation.

Good as NDSolve is, it isn't perfect.  Consider the following
example (taken from the Mathematica supplement).
		y'(x) = y(x) - 3(e^x),    y(0) = 1
This is linear and can be solved exactly:

\(actual[x_]\  = \ 
      First[y[x] /. 
          DSolve[
			\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {
              \(y'\)[x]\  == y[x]\  - \ 3  Exp[\(-2\) x], 

              				\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
              y[0]\  == \ 1}, \ \ 

            			\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ y[x], x]]
        \)

If we solve numerically and compare, we get

\(\(
approx[x_]\  = \ 
      First[y[x] /. 
          NDSolve[{\(y'\)[x]\  == y[x]\  - \ 3  Exp[\(-2\) x], 

              				\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ y[0]\  == \ 1}, 
            
				\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
            y[x], {x, 0, 11}]]\)\)

Plot[approx[x] - actual[x], {x, 0, 11}, PlotRange -> {-1, 1}]

Something definitely goes wrong when x reaches 8. Asking for
more accuracy does little to help.

\(approx[x_]\  = \ 
      First[y[x] /. 
          NDSolve[{\(y'\)[x]\  == y[x]\  - \ 3  Exp[\(-2\) x], 

              				\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ y[0]\  == \ 1}, 
            
				\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
            y[x], {x, 0, 15}, 

            				\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
            AccuracyGoal -> 8, 

            \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
            PrecisionGoal\  -> \ 8, \ 

            				\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
            WorkingPrecision -> \ 18]]\)

Plot[approx[x] - actual[x], {x, 0, 15}, PlotRange -> {-1, 1}]

We've effectively asked for a hundred times more precision and
gotten an approximation that only holds out for another two
units along the x axis.  The problem is that the exact solution
of our differential equation is a sort of unstable equilibrium

Needs["Graphics`PlotField`"]

\(\(PlotVectorField[{1, y - 3  Exp[\(-2\) x]}, {x, 0, 10}, {y, \(- .5\), 
        1.5}, \ Axes -> \ True, 
	ScaleFunction -> \((1&)\), 
      Ticks -> None, Frame -> True, AspectRatio -> 1, \ 
      PlotPoints\  -> \ 10]; \)\)

As soon as the approximate solution is even a tiny bit off, the
direction field will rapidly carry it further away from the
actual solution.