Skip to main content

PreCourse 4

NERC Numerical Earth Science Modelling Short Course: Lecture 4 – Approximate methods for differentiation and integration

Dr Simon A Mathias

Department of Earth Sciences

Durham University

Contents

This week’s leaning objectives

At the end of the week you should be able to:

  • Describe how derivatives can be approximated using finite differences.
  • Approximate the derivative of a variable using diff.
  • Show how finite differences link with trapezoidal integration.
  • Approximate the integral of a variable using cumtrapz.

The main focus of this module, Modelling Earth Processes, is to develop numerical models of various earth processes. Generally, these models involve the numerical solution of sets of partial differential equations. You will learn more about these in subsequent lectures. In this lecture we will look at the more fundamental process of approximating derivatives and integrals using finite differences and trapezoidal integration, respectively.

Finite differences

The basic idea about finite differences is as follows:

$y=f(x)$

Consider a continuous function .

$N$

Now consider a set of number of discrete points along the $x$ axis $x=x_1,x_2,x_3,...,x_N$ .

$$ y $$

The corresponding values of can be written as $y_1,y_2,y_3,...,y_N$ .

$y(x=x_1)=y_1$

Understand that, $y(x=x_2)=y_2$ and so on.

$y$

An approximation of the derivative of with respect to $x$ can be obtained as follows:

$\left. \frac{dy}{dx} \right|_{x=x_{i-1/2}} \approx \frac{y_i-y_{i-1}}{x_i-x_{i-1}}$

where $i = 2,3,... N$

$x=x_{i-1/2}$

The value of the derivative corresponds to the location, where $x_{i-1/2}=(x_i+x_{i-1})/2$ .

$y$

Note that calculus is often referred to as “infinitesimal calculus”. The infinitesimal term refers to the focus on how changes over infinitesimal distances. A definition of a derivative can be written as follows:

$\lim_{x_i\rightarrow x_{i-1}}\left. \frac{dy}{dx} \right|_{x=x_{i-1/2}} = \frac{y_i-y_{i-1}}{x_i-x_{i-1}}$

Approximating derivatives using finite differences

Consider the trigonometric equation

$y = \cos 2x$
$x$

Differentiating with respect to gives

$\frac{dy}{dx} = -2\sin 2x$

Creat a new script file and type:

function ExamplesFromLecture4

% Some examples demonstrating the principal of differentiation and
% integration using finite  differences and trapezoidal integration,
% respectively.

at the top and save in an appropriate place.

$N$

First we will make a series of number of equally spaced points from zero to five using the linspace command (type “help linspace” in the Command Window to find out more). Following from this will calculate some corresponding values of $y$ (where $y=cos 2x$).

%Define number of points for the discretisation of x
N=7;
x=linspace(0,5,N)
%Calculate corresponding values of y
y=cos(2*x)

x =

  Columns 1 through 3

         0    0.8333    1.6667

  Columns 4 through 6

    2.5000    3.3333    4.1667

  Column 7

    5.0000


y =

  Columns 1 through 3

    1.0000   -0.0957   -0.9817

  Columns 4 through 6

    0.2837    0.9274   -0.4612

  Column 7

   -0.8391

Now we will estimate the derivative using finite differences. Note that you can use end in the code to automatically identify the number of elements for a given dimension of an array of concern.

%Calculate the change in x, i.e. x(i)-x(i-1)
dx=x(1,2:end)-x(1,1:end-1)
%Calculate the change in y, i.e. y(i)-y(i-1)
dy=y(1,2:end)-y(1,1:end-1)
%Estimate the derivative with respect to x using finite difference
dydxFD=dy./dx

dx =

  Columns 1 through 3

    0.8333    0.8333    0.8333

  Columns 4 through 6

    0.8333    0.8333    0.8333


dy =

  Columns 1 through 3

   -1.0957   -0.8860    1.2653

  Columns 4 through 6

    0.6437   -1.3886   -0.3779


dydxFD =

  Columns 1 through 3

   -1.3149   -1.0631    1.5184

  Columns 4 through 6

    0.7724   -1.6663   -0.4534

Interestingly, we have only six values of the derivative although we started with seven points. It is important to realise that these derivatives do not correspond directly to the points:

$x_1,x_2,x_3,...,x_N$

Rather they correspond to a different set of points:

$x_{2-1/2},x_{3-1/2},x_{4-1/2},...,x_{N-1/2}$
$(N-1)$

