Back to November Calendar            ←Previous Entry                             Next Entry→

November 14, 2005

Spline Representations VI

Let’s repeat the same thing using MatLab, just to see how to form these notion/motions in that environment.

I start with the coefficient matrix:

>> A = [1 2 0 0 0;1 4 1 0 0;0 1 4 1 0;0 0 1 4 1;0 0 0 2 1];

Then the vector of the knot’s y-coordinates:

>> y = [0 2 1 3 1];

Next I dimension a row vector b of 5 zeros and load it up with the weighted average of slopes we need to compute the slopes at the knots and ends of the not knot spline:

>> b=zeros(1,5);

>> for i = 1:4 mu(i) = (y(i+1)-y(i))/1; end;

>> b(1)=(5*mu(1)+mu(2))/6

b =

    1.5000         0         0         0         0

>> b(5)=(5*mu(4)+mu(3))/6

b =

    1.5000         0         0         0   -1.3333

>> for i = 2:4 b(i) = mu(i-1)+mu(i); end;

>> b=3.*b

b =

    4.5000    3.0000    3.0000         0   -4.0000

>> b=b';

I transpose the b vector to a column vector so that we can solve for the tangent line slopes like so:

>> m = A\b

m =

    6.5833

   -1.0417

    0.5833

    1.7083

   -7.4167

Next, I use the formulas described above to solve for the 3^rd and 2^nd order coefficients, and then, just to be consistent, I use the same “add a k” naming scheme to name the 1^st order and constant coefficients:

>> for k=1:4 ak(k) = (m(k)+m(k+1))/h^2-2*(y(k+1)-y(k))/h^3; end;

>> for k=1:4 bk(k) = 3*(y(k+1)-y(k))/h-(2*m(k)+m(k+1))/h; end;

>> ck=m;

>> dk=y;

Now, there’s got to be a better way to do this, but I resorted to giving each domain section a separate name and then plotted each as a separate x,y pair as shown:

>> x1=[0:.05:1];

>> x2=[1:.05:2];

>> x3=[2:.05:3];

>> x4=[3:.05:4];

>> plot(x1,splines03(x1,ak,bk,ck,dk,1),x2,splines03(x2,ak,bk,ck,dk,2),x3,splines03(x3,ak,bk,ck,dk,3),x4,splines03(x4,ak,bk,ck,dk,4) ,'LineWidth',2)

The spline03 function is the following m-file:

function y = splines03(x,a,b,c,d,k)

y = a(k).*(x-k+1).^3.+b(k).*(x-k+1).^2.+c(k).*(x-k+1)+d(k);

So the plots are neatly spliced together in a spline:

What is h is not 1? I don’t know.  I tried it and got non-sense!

How about a bigger example?  So I return to National Oceanic Atmospheric Administration for Mt. San Jacinto elevation data, like I did back on July 20. 

Entering latitude/longitude data and checking the box for “” I get header file with this info:

file_title               = My Selection

data_type                = raster

grid_cell_registration   = center

map_projection           = Lat/Lon

left_map_x               = -117.00000000000

right_map_x              = -116.90000000000

upper_map_y              = 33.00000000000

lower_map_y              = 34.00000000000

number_of_rows           = 120

number_of_columns        = 12

grid_size                = 0.00833333333

elev_m_unit              = meters

minimum_value            = 122

maximum_value            = 1679

elev_m_min               = 122

elev_m_max               = 1679

elev_m_missing_flag      = -500

number_of_display_colors = 256

data_value_unit          = elev_m

data_byte_order          = little_endian

 So, evidently the “grid size” is 1/120 = 0.00833333 of a degree, which, if you take the radius of Earth as 6370 km, is a distance of about 1/120*(Pi/180)*6.37e6 = 926m.  While this works along a Meridian, but the longitudinal distances would be somewhat less because such a circular cross-section of Earth would have a smaller radius. 

The data returned are the elevations at lattice points on a 120x12 grid.  What I propose is to take 1/12 (every twelfth, including the first, 13^th, etc.) of the latitude data for some particular longitude, fit a cubic spline and then compare it with the line plot for the larger data set.  Sound like fun?  Here we go!

As for the right side, I’ve got these data copied into Excel, so it’s easy enough to do the computations there.  Wait a moment, would you?  Hold on…

First enter the coefficient matrix in, say, K2:U12 like so:

Then enter MINVERSE(K2:U12) in, say, W2 and highlight an 11x11 block of cells, say, W2:AG12, and hold down ctrl + shift + enter simultaneously to produce the inverse matrix:

In MatLab, you’d enter

>> A = [1 2 0 0 0 0 0 0 0 0 0 ; 1 4 1 0 0 0 0 0 0 0 0 ; 0 1 4 1 0 0 0 0 0 0 0 ; 0 0 1 4 1 0 0 0 0 0 0 ; 0 0 0 1 4 1 0 0 0 0 0 ; 0 0 0 0 1 4 1 0 0 0 0 ; 0 0 0 0 0 1 4 1 0 0 0 ; 0 0 0 0 0 0 1 4 1 0 0 ; 0 0 0 0 0 0 0 1 4 1 0 ; 0 0 0 0 0 0 0 0 1 4 1 ; 0 0 0 0 0 0 0 0 0 2 1 ]

and then, to find the inverse,

>> inv(A)

