|
Spline Representations IV
Any road, the NLib example program 4.4.2: Cubic Spline Fit finally produces the spline plot shown here:
Let’s compare this curve with the more well-known 8th degree interpolating polynomial whose coefficients c are obtained by solving the system Where V is the Vandermonde matrix, whose elements are powers of xi. In MatLab, you can then use the backslash operator, \, to solve the least squares problem as c = y\V.
More simply, you can use the m-file polyfit. This is shown here:
Define the domain of the control points: >> t = [-4:1:4] t = -4 -3 -2 -1 0 1 2 3 4
Define the range of the control points: >> y=1./(1+t.*t) y = Columns 1 through 6 0.0588 0.1000 0.2000 0.5000 1.0000 0.5000 Columns 7 through 9 0.2000 0.1000 0.0588
Define the domain x and range f of the polynomial fit points and plot the control points and the interpolating polynomial together. Note that p will contain just the coefficients of the fitting polynomial: >> x = (-4: 0.1: 4)'; >> p = polyfit(t,y,8); >> f = polyval(p,x); >> plot(t,y,'o',x,f,'-')
Happily, this is corroborated by the Scientific Notebook result (using its Maple engine) to enter
How about computing and comparing the measure of total curvature for the interpolating polynomial and the spline?
Here’s some
nice commentary from the MatLab
help on integration: The area beneath a section of a function F(x) can be determined by numerically integrating F(x), a process referred to as quadrature. The MATLAB quadrature functions are:
quad
Use adaptive Simpson quadrature dblquad Numerically
evaluate double integral
To integrate the function defined by humps.m from 0 to 1, use q = quad(@humps,0,1)
q = 29.8583
Both quad and quadl operate recursively. If either method detects a possible singularity, it prints a warning. You can include a fourth argument for quad or quadl that specifies a relative error tolerance for the integration. If a nonzero fifth argument is passed to quad or quadl, the function evaluations are traced.
I should note that I tried to use an m-file to define a polynomial function from its coefficient array, which I succeeded in doing (just a for loop that does pretty much what polyval does) but quad rejected it…I don’t know why. Perhaps it’ll die.
Anyroad, it
got under my hat and between my ears to try using MatLab’s quad for
integrating (aka quadrature) to compute the total curvature where r(s) is the curve parameterized by the arc length, s, and T is the tangent vector.
For another parameterization parameter, say, t, (arc length is often tricky, if not completely impracticable) the formula can be worked this way:
Differentiating dr/dt we get This allows the following formation of the cross product of first and second derivatives of the position vector per parameter t: Equating the norms, I wonder if this derivation for k came about through mere luck and that’s why the justification seems to start with an object out of the blue and end in a formula for curvature we can savor:
This formula works for 3-D space curves with x, y and z parameterized by t, but a planar (2-D) curve with z=0 (the xy plane) is a special case worthy of simplification. In this case, My worry at this juncture, is that, if I am to develop an integral for the accumulated curvature over some length of curve I’ll need to integrate something like which strikes me as an intractable unwieldliness of no low order. However, it would seem to suffice to measure the total curvature as the area under the curvature curve (heh) and since
a more
convenient measure might be function y = polyCurveFunk(a,b) % [a,b] is the interval over which we'll sample % x at integer inputs and % y as the corresponding 1/(1+x^2) function values x = [a:1:b]; y = 1./(1.+x.^2); % n is the number of subintervals in this sampling n = length(x)-1; % p is the interpolating polynomial p = polyfit(x,y,n); % resample x for high resolution x = [a:(b-a)/10000:b]; % generate interpolating polynomial y-values for this higher resolution pv = polyval(p,x); % p1 has coefficients for the derivative function of p p1 = polyder(p); % these are the coefficients of the polynomial for p1^2 p1square = conv(p1,p1); % p2 has the coefficients for the derivative of p1 p2 = polyder(p1); % generate values for terms in curvature function and generate p1squarevals = polyval(p1square,x); p2absvals = abs(polyval(p2,x)); curvature = p2absvals./(1.+p1squarevals).^1.5; % plots plot(x,pv,x,curvature) % zero out left endpoint to get right-endpoint sum of integral curvature(1) = 0; sum(curvature)*(b-a)/10000 % compute simpler measure for comparison sum(p2absvals.*p2absvals)*(b-a)/10000
Here’s calling the function:
>> polyCurveFunk(-4,4); ans = 8.54760869196934 ans = 2.404317873342450e+002
So the
compare that with the sum of the total curvatures of the 8 cubic splines.
To compute the equivalent measures of total curvature for the spline, I think we’ll somehow have to get the actual cubics, rather than just do the magic function call! For me, this has meant commandeering the whole MatLab/NLib into the land of c++. I see three good reasons for this approach:
|
