Numeric Integration

The distinction between numeric quadrature and numeric integration has to do with the nature of the differential equation that is to be solved. Occasionally, a differential equation is explicitly expressed in terms of the independent variable x, such as the following.

(dy/dx) = f (x), (Subscript[x, 0]) ≤ x ≤ (Subscript[x, n])

In this case the solution is given by the following definite integral.

\(∫\_\(x\_0\)\%\(x\_\(\(\ \)\(n\)\)\)f \((x)\) dx\)

Numeric quadrature methods such as the Trapezoid Method or Simpson's Method yield a solution.  Otherwise, the differential equation is not solved explicitly in terms of the independent variable but is defined implicitly in terms of both the independent and dependent variables, such as the following.

(dy/dx) = f (x,  y), (Subscript[x, 0]) ≤ x ≤ (Subscript[x, n])

In this case,we must use numeric integration to solve the differential equation.  Finding a function F(x) that satisfies the equation (dy/dx) = f (x, y) would complete the problem.  However, this is typically not possible. The goal then becomes to find a sequence of values (Subscript[y, 0]),(Subscript[y, 1]),..., (Subscript[y, n]) with each (Subscript[y, i]) ≐ F((Subscript[x, i])).

Euler's Method

With Euler' s Method, we make a forward difference approximation to the derivative at the points(Subscript[x, 0]),(Subscript[x, 1]),...,(Subscript[x, n - 1]).  Although it is not necessary that the sequence  of points{(Subscript[x, i])}_(i  =  0)^nbe equidistant from each other, it will make our  calculations easier if we assume that this is the case . So at each point(Subscript[x, i]),for 0 ≤ i<n,let Δ x=(Subscript[x, i + 1])- (Subscript[x, i]).In this instance, \(the\ forward\ difference\ looks\ like\ the\ \(\(following\)\(.\)\)\)

(((Subscript[y, i + 1]) -  (Subscript[y, i]))/((Subscript[x, i + 1]) -  (Subscript[x, i])))= (((Subscript[y, i + 1]) -  (Subscript[y, i]))/(Δ x))≐ f ((Subscript[x, i]),  (Subscript[y, i]))

This implies the following.

(Subscript[y, i + 1]) ≐ (Subscript[y, i])+ f( (Subscript[x, i]),  (Subscript[y, i])) Δ x

Figure 4.1 shows a single step of the Euler Method applied to the point ((Subscript[x, n]),(Subscript[y, n])) in an attempt to approximate ((Subscript[x, n + 1]),F((Subscript[x, n + 1])) ) of the function y = F(x).The red line that is tangent to the curve determines the approximate point ((Subscript[x, n + 1]),(Subscript[y, n + 1])).

Therefore,our iteration steps will be based upon the following equations.

(Subscript[y, 1]) ≐ (Subscript[y, 0]) + f((Subscript[x, 0]),  (Subscript[y, 0])) Δ x
(Subscript[y, 2]) ≐ (Subscript[y, 1]) + f((Subscript[x, 1]),  (Subscript[y, 1])) Δ x
...
(Subscript[y, n]) ≐ (Subscript[y, n - 1]) + f((Subscript[x, n - 1]),  (Subscript[y, n - 1])) Δ x

Figure 4.2 shows several iterations of the Euler Method to approximate the function y = F(x) that satisfies the equation y' = f ( x, y).  Each of the blue points ( (Subscript[x, k]), (Subscript[y, k])) is an approximation to the corresponding blue point ( (Subscript[x, k]), F( (Subscript[x, k]))).

For this method,the discretization error  is bounded above by a constant times Δ x which is referred to as O(Δ x).   The (global) discretization error in using an iterative numeric method is defined to be (max_(0 ≤  i  ≤  n))|(Subscript[y, approx]) -  (Subscript[y, actual])|.

Examples

The following Mathematica excerpt shows n = 5 iterations of the Euler Method starting at the point ( (Subscript[x, 0]),(Subscript[y, 0])) = (0, 1) until (Subscript[x, n]) = 0.1 with Δx = 0.02.

Clear[f, x, 
        h];
