Review: Numerical Methods for Differential Equations

ECE Generic and 2331

Han Q. Le(c)

Segment 1 - basic concept

The app below is used throughout this lecture. It is best to run it outside this notebook (cut and paste into a new notebook) or just download it from the website and run it in parallel with this notebook.

1. Introductory example

1.1 Discussion

Suppose we a problem like this:   ((dy)/(dx))=f[x;y].
Let’s start out with one point: (Subscript[x, 0]) and a known initial or boundary value (Subscript[y, 0]). Then:
                                                               (Subscript[((dy)/(dx))|, (Subscript[x, 0])])=f[(Subscript[x, 0]);(Subscript[y, 0])]  
How do we obtain a numerical estimate for the next point, and eventually for all other points over the range x of  interest? For the next point (Subscript[x, 1]), we can do this: (see Taylor’sseries below)
                                           (Subscript[(Subscript[y, 1])~ (Subscript[y, 0]) +((dy)/(dx))|, (Subscript[x, 0])]) ((Subscript[x, 1])-(Subscript[x, 0]) )=(Subscript[y, 0]) +h  f[(Subscript[x, 0]);(Subscript[y, 0])] ; where  (Subscript[x, 1])=(Subscript[x, 0]) +h
then, this can be repeated for the next point and so on.

This simplistic approach is actually known as the Euler’s method. Let’s consider a simple example:
                                                                ((dy)/(dx))=f(x,y)=((x^2)/(Sqrt[y])) 
Clearly, this is solvable because we can write:    (Sqrt[y])dy=(x^2) dx 
and integrate both side:

∫(Sqrt[y])\[DifferentialD]y == ∫(x^2) \[DifferentialD]x

((2 (y^(3/2)))/3)==((x^3)/3)

In fact, we can just ask Mathematica to do it, so that it includes an arbitrary constant:

