In practical real life engineering designs, we often come across problems which consist of numerical operations such as integration or differentiation. For example, in a standard tracing of a robotic arm, the computation of the total error involves integration of instantaneous errors. While mathematical problems require analytic solutions which are sometimes very hard or impossible to obtain, most engineering problems require only numerical solutions. Sometimes the variation of values with parameters is considered sufficient instead of actually finding the closed form relation in terms of parameters. Such a method which involves the numerical computation of the solution is known as a numerical method. These methods form the backbone of engineering and are widely available in MATLAB. In the present chapter, we will study a few basic numerical techniques with which almost all numerical analysis can be done.

One of the prime concerns in numerical methods is that computation platforms only support discrete data. There is no sinusoid signal in numerical software, but only time sampled data of a sinusoidal signal. You can increase the sampling rate to make it closer to a continuous signal, but true continuity can never be achieved. The most important thing in discrete representations is step size and sampling with such a step size gives a sequence of numbers with which we represent the data in MATLAB.

Differentiation

Recall that differentiation of any continuous function f(x) with respect to x is defined as

and f’(x) is known as the derivative of f. For discrete level methods, h can be small but can never be too close to 0. Instead, we choose a very small value (sufficiently small) and compute the derivative with that h. The accuracy of computation is dependent on the value of h. However, a small h size will generate large sized vectors when sampling. It is generally found that with sufficiently small h, we can achieve an accuracy which is accurate enough for all practical purposes. However, what the sufficient value of h is must be determined by taking care of the nonlinearity of the function and sample points.

Let us compute the derivative of the function . We will first choose a value of h and then compute as the derivative.

We will first make a function func which defines the objective function:

function y=func(x)
y=sin(x).*exp(-x.^2);

You can easily guess the implementation for the discrete time differentiation of func:

x=4;h=0.01;
ddxofy=1/h*(func(x+h)-func(x));

We can repeat the above computation to compute the derivative for a whole range, say [0:.01:1]. Let us loop over the x range and calculate the derivative over all points using the above computation,

i=1;h=0.01;xrange=0:h:1;
for x=xrange
       ddxofy(i)= 1/h*(func(x+h)-func(x));
       y(i)=func(x);i=i+1;
end
plot(xrange,y);
hold on;
plot(xrange,ddxofy);

Figure 3-1 shows the output MATLAB displays. We can see that even if we don’t know the analytic form of the derivative, we can see the behavior of the derivative and can easily find its value at any x.

The following example attempts to perform the same operation via vectorization. We first construct a vector x with elements distance h apart, then we compute the function over the whole vector and then subtract the adjacent element’s function values to get f(x+h)-f(x) and finally divide by h. We can use the diff function, which computes the adjacent differences. For example,

x=[1 3 4 2 4];
y=diff(x);

will give y=[2 1 -2 2], which has one element less than the original vector. Using diff, we can implement the following

x=[0:h:1+h];
y=func(x);
ddxofy=diff(y)/h;
x=x(1:end-1);

Since diff reduces the length by 1, we need to reconstruct x such that x and ddxofy have the same size and maintain the element-to-element correspondence. One approach is to remove the last element as performed in the above implementation. Another approach can be to center the interval via the following code

x=(x(1:end-1)+x(2:end))/2;

Another important aspect in all such computations is the selection of the step size. h should be small enough to give a good approximate solution, but if h is smaller than such “good” h it will cause an unnecessarily high computation time due to the increased size of the vector, without any significant improvement in accuracy.

Partial Differentiation

We can perform the same steps for multivariable functions to compute the partial derivative with respect to each variable. Consider the objective function

We will first define an anonymous function to represent the objective function

f=@(u,v) exp(-u.^2-v.^2)

Notice the use of dots to make the function vectorized so that it can be directly applied to a whole matrix or vector to compute the function value elementwise. Now, since the function has two inputs, its representation will be two-dimensional. Instead of constructing the vector, we will construct a 2D grid with one dimension representing u and the other v, with the help of meshgrid

urange=[0:.01:1];
vrange=xrange;
[U V]=meshgrid(urange,vrange);

U and V are two-dimensional matrices such that (U(i, j), V(i, j)) denotes the x, y coordinate of the (i, j)th point in the Cartesian plane. We can evaluate the whole function over the grid

fvalue=f(U,V);

Now we can call diff in the row direction (second dimension) (i.e. the difference of adjacent elements in the same row) to compute and divide by du to get