f[x_, y_] := x + y;
For[
x = 0; h =  .02; 
      y = 1, 
x <  .1, 
x += 
        h, 
y = y + f[x, y]*h; 
Print[
        y] ]

Function euler[] accepts a function of independent variable x, a domain over which to integrate, (Subscript[x, 0])≤ x ≤ (Subscript[x, n]), two y values over which to sketch the graph (Subscript[y, min]) ≤ y ≤ (Subscript[y, max]) and a number of subintervals, n.  Tangent lines are graphed on the curve to see the value of the derivative that is being used.

Clear[x];euler[f_, {x0_, xn_}, {ymin_, ymax_}, n_]:=
Module[{deltax = (xn - x0)/n, xxvals,  yyvals, l, ll, p1, 
              p2, p3, p4, p5, lin, disc},
  
xxvals  =  Table[N[x0  +  i*deltax], {i, 0, n}];
yyvals   =  
                f /. x -> xxvals;
  
lin  =  
                Table[Line[{{xxvals[[i]], 
                        yyvals[[i]]}, {xxvals[[i + 1]],  
                        yyvals[[i]] + (D[f, x] /. 
                                x -> xxvals[[i]] )*deltax}}], {i, 
                    1, n}];
  
ll = {{xxvals[[1]], 
                    yyvals[[1]]}};
Do[
                ll = Append[
                    ll, {xxvals[[i]], 
                      ll[[i - 1, 
                            2]] + (D[f, x] /. 
                              x -> xxvals[[i - 1]])*deltax}], {i, 
                  2, n + 1}];
l  = Transpose[{xxvals, yyvals}];
disc = 
                Max[Abs[Transpose[ll][[2]] - yyvals]];
  
p1  =  
                Plot[f, {x, x0, xn}, PlotRange -> {ymin, ymax}, 
                  DisplayFunction -> Identity, 
                  PlotStyle -> RGBColor[1, 0, 0], 
                  PlotLabel -> "Euler (N, Max Error) = " 
                      StandardForm[{n, N[disc]}]];
 p2  =  
                ListPlot[l, PlotRange -> {ymin, ymax}, 
                  PlotStyle -> {PointSize[0.02], RGBColor[0, 0, 1]}, 
                  DisplayFunction -> Identity];
p3  = 
                Graphics[lin, DisplayFunction -> Identity ];
p4 = 
                ListPlot[ll, PlotRange -> {ymin, ymax}, 
                  PlotStyle -> {PointSize[0.02], 
                      RGBColor[0, 0, 
                        0]}, 
DisplayFunction -> 
                    Identity];
p5 = 
                ListPlot[ll, PlotRange -> {ymin, ymax}, 
                  PlotStyle -> {PointSize[0.02], 
                      RGBColor[ .25, 1,  .25]}, 
                  PlotJoined -> 
                    True, 
DisplayFunction -> 
                    Identity];
Show[p1, p2, p3, p4, p5, 
                DisplayFunction -> $DisplayFunction];]

Evaluate function euler[] over 0  ≤  x ≤ 2 π with -1 ≤ Cos[x] ≤ 2 using n = 8 panels.  Changing the value of n will show how small numbers of iterations do a poor job of approximating the function.  Larger n values will demonstrate a reduction in  the discretization error.

euler[Cos[x] , {0, 2  Pi}, {-1, 2}, 4];

Evaluation of the following Table[] command will generate a series of graphs that when animated,will demonstrate how the error and approximation improve as n grows.See if you can discover why for integer values of k,the Euler approximation to Cos[k x] or Sin[k x] will always agree with the function after one period.

Table[
        euler[Cos[x], {0, 2  Pi}, {-1, 2}, i], {i, 2, 10}];

Modified Euler's Method

