Analysis of Functions, Interpolation, Curve Fitting, Integrals and Differential Equations 2

Localize minima and maxima of functions

Let us try to find the local minima and maxima for the function func(x). The interval of interest is [-6 0]. The algorithms are iterative. There are 2 methods to use. The first one decides x in a given interval, and the second one looks for x around an initial guess. To decide the maxima we are looking for an x that minimizes the negative function: -func(x).

fminbnd(func,x1,x2): Gives the x-value for a local minima.

[x,y,flag,info]=fminbnd(func,x1,x2)

As above, but we also get the y-value. Flag gives a positive number, if minima is localized and a negative number otherwise. We store the information about the iterations in the variable info.
fminsearch(func, x0,): Gives the x-value for a local minima somewhere close to the initial guess x0.
Decide the global minima and maxima that exist on the interval -8 to 0 for the function func(x). This gives:

x =
-1.7326
y =
-5.9511
flag =
1

This seems to be true for our function. If we want to find the maxima values, the same command can be used. The only difference will be a minus sign in front of the function. Whenever we call the function fzero. The command looks for a minima but will in fact localize the maxima due to the minus sign. See below!

function f=func(x)
f=-(3*sin(x)-2+cos(x.^2)); % Note the minus sign in front of parenthesis.
>> fminbnd(@func,-8,0)

Look below for the answer. Note that the y-value is negated. Compare the plot func(x to our result below.

x = y =
-1.8691 -5.0046

Our maxima occurs at x=-1.8691, and the value of y becomes y=5.0046. Is this an accurate result according to your opinion?

Interpolation of data sets

Assume that we have made a number of measurements and thereby guessed a function relating the input and output. Then we also have some knowledge of the values between the measured points. This is very common in sampled applications. Think of a sampled continuous sine wave function with a very low sampling frequency. This can be done very simply, if we use a vector x with low resolution. Let us create one. Use the x vector values and calculate the corresponding sine function values and plot it in bellow Figure.

x=0:20;
y=sin(x), plot(x,y),grid



Low resolution sine plot

We could think of it as a rather badly shaped sine wave curve probably due to a too long sampling period. By interpolation one decides on a function P(x) that passes through certain known data points. The interpolated function can be used to calculate the function value approximately for any x. Normally one uses polynomials or functions combined from polynomials (spline functions) to make interpolations, but also sine and cosine waves can be useful.
Now, let us see what happens when we try to make a curve fit 10 random numbers.

y=rand(1,10) % 10 random numbers in a row vector
plot(y,'*')

See the figure below. Investigate whether it is possible to find a polynomial that can be adjusted to the random data.
Enter the figure window under Tools-> Basic Fitting. Mark the ninth-degree polynomial. This gives a polynomial that passes through all random data points. The basic fitting uses the least square method. Even if we use a lower degree it is still the best solution according to the least square method.

Mark the small box show equation. We can then explicitly see the polynomial in the figure window. Remove the ninth-degree polynomial and choose a linear equation instead. Matlab will now try to find the best linear equation. Mark the small box plot residuals. The error residual between each data point and the linear equation is calculated and shown. Despite large discrepancies this is the best first-order polynomial.

The error residual between each data point and the linear equation

We will now try to accomplish the same thing with Matlab commands that we did with the interactive basic fitting in the figure window.

y=rand(1,10), x=1:10;
plot(x,y,'*') ,hold on
p1=polyfit(x,y,1) % Best coefficients for a first order polynomial
% according to the least square method. Stored in
% a ector p1.
p9=polyfit(x,y,9) % Ninth-order coefficients.
xx=1:0.01:10; % Improves the resolution in the x-vector.
p1_func=polyval(p1,xx);% Calculates polynomial values for x-values.
p9_func=polyval(p9,xx);% As above.
plot(xx,p1_func,xx,'',p9_func,'-.')
legend('y=randn(1,10)','linear','9th degree')

The result of the commands can be seen in the figure below. The higher the order of the polynomial, the higher the oscillation between the data points. Notice the huge peaks between values 1 and 2 and 9 and 10. You should have been alarmed if this had been real data and your mission was to find an interpolation function that nicely fits the measurements. It is definitely not common sense to look for a ninth-order polynomial for a data set of 10 points.

The higher the order of the polynomial, the higher the oscillation between the data points

Integrals

Matlab can solve integrals and differential equations with numerical methods. Later on in this course we will show how integrals and differential equations can be solved in Symbolic Math toolbox. Numerical integration in Matlab is performed by the command quad. This stands for numerical quadrature. We can solve integrals only with defined limits. This can be done for single, double and triple integrals. To solve generalized integrals or integrals expressed with symbols, we must use the Symbolic Math toolbox. Let us exemplify with some numerical examples.

Example 1: Single integrals

Decide the integral below.

First we put the integrand in a function file: integrand.m.

function y=integrand(x)
y=sqrt(exp(x));

Write!

I1=quad('integrand',0,10)

Matlab gives:

I1 = 294.8263

It seems reasonable enough, if we plot the function. We can also use quad with a fourth argument, namely the tolerance for the calculation. This decides when the numerical integration will stop. The two following examples can be hard to understand, if one lacks the background of multivariable calculus.

Example 2: Double integrals

Calculate the double integral below.

Make as in the previous example a functionwhere we put the primitive function. We name the m-file as integrand2.m

function z=integrand2(x,y)
z=sqrt(exp(x)).*sqrt(exp(y));

Let us now call the function dblquad instead of it.

I2=dblquad('Integrand2',0,1,0,1)

and the answer is:

I2=1.6834

We will now take a look of the integrated area.

[X,Y]=meshgrid(0:0.1:1,0:0.1:1); % Makes a grid of calculation points.
Z=Integrand2(X,Y); % Calculates the integral.
mesh(X,Y,Z) ; % Makes a plot that connects these points.
% punkter.

The last three lines achieves Figure below.

Double integral

Matlab can also perform triple integrals. This example can be found by help triplequad.

Example 3: Triple integrals

A function f = y*sin(x)+z*cos(x) is integrated for the intervals 0 < x < pi, 0 < y < 1 and -1 < z < 1. This is achieved by:
First we create a function file integrnd.m.

function f = integrnd(x,y,z) % 3 input arguments and 1 output argument.
f = y*sin(x)+z*cos(x);

We refer to this m-file by using its handle. In the matlab command window we simply write:

triplequad(@integrnd,0,pi,0,1,-1,1);

You may also like...