ans =

  Columns 1 through 7

    2.1547   -1.1547    0.3094   -0.0829    0.0222   -0.0060    0.0016

   -0.5774    0.5774   -0.1547    0.0415   -0.0111    0.0030   -0.0008

    0.1547   -0.1547    0.3094   -0.0829    0.0222   -0.0060    0.0016

   -0.0415    0.0415   -0.0829    0.2902   -0.0777    0.0208   -0.0056

    0.0111   -0.0111    0.0222   -0.0777    0.2888   -0.0774    0.0207

   -0.0030    0.0030   -0.0060    0.0208   -0.0774    0.2887   -0.0774

    0.0008   -0.0008    0.0016   -0.0056    0.0207   -0.0774    0.2888

   -0.0002    0.0002   -0.0004    0.0015   -0.0056    0.0208   -0.0777

    0.0001   -0.0001    0.0001   -0.0004    0.0016   -0.0060    0.0222

   -0.0000    0.0000   -0.0001    0.0002   -0.0008    0.0030   -0.0111

    0.0001   -0.0001    0.0001   -0.0004    0.0016   -0.0060    0.0222

  Columns 8 through 11

   -0.0004    0.0001   -0.0001    0.0001

    0.0002   -0.0001    0.0000   -0.0000

   -0.0004    0.0001   -0.0001    0.0001

    0.0015   -0.0004    0.0002   -0.0002

   -0.0056    0.0016   -0.0008    0.0008

    0.0208   -0.0060    0.0030   -0.0030

   -0.0777    0.0222   -0.0111    0.0111

    0.2902   -0.0829    0.0415   -0.0415

   -0.0829    0.3094   -0.1547    0.1547

    0.0415   -0.1547    0.5774   -0.5774

   -0.0829    0.3094   -1.1547    2.1547

Though, computationally, it’s usually a waste (and often not possible) to compute the inverse matrix for the purpose of solving a system of equations.  A tridiagonal matrix with loads of symmetry like:

A =

     1     2     0     0     0     0     0     0     0     0     0

     1     4     1     0     0     0     0     0     0     0     0

     0     1     4     1     0     0     0     0     0     0     0

     0     0     1     4     1     0     0     0     0     0     0

     0     0     0     1     4     1     0     0     0     0     0

     0     0     0     0     1     4     1     0     0     0     0

     0     0     0     0     0     1     4     1     0     0     0

     0     0     0     0     0     0     1     4     1     0     0

     0     0     0     0     0     0     0     1     4     1     0

     0     0     0     0     0     0     0     0     1     4     1

     0     0     0     0     0     0     0     0     0     2     1

ought to lend itself to efficiencies.  Generally, a good computing library such as MatLab’s has built in these efficient methods, so when you invoke the command x=A\b, a more efficient method is used to find the solution to Ax = b.

So let’s set about to find the right hand side (RHS) for this system.  We start with elevation values

>> y = [564 163 604 556 1002 615 714 646 454 799 1024];

Then, as before

>> for i = 1:10 mu(i) = (y(i+1)-y(i))/1; end;

>> h=1;

>> b(1)=(5*mu(1)+mu(2))/(6*h);

>> b(11)=(5*mu(10)+mu(9))/(6*h);

>> for i = 2:10 b(i) = mu(i-1)+mu(i); end;

>> b=3.*b

Thus, the right hand side is

b =

        -782

         120

        1179

        1194

         177

        -864

          93

        -780

         459

        1710

         735

And the tangent line slopes at the knots are (note the 1000 factor that applies to all):

>> m=A\b

m =

  1.0e+003 *

   -1.5482

    0.3831

    0.1358

    0.2527

    0.0476

   -0.2659

    0.1519

   -0.2487

    0.0629

    0.4560

   -0.1771

Implement the formulas for the coefficients.  I call them ak, bk, ck and dk:

>> for k = 1:10 ak(k) = (m(k)+m(k+1))/h^2-2*(y(k+1)-y(k))/h^3; end;

>> for k = 1:10 bk(k) = 3*(y(k+1)-y(k))/h-(2*m(k)+m(k+1))/h; end;

>> ck=m;

>> dk=y;

The domain variable is set up in an 11x10 matrix where each row is the domain of cubic:

>> for i=1:11 for j=1:10 xi(i,j) = (i-1)*h + (j-1)*0.1; end; end;

The “hold on” command is used so that the various plots don’t erase the previous.  The splines04 m-file is listed too:

>> hold on

>> for i=1:10 plot(xi(i,:),splines04(xi(i,:),ak,bk,ck,dk,i,h)); end;

function y = splines04(x,a,b,c,d,k,h)

y = a(k).*(x-h*(k-1)).^3.+b(k).*(x-h*(k-1)).^2.+c(k).*(x-h*(k-1))+d(k);

Then I add dots to the knots...dots’ knots’ nice!

>> hold on

>> xy=[0,1,2,3,4,5,6,7,8,9,10];

>> plot(xy,y,'o');

Now, if you recall, The 11 control points we choose were actually just (about) 1/12 of the elevation data we obtained for longitude 116.96˚.  Here’s the whole set of 120 data points:

>> msj = [564 517 487 454 471 445 401 385 276 182 138 143 163 291 300 379 446 522 609 589 652 492 599 592 604 554 544 373 306 286 303 270 263 298 330 470 556 551 570 736 654 1065 1185 1225 1479 1250 1397 1256 1002 935 916 953 928 630 710 511 447 436 430 480 615 682 674 709 728 650 645 647 593 594 675 787 714 737 790 770 659 633 556 527 517 504 505 516 646 634 503 493 497 491 485 480 471 475 475 465 454 478 804 1089 782 693 601 612 657 747 753 782 799 769 779 791 821 819 837 856 894 930 988 1024]

These are plotted by setting up the appropriate 120 domain values (had to delete the last one to match the dimensions) and then the line connecting these (the better approximation to the mountain’s cross-section at that longitude) are plotted in red:

>> xmsj=[0:1/12:10]

>> xmsj(121)=[]

>> plot(xmsj,msj,'r')