Although Euler's method is easy to understand and implement, it is not very accurate.  One method that improves upon its accuracy is known as the Modified Euler Method. It is based upon a weighted averaging of two adjacent derivatives, f( (Subscript[x, n]), (Subscript[y, n])) and f( (Subscript[x, n + 1]) ,(Subscript[y, n + 1])).  The challenge with this is that f( (Subscript[x, n + 1]) ,(Subscript[y, n + 1])) is dependent upon (Subscript[y, n + 1]) which by Euler is dependent upon f( (Subscript[x, n]), (Subscript[y, n])).  So the modification will use a predictor-correct mechanism that predicts the next value of y, call it (Subscript[Y, n + 1]) and then uses this to compute f ( x, Y ) so that the weighted average can then correct the value to (Subscript[y, n + 1]).  Our prediction strategy will be a regular Euler Method step at (Subscript[x, i]).

(Subscript[Y, i + 1]) = (Subscript[y, i])+ f( (Subscript[x, i]),  (Subscript[y, i])) Δ x

For correction, we will use the point ( (Subscript[x, i]) + Δx / 2, ((Subscript[y, i]) + (Subscript[Y, i + 1])) / 2 ) to then perform the corrector step.  Geometrically, this point is the center of the rectangle formed by the points ( (Subscript[x, i]), (Subscript[y, i]) ), ( (Subscript[x, i + 1]), (Subscript[y, i])), ( (Subscript[x, i + 1]),(Subscript[Y, i + 1])) and ( (Subscript[x, i]), (Subscript[Y, i + 1])).  The correction will involve a the repetition of a regular Euler Method step using ( (Subscript[x, i]) + Δx / 2, ( (Subscript[y, i]) + (Subscript[Y, i + 1])) / 2) to evaluate the derivative.  The equation that is
used is the following.

(Subscript[y, i + 1]) ≐ (Subscript[y, i]) + f ( (Subscript[x, i]), Δx / 2, (Subscript[y, i]) + f ( (Subscript[x, i]), (Subscript[y, i])) Δx / 2 ) Δx

In this representation,what we have been calling (Subscript[Y, i + 1]) is equal to (Subscript[y, i]) + f((Subscript[x, i]), (Subscript[y, i])) Δx.

Figure 4.3 shows the geometry involved with the Modified Euler Method. The point ((Subscript[x, n + 1]), (Subscript[y, Euler])) has been predicted by
a single Euler iteration from the point ((Subscript[x, n]), (Subscript[y, n])).  The midpoint is used to generate the tangent point that is shown in red which is used to generate the derivative that will correct the prediction.  This is displayed by transposing the blue line to the parallel green line through ((Subscript[x, n]), (Subscript[y, n])) that determines (Subscript[y, Mod - Eul]).

Thus, an iteration of the Modified Euler Method would give the following.

(Subscript[y, 1]) ≐ (Subscript[y, 0]) + f ( (Subscript[x, 0]) + Δ x /2, (Subscript[y, 0]) + f( (Subscript[x, 0]),  (Subscript[y, 0])) Δx / 2) Δ x
(Subscript[y, 2]) ≐ (Subscript[y, 1]) + f ( (Subscript[x, 1]) + Δ x /2, (Subscript[y, 1]) + f( (Subscript[x, 1]),  (Subscript[y, 1])) Δx / 2) Δ x
...
(Subscript[y, n]) ≐ (Subscript[y, n - 1]) + f ( (Subscript[x, n - 1]) + Δ x /2, (Subscript[y, n - 1]) + f( (Subscript[x, n - 1]),  (Subscript[y, n - 1])) Δx / 2) Δ x

Figure 4.4 shows the red tangent line that was used to predict ((Subscript[x, n + 1]),(Subscript[y, Euler])) and the green derivative that determines the
corrected value ((Subscript[x, n + 1]), (Subscript[y, Mod - Eul])).  In this particular case the corrected value ((Subscript[x, n + 1]), (Subscript[y, Mod - Eul])) is virtually indistinguishable from the actual value ((Subscript[x, n + 1]), F((Subscript[x, n + 1]))).

The error term in using the Modified Euler's Method is bounded above by a constant times ((Δx)^( 2)) which is referred to as O(((Δx)^( 2))).  Note that this improved accuracy is at a cost of evaluating the function f ( x, y) twice per iteration.

Runge-Kutta 4