dduofz=diff(z,1,2);

Similarly, for we need to evaluate diff in the first dimension:

ddvofdz=diff(z,1,1);

or just

ddvofdz=diff(z);

MATLAB also provides a function called gradient to evaluate the same

[px,py] = gradient(z,ustepsize,vstepsize);

Similarly, we can compute the derivative of higher dimensional matrices and higher order derivatives. Let us consider, for example, the computation of a Jacobian. We just need to do the previous operation for each function output component. Recall that the Jacobian is defined as

Let’s consider the function

We can compute the Jacobian via the following code

range = -2:0.2:2;
[x,y] = meshgrid(range);
Z1 = (x.^2 + y.^2);
[px1,py1] = gradient(z1,.2,.2);
Z2 = (x.^2 - y.^2);
[px2,py2] = gradient(z2,.2,.2);
J(:,:,1)=[px1];
J(:,:,2)=[px2];
J(:,:,3)=[py1];
J(:,:,4)=[py2];

Computing Higher Derivatives

Applying the same steps multiple times will result in higher order derivatives. However, we should remember that diff decreases the vector size by 1, so we should rectify the length at each differentiation step by either of the two methods discussed earlier. The following example describes this process

x=[0:h:1+h];
y=func(x);
ddxofy=diff(y)/h;
x1=1/ 2 (x(1:end-1)+x(2:end));
d2dx2y=diff(ddxofy)/h;
x2=1/ 2 (x1(1:end-1)+x1(2:end));

Integration

In a discrete simulation environment, definite integration is just a fine summation. For example, consider the following function,

and the integral , which gives the area under the curve in Figure 3-2(a), and which can be approximated by the shaded area in Figure 2(b). Tuning the value of h, this approximation can be made more accurate. In mathematical terms, the same approximation can be represented as

In MATLAB or any numerical software, we cannot implement the continuous representation but we can always implement the discrete representation and evaluate the approximation with some error. Let us first define the x vector and evaluate the function over it

h=1/100;
x=[0:h:1];
z=sin(2x);

and then compute the summation via

intz=sum(z)*h;

which equals  0.4639, which is very close to the actual value of 0.4597. Since trapezoidal division-based approximation is more accurate, MATLAB provides the trapz function, which computes the summation by trapezoidal area approximation instead of rectangle area approximation.

intz=trapz(x,z);

This results in 0.4596, which is far closer to the actual value. Indefinite integrals can also be computed in numerical form. For other values, you cannot compute the exact analytical form of an indefinite integral, but you can find its numerical values at different x and plot its variation, which may be sufficient from an engineering perspective.

h=1/100;
x=[0:h:1];
z=sin(2x);
intz=cumsum(z)*h;
x1=x(1:end-1);

Multi-dimensional Integration

Let us consider a two-dimensional function and integration over both the variables. This can be extended to higher dimensional functions. Consider the function

We are interested in evaluating

We construct a meshgrid to evaluate the function over the full range and then use trapz twice to compute the double integral. When used on 2D matrices, trapz computes the integration for each column and returns a row vector which is equivalent to integration with respect to v to get

We can apply trapz again on I(u) to compute the integral:

[U V]=meshgrid(range);
f=func(U,V);
Iu=trapz(V,f);
I=trapz(U(1,:),Iu);

We see that, since each of the rows of U are the same and Iu is just a vector, we need only the first row of U to serve as x in trapz.

Integration Over an Infinite Range

If the function is such that its integral over the infinite range can be approximated by an integral over a finite range, we can compute the integration using MATLAB. One example of such a function is any function which decreases quickly as . The crucial point here is to find the finite range which can approximate the integral with small errors.

Consider the following function

Suppose we are interested in computing the integral over the range 0 to infinity. Since the function decreases fast enough, let us compute the integral over the ranges [0 2], [0 10] and [0 100].

upperlimit=[2 10 100];
for i=1:3
        x=0:.01:upperlimit(i);
        y=exp(-2*x);
        I(i)=trapz(x,y);
end

You can see that I(2) and I(3) are almost the same and are different from I(1). So [0 10] is a good enough range for this particular example of integration. Moreover, as usual, you don’t want to take a bigger range (for example, [0 100] here) if the same results can be obtained with a lower range (for example, [0 10] here) because a bigger range will result in a bigger vector x which will slow down the computation without any significant improvement in accuracy.

