|
November
16, 2005
Spline
Representations VIII
Here’s a run-through
of the c++ code I finagled for fitting cubic splines. First, here’s the main routine which sets
up the 9 control points and shows their coordinates on the console with showvec calls. The first call to spline with ic
= 1 computes the coefficients c in  while subsequent calls (in the for
loop) with ic = 0 compute actual
spline y-values.
Note:
- Variable
y is a function of x here, instead of x a function of t, as previously formulated. We may want to go back to parameterizing
x, y and z in terms
of t, but, for this instance
of code, we’ve got y = f(x).
- Lower
case x and y contain control point
data while X, Y contain both control point
data and spline data.
- I
commented out the graphmm call to the MatLab
m-file because, well, this programs runs independent of MatLab. It would be fairly simple to plot the
curve with OpenGL.
:
|
// do the
spline thing
#include <vector>
#include <iostream>
#include <cmath>
using namespace std;
void showvec(char *, vector<double>,
int);
double spline(double, vector<double>
&, vector<double> &,
vector<double> &, int, int);
double
tridec(vector<vector<double> >
&, vector<double> &, int, vector<double>
&);
double epsilon(void);
int MIN(int, int);
int main (void)
{
/* This is hardwired to fit the function 1/(1+x^2)
to
samples corresponding to x = -4, -3, -2, ... , 4.
n = number of sample to fit polynomial on
m = number of samples in the finer
resolution of fitting polynomial
i and j are just counters--can't use 'em
in MatLab
d = the distance from 0 of the extreme x
values
a = arbitrary x value */
int n = 9, m = 10*n, i, j;
double d = 4, a;
/* x and y are
intial samples. c are linear term
coefficients X and Y are
2x90 vectors containing all info
about coordinates to be plotted
create x, y and c zero-vectors of
dimension n */
vector<double> x(n,0.), y(n,0.), c(n,0.);
vector<vector<double> > X, Y;
vector<double>
dummy(m,0.);
X.push_back(dummy); X.push_back(dummy);
Y.push_back(dummy); Y.push_back(dummy);
cout << "\nCubic Spline
Fit\n";
for (i =
0; i < n; i++) {
a = 1-d + 2*d*(i-1)/(n-1);
x[i] = a;
y[i] = 1/(1 + a*a);
}
cout << "\n\nInput/output
control points coordinates are:";
showvec ("x",x,n);
showvec ("y",y,n);
spline (0,x,y,c,n,1);
for (i =
0; i < m; i++) {
a = x[0] + (i)*(x[n-1] - x[0])/(m-1);
X[0][i] = a;
Y[0][i] = spline
(a,x,y,c,n,0);
j = MIN (i,n-1);
X[1][i] = x[j];
Y[1][i] = y[j];
}
//graphmm
(X,Y,2,m,"Cubic Spline
Fit","x","y","fig442",1);
char ch;
cin >> ch;
return
0;
}
|
Then the
spline function, changed from the previous incarnation to be consistent with STL
vectors instead of NLib/MatLab vectors, is shown below. The main difference is that the index of
the first element of a vector is now 0, instead of 1. The constructors are also quite different
in syntax and I changed the vector variables to pass-by-reference.
Also, there is
code to print out the coefficient parameters for the fitting cubics which
we'll look at again in a moment.
|
/*
-------------------------------------------------------------------- */
double spline(double a, vector<double>
&x, vector<double> &y,
vector<double> &c, int n, int ic)
{
/*
-------------------------------------------------------------------- */
/*
Description: Evaluate a cubic spline function at x.
On entry: a
= evaluation point
x = n by 1 vector of strictly increasing
independent
variables
y = n by 1 vector of dependent variables
c = n by 1 vector containing the spline
coefficients,
n = number of sample points (n >= 3)
ic = initial condition code
defined as follows:
0 =>
use previously computed c
1 =>
compute vector c
On exit: spline = natural cubic spline function
evaluated at a
Notes: 1. spline is initialized using ic
<> 0 and then
subsequent calls use ic
= 1.
2. spline uses the function
tridec in linear.c to
compute the cubic spline
coefficients. */
int k =
0, i;
double s, b0, b1, p;
vector<double> h, w, v;
vector<vector<double> > T;
if (ic)
{
/* Compute
spline coefficients */
h.assign( n-1, 0. );
w.assign( n-1, 0. );
v.assign( n, 0. );
T.push_back(v);
T.push_back(v);
T.push_back(v);
/* Initialize right-hand side v
*/
for
(i = 0; i < n-1; i++) {
h[i] = x[i+1] - x[i];
if
(h[i] > 0)
w[i] = (y[i+1] - y[i])/h[i];
else {
cout << "\nThe independent variables in
spline must be";
cout << "\nstrictly increasing.\n";
return
0;
}
}
v[0] = 0.;
v[n-1] = 0.;
for
(i = 1; i < n-1; i++)
v[i] = 6*(w[i] - w[i-1]);
/* Solve
equations */
for
(i = 0; i < n-1; i++)
{
if (i >
0)
{
T[0][i] = h[i];
T[1][i] = 2*(h[i-1] + h[i]);
}
if
(i < n-2)
T[2][i] = h[i];
}
T[1][0] = T[1][n-1] = 1;
/* TRIDEC !!!!!!!!!!!!!!!!!!!!!! */
tridec (T,v,n,c);
showvec("c",c,n);
return 0.;
}
/* Evaluate
spline function at a */
if (ic
== 0) {
while ((a > x[k+1])
&& (k < n-1)) k++; /* go to segment k
*/
p = x[k+1] - x[k];
b0 = (6*y[k+1] - p*p*c[k+1])/(6*p); // coefficients
of linear terms
b1 = (6*y[k] -
p*p*c[k])/(6*p);
s = b0*(a - x[k]) + b1*(x[k+1] - a);
s += (double)
(pow(a-x[k],3)*c[k+1] + pow(x[k+1]-a,3)*c[k])/(6*p);
cout << "\nb0 = " << b0 <<
"\tb1 = " << b1
<< "\tc[" << k+1
<< "] = " << c[k+1]
<<
"\tc[" << k << "] = " << c[k];
return s;
}
}
|
The call to tridec is used to solve the big
matrix equation for the coefficients c. Tridec
is optimized for solving a tridiagonal matrix system.
|
/*
------------------------------------------------------------------- */
double tridec
(vector<vector<double> >
&T, vector<double> &b, int n, vector<double>
&x)
{
/*
------------------------------------------------------------------- */
/*
Descripton: Solve the tridiagonal linear algebraic system Ax = b
using the LU decomposition
method.
On entry: T = 3 by n matrix containing the nonzero
entries of A:
row contents
---------------------------------------------------
0 A[0][1], A[1][2], ... ,
A[n-2][n-1], ---
1 A[0][0], A[1][1], ... , A[n-2][n-2],
A[n-1][n-1]
2 A[1][0], A[2][1], ... ,
A[n-1][n-2], ---
---------------------------------------------------
b = n by 1 right-hand side
vector
n = number of unknowns (n
>= 3)
On exit: x = n by 1 solution vector (cubic
coefficients for splines)
Returned: The determinant of A. If it is nonzero, then
x has been successfully
computed. */
int k;
double eps = epsilon()/100, delta = 1;
vector<double> y;
vector<vector<double> >
Q;
/*
Initilialize */
y.assign( n, 0. );
Q.push_back(y);
Q.push_back(y);
Q.push_back(y);
/* Factor T
into L and U */
Q[2].assign( T[2].begin(), T[2].end() );
Q[1][0] = T[1][0];
for (k =
0; k < n-1; k++) {
delta *= Q[1][k];
if
(fabs(Q[1][k]) < eps)
return
0;
else {
Q[0][k] = T[0][k]/Q[1][k];
Q[1][k+1] = T[1][k+1] -
T[2][k]*Q[0][k];
}
}
delta *= Q[1][n];
/* Perform forward substitution to find
y */
for (k =
1; k < n; k++)
y[k] = (b[k] -
Q[2][k-1]*y[k-1])/Q[1][k];
y[0] = b[0]/Q[1][0];
showvec("tridec y",y,n); //cin >> ch;
/* Perform back substitution to find x
*/
x[n-1] = y[n-1];
for (k =
n-2; k >= 0; k--) {
x[k] = y[k] - Q[0][k]*x[k+1];
}
showvec("x",x,n);
return
delta;
}
|
Assuming that
tridec has returned a non-zero determinant, the vector x (named c from the
calling function) will contain the cubic coefficients for the splines. The calls to spline in the for-loop of the
main function then compute the 90 y-coordinates
of the splines.
There are a
couple of other little functions involved, which may be of interest: MIN, epsilon
and showvec.
|
int MIN(int a, int b) {
return ((a < b) ? a :
b);
}
/*
------------------------------------------------------------------- */
double epsilon (void) {
/*
------------------------------------------------------------------- */
/* Description:
Compute the machine epsilon for type float.
Returns:
The smallest number e such that 1 + e > 1. */
double e
= 1, x;
do {
e /= 2;
x = 1 + e;
}
while (x
> 1);
return
2*e;
}
/*
-------------------------------------------------------------- */
void showvec(char * name, vector<double>
v, int n) {
cout << endl << name << " = [";
for(int i = 0; i < n-1; ++i)
cout << v[i] << ", ";
cout << v[n-1] << "]" << endl;
}
|
Here’s the
output of the program:
|
Cubic
Spline Fit
Input/output
control points coordinates are:
x = [-4,
-3, -2, -1, 0, 1, 2, 3, 4]
y =
[0.0588235, 0.1, 0.2, 0.5, 1, 0.5, 0.2, 0.1, 0.0588235]
tridec y
= [0, 0.0882353, 0.296471, 0.242017, -1.6725, 0.769683, 0.115303, 0.063
675, 0]
x = [0,
0.063675, 0.0982414, 0.74336, -1.87168, 0.74336, 0.0982414, 0.063675, 0]
c = [0,
0.063675, 0.0982414, 0.74336, -1.87168, 0.74336, 0.0982414, 0.063675, 0]
b0 =
0.0893875 b1 = 0.0588235 c[1] = 0.063675 c[0] = 0
b0 =
0.183626 b1 = 0.0893875 c[2] = 0.0982414 c[1] = 0.063675
b0 =
0.376107 b1 = 0.183626 c[3] = 0.74336 c[2]
= 0.0982414
b0 =
1.31195 b1 = 0.376107 c[4] = -1.87168 c[3] = 0.74336
b0 =
0.376107 b1 = 1.31195 c[5] = 0.74336 c[4]
= -1.87168
b0 =
0.183626 b1 = 0.376107 c[6] = 0.0982414 c[5] = 0.74336
b0 =
0.0893875 b1 = 0.183626 c[7] = 0.063675 c[6] = 0.0982414
b0 =
0.0588235 b1 = 0.0893875 c[8] = 0 c[7] = 0.063675
|
It was a bit
tedious, but I then entered these into a graphing program as
Y(x) = 0.0893875(x+4)+0.0588235(-3-x)+(0.063675(x+4)^3)/6
Y(x) =
0.183626(x+3)+0.0893875(-2-x)+ (0.0982414(x+3)^3+0.063675(-2-x)^3)/6
Y(x) =
0.376107(x+2)+0.183626(-1-x)+(0.74336(x+2)^3+.0982414(-1-x)^3)/6
Y(x) =
1.31195(x+1)+0.376107(-x)+(-1.87168(x+1)^3+.74336(-x)^3)/6
Y(x) = 1.31195(-x+1)+0.376107(x)+(-1.87168(-x+1)^3+.74336(x)^3)/6
Y(x) =
0.376107(-x+2)+0.183626(-1+x)+(0.74336(-x+2)^3+.0982414(-1+x)^3)/6
Y(x) =
0.183626(-x+3)+0.0893875(-2+x)+ (0.063675(-2+x)^3+.0982414(-x+3)^3)/6
Or, in pretty
type, and rounding coefficients to 3-digits:
When you graph
these without restricting to the spline domains, you get a graph like this:
Now to compute
the sum of the spline curvatures
Using MatLab, I start by making an m-file
for the cubic spline function in terms of the horizontal coordinates array
(named u, here) and the linear
and cubic coefficients array (named c.)
function y = spline02(a,c,k,u,x)
y =
a(k+1).*(x-u(k))+a(k).*(u(k+1)-x)+(c(k+1).*(x-u(k)).^3+c(k).*(u(k+1)-x).^3)./6;
To be concise:
a = linear coefficients array
c = cubic coefficients array
k = index
u = array of control point inputs
x = domain variable
y =
range (output) variable
In the command
window then you set everything up this way:
>> a = [0.0588235
0.0893875 0.183626 0.376107 1.31195]
a =
0.0588
0.0894 0.1836 0.3761
1.3119
>> c=[0
0.063675 0.0982414 0.74336 -1.87168]
c =
0
0.0637 0.0982 0.7434
-1.8717
>> u = -4:1:4
u =
-4
-3 -2 -1
0 1 2
3 4
>> x1 = u(1):.05:u(2);
>> x2 =
u(2):.05:u(3);
>> x3 =
u(3):.05:u(4);
>> x4 =
u(4):.05:u(5);
>> y1 =
spline02(a,c,1,u,x1);
>> y2 =
spline02(a,c,2,u,x2);
>> y3 =
spline02(a,c,3,u,x3);
>> y4 =
spline02(a,c,4,u,x4);
>>
plot(x1,y1,x2,y2,x3,y3,x4,y4,'LineWidth',2)
To produce this
graph:
Note that we
produced these splines before by NLib’s spline example, but we didn’t know
the actual coefficients.
Now that we
have coefficient arrays, let’s see if we can develop a nice function for the
total curvature
Substituting
,
this becomes
The indices
are off by 1, as seems to be plaguing my work with MatLab…oh well. In MatLab,
we’d maybe write a for loop, like so:
>> splineCurvTot =
0;
>> for k = 1:4
splineCurvTot =
splineCurvTot + 2*(c(k)*c(k)+c(k)*c(k+1)+c(k+1)*c(k+1))/3
end
splineCurvTot =
0.0027
splineCurvTot =
0.0160
splineCurvTot =
0.4395
splineCurvTot =
2.2158
The factor of
2 doubles the curvature for the left side of the symmetric curve. Provided all these computations are good, and I think they
are, the curvature of the spline is about 1% of the curvature of the single
interpolating polynomial.
I haven’t
computed the more complicated formula measuring (perhaps more precisely) the
curvature. That’s a good exercise! Oh, ok, I’ll do it. Here’s my attempt, done all in the command
window
>> a = [0.0588235
0.0893875 0.183626 0.376107 1.31195];
>> c = [0 0.063675 0.0982414 0.74336 -1.87168];
>> u = -4:1:4;
>> x1 = u(1):.05:u(2);
>> x2 =
u(2):.05:u(3);
>> x3 =
u(3):.05:u(4);
>> x4 =
u(4):.05:u(5);
>> y1 =
spline02(a,c,1,u,x1);
>> y2 =
spline02(a,c,2,u,x2);
>> y3 =
spline02(a,c,3,u,x3);
>> y4 =
spline02(a,c,4,u,x4);
>> dydx1 =
diff(y1)./diff(x1);
>> dydx2 =
diff(y2)./diff(x2);
>> dydx3 =
diff(y3)./diff(x3);
>> dydx4 =
diff(y4)./diff(x4);
>> for k=1:20 x1mids(k) =
(x1(k)+x1(k+1))/2; end
>> dy2dx21 =
diff(dydx1)./diff(x1mids);
>> for k = 1:19 dydx1mids(k) =
(dydx1(k)+dydx1(k+1))/2; end
>> for k=1:20 x2mids(k) = (x2(k)+x2(k+1))/2;
end
>> for k=1:20 x3mids(k) =
(x3(k)+x3(k+1))/2; end
>> for k=1:20 x4mids(k) =
(x4(k)+x4(k+1))/2; end
>> for k = 1:19 dydx2mids(k) =
(dydx2(k)+dydx2(k+1))/2; end
>> for k = 1:19 dydx3mids(k) =
(dydx3(k)+dydx3(k+1))/2; end
>> for k = 1:19 dydx4mids(k) =
(dydx4(k)+dydx4(k+1))/2; end
>> dy2dx22 =
diff(dydx2)./diff(x2mids);
>> dy2dx23 =
diff(dydx3)./diff(x3mids);
>> dy2dx24 =
diff(dydx4)./diff(x4mids);
>> int1 =
abs(dy2dx21)./(1+dydx1mids.*dydx1mids).^1.5;
>> int2 =
abs(dy2dx22)./(1+dydx2mids.*dydx2mids).^1.5;
>> int3 =
abs(dy2dx23)./(1+dydx3mids.*dydx3mids).^1.5;
>> int4 =
abs(dy2dx24)./(1+dydx4mids.*dydx4mids).^1.5;
>> curve(1) =
.05*sum(int1)+.025*(int1(1)+int1(19));
>> curve(2) =
.05*sum(int2)+.025*(int2(1)+int2(19));
>> curve(3) =
.05*sum(int3)+.025*(int3(1)+int3(19));
>> curve(4) =
.05*sum(int4)+.025*(int4(1)+int4(19));
>> 2*sum(curve)
ans =
2.1568
This is about the curvature found using the same measure
on the single interpolating polynomial: 8.54760869196934,
as found yesterday.
http://www.utexas.edu/its/rc/tutorials/matlab/
http://dept.seattlecolleges.com/southengineering/introductiontomatlab.pdf
enamel
http://www.applet-magic.com/gaussbf.htm
http://www.cs.berkeley.edu/~sequin/CS184/GuestLecture_05/SplineLect1.html
|