The Runge-Kutta Methods use a weighted average of derivatives similar to the Modified Euler Method.  Both of these methods can be justified by means of algebraic manipulations of Taylor series that we are not at liberty to discuss at this point.  However, with the RK4 Method, there are 3 levels of correction after an initial Euler Method prediction.  In general,we have the
following.

(Subscript[y, i + 1]) ≐  (Subscript[y, i]) + ( (Subscript[r, 1]) + 2 (Subscript[r, 2]) + 2 (Subscript[r, 3]) + (Subscript[r, 4])) / 6   where
(Subscript[r, 1]) = f  ( (Subscript[x, i]), (Subscript[y, i])) Δx
(Subscript[r, 2]) = f ( (Subscript[x, i]) + Δx / 2, (Subscript[y, i]) + (Subscript[r, 1]) / 2) Δx
(Subscript[r, 3]) = f ( (Subscript[x, i]) + Δx / 2, (Subscript[y, i]) + (Subscript[r, 2]) / 2) Δx
(Subscript[r, 4]) = f ( (Subscript[x, i]) + Δx, (Subscript[y, i]) + (Subscript[r, 3])) Δx

Note that the calculation of (Subscript[r, 1]) and (Subscript[r, 2]) is nothing more than the Modified Euler Method described above.  The addition of the terms (Subscript[r, 3]) and (Subscript[r, 4]) will provide additional accuracy.  Thus, when we iterate the RK4 Method, we find ourselves with the following sets of calculations.

(Subscript[y, 1]) ≐  (Subscript[y, 0]) + ( (Subscript[r, 1]) + 2 (Subscript[r, 2]) + 2 (Subscript[r, 3]) + (Subscript[r, 4])) / 6   where
(Subscript[r, 1]) = f ((Subscript[x, 0]), (Subscript[y, 0])) Δx
(Subscript[r, 2]) = f ((Subscript[x, 0]) + Δx / 2, (Subscript[y, 0]) + (Subscript[r, 1]) / 2) Δx
(Subscript[r, 3]) = f ((Subscript[x, 0]) + Δx / 2, (Subscript[y, 0]) + (Subscript[r, 2]) / 2) Δx
(Subscript[r, 4]) = f ((Subscript[x, 0]) + Δx, (Subscript[y, 0]) + (Subscript[r, 3])) Δx

(Subscript[y, 2]) ≐  (Subscript[y, 1]) + ( (Subscript[r, 1]) + 2 (Subscript[r, 2]) + 2 (Subscript[r, 3]) + (Subscript[r, 4])) / 6   where
(Subscript[r, 1]) = f ((Subscript[x, 1]), (Subscript[y, 1])) Δx
(Subscript[r, 2]) = f ((Subscript[x, 1]) + Δx / 2, (Subscript[y, 1]) + (Subscript[r, 1]) / 2) Δx
(Subscript[r, 3]) = f ((Subscript[x, 1]) + Δx / 2, (Subscript[y, 1]) + (Subscript[r, 2]) / 2) Δx
(Subscript[r, 4]) = f ((Subscript[x, 1]) + Δx, (Subscript[y, 1]) + (Subscript[r, 3])) Δx
...
(Subscript[y, n]) ≐  (Subscript[y, n - 1]) + ( (Subscript[r, 1]) + 2 (Subscript[r, 2]) + 2 (Subscript[r, 3]) + (Subscript[r, 4])) /6   where
(Subscript[r, 1]) = f ((Subscript[x, n - 1]), (Subscript[y, n - 1])) Δx
(Subscript[r, 2]) = f ((Subscript[x, n - 1]) + Δx / 2, (Subscript[y, n - 1]) + (Subscript[r, 1]) / 2) Δx
(Subscript[r, 3]) = f ((Subscript[x, n - 1]) + Δx / 2, (Subscript[y, n - 1]) + (Subscript[r, 2]) / 2) Δx 
(Subscript[r, 4]) = f ((Subscript[x, n - 1]) + Δx, (Subscript[y, n - 1]) + (Subscript[r, 3])) Δx

The error term in using the RK4 Method is bounded above by a constant times ((Δx)^( 4)) which is referred to as O(((Δx)^( 4))).  Note that this improved accuracy is at a cost of evaluating the function f ( x, y) four times per iteration.

Examples

The following Mathematica excerpt shows n = 5 iterations of the RK4 Method starting at the point ((Subscript[x, 0]),(Subscript[y, 0])) = (0, 1) until (Subscript[x, n]) = 0.1 with Δx = 0.02.

Clear[f, x, y, h,  r1, 
        r2, r3, r4];
f[x_,  y_]  :=  x  +  y;
For[
x = 0; h =  .02; 
      y = 1, 
x <  .1, 
x +=  
        h, 
r1 = h*f[x, y];  
r2 = 
        h*f[x + h/2, y + r1/2];  
r3 = 
        h*f[x + h/2, y + r2/2];  
r4 = 
        h*f[x + h, y + r3];  
y = 
        y + (r1 + 2*r2 + 2*r3 + r4)/6; 
Print[
        y]]

Function rk4[] accepts a function of independent variable x, a domain over which to integrate, (Subscript[x, 0])≤ x ≤ (Subscript[x, n]), two y values over which to sketch the graph (Subscript[y, min]) ≤ y ≤ (Subscript[y, max]) and a number of subintervals, n.

Clear[x];rk4[f_, {x0_, xn_}, {ymin_, ymax_}, n_]:=
Module[{deltax = (xn - x0)/n, xxvals,  yyvals, l, ll, p1, 
              p2, p4, p5, lin, disc},
  
xxvals  =  Table[N[x0  +  i*deltax], {i, 0, n}];
yyvals   =  
                f /. x -> xxvals;
  
ll = {{xxvals[[1]], 
                    yyvals[[1]]}};
Do[
                r1  =  (D[f, x] /. x -> xxvals[[i - 1]])*
                    deltax; 
r3  =  r2  =  (D[f, 
                            x] /. x -> (xxvals[[i - 1]] + 
                                deltax/2))*
                      deltax; 
r4  =  (D[f, x] /. 
                        x -> (xxvals[[i - 1]] + deltax))*
                    deltax; 
ll  =  
                  Append[ll, 
{xxvals[[i]], 
                      ll[[i - 1, 
                            2]] + (r1 + 2 r2  +  2 r3  +  r4)/
                          6}], {i, 2, n + 1}];
disc = 
                Max[Abs[Transpose[ll][[2]] - yyvals]];
l  = Transpose[{xxvals, yyvals}];
  
p1  =  
                Plot[f, {x, x0, xn}, PlotRange -> {ymin, ymax}, 
                  DisplayFunction -> Identity, 
                  PlotStyle -> RGBColor[1, 0, 0], 
                  PlotLabel -> "RK4 (N, Max Error) = " 
                      StandardForm[{n, N[disc]}]];
 p2  =  
                ListPlot[l, PlotRange -> {ymin, ymax}, 
                  PlotStyle -> {PointSize[0.02], RGBColor[0, 0, 1]}, 
                  DisplayFunction -> Identity];
p4 = 
                ListPlot[ll, PlotRange -> {ymin, ymax}, 
                  PlotStyle -> {PointSize[0.02], RGBColor[0, 0, 0]}, 
                  DisplayFunction -> Identity];
p5 = 
                ListPlot[ll, PlotRange -> {ymin, ymax}, 
                  PlotStyle -> {PointSize[0.02], 
                      RGBColor[ .25, 1,  .25]}, PlotJoined -> True, 
                  DisplayFunction -> Identity];
Show[p1, p2, p4, p5, 
                DisplayFunction -> $DisplayFunction];]

Evaluate function rk4[] over 0  ≤  x ≤ 2 π with -1 ≤ Cos[x] ≤ 2 using n = 8 panels.

rk4[Cos[x] , {0, 2  Pi}, {-1, 2}, 
        8];

Evaluation of the following Table[] command will generate a series of graphs that when animated, will demonstrate how the error and approximation improve as n grows.

Table[
        rk4[Cos[x], {0, 2  Pi}, {-1, 2}, i], {i, 2, 12}];