fspline6.m

Incorporating Prior Knowledge in Cubic Spline Approximation - Application to the Identification of R

function  [SP1,SP2,SP3] = fspline6(x,y,knots,kin,s,lambda1,AA,bb,AAne,bbne)
%Iterativ 6 spline-t identifikalo rutin
%

%Start
n = length(knots);
m = size(x,1);
ns = length(s);
nc = 6;

%--------------------------------------------------------
%Destination vector
%--------------------------------------------------------
Y = [ y(:,1); y(:,2); y(:,3); zeros(ns,1) ];
X = zeros(m+m+m+ns,2*n+2*n+2*n);

%--------------------------------------------------------
%First part of X matrix: A,B,C points
%--------------------------------------------------------
for k=1:3,
  j = 1;
  while (j<=m),
    %search index of interval
    i = 1;
    while (i<n) & (x(j,k)>knots(i+1)),
      i = i+1;
    end
    if (i>=n) | (x(j,k)<knots(i)),
      error('Error: Out of intervalls!');
    end
    %calculate linear parameters
    k1 = knots(i); k2 = knots(i+1); xx = x(j,k); h = k2-k1;
    a = (k2-xx)^2*(xx-k1)/h^2; b = -(xx-k1)^2*(k2-xx)/h^2;
    c = (k2-xx)^2*(2*(xx-k1)+h)/h^3; d = (xx-k1)^2*(2*(k2-xx)+h)/h^3;
    %pack them into X matrix
    X((k-1)*m+j,(k-1)*2*n+i*2-1) = c;
    X((k-1)*m+j,(k-1)*2*n+i*2+0) = a;
    X((k-1)*m+j,(k-1)*2*n+i*2+1) = d;
    X((k-1)*m+j,(k-1)*2*n+i*2+2) = b;
    % next point in x()
    j = j+1;
  end
end

%--------------------------------------------------------
% Second part of X matrix: dca/dt+..+dca*/dt+dcc*/dt = 0
%--------------------------------------------------------
j = 1;
while (j<=ns),
  for k=1:3,
    %search index of interval
    i = 1;
    while (i<n) & (s(j)>knots(i+1)),
      i = i+1;
    end
    if (i>=n) | (s(j)<knots(i)),
      error('Error: Out of intervalls!');
    end
    %calculate linear parameters
    xx = s(j);
    k1 = knots(i); k2 = knots(i+1); h = k2-k1;
    a = (k2-xx)^2*(xx-k1)/h^2;
    b = -(xx-k1)^2*(k2-xx)/h^2;
    c = (k2-xx)^2*(2*(xx-k1)+h)/h^3;
    d = (xx-k1)^2*(2*(k2-xx)+h)/h^3;
    da = (-2*(k2-xx)*(xx-k1)+(k2-xx)^2)/h^2;
    db = -(2*(xx-k1)*(k2-xx)-(xx-k1)^2)/h^2;
    dc = (-2*(k2-xx)*(2*(xx-k1)+h)+(2*(k2-xx)^2))/h^3;
    dd = (2*(xx-k1)*(2*(k2-xx)+h)-(2*(xx-k1)^2))/h^3;
    %add to X
    rw = 3*m+j;  cl = (k-1)*2*n+i*2-1;
    X(rw,cl+0) = X(rw,cl+0) + lambda1*(dc+kin(k)*c);
    X(rw,cl+1) = X(rw,cl+1) + lambda1*(da+kin(k)*a);
    X(rw,cl+2) = X(rw,cl+2) + lambda1*(dd+kin(k)*d);
    X(rw,cl+3) = X(rw,cl+3) + lambda1*(db+kin(k)*b);
  end
  % next point in s()
  j = j+1;
end

%--------------------------------------------------------
% Solve least-squares problem
%--------------------------------------------------------
T = lsqlin(X,Y,AAne,bbne,AA,bb);

%Spline - A
SP1.x = knots;
SP1.y = T([1:2:2*n-1]);
SP1.dy = T([2:2:2*n]);
%Spline - B
SP2.x = knots;
SP2.y = T([2*n+1:2:4*n-1]);
SP2.dy = T([2*n+2:2:4*n]);
%Spline - C
SP3.x = knots;
SP3.y = T([4*n+1:2:6*n-1]);
SP3.dy = T([4*n+2:2:6*n]);