Matlab Demo 4
Polynomial and Spline Interpolation
- Working With Polynomials
- Polynomial Interpolation and Approximation
- Spline Interpolation
- Fourier Series
- Nonlinear Least Squares
Working With Polynomials
Polynomials in Matlab are represented as vectors of coefficients, highest power first. There are a few commands for manipulating them, unfortunately no polyadd.
>> help polyfun Interpolation and polynomials.
...
Polynomials.
roots - Find polynomial roots.
poly - Convert roots to polynomial.
polyval - Evaluate polynomial.
polyvalm - Evaluate polynomial with matrix argument.
residue - Partial-fraction expansion (residues).
polyfit - Fit polynomial to data.
polyder - Differentiate polynomial.
polyint - Integrate polynomial analytically.
conv - Multiply polynomials.
deconv - Divide polynomial
Examples:
>> p = poly([1,2,3])
p =
1 -6 11 -6
>> roots(p)
ans =
3.000000000000001
1.999999999999998
1.000000000000001
>> p2 = conv(p,p)
p2 =
1 -12 58 -144 193 -132 36
>> roots(p2)
ans =
3.000000000000197 + 0.000000542678441i
3.000000000000197 - 0.000000542678441i
1.999999999999772 + 0.000000308732273i
1.999999999999772 - 0.000000308732273i
1.000000129717776
0.999999870282281
>> deconv(p2,p)
ans =
1 -6 11 -6
>> [n,d,pp] = residue([1,0,0,0,0],p)
n =
40.499999999999943
-15.999999999999936
0.499999999999998
d =
3.000000000000001
1.999999999999998
1.000000000000001
pp =
1 6
>> polyder(p)
ans =
3 -12 11
>> y = polyval(p,0:0.1:4);
>> x = 0:0.1:4;
>> y = polyval(p,x);
>> plot(x,y)
>> roots(polyder(p))
ans =
2.577350269189626
1.422649730810374
>> P = polyint(p)
P =
0.2500 -2.0000 5.5000 -6.0000 0
Polynomial Interpolation and Approximation
The polyfit command can do both interpolation and approximation.
Example 2: US Population
Here is another example, based on the US population during the 20th century:
>> load('population.dat')
>> population
population =
1900 75994575
1910 91972266
1920 105710620
1930 122775046
1940 131669275
1950 150697361
1960 179323175
1970 203235298
1980 226547082
1990 248709873
2000 281421906
>> t = population(:,1);
>> pop = population(:,2);
>> p = polyfit(t,pop,10)
Warning: Polynomial is badly conditioned. Add points with distinct X
values, reduce the degree of the polynomial, or try centering
and scaling as described in HELP POLYFIT.
> In polyfit at 80
p =
1.0e+021 *
Columns 1 through 8
-0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0000
Columns 9 through 11
0.0000 -0.0053 1.6020
>> tt = 1890:2010;
>> popp = polyval(p,tt);
>> plot(tt,popp,t,pop,'o')
Let's zoom in a little:
>> axis([1890,2010,0,300000000])
What is the problem? The problem is that the x values run between 1900 and 2000. There are large differences between x, x5 and x10, so the coefficients differ by many orders of magnitude. The situation improves if you just change the x axis:
>> t = (0:0.1:1)';
>> tt = (-.1:0.01:1.1)';
>> p = polyfit(t,pop,10)
p =
1.0e+012 *
Columns 1 through 8
0.3131 -1.0786 1.0986 0.3621 -1.6759 1.5091 -0.6674 0.1564
Columns 9 through 11
-0.0181 0.0009 0.0001
>> popp = polyval(p,tt);
>> plot(tt,popp,t,pop,'o')
Matlab has a builtin way of dealing with this: Using polyfit with 3 output arguments centers and scales the data.
>> [p,s,mu] = polyfit(t,pop,10)
p =
1.0e+008 *
Columns 1 through 8
0.0504 0.2365 -0.3403 -1.1223 0.8735 1.7642 -1.0357 -1.0736
Columns 9 through 11
0.6151 0.8741 1.5070
s =
R: [11x11 double]
df: 0
normr: 1.6687e-006
mu =
0.5000
0.3317
>> popp = polyval(p,(tt - mu(1))/mu(2));
>> plot(tt,popp,t,pop,'o')
Even better is to give up on high degree polynomials, and use approximation instead of interpolation:
>> p = polyfit(t,pop,2)
p =
1.0e+008 *
0.9454 1.0775 0.7831
>> popp = polyval(p,tt);
>> plot(tt,popp,t,pop,'o')
Spline Interpolation
The Spline Toolbox
The spline toolbox is not for the faint of heart. It is more for the experts, and not needed for simple spline interpolations. Type help splines for more info.
Standard Matlab (without the toolbox) can handle not-a-knot, clamped, and cubic Hermite splines.
Not-A-Knot Spline
>> x = -5:5;
>> y = 1./(1 + x.^2);
>> xx = -5:0.1:5;
>> yy = spline(x,y,xx);
>> plot(xx,yy,x,y,'o')
Let's look at the inside of a spline structure:
>> pp = spline(x,y);
>> yy = ppval(pp,xx);
>> pp
pp =
form: 'pp'
breaks: [-5 -4 -3 -2 -1 0 1 2 3 4 5]
coefs: [10x4 double]
pieces: 10
order: 4
dim: 1
>> [breaks,coefs,l,k,d] = unmkpp(pp)
breaks =
-5 -4 -3 -2 -1 0 1 2 3 4 5
coefs =
0.0062 -0.0082 0.0224 0.0385
0.0062 0.0104 0.0246 0.0588
0.0069 0.0290 0.0640 0.1000
0.1072 0.0499 0.1429 0.2000
-0.4357 0.3715 0.5643 0.5000
0.4357 -0.9357 0.0000 1.0000
-0.1072 0.3715 -0.5643 0.5000
-0.0069 0.0499 -0.1429 0.2000
-0.0062 0.0290 -0.0640 0.1000
-0.0062 0.0104 -0.0246 0.0588
l =
10
k =
4
d =
1
Clamped Spline
For the clamped spline, put in the values of the endpoint derivatives before and after the y data. In the example below, we want a derivative of 0 at the left end, and -1 at the right end:
>> pp = spline(x,[0,y,-1]);
>> [breaks,coefs,l,k,d] = unmkpp(pp)
breaks =
-5 -4 -3 -2 -1 0 1 2 3 4 5
coefs =
-0.0102 0.0305 0 0.0385
0.0106 0.0000 0.0306 0.0588
0.0058 0.0318 0.0624 0.1000
0.1072 0.0493 0.1435 0.2000
-0.4348 0.3710 0.5638 0.5000
0.4321 -0.9334 0.0014 1.0000
-0.0934 0.3628 -0.5693 0.5000
-0.0583 0.0824 -0.1241 0.2000
0.1855 -0.0925 -0.1342 0.1000
-0.7219 0.4641 0.2374 0.0588
l =
10
k =
4
d =
1
>> yy = ppval(pp,xx);
>> plot(xx,yy,x,y,'o')
Cubic Hermite Spline
In general, Hermite interpolation refers to interpolation with point values and derivatives given. A cubic Hermite spline is a piecewise cubic function with given point values and derivatives at the knots. It has one continuous derivative.
Point values and derivatives give 4 pieces of information per subinterval, which determines the spline uniquely. The derivates are estimated from the point values. Basically, you look at the nearest few points, interpolate a low order polynomial, and differentiate that.
>> pp = pchip(x,y);
>> [breaks,coefs,l,k,d] = unmkpp(pp)
breaks =
-5 -4 -3 -2 -1 0 1 2 3 4 5
coefs =
-0.0035 0.0139 0.0100 0.0385
0.0032 0.0107 0.0272 0.0588
0.0083 0.0333 0.0583 0.1000
-0.0750 0.2250 0.1500 0.2000
-0.6250 0.7500 0.3750 0.5000
0.6250 -1.1250 0 1.0000
0.0750 -0.0000 -0.3750 0.5000
-0.0083 0.0583 -0.1500 0.2000
-0.0032 0.0204 -0.0583 0.1000
0.0035 0.0034 -0.0272 0.0588
l =
10
k =
4
d =
1
>> yy = ppval(pp,xx);
>> plot(xx,yy,x,y,'o')
Here is a comparison of cubic and Hermite splines, side by side:
>> hermite = pchip(x,y);
>> cubic = spline(x,y);
>> yh = ppval(hermite,xx);
>> ys = ppval(cubic,xx);
>> plot(xx,yh,'g',xx,ys,'r',x,y,'o')
Fourier Series
The commands for FFT and inverse FFT are fft and ifft.
Look at the Matlab built-in demos sunspots and fftdemo for examples.
Nonlinear Least Squares
Look at the Matlab built-in demo fitdemo.
This demo does attempts to fit a function of the form
f(x) = c1 exp(λ1x) + c2 exp(λ2x)
to given data. This is done as follows:
Write a function fitfun which takes 3 arguments: lambda, x, y. Lambda has to go first, because we will be optimizing with respect to it.
For a given fixed lambda, this is a linear least squares problem in the c's. Find the optimal c1, c2 and return the residual. Then use fitfun as the input for an optimizing routine (fminsearch, in the demo), and find the optimal lambdas.
Last Updated: February 12, 2009