DSolve[y'[x]== ((x^2)/(Sqrt[y[x]])),y[x],x]

But we will use this example so that we can compare the various numerical methods to the exact result.

1.2 Numerical illustration

We will choose an initial condition so that we can test the simple approach above. W

y(x)->((((3 (Subscript[c, 1])+(x^3)))^(2/3))/(2^(2/3))) /. x-> 0

We will let y(0)=1, so that (Subscript[c, 1])=2/3. We have the APP below (click on Euler).

Observe for discussion (next section):
- How does the error change vs. step size h?
- What are the causes of error?

Concept review: types of error

"To err is human"

"... and so is computer!"

Error: how does it happen?

Classification of error types

1.	Truncation error: Technically, this name refers to errors that are caused by truncating or stopping a numerical procedure before it has really converged or stopped changing.  Many numerical procedures would have to be run forever in order to truly converge if you could record intermediate values with infinite precision.  However, in practice there is always some sort of stopping criterion, and the difference between the final intermediate value and the true value represents the truncation error associated with the procedure.

2.	Inherent errors: Errors resulting from the inaccuracies of the measurement process itself, or possibly from inaccuracies in the mathematical model that the analyst uses to represent the physical world, are lumped together and called inherent errors.

3. 	Round-off Error: Since all computers have a finite word length and use the binary number system, most decimal numbers cannot be represented with complete accuracy in a computer.  Note that this is true whether the numbers are stored in fixed-point format or floating-point format. Whether the infinite binary string that represents the true value of the decimal number is truncated (i.e., chopped off) or rounded somehow to produce an internal representation, this type of error is simply called round-off error.

4. 	Propagated error: This is error in succeeding steps of a process that is propagated forward from a previous inaccurate result.  This can result in a procedure that never converges or one which “blows up.”

2. Euler's method and error

2.1 Error consideration

The method we use above is actually the basic Euler’s method. We also see that it can have significant error. The larger step size h is, the larger is the error. The issue is how to determine the magnitude of this error and how it varies vs. h. In order to do this, we will review Taylor’s series theorem

2.1.1 Review of Taylor’s series

Taylor’s series theorem states that if a function is continuous and differentiable of all orders in an interval [x,(Subscript[x, 0])], then:
                  f[x]=f[(Subscript[x, 0])]+ f'[(Subscript[x, 0])] (x-(Subscript[x, 0]))+(1/(2!)) f''[(Subscript[x, 0])](((x-(Subscript[x, 0])))^2)+(1/(3!)) f'''[(Subscript[x, 0])](((x-(Subscript[x, 0])))^3) +...+(1/(n!)) (f^((n)))[(Subscript[x, 0])](((x-(Subscript[x, 0])))^n)+...  
or, we can write,    let x-(Subscript[x, 0])≡ h
                  f[x]=f[(Subscript[x, 0])]+ f'[(Subscript[x, 0])] h+(1/(2!)) f''[(Subscript[x, 0])] (h^2)+...+(1/(n!)) (f^((n)))[(Subscript[x, 0])](h^n)+O[(h^(n+1))]

Here we introduce an important notation: O[(h^(n+1))], which is used to indicate that the sum of remaining terms (for n-> infinity), varies on the order of (h^(n+1)).
This is important, because if we stop at the nth term when computing f[x], we can set an upper limit how the remaining error scales as a function of step size. For example, let n+1=5, then if we reduce step size a factor of 2, the error will be reduce by (1/(2^5))=1/32.

2.1.2 Apply to Euler’s method

If we go back to the APP above and click of error check, observe how the error at one point (x=1) varies as a function of step size. What can be seen is that it scales as O[h].

We can see why. here is what is computed:           (Subscript[y, i+1])= (Subscript[y, i]) +h  f[(Subscript[x, i]);(Subscript[y, i])] 
where as                                      (Subscript[(OverHat[y]), i+1])= (Subscript[(OverHat[y]), i]) +h  f[(Subscript[x, i]);(Subscript[(OverHat[y]), i])] + O[(h^2)]
The O[(h^2)] principally includes the term of second order:
                                                    O[(h^2)]=(1/(2!)) y''[(Subscript[x, i])](h^2) +O[(h^3)] 
This error is at (Subscript[x, i+1]), but it will propagate to the next and the next and so on...
Let (Subscript[(OverHat[y]), i]) indicate the exact value of y, to be distinguished from the estimate value (Subscript[y, i]). We can expand:
                                          (Subscript[y, i+1])-(Subscript[(OverHat[y]), i+1])= (Subscript[y, i])-(Subscript[(OverHat[y]), i])  +h  (f[(Subscript[x, i]);(Subscript[y, i])]-f[(Subscript[x, i]);(Subscript[(OverHat[y]), i])])-O[(h^2)]
                                          (Subscript[y, i+1])-(Subscript[(OverHat[y]), i+1])= (Subscript[Δy, i]) +h  ((∂/(∂y))f[(Subscript[x, i]);(Subscript[y, i])](Subscript[Δy, i]) +O'[((Δy)^2)])-O[(h^2)] 
We see that the error (Subscript[Δy, i]) becomes a sum whose length increases proportional with inverse h: the smaller the step size is, the more points there are and the larger the accumulation. Hence, the error doesn’t scale as O[(h^2)], but worse, only as O[h]

2.2 More general consideration of error

From the APP, we also see that: 
- smaller step size: better approximation initially, but 
- too small step size: computer round-off errors (truncation, loss of precision)

The general principle is that the algorithm should be designed such that the order of error, i. e.  O[(h^(n+1))] should be known, which means that all other terms of higher order should be such that they are canceled out. This way, regardless of the error magnitude, we can “cap” the error with a scaling limit. This leads to improved and generalized method such as Runge-Kutta.

3. Improved Euler's method

3.1 Overview of approaches to improve

See the APP on method illustration.

3.2 Heun's method

To improve on the Euler's method, the issue is the error for calculating the derivative at just one side. This leads to an accumulation that can scale as O[h] instead of second-order derivative error of O[(h^2)].
An improvement is to use the derivative as the average of both end points. In other words:
                                (Subscript[y, i+1])  = y((Subscript[x, i]))+h ((f((Subscript[x, i]),(Subscript[y, i])) + f((Subscript[x, i+1]),(Subscript[y, i+1])))/2)
But if this is implemented exactly, it means that (Subscript[y, i+1]) has to be solved, which, while feasible with find root method, will cost extra computing processes. Heun’s method is that there is no need for (Subscript[y, i+1]) in f((Subscript[x, i+1]),(Subscript[y, i+1])) to be self-consistent, but just an approximation:
                                      ((Subscript[y, i+1])^((0))) =y((Subscript[x, i]))+hf((Subscript[x, i]),(Subscript[y, i])) ; and then:
                                    (Subscript[y, i+1])  = y((Subscript[x, i]))+h ((f((Subscript[x, i]),(Subscript[y, i])) + f((Subscript[x, i+1]),((Subscript[y, i+1])^((0)))))/2)     
This algorithm costs only one extra-step (see APP illustration above), hence Heun's method is also known as two-step improved Euler’s method. This improvement can bring the residual errors to the order of O[(h^2)] as shown below.

3.3 Heun’s method error scaling

3.4 Midpoint method (modified Euler)

Similar to Heun's method, we can also use two steps instead of just one step, but this time, instead of averaging the derivatives on both ends, we just guess a value at midpoint, and use the derivative at that point. The first step is to get "mid-point" value
                  (Subscript[y, i+1/2])=y((Subscript[x, i]))+(h/2)f((Subscript[x, i]),(Subscript[y, i])) 
The next step: use this to get "mid-point" estimate of derivative and then the next step:
                  (Subscript[y, i+1])  =y((Subscript[x, i]))+h f((Subscript[x, i+1/2]),(Subscript[y, i+1/2])).
Use the APP above to see how it works.

3.5 Error comparison

Below is an error comparison of all three methods discussed previously.

Note that both Heun and mid-point methods indeed achieve error scaling of the second order O[(h^2)].

4. Runge-Kutta method concept (summary only)

4.1 Basic idea

Use the APP for the discussion.

Theoretically, if we use Taylor’s series theorem, we can, in principle keep on including higher and higher order derivatives in the estimate such that  
                                  (Subscript[y, i+1])-(Subscript[(OverHat[y]), i+1])-> O[(h^n)]
as high n as desired.
This means to estimate 1st order derivative, use it to estimate the second, and so forth, which is a recursive process. The process can be parametrized (as illustrated in the APP for second order) with the constraint and criterion such that the residual error is minimum.
This is a complex optimization problem but feasible, and the essence is only to achieve a “recipe” with optimal parameters. The implementation is then simple in concept: Every order of derivative is a function of the previous, which can be traced all the ways back to f((Subscript[x, i])+p h,f[(Subscript[x, i])+p h,Q[(Subscript[y, i])]]), which represents an estimate at point (Subscript[x, i])+p h where p is a parameter, and Q[(Subscript[y, i])] represents a corresponding value.

Thus, the general expression is:       (Subscript[y, i+1])=(Subscript[y, i])+ϕ((Subscript[x, i]),(Subscript[y, i]),h) h
where ϕ((Subscript[x, i]), (Subscript[y, i]), h) is a function whose form is to be determined for minimal error.
- For Heun's method, note that:
                 ϕ((Subscript[x, i]),(Subscript[y, i]),h)= ((f((Subscript[x, i]),(Subscript[y, i])) + f((Subscript[x, i+1]),((Subscript[y, i+1])^((0)))))/2)=((f((Subscript[x, i]),(Subscript[y, i])) + f[  (Subscript[x, i+1]), (Subscript[y, i]) + h f((Subscript[x, i]),(Subscript[y, i]))   ])/2)
- For mid-point method:
                 ϕ((Subscript[x, i]),(Subscript[y, i]),h)=f((Subscript[x, i+1/2]),(Subscript[y, i+1/2])) =f((Subscript[x, i]) + h/2 ,(Subscript[y, i]) + (h/2) f((Subscript[x, i]),(Subscript[y, i])))

Hence, Runge-Kutta’s method generalizes the concept to create an algorithm (recipe) to obtain ϕ((Subscript[x, i]), (Subscript[y, i]), h) for decreasing error scaling  (Subscript[y, i+1])-(Subscript[(OverHat[y]), i+1])-> O[(h^n)], which is to increase the accuracy as needed. The algorithm is recursive in a manner suitable for computer implementation.

4.2 Example: R-K for second order

4.2.1 Mid-point method

Here is an implementation of R-K method to the second order:
(Subscript[y, i+1]) = (Subscript[y, i]) + ϕ((Subscript[x, i]), (Subscript[y, i]), h) h; 
where 
ϕ((Subscript[x, i]), (Subscript[y, i]), h)  = (Subscript[a, 1]) (Subscript[k, 1])  + (Subscript[a, 2]) (Subscript[k, 2])  ; where
(Subscript[k, 1]) = f((Subscript[x, i]),(Subscript[y, i])) ; and
(Subscript[k, 2]) = f((Subscript[x, i]) + (Subscript[p, 1]) h,(Subscript[y, intermediate]))  =f((Subscript[x, i]) + (Subscript[p, 1]) h,(Subscript[y, i]) + (Subscript[q, 11]) f((Subscript[x, i]),(Subscript[y, i])) h) ;
But what are the values (Subscript[a, 1]) , (Subscript[a, 2]), (Subscript[p, 1]) , (Subscript[q, 11])?
Suppose we are doing mid-point method, what are those?
(Subscript[y, i+1])  =  y((Subscript[x, i])) + h f((Subscript[x, i]) +(1/2) h,(Subscript[y, intermediate])).
(Subscript[y, intermediate])  = (Subscript[y, i]) + (h/2) f((Subscript[x, i]),(Subscript[y, i])) .
We notice that:
(Subscript[a, 1])=0 ;
(Subscript[a, 2]) = 1;
(Subscript[k, 2]) = f((Subscript[x, i]) + (Subscript[p, 1]) h,(Subscript[y, intermediate]))  = f((Subscript[x, i]) +(1/2) h,(Subscript[y, i]) + (h/2) f((Subscript[x, i]),(Subscript[y, i])) )
so: 
(Subscript[p, 1]) = 1/2 ;  (Subscript[q, 11]) = 1/2.
Thus, we say that mid-point method is a special case of 2nd order R-K with:
(Subscript[a, 1])=0 ;
(Subscript[a, 2]) = 1;
(Subscript[p, 1]) = 1/2 ; and  (Subscript[q, 11]) = 1/2.

4.2.2 Heun’s method

Are the above parameters the only choices?

No, in general, they are chosen to minimize the residual errors: this is the fundamental concept in the R-K's method.
For example:
(Subscript[a, 1]) = (Subscript[a, 2]) = 1/2; then:  (Subscript[p, 1]) = 1 and (Subscript[q, 11]) =1:
Here:
(Subscript[y, i+1]) = (Subscript[y, i]) + ((Subscript[a, 1]) f((Subscript[x, i]),(Subscript[y, i]))   + (Subscript[a, 2])  f((Subscript[x, i]) + (Subscript[p, 1]) h,(Subscript[y, i]) + (Subscript[q, 11]) f((Subscript[x, i]),(Subscript[y, i])) h)   )  h;

Let's compare this expression with Heun's method:
(Subscript[y, i+1])  =  y((Subscript[x, i])) + h((f((Subscript[x, i]),(Subscript[y, i])) + f((Subscript[x, i+1]),((Subscript[y, i+1])^((0)))))/2)
         = (Subscript[y, i]) + h((f((Subscript[x, i]),(Subscript[y, i])) + f((Subscript[x, i]) + h ,  (Subscript[y, i]) + h f((Subscript[x, i]),(Subscript[y, i])) ))/2)
we see that it is exactly the same. Heu’s method is one choice for 2nd order R-K.

4.2.3 Another example: Ralston's method

If we choose  (Subscript[a, 2]) = 2/3, then we have so-called Ralston's method.
Here (Subscript[a, 1]) = 1/3, (Subscript[p, 1]) = 3/4 and (Subscript[q, 11]) = 3/4.

5. Recipes for R-K for higher order

5.1 Third order

A third order R-K expression is:
  (Subscript[y, i+1]) =(Subscript[y, i])+(1/6)h ((Subscript[k, 1])+4(Subscript[k, 2])+(Subscript[k, 3]));
where:
(Subscript[k, 1]) =  f((Subscript[x, i]),(Subscript[y, i])) ;
(Subscript[k, 2]) =  f((Subscript[x, i]) + (1/2)h,(Subscript[y, i]) +  (1/2)(Subscript[k, 1]) h) ;
(Subscript[k, 3]) =  f((Subscript[x, i]) + h,(Subscript[y, i]) +(2 (Subscript[k, 2]) -  (Subscript[k, 1]) ) h  ) ;

Or: (Subscript[y, i+1]) =(Subscript[y, i])+(h/6) (({{1, 4, 1}})) . (({{(Subscript[k, 1])}, {(Subscript[k, 2])}, {(Subscript[k, 3])}}))

5.2 Fourth order

A fourth order expression:
 (Subscript[y, i+1])=(Subscript[y, i])+(1/6)h ((Subscript[k, 1])+2(Subscript[k, 2])+2(Subscript[k, 3])+(Subscript[k, 4])); 
where:
(Subscript[k, 1]) =  f((Subscript[x, i]),(Subscript[y, i])) ;
(Subscript[k, 2]) =  f((Subscript[x, i]) + (1/2)h,(Subscript[y, i]) +  (1/2)(Subscript[k, 1]) h) ;
(Subscript[k, 3]) =  f((Subscript[x, i]) +(1/2) h,(Subscript[y, i]) + (1/2)(Subscript[k, 2]) h  ) ;
(Subscript[k, 4]) =  f((Subscript[x, i]) + h,(Subscript[y, i]) + (Subscript[k, 3]) h  ) ;
Or:
(Subscript[y, i+1]) =(Subscript[y, i])+(h/6) (({{1, 2, 2, 1}})) . (({{(Subscript[k, 1])}, {(Subscript[k, 2])}, {(Subscript[k, 3])}, {(Subscript[k, 4])}}))

5.3 Fifth order R-K (or Butcher's method)

(Subscript[y, i+1]) =(Subscript[y, i])+(h/90) (({{7, 32, 12, 32, 7}})) . (({{(Subscript[k, 1])}, {(Subscript[k, 3])}, {(Subscript[k, 4])}, {(Subscript[k, 5])}, {(Subscript[k, 6])}})) 
(note that (Subscript[k, 2]) is not included as it is only an intermediate value for recursion)
where:
(Subscript[k, 1]) =  f((Subscript[x, i]),(Subscript[y, i])) ;
(Subscript[k, 2]) =  f((Subscript[x, i]) + (1/4)h,(Subscript[y, i]) +  (1/4)(Subscript[k, 1]) h) ;
(Subscript[k, 3]) =  f((Subscript[x, i]) +(1/4) h,(Subscript[y, i]) + (1/8)(Subscript[k, 1]) h + (1/8)(Subscript[k, 2]) h  ) ;
(Subscript[k, 4]) =  f((Subscript[x, i]) + (1/2)h,(Subscript[y, i]) - (1/2)(Subscript[k, 2]) h  +  (Subscript[k, 3]) h  ) ;
(Subscript[k, 5]) =  f((Subscript[x, i]) + (3/4)h,(Subscript[y, i]) + (3/16)(Subscript[k, 1]) h  +(9/16) (Subscript[k, 4]) h  ) ;
(Subscript[k, 6]) =  f((Subscript[x, i]) + h,(Subscript[y, i]) - (3/7)(Subscript[k, 1]) h  +(2/7)(Subscript[k, 2]) h +(12/7)(Subscript[k, 3]) h-(12/7) (Subscript[k, 4]) h  + (8/7)  (Subscript[k, 5]) h) ;

5.4 Additional note:

Runge-Kutta methods are useful for numerically solving certain types of ordinary differential equations. Deriving high-order Runge-Kutta methods is no easy task, however. There are several reasons for this. The first difficulty is in finding the so-called order conditions. These are nonlinear equations in the coefficients for the method that must be satisfied to make the error in the method of order O((h^n)) for some integer n where h is the step size. The second difficulty is in solving these equations. Besides being nonlinear, there is generally no unique solution, and many heuristics and simplifying assumptions are usually made. Finally, there is the problem of combinatorial explosion. For a twelfth-order method there are 7813 order conditions! 
For Mathematica, specify “method” option in function NDSolve allows one to choose specific algorithm, and newer version of Mathematica also does automatic selection of R-K order.

6. Error comparison

Use the APP to study error. Verify that the scaling is correct for 3rd, 4th, and 5th order R-K.

Note about error near machine epsilon. Mathematica gives option for much higher precision if needed.