Multidimensional Integration over Non-rectangular Intervals

The integration need not be over a rectangular area. Consider the following integral

which is equivalent to the integral

The procedure is almost the same, but we don’t have non-rectangular matrices in MATLAB. So we define an auxiliary function g(x, y). We first find a rectangular area which covers the given region R and construct a meshgrid over this rectangular area.

range=0:.01:10;
[U V]=meshgrid(range);

Then we compute the value of the function f for all grid points in the same fashion.

f=func(U,V);

To get the value of the function g from the function f,we just make all the values corresponding to grid points outside the given interval equal to zero.

g=f;
g(U<V)=0;

Now we can compute the integral of this modified function g over the rectangular region which will be equal to the integral of the original function f over the original non-rectangular region.

Iu=trapz(V,g);
I=trapz(U(1,:),Iu);

Solving Equations

In this section, we will talk about solving a single variable equation which may contain linear or nonlinear functions. For example, suppose we want to solve

In other words, we are interested in finding the zeros of a function. This problem is one of the classical problems of numerical analysis. It has variety of applications, such as finding the minimum of a function by finding the zeros of its first derivative, finding the equilibrium of a system, etc.

Polynomial Functions

If the function is a polynomial, its zeros can be found using the roots function. Refer to the second chapter or a MATLAB reference on polynomials for how polynomials can be defined. For example, if we want to find the roots of , we do the following:

P=[1 2 2 0 1];
x=roots(P);

Zeros of a General Function

MATLAB provides a very wide set of functions to solve numerical problems, especially for finding zeros and optimization. We will explore more of them throughout the book. Here, we will focus on a special function called fzero, which is used to find the solution of a one variable function.

Consider the simple function

The important thing to remember is that MATLAB will attempt to solve this equation by numerical techniques. Numerical techniques are very different to analytic techniques. To illustrate the difference, let us consider the following simple example

An analytic approach is to transfer -5 to the RHS to get .

However, a numerical approach is iterative. The simple procedure will start with an initial guess x = 1 and compute its value which is less than 2. It will first increase the value of x to 5 and compute . It will observe that f(x) is getting closer to the RHS, so it will increase further to 9 and observe that the function value f(x) = 4 has exceeded the RHS. It would then decrease the value to the mid-value of x = 7 and find that f(x) = 2 matches the RHS and output the answer. This was a linear approach and therefore the solution algorithm was very easy. Even in such a case, the number of steps to reach the answer depends on the initial guess. As the equation degree or complexity increases, the need for a more sophisticated algorithm becomes evident. Also, for some cases, the choice of initial guess may affect the output you get.

Now let us concentrate on solving the original problem with fzero. Like any numerical method, to call fzero, we need an initial guess x0,

X0=1;
x=fzero('cos(x)-x',x0);

The underlying algorithm behind fzero is very simple and consists of two steps.

Interval search step: It first searches for a value x1 for which the function value differs in sign from the function value at x0. When it finds such a value, it sets the search interval as [x0 x1].

Zero search step: It now uses interpolation to find an approximate zero x2. If f(x2) has the same sign as f(x1), the next search interval is set to [x0 x2] otherwise the new search interval is [x2 x1]. It then repeats the process. It may need many iterations to reach the desired level of accuracy.

If we know an interval for which the function values differs in sign at its boundaries, we can directly specify this interval rather than an initial guess at the solution

x=fzero('cos(x)-x',[0 1]);

A numerical method may have multiple options. For example, termination condition specifications, such as when to terminate or what the tolerance error is after which the iteration should terminate, or display specifications. You can specify these specifications using optimset. For example, you can specify the tolerance by writing

x=fzero('cos(x)-x',x0,optimset('TolX',1e-3));

The optimset function returns an option structure which contains control parameters for the fzero function. More information can be seen by typing help optimset and fzero. Similarly, we can instruct fzero to display each iteration result in the command window by setting display in optimset

x=fzero('cos(x)-x',x0,optimset('TolX',1e-3,'Display','iter'));

The result would be something like

