In this tutorial we will deal with analysis of functions, interpolation, curve fitting, integrals and differential equations. Firstly, we will need to use polynomials and therefore we have to be familiar with the representation of these. A general polynomial looks like: p(x)=anxn + an-1xn-1 +……….+ a1x + a0 and is represented by a vector in Matlab [1]:
p=[ an an-1 ……. a1 a0 ]
Here we have a list of basic commands dealing with polynomials.
polyval(p,x): Calculates the value of polynomial p for different x. If x is a vector then the polynomial is evaluated for each element in the vector x.
poly(A): Gives a vector that represents the characteristic polynomial for the matrix [2] A.
roots [3](p): Gives a vector with the zeros for the polynomial p(x)=0.
polyder(p): Gives a vector that represents the time-derivative of the polynomial p(x). The coefficients are sored in the vector p.
conv(p,q): Multiplies the polynomials p and q with each other. Returns a coefficient vector.
polyint(p): Integrates the polynomial p analytically and uses the constant of the integration c. The constant c is assigned to 0, if it is not explicitly given.
residue(p,q): Makes a partial fraction expansion of p(x)/q(x).
Example 1: Zeros of a polynomial
Represent the polynomial p(x)=3x3 + 2x2 -2x + 4 in Matlab and find its zeros. Let’s plot the function and check the zeros. This gives a quick idea of what the function looks like. See the resulting figure below.
>> x=-10:0.1:10;
>> plot(x,3*x.^3+2*x.^2-2*x+4), title('p(x)=3*x^3+2*x^2-2*x+4')
>> xlabel('x'), grid
Define the polynomial. The coefficients in the polynomial are arranged in descending order in the vector p. The orders that are nonzero in the polynomial will be represented by zeros in the vector p.
>> p=[3 2 -2 4] % Represents the coefficients in p(x)
With polyval we can easily calculate the value of the polynomial in different x-values.
>> polyval(p,6), polyval(p,7), polyval(p, -5)
ans= 712 , ans = 1117 , ans = -311
What do you think? Are these values correct, if we use the plot below? Make some thinking and check your result.
Plot of a polynomial
Let us try some of the other functions that can be applied to polynomials, like polyder and polyint. They perform the time-derivative and integrate the polynomials p(x)=3x3 + 2x2 -2x + 4. The time-derivative becomes: p'(x)= 9x2 + 4x -2 and integration gives: P(x)= (3/4)x4 + (2/3)x3 – x2 + 4x+C
Now compare what Matlab gives us:
>> C=1 % C is a integration constant.
>> Pder=polyder(p), Pint=polyint(p,C)
Pder =
9 4 -2
Pint =
0.7500 0.6667 -1.0000 4.0000 1.0000
Note that we only obtain the coefficients in the new polynomials. Introduce another polynomial q(x)=x.
>> q=[1 0]
Multiply the polynomial q(x) with p(x), it becomes: pq(x)= 3x4 + 2x3 – 2x2 + 4x and Matlab gives:
>> conv(p,q) % Performs a polynomial multiplication of p and q.
ans= 3 2 -2 4 0
Let us continue with other functions. Now, check the zeros in the polynomials. This is done with the Matlab command root.
roots(p) % Gives a vector with zeros to the polynomial p(x).
ans =
-1.6022
0.4678 + 0.7832i
0.4678 - 0.7832i
Above we can see something quite obvious. There are 3 zeros to a third order polynomial. It is nothing to be astounded by, but only one of these zeros is real. Can we foretell this by looking at the plot in first figure in the tutorial. I would say yes, because if we zoom the curve, we can find the zero-crossing. This gives us a real-valued zero. In our example there is only one, but what happens to the other two zeros? Since they are complex-conjugated, they are not visible.
Finally we will also take a look at the residue command. It makes a partial fraction expansion of the polynomials p(x) and q(x). Look at the ratio q(x)/p(x)!
[ t,n,the_rest]=residue(q,p) % There are 3 output arguments from residue.
t =
-0.1090
0.0545 - 0.0687i
0.0545 + 0.0687i
n =
-1.6022
0.4678 + 0.7832i
0.4678 - 0.7832i
the_rest =
[]
A partial fraction expansion looks like: R(x)= t1/ ( x-n1) ) + t2/ ( x-n2) + t3/ ( x-n3) + the_rest
Now let us define a function in Matlab. As you hopefully remember this is nothing more than a m-file. We will call it func.m and it should be used together with an input argument x. X is a real-valued argument.
% The function func.m created by MatlabCorner.Com
function f=func(x)
f=3*sin(x)-2+cos(x.^2);
We will now take a look at a plot of the function, but first we must decide what region we are interested in. See the plot below.
>> x=-10:0.1:10;
>> plot(x,func(x)), grid on % Note that: fplot(@func,[-10,10])
>> title( 'func(x)') % works equally well.
The command fplot accepts a string as an input argument, but also a handle, @func. We will also use handles later when dealing with figure windows for more advanced plotting purposes, as well as when we work with GUI.
Regional Polynomial Plot
In the previous example the zeros for a polynomial were decided with the command roots. Here on the other hand we have no polynomial, but only a function, and we can instead use the command fzero. The fzero command uses a repetetive algorithm. We must always add an initial guess and then the fzero tries to localize a zero closest to the guess.
Assume we would like to find the zero in the interval 0-1 in the example above. The function has an infinite number of zeros outside this interval 0-1. Our guess becomes:
>> fzero(@funk, 0.5), fzero(@funk, 0.9), fzero(@funk, -1.5)
They all produce the same zero!
ans= 0.3423
Our three guesses seems to use the fact that all of them have the same sign of the time-derivative, which and is why the algorithm converges toward 0.3423. If we change the initial guess to be on the other side of the peak, fzero will give us a new answer (zero).
>> fzero(@funk, 1.5 )
ans = 1.7053
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 Figure 3.
>> 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);