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

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)=a_nxⁿ + a_n-1x^n-1 +……….+ a₁x + a₀ and is represented by a vector in Matlab:
p=[ aⁿ a^n-1 ……. a¹ a⁰ ]

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 A.
roots(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)=3x³ + 2x² -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

[help]What do you think? Are these values correct, if we use the plot below? Make some thinking and check your result.[/help]

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)=3x³ + 2x² -2x + 4. The time-derivative becomes: p’(x)= 9x² + 4x -2 and integration gives: P(x)= (3/4)x⁴ + (2/3)x³ – x² + 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

[warning]That we only obtain the coefficients in the new polynomials. Introduce another polynomial q(x)=x.[/warning]

q=[1 0]

Multiply the polynomial q(x) with p(x), it becomes: pq(x)= 3x⁴ + 2x³ – 2x² + 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)= t₁/ ( x-n₁) ) + t₂/ ( x-n₂) + t₃/ ( x-n₃) + 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