fspline3.m

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

function  [SP1,SP2,SP3] = fspline3(x,y,knots,s,lambda1,lambda2,w,sum1,sum2,AA,bb,AAne,bbne)
% FSPLINE3 Splines identification function for three cubic spline
%
%  [sp1,sp2,sp3] = fspline3(x,y,knots,s,lambda1,lambda2,w,sum1,sum2,A,b,Anoeq,bnoeq)
%
%   x: (n,3) matrix, y: (n,3) matrix / x -> y /
%   knots: vector /common/
%   s: points for soft constr. /vector/
%   lambda1,lambda2: weight of value and deriv. constr.
%   w: stoichiometric of A,B,C (value and deriv.) /vector/
%   sum1,sum2: desired values without lambda /vector/
%   A,b: equality constr.

if (nargin<11),
  error('Error: Not enough input arguments');
end

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

%--------------------------------------------------------
%Destination vector
%--------------------------------------------------------
if (ns>0)
  Y = [ y(:,1); y(:,2); y(:,3); ones(ns,1)*sum1*lambda1 ];
  X = zeros(m+m+m+ns,2*n+2*n+2*n);
else
  Y = [ y(:,1); y(:,2); y(:,3) ];
  X = zeros(m+m+m,2*n+2*n+2*n);
end

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

%--------------------------------------------------------
% Second part of X matrix: amount of spline-values
%--------------------------------------------------------
j = 1;
while (j<=ns),
  for l=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
    k1 = knots(i); k2 = knots(i+1); xx = s(j); 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(3*m+j,(l-1)*2*n+i*2-1) = lambda1*w(l)*c;
    X(3*m+j,(l-1)*2*n+i*2+0) = lambda1*w(l)*a;
    X(3*m+j,(l-1)*2*n+i*2+1) = lambda1*w(l)*d;
    X(3*m+j,(l-1)*2*n+i*2+2) = lambda1*w(l)*b;
  end
  % next point in s()
  j = j+1;
end

%--------------------------------------------------------
% Solve least-squares problem
%--------------------------------------------------------
warning off
if nargin<10,
  T = lsqlin(X,Y,[],[],[],[]);
else
  if nargin<12,
    T = lsqlin(X,Y,[],[],AA,bb);
  else
    T = lsqlin(X,Y,AAne,bbne,AA,bb);
  end
end
warning on
%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]);