Search for an interval around 0 containing a sign change:
Func-count a f(a) b f(b) Procedure
      1 0 1 0 1 initial interval
      3 -0.0282843 1.0278  0.0282  0.971316 search
      5 -0.04 1.0392  0.04 0.9592 search
      7 -0.0565685   1.0549  0.0565  0.941832 search
      9 -0.08 1.0768 0.08    0.916802 search
 11 -0.113137   1.1067  0.1131  0.88047      search
 13 -0.16 1.14723 0.16    0.827227 search
 15 -0.226274 1.20078 0.2262  0.748235 search
 17 -0.32 1.26924 0.32 0.629235 search
 19 -0.452548 1.3518  0.45254 0.446787 search
 21 -0.64 1.4421  0.64  0.162096 search
 23 -0.905097 1.5227  0.9050  -0.287487 search
Search for a zero in the interval [-0.905097, 0.905097]: 
Func-count x f(x) Procedure
 23 0.905097      -0.287487 initial
 24 0.61761 0.197656 interpolation
 25 0.734737       0.00727045 interpolation
 26 0.739099      -2.37581e-005 interpolation
 27 0.739099      -2.37581e-005 interpolation

fzero also accepts function handles

f=@(x) sin(x)-.4;
x=fzero(f,0);

or

x=fzero(@myfunc1,[0 1]);

where myfunc1 is defined in myfunc1.m

If the function is dependent not only on x but also on a few more inputs, say c, and you want to know the value of x such that , you can define an anonymous function g(x) = f(x, c) for fixed x and fzero must be called over this g(x). The following example describes the same method for the variable c.

for c=0:1:10
g=@(x) sin(x)-c/10;
x(c+1)=fzero(g,[0 3.14]);
end

Interpolation

When we have discrete data instead of continuous data, the need to reconstruct the continuous data is evident. For example, with the knowledge of output voltage sampled at .1 seconds starting from t=0, can we compute (approximate) the value at 2.314 sec? In many cases we don’t have data for all the points, but only at some. In that case, the common problem is to generate the data for other points. This procedure is known as interpolation. In other words, given x_n for n = 1,2, . . . k+n and y_n values for n=1,2,…k, we need y_n for where .

One-dimensional Interpolation

One-dimensional interpolation can be done using the interp1 function. Suppose we know the vector Y containing the function values y=f(x) at x values contained in X and we wish to evaluate the function over a finer grid XI. We can do so by using the interp1 function in the following format:

YI = interp1(X, Y, XI, method),

where the vectors X and Y are the vectors holding the x- and the y- coordinates of points to be interpolated, respectively, XI is a vector holding the points of evaluation, i.e., new points, and method is an optional string specifying an interpolation method.

The following methods work with the function interp1.

Nearest Neighbor Interpolation

Method = ‘nearest’. Produces a locally piecewise constant interpolant.

Linear Interpolation

Method = ‘linear’. Produces a piecewise linear interpolant.

Cubic Spline Interpolation

Method = ‘spline’. Produces a cubic spline interpolant.

Cubic Interpolation

Method = ‘cubic’. Produces a piecewise cubic polynomial.

In this example, the y data for the x points are given as

x = 0:pi/5:pi;
y = sin(2.*x);

which need to be interpolated using the method of interpolation ‘nearest’ . The interpolant is evaluated at the following points

xi = 0:pi/100:pi;
yi = interp1(x, y, xi, 'nearest');

Spline Interpolation

The MATLAB function spline is designed for computations with cubic splines (n = 3) that are twice continuously differentiable (k = 2) on the interval. The MATLAB command

yi = spline(x, y, xi)

evaluates a cubic spline s(x) at points stored in the array xi.

Data Fitting and Polynomial Interpolation

One method of interpolation is to fit an nth degree polynomial to the data and then calculate polynomial values at interpolated points. This can be done easily as

p=polyfit(x,y,n);
yi=polyval(p,xi);

For example, suppose we are given data x and y as

x=0:.01:1;
y=x.^3+2*x+0.02*rand(size(x));

Let us not worry about how x and y are generated and assume the vectors are given to us. We can fit a third degree polynomial to this data as

p=polyfit(x,y,3)

which gives us p = [0.9686 0.0489 1.9758 0.014]. Once we get the polynomial coefficients, the values at any x or range can be computed as

xi=0:.005:2;
yi=polyval(p,xi)

Figure 3-3 shows the MATLAB output when x,y and xi,yi are plotted on the same graph. Since the limits of xi are beyond x, it is in fact an extrapolation.

Arbitrary Interpolation

You may need to interpolate your data according to your model. This problem comes under model parameter analysis and can be solved using optimization. We will revisit this method in later chapters.