fk3.m

来自「Incorporating Prior Knowledge in Cubic S」· M 代码 · 共 61 行

function  [k] = fk3(spa,spb,spc,te,M0)
% FK3 Megbecsuli a k1,k2,k3-at a spline gorbebol
%
%

splines{1} = spa;
splines{2} = spb;
splines{3} = spc;

Y = M0*ones(length(te),1);
X = zeros(length(te),3);

%--------------------------------------------------------
%X matrix: Integrals of A,B,C
%--------------------------------------------------------
for k=1:3,
 m = length(te);
 for j = 1:m,
  X(j,k) = intspline1(splines{k},te(j));
 end
end

%--------------------------------------------------------
%Y vector: M0-(A+B+C)
%--------------------------------------------------------
for k=1:3,
 knots = splines{k}.x;
 n = length(knots);
 m = length(te);
 for j = 1:m,
  %search index of interval
  i = 1;
  while (i<n) & (te(j)>knots(i+1)),
   i = i+1;
  end
  if (i>=n) | (te(j)<knots(i)),
   error('Error: Out of intervalls!');
  end
  %Calculate the spline value in te(j) 
  k1 = knots(i); k2 = knots(i+1);
  xx = te(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;
  yy = c*splines{k}.y(i) + a*splines{k}.dy(i) + d*splines{k}.y(i+1) + b*splines{k}.dy(i+1);
  %Add to Y
  Y(j) = Y(j)-yy;
 end
end

%--------------------------------------------------------
% Solve least-squares problem
%--------------------------------------------------------

Ane = [];
bne = [];
Ane = -eye(3,3);
bne = zeros(3,1); 
T = lsqlin(X,Y,Ane,bne,[],[]);

k = T;