of which there are only in number.

As discussed earlier, it is possible to locate this latter set of points from

$x_{i-1/2}=(x_i+x_{i-1})/2$

, $i=2,3,...,N-1$

Add the following to your code:

%Calculate the locations at which the derivatives apply
xMP=(x(1,2:end)+x(1,1:end-1))/2

xMP =

  Columns 1 through 3

    0.4167    1.2500    2.0833

  Columns 4 through 6

    2.9167    3.7500    4.5833
$-2\sin 2x$

Recall that we know the derivative is , in this case,. So lets compare:

%Calculate derivative analytically at the xMP points
dydz=-2*sin(2*xMP)

dydz =

  Columns 1 through 3

   -1.4804   -1.1969    1.7095

  Columns 4 through 6

    0.8697   -1.8760   -0.5105

Now lets make a plot to compare graphically

%Plot results
figure(1)
hold off %Destroy previous plots in figure(1)
plot(x,y,'k',xMP,dydz,'b',xMP,dydxFD,'ro')
legend('y = cos 2x','dy/dx Analytical','dy/dx Finite difference')
xlabel('x (-)')
ylabel('see legend')
title('Approximate calculus')

The accuracy is quite poor because we are using only seven points for our discretisation. Try running your code again with 20 points.

Although using finite difference techniques to differentiate trigonometric functions seems slightly trivial, we are just using these as a demonstration. The main point to understand is that we can use similar code to approximate the derivative of any function of interest.

Approximating derivatives using the diff command

The diff command is a MATLAB command, which can be used to perform the above analysis in a computationally more efficient manner. Type “help diff” in the command window to find out more.

Note that diff(X,N,DIM), where X is the name of the array to be differenced, N is the number of times the array us to be differences, and DIM is the dimension of the array along which the differencing is to take place.

To aid our demonstration of diff, type the following 2D array in the command window:

A=[1 6 5 4 3 2; 5 6 4 4 2 1; 7 8 5 1 1 2]

A =

     1     6     5     4     3     2
     5     6     4     4     2     1
     7     8     5     1     1     2

Now we are going to look at the difference between each column

A(:,2:end)-A(:,1:end-1)

ans =

     5    -1    -1    -1    -1
     1    -2     0    -2    -1
     1    -3    -4     0     1

The diff command can be used to do exactly the same thing by typing:

diff(A,1,2)

ans =

     5    -1    -1    -1    -1
     1    -2     0    -2    -1
     1    -3    -4     0     1

If instead we wanted to look at the difference between each row, we would type:

diff(A,1,1)

ans =

     4     0    -1     0    -1    -1
     2     2     1    -3    -1     1

If we wanted to look at the difference of each difference between each row, we would type:

diff(A,2,1)

ans =

    -2     2     2    -3     0     2

Go back to your script file that you created and add the following code

dydxFD=diff(y,1,2)./diff(x,1,2)

dydxFD =

  Columns 1 through 3

   -1.3149   -1.0631    1.5184

  Columns 4 through 6

    0.7724   -1.6663   -0.4534

Above we have calculated the derivative in one line. Note that the 1, in the diff argument, is because we are only taking one difference and the 2 is because we want to difference the columns as opposed to the rows. Recall that x and y are row vectors.

Trapezoidal rule and its relation to finite difference

Recall that integration is the reverse of differentiation. Whereas differentiation corresponds to finding the gradient of a line on a graph, integration corresponds to finding the area under a line.

A simple way to approximate an integral is to apply the so-called trapezoidal rule. Click on this link http://mathworld.wolfram.com/TrapezoidalRule.html TrapezoidalRule for an illustration.

Recalling that integration is the reverse of differentiation. Consider again the finite difference approximation

$\left. \frac{dy}{dx} \right|_{x=x_{i-1/2}} \approx \frac{y_i-y_{i-1}}{x_i-x_{i-1}}$

where $i = 2,3,... N$

$dy/dx=J$

For convenience, let us say such that

$J_{i-1/2} \approx \frac{y_i-y_{i-1}}{x_i-x_{i-1}}$

It can also be said

$y=\int J dx$
$y_i$

Rearranging the finite difference equation above such that is the subject of the formula we get

$y_i \approx J_{i-1/2}(x_i-x_{i-1})+y_{i-1}$
$x_{i-1/2}=(x_i+x_{i-1})/2$

In the same way that, it can also be said that $J_{i-1/2}=(J_i+J_{i-1})/2$, from which it can be seen

$y_i \approx \left(\frac{J_{i}+J_{i-1}}{2}\right)(x_i-x_{i-1})+y_{i-1}$

which is essentially a mathematical expression of the trapezoidal rule.

Approximating integrals using the trapezoidal rule

Consider again the trigonometric equation

$y = \cos 2x$
$x$

Integrating with respect to gives

$F = \int y dx = \frac{1}{2}\sin 2x + C$
$C$

where is an integration constant.

$F=3$

Let us now impose the constraint that when $x =0$ . It follows that $C = 3$ . Therefore

$F = 3+\frac{1}{2}\sin 2 x$
$F$

We now wish to approximate using the trapezoidal rule.

From the previous discussion we can say that

$F_i \approx \left(\frac{y_{i}+y_{i-1}}{2}\right)(x_i-x_{i-1})+F_{i-1}$

To explore this further let us add the following code to our script file:

%Initialise a vector of zeros called Ftrapz
Ftrapz=zeros(1,N);
%Estimate the integral of y with respect to x using the trapezoidal rule
for i=2:N
    Ftrapz(1,i)=(y(1,i)+y(1,i-1))/2*(x(1,i)-x(1,i-1))+Ftrapz(1,i-1);
end
%Display result in the command window
disp(Ftrapz)

  Columns 1 through 3

         0    0.3768   -0.0721

  Columns 4 through 6

   -0.3630    0.1416    0.3359

  Column 7

   -0.2059

Now lets add some code to compare the result to the analytical solution:

%Calculate values F using the analytical function from notes
F=3+sin(2*x)/2

F =

  Columns 1 through 3

    3.0000    3.4977    2.9047

  Columns 4 through 6

    2.5205    3.1871    3.4436

  Column 7

    2.7280
$F(x=0)=3$

Note that we also need to apply our constraint to the trapezoidal integration result. This can be done as follows:

%Apply constraint to trapezoidal integration
Ftrapz=Ftrapz-Ftrapz(1,1)+3

Ftrapz =

  Columns 1 through 3

    3.0000    3.3768    2.9279

  Columns 4 through 6

    2.6370    3.1416    3.3359

  Column 7

    2.7941
$x = 0$

What we have done is subtracted the first value of Ftrapz, which corresponds to in this case, from all the values of F. We have then added on the correct value at $x = 0$ , which according to our constraint, is 3.

%Plot results
figure(1)
hold on %Plot new results in figure(1) alongside previous results
plot(x,F,'g',x,Ftrapz,'k:s')
legend('y = cos 2x','dy/dx Analytical','dy/dx Finite difference',...
    '\int y dx Analytical','\int y dx Trapezoidal')
legend('location','eastoutside')

Approximating integrals using the cumtrapz command

The cumtrapz command is a MATLAB command, which can be used to perform exactly the same analysis above in a computationally more efficient manner. Type “help cumtrapz” in the command window to find out more.

Note that cumtrapz(X,Y,DIM), where X is the name of the array to be integrated with respect to, Y is the array to be integrated, and DIM is the dimension of the array along which the integrating is to take place.

Add the following code to your script file to repeat the integration using cumtrapz:

%Use cumtrapz to integrate y with respect to x
Fcumtrapz=cumtrapz(x,y,2);
%Apply constraint
Fcumtrapz=Fcumtrapz-Fcumtrapz(1,1)+3

Fcumtrapz =

  Columns 1 through 3

    3.0000    3.3768    2.9279

  Columns 4 through 6

    2.6370    3.1416    3.3359

  Column 7

    2.7941

Assignment

$y$

Make a new script file and write a subfunction to both differentiate and integrate a given set of dependent and independent variables using diff and cumtrapz (if one differentiates with respect to $x$ , $y$ and $x$ are known as the dependent and independent variables, respectively).

$0 < x < 2$

Then write some additional script to use your subfunction to both differentiate and integrate the following functions for :

$y = 5x^2 + 2x +2$

1)

$y = e^x$

2)

$y = (1+x^2)^{-1}$

3)

$x = 0$

For each equation apply a constraint such that the integrals all equal 1 when.

Write further script to plot your results appropriately on a single figure with three subplots. Compare your results with analytical results where possible. You may find www.wolframalpha.com useful in this respect.


Published with MATLAB® 7.13