Chapter 5: Numerical Differentiation and Integration

5.1 Numerical differentiation: forward, backward and centered schemes

It was shown in class how approximations to a derivative can be obtained by expanding a function using Taylor series. We saw three schemes for calculating numerically the derivative of a continuous functions. These are called the forward, backward and centered difference schemes.

((df((Subscript[x, 0])))/dx)≈((f((Subscript[x, 0])+Δx)-f((Subscript[x, 0])))/(Δx))   	(Forward scheme)

((df((Subscript[x, 0])))/dx)≈((f((Subscript[x, 0]))-f((Subscript[x, 0])-Δx))/(Δx))   	(Backward scheme)

((df((Subscript[x, 0])))/dx)≈((f((Subscript[x, 0])+Δx)-f((Subscript[x, 0])-Δx))/(2Δx))	(Centered scheme)

Here's how these schemes can be coded

Remove["Global`*"]

NDf[f_,x0_,Δx_]:=((f[x0+Δx]-f[x0])/(Δx))   (* Forward difference scheme *)

NDb[f_,x0_,Δx_]:=((f[x0]-f[x0-Δx])/(Δx)) (* Backward difference scheme *)

NDc[f_,x0_,Δx_]:=((f[x0+Δx]-f[x0-Δx])/(2Δx)) (* Centered difference scheme *)

Let's define a function that can be derived analytically, for example sin(x),

f[x_]:=Sin[x]

Plot[f[x],{x,0,2π},AxesLabel->{"x","f(x)"}]

and let's compare the derivative at π/4 of the 3 schemes, using say Δx=0.1, with the exact solution,

x0=π/4;Δx=0.1;
{NDf[f,x0,Δx],NDb[f,x0,Δx],NDc[f,x0,Δx],f'[x0]//N}

Now let's try with Δx=0.01

x0=π/4;Δx=0.01;
{NDf[f,x0,Δx],NDb[f,x0,Δx],NDc[f,x0,Δx],f'[x0]//N}

Would be more illuminating to look at the absolute relative error. Let's do this

x0=π/4;Δx=0.01;
{Abs[(NDf[f,x0,Δx]-f'[x0])/f'[x0]],
Abs[(NDb[f,x0,Δx]-f'[x0])/f'[x0]],
Abs[(NDc[f,x0,Δx]-f'[x0])/f'[x0]]}

We see that the forward and backward schemes have a similar error but that the centered scheme is much more accurate.

The second order derivative can be constructed from the forward and backward scheme as seen in class and is given by

(((d^2)f((Subscript[x, 0])))/(dx^2))≈((f((Subscript[x, 0])+Δx)-2f((Subscript[x, 0]))+f((Subscript[x, 0])-Δx))/((Δx)^2))   	(Centered scheme for 2nd order derivative)

Which can be coded like this:

NDc2[f_,x0_,Δx_]:=((f[x0+Δx]-2f[x0]+f[x0-Δx])/((Δx)^2))

And compared with the exact solution

x0=π/4;Δx=0.01;
NDc2[f,x0,Δx]
f''[x0]//N

5.2 Numerical differentiation of discrete functions

Remove["Global`*"]

Numerical differientiation is largely used for estimating derivatives from experimental data. Consider the following set of measurements of a ball thrown in the air. The list of pairs give the time-position of a ball.

BallData={{0,1.302},{0.03333,1.411},{0.06667,1.5},{0.1,1.578},{0.1333,1.6456},{0.1667,1.703},{0.2,1.745},{0.2333,1.781},{0.2667,1.807},{0.3,1.828},{0.3333,1.818},{0.3667,1.818},{0.4,1.807},{0.4333,1.776},{0.4667,1.734},{0.5,1.682},{0.5333,1.63},{0.567,1.552},{0.6,1.469},{0.6333,1.37},{0.667,1.266},{0.7,1.151},{0.733,1.026},{0.7667,0.875},{0.8,0.719},{0.8333,0.5573},{0.867,0.3854},{0.9,0.193},{0.9333,0.0052}};

ListPlot[BallData,AxesLabel->{"t (s)","y (m)"}]

What is the acceleration a(t)? The acceleration is 
	a=(((d^2)y)/(dt^2)). 
Numerically this can be coded as

a[data_,i_]:=(((data[[i+1,2]]-2data[[i,2]]+data[[i-1,2]]))/(((data[[i+1,1]]-data[[i,1]]))^2))

We can make a list of pair to plot the acceleration as a function of time. Note that we do not attempt to calculate the acceleration at i = 1 and i = 29 because this would require data outside the list.

accel=Table[{BallData[[i,1]],a[BallData,i]},{i,2,28}];

We can plot the results

ListPlot[accel,Joined->True,AxesLabel->{"t (s)","a (m (s^(-2)))"}]

We can get the mean acceleration over the experiment with the command Mean, i.e.

Mean[accel]

5.3 Numerical integration: Rectangle rule

The simplest algorithm to estimate numerically a definite integral of this sort

(∫_a^b)f(x)dx		(1)

is to divide the area between f(a) and f(b) into n rectangles of constant width

Δx=(((b-a))/n)

and heights

f((Subscript[x, 1])), f((Subscript[x, 2])), f((Subscript[x, 3])),⋯,f((Subscript[x, n])),

and to sum the areas of these rectangles. 
This leads to the following approximation called the rectangle rule

(∫_a^b)f(x)dx≈[f((Subscript[x, 1]))+f((Subscript[x, 2]))+f((Subscript[x, 3]))+⋯+f((Subscript[x, n]))] Δx=(∑_(i=1)^n)f((Subscript[x, i]))Δx

where

(Subscript[x, i])=a+(i-1)Δx

Let's create our own function that will achieve this type of numerical integration. I will call this function NIRect for Numerical Integration using the Rectangle rule. I want to create this function with the following four arguments
	f: the function that we want to integrate numerically
	a: the lower limit of integration
	b: the upper limit of integration
	n: the number of elements (i.e. rectangles) to use for the integration

Your function may look like this:

Remove["Global`*"]

NIRect[f_,a_,b_,n_]:=(((b-a)/n))(∑_(i=1)^n)f[a+(i-1)(((b-a)/n))]

Here's how you can use this function. First define the function that you need to integrate, say sin(x) from 0 to 2π, and then call NIRect with the appropriate arguments. Let's use only 10 rectamgles, i.e.,

f[x_]:=Sin[x]
a=0;b=π/2;n=10;
NIRect[f,a,b,n]//N

The integral

(∫_0^(π/2))sin(x)dx

can in fact be integrated anaytically and we know that the exact answer should be

(∫_0^(π/2))Sin[x]\[DifferentialD]x

We see that using the rectangle rule with 10 rectangles returns 0.919403, which is an error of roughly 8% compared to the exact solution which is exactly 1. Let's try with 100 rectangles

f[x_]:=Sin[x]
a=0;b=π/2;n=100;
NIRect[f,a,b,n]//N

We see that increasing the number of rectangles produces a more accurate results. Yet, it is not 1. In fact, it will never be 1 unless you use an infinite number of rectangles. So, when doing numerical calculations like this there is always the issue of determining when to stop the integration. For example, if we use 10 billion rectangles the computer may take days to perform the integration. There is always a trade off between limits of computational power and accuracy needed for your application.

Note that we don't need to use the variable f for the input function. Here's another example with another function that cannot be integrated analytically

g[x_]:=(E^(-(Sqrt[(x^2)+x])))
a=0;b=2;n=10;
NIRect[g,a,b,n]//N

Notice how I called my function NIRect this time with the variable name g even though the function is originally defined with the variable f. This illustrates the power and flexibility of defining functions. Same can be said with the other arguments. For example a and b can really be anything else. For example

MyFunction[x_]:=3(x^2)+x
LowerLimit=2.3;UpperLimit=10.2;NumberOfRectangle=100;
NIRect[MyFunction,LowerLimit,UpperLimit,NumberOfRectangle]//N

Sometime you may want to be more explicit and clearer in the way you code functions. For example, it is not immediately clear what (((b-a)/n)) represents in the definition of NIRect. You could be clearer using the following syntax:

NIRect[f_,a_,b_,n_]:=(Δx=(((b-a)/n));Δx(∑_(i=1)^n)f[a+(i-1)Δx])

The first term in the parenteses defines Δx and the second part is the actual integration using the definition of Δx. You can check that this will return the same answer as in the previous method. For example,

g[x_]:=(E^(-(Sqrt[(x^2)+x])))
a=0;b=2;n=10;
NIRect[g,a,b,n]//N

5.4 Numerical integration: Trapezoidal rule

A more accurate method for solving an integral like (1) is to use n trapezoids instead of n rectangles (see notes given in class for more details). You must understand where this formulation comes from.

(∫_a^b)f(x)dx≈[(1/2)f((Subscript[x, 1]))+(1/2)f((Subscript[x, n+1]))+(∑_(i=2)^n)f((Subscript[x, i]))] Δx,

where as in the rectangle rule

Δx=(((b-a))/n),

and

(Subscript[x, i])=a+(i-1)Δx.

Note how although we use n trapezoids the function needs to be evaluated up to f((Subscript[x, n+1])) whereas in the rectangle rule the function only needed to be evaluated up to f((Subscript[x, n])). Can you tell why?

Here's how this can be coded within a new function that I will call NITrap

NITrap[f_,a_,b_,n_]:=(Δx=(((b-a)/n));Δx((1/2)f[a]+(1/2)f[b]+(∑_(i=2)^n)f[a+(i-1)Δx]))

f[x_]:=Sin[x]
a=0;b=π/2;n=10;
NIRect[f,a,b,n]//N
NITrap[f,a,b,n] //N

From these outputs it is quite clear that the trapezoid rule is more accurate using the same number of trapezoids as the rectangle rule

5.5 Numerical integration: Simpson's rule

The derivation for this rule was done in class. See your notes for details. Briefly, this method fits a quadratic function to blocks of three points of the function to integrate. Mathematically, the rule is,

(∫_a^b)f(x)dx≈[f((Subscript[x, 1]))+4f((Subscript[x, 2]))+2f((Subscript[x, 3]))+4f((Subscript[x, 4]))+⋯+2f((Subscript[x, n-1]))+4f((Subscript[x, n]))+f((Subscript[x, n+1]))] ((Δx)/3)

Or more compactly,

(∫_a^b)f(x)dx≈[f((Subscript[x, 1]))+f((Subscript[x, n+1]))+2(∑_(i=1)^(n/2-1))f((Subscript[x, 2i+1]))+4(∑_(i=1)^(n/2))f((Subscript[x, 2i]))] ((Δx)/3)

Be careful with the index when you code this method. 
Important note: Simpson's rule can only be applied by dividing the area under the curve  into an even number of elements.  In other words, n must be even.

5.6 Numerical integration of discrete functions

Often in Physics we do not deal with continuous functions. Often we deal with a discrete function, e.g. a set of measurements, that also needs to be integrated to find the solution to a problem.
Here's an illustrative example of this. Consider a set of measurements of the velocity of a particle as a function of time sampled at 0.5 Hz. This gives the following table of data, where the first column is the time in seconds and the second column is the velocity in m/s.

({{0.0,    0.050}, {0.5`,    0.136}, {1.0,    0.366}, {1.5`,    0.967}, {2.0,    2.431}, {2.5`,    5.486}, {3.0, 10.200}, {3.5`, 14.914}, {4.0, 17.968}, {4.5`, 19.432}, {5.0, 20.033}, {5.5`, 20.263}, {6.0, 20.350}, {6.5`, 20.381}, {7.0, 20.393}, {7.5`, 20.397}, {8.0, 20.399}, {8.5`, 20.399}, {9.0, 20.399}, {9.5`, 20.400}, {10.0, 20.400}})

Here's the question: Given this set of measurements what is the total distance traveled by the particle from t=0 to t=10 s? If the velocity v(t) were a continuous function we would write that

v=(dx/dt)

from which we would solve for position as

Total distance=∫dx=(∫_0^10)v dt  	(2)

But this is not the case here because we have a discrete functions. Instead we need to make a numerical integration.
Let's enter these data into a two-column list. I'll call this list vt:

Remove["Global`*"]

vt={{0,0.050},{0.5,0.136},{1.0,0.366},{1.5,0.967},{2.0,2.431},{2.5,5.486},{3.0,10.200},{3.5,14.914},{4.0,17.968},{4.5,19.432},{5.0,20.033},{5.5,20.263},{6.0,20.350},{6.5,20.381},{7.0,20.393},{7.5,20.397},{8.0,20.399},{8.5,20.399},{9.0,20.399},{9.5,20.40},{10.0,20.400}};

Let's plot the data. Always a good idea.

ListPlot[vt,AxesLabel->{"t (s)","v (m/s)"}]

I am purposely not joining the dots to make it clear that this is a discrete function.
Below are two new functions DNIRect and DNITrap that I created to make the integration of a discrete function when tabulated in a 2 column list. The D in DNIRect is for "discrete".

DNIRect[list_]:=(Δx=list[[2,1]]-list[[1,1]];n=Length[list]-1;Δx((∑_(i=1)^n)list[[i,2]]))

DNITrap[list_]:=(Δx=list[[2,1]]-list[[1,1]];n=Length[list]-1;Δx((1/2)list[[1,2]]+(1/2)list[[n+1,2]]+(∑_(i=2)^n)list[[i,2]]))

These two functions I created will always integrate the full list of discrete measurements provided. I haven't made it flexible enough to allow for performing an integral over a sub-range of the list of measurements. This would be better and more general but would make the example more complicated. We'll leave these functions like that for now.
Note that here, since we are dealing with a discrete function we do not specify the number of element n to use for the integration, as we did before when performing a numerical integration of a continuous function. The number n is determined by the number of data points in the tabulated function. You always want to be as accurate as possible so you use the maximum number of data points available. Note that for n elements, i.e. n rectangles for the rectangle rule, there are n+1data points. 
We can then answer the question by doing the numerical integration with the rectangle and the trapezoid rule like this (i.e. solving equation (2) above):

DNIRect[vt]
DNITrap[vt]

Here we have no real means to tell which integration gives the more accurate result. We know however, from previous tests with continuous functions, that the trapezoid rule is more accurate. Be should therefore believe that the answer returned by the trapezoid rule is more accurate, unless of course we made an error in coding it which is something you always have to be VERY careful about.

5.7 Comparing rectangle rule and trapezoid rule

The goal of this section is to compare visually the error associated with the rectangle and the trapezoidal rule as a function of the number of elements (rectangles or trapezoids) used. This will require some coding.  
The integral to evaluate is

(∫_a^b)f(x)dx   (1)

where

f(x)=sin(x)+(x^3)
a=0
b=5

1. Define f(x) so that for later manipulation you want have to deal with the full expression sin(x)+(x^3).

Remove["Global`*"]

f[x_]:=Sin[x]+(x^3)

2. Plot this function between x=-10 and x=10.

x1=-10;x2=10;
Plot[f[x],{x,x1,x2},AxesLabel->{"x","f(x)"}]

3. Calculate the exact value of the integral (1) above.

a=0;b=5;
(∫_a^b)f[x]\[DifferentialD]x//N

4. Make your own code to evaluate numerically the integral (1) using the rectangle rule with 100 rectangles.

NIRect[f_,a_,b_,n_]:=(Δx=(((b-a)/n));Δx(∑_(i=1)^n)f[a+(i-1)Δx])

a=0;b=5;n=100;
NIRect[f,a,b,n]//N

5. make your own code to evaluate numerically the integral (1) using the trapezoid rule with 100 trapezoids.

NITrap[f_,a_,b_,n_]:=(Δx=(((b-a)/n));Δx((1/2)f[a]+(1/2)f[b]+(∑_(i=2)^n)f[a+(i-1)Δx]))

a=0;b=5;n=100;
NITrap[f,a,b,n]//N

6. Briefly comment on the answers you obtained in 4. and 5.

Comment: With the same number of elements the trapezoid rule is more accurate.

7. The absolute relative error of the numerical integration can be calculated as follow:

ϵ=((|Numerical Solution - Exact Solution|)/(|Exact Solution|))

Code this expression for the relative error and calculate ϵ for the results obtained in 4. and 5. above.

a=0;b=5;n=10;
ϵ1=((Abs[NIRect[f,a,b,n]-(∫_a^b)f[x]\[DifferentialD]x])/(Abs[(∫_a^b)f[x]\[DifferentialD]x]))//N
ϵ2=((Abs[NITrap[f,a,b,n]-(∫_a^b)f[x]\[DifferentialD]x])/(Abs[(∫_a^b)f[x]\[DifferentialD]x]))//N

8. Redo the same as in 7. but this time calculate the error when using 100 rectangles or trapezoids to evaluate (1).

a=0;b=5;n=100;
ϵ1=((Abs[NIRect[f,a,b,n]-(∫_a^b)f[x]\[DifferentialD]x])/(Abs[(∫_a^b)f[x]\[DifferentialD]x]))//N
ϵ2=((Abs[NITrap[f,a,b,n]-(∫_a^b)f[x]\[DifferentialD]x])/(Abs[(∫_a^b)f[x]\[DifferentialD]x]))//N

9. It should be clear from your previous calculations that the error reduces significantly by increasing the number of elements and that the trapezoid rule is always more accurate than the rectangle rule.  
What would be illuminating would be to plot the relative error for each algorithms (i.e. rectangle and trapezoid rules) as a function of the number n of elements used to approximate numerically the integral. Before coding this you need to become a little more familiar with the command Table.
First, use the command Table to automatically create the following list: {1,2,3,4,5,6,7,8,9,10}

Table[i,{i,1,10}]

Now, use the command Table to automatically create the following list: {2,4,6,8,10,12,14,16,18,20}

Table[i,{i,2,20,2}]

Try again to create the following list {2,4,8,16,32,64,128}

Table[(2^x),{x,1,7}]

With the command Table you can also create a matrix, i.e. a list of lists. This will become useful for plotting later because we will need to plot the error ϵ as a function of the number n of elements used. We will therefore need to create a 2 column matrix. 
The following matrix  {{1,2},{2,4},{3,8},{4,16},{5,32},{6,64},{7,128}} can be created like this:

Table[{x,(2^x)},{x,1,7}]

Make sure you understand this syntax of the command Table. We will use this command a lot in the remaining of this course.

The advantage of creating a two column matrix like we did is that we can plot one column as a function of the other column, like a tabulated function, with the command ListPlot. For example try these two instructions:

MyList=Table[{x,(2^x)},{x,1,7}];

ListPlot[MyList]

You should see that this has created a graph of the tabulated function defined with the command Table.

10. So now the challenge is to create a two column matrix with the first column containing the number n of elements used in the numerical integration and in the second column the relative error ϵ. 
Here's how I would code this. Here I am evaluating the error of the rectangle rule with 1,2,3,...,100 rectangles

Make sure you understand this syntax. Here NIRect is the name of the function I have created to perform a numerical integration with the rectangle rule.

errRect=Table[{n,((Abs[NIRect[f,a,b,n]-(∫_a^b)f[x]\[DifferentialD]x])/(Abs[(∫_a^b)f[x]\[DifferentialD]x]))},{n,1,100}];

Now you can plot the results of this analysis like this:

ListPlot[errRect,Joined->True,AxesLabel->{"n","ϵ"}]

11. Redo the same analysis with a graph for the trapezoid rule.

errTrap=Table[{n,((Abs[NITrap[f,a,b,n]-(∫_a^b)f[x]\[DifferentialD]x])/(Abs[(∫_a^b)f[x]\[DifferentialD]x]))},{n,1,100}];

ListPlot[errTrap,Joined->True,AxesLabel->{"n","∈"}]

12. Finally compare the results obtained in 11 for the rectangle rule and in 12 for the trapezpoid rule on a single graph. Draw a brief conclusion to the results obtained graphically.

ListPlot[{errRect,errTrap},Joined->{True,True},AxesLabel->{"n","∈"},PlotStyle->{{Black,Thickness[0.01]},{Red,Thickness[0.01]}},PlotLegends->Placed[{Style["Rect.",Black,FontFamily->"Times",18],Style["Trap.",Red,FontFamily->"Times",18]},{1,0.5}]]

Comment: The trapezoidal rule converges much faster than the rectangle rule to the exact solution.