Curve
Fitting with Scilab
Neither
Scilab nor Scicoslab have a function for straight curve fitting, such
as the
polyfit function that we can find in Matlab. However, it’s not that
difficult
to develop (or find?) a custom made function for that purpose.
In [1]
we can find a suggestion for this task. This is an edited version of
the transcription
of the 15line code.
function
p = polyfit(x, y, n)
if
length(x) ~= length(y)
error('x
and y vectors must be the same
size')
end
x = x( : );
y = y(
: );
V =
ones(length(x), n+1);
for j = n : 1 : 1
V(:,
j) = x .* V (:, j+1);
end
[Q, R] = qr(V);
QTy = Q' * y;
p = R(1 : n+1, 1 : n+1) \ QTy(1 : n+1);
p = p.';
endfunction
There
are three inputs: vectors x and y, and a scalar n. x and y are the
vectors
defining our points to be fitted, and n is the order of the polynomial
to be
found.
On the
other hand, we also need an equivalent of the polyval function given in
Matlab.
Naturally there are many possibilities here. The following code is much
simpler
than the previous one and we’re showing our suggestion.
function
y = polyval(p, x)
y = 0*x;
p = mtlb_fliplr(p);
for ix = 1 : length(p)
y =
y + p(ix) * x.^(ix1);
end
endfunction
The
goal of this function is to use the polynomial found by
polyfit, receive an
xvalue and release the corresponding yvalue produced by that
polynomial. There
are two inputs: we need the coefficients of the polynomial, in p, and
the
xvalue, in x. We just do the required math in a simple and straight
way. We
use the function mtlb_fliplr (provided in the M2SCI toolbox) just to
sort the
order of the exponent of the variable.
Now,
let’s see in action both functions.
//
Prepare our environment
clear;
clc
//
Declare our two functions
getf('polyfit.sci')
getf('polyval.sci')
//
Define (arbitrarily) the number of points to take into account
np =
10;
//
Define the xvector and two functions, the second function
// is a
noised version of the first one
x =
linspace(0, 1, np);
y1 =
x.^3  5*x.^2  3*x  7;
y2 = y1
+ .1*rand(1, np);
//
Enter the x and y vectors, and the order of the polynomial
// that
we want to obtain
p =
polyfit(x, y2, 3)
//
Define other xvalues and find the original function
x = 0 :
.4 : 1
y1 =
x.^3  5*x.^2  3*x  7
//
Use
polyval to find the equivalent values in the
//
noised function
y =
polyval(p, x)
//
Divide the values just for comparison purposes
ratio =
y1./y
Now, the
(commented) results shown by Scilab are:
//
This
is our 3rd. order polynomial from the noised function
p =
1.1103879

5.0400213  3.052195
 6.9291697
//
Another
set of points generated
x =
0.
0.4
0.8
y1 =  7. 
8.936  12.088
//
Find
the yvalues of the noised function
y =

6.9291697 
8.8853863 
12.028021
//
Compare the original vs the noised results
ratio =
1.0102221
1.0056963
1.0049866
We can
see that the polyfit and polyval functions coded in Scilab are
working pretty
fine.
Reference:
[1]
Medrano Sanchez, C., Plaza Garcia, I.; Software libre para Calculo
Numerico;
Alfaomega; Mexico, DF, 2009.
From
'Curve Fitting Scilab' to Matlab home
From
'Curve Fitting' to Scilab menu
