home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 5.ddi / SPLINES.DI$ / FNINT.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  2.0 KB  |  56 lines
  1. function intgrf=fnint(f)
  2. %FNINT    Indefinite integral of a function represented by splines.
  3. %        intgrf=fnint(f)
  4. %
  5. %  returns the (representation of the) indefinite integral of the function 
  6. %  contained in  f  (and in the same form). 
  7.  
  8. % C. de Boor /  Feb.25, 1989
  9. % C. de Boor /  Mar.21, 1991 (correct misprint in line 15)
  10. % C. de Boor /  Sep.23, 1991 (add more knots at right end for B-form)
  11. % C. de Boor /   6 Feb  1992  (finish the last correction)
  12. % Copyright (c) 1990-92 by Carl de Boor and The MathWorks, Inc.
  13.  
  14. if (f(1)==10), % the function is in pp-representation:
  15.    [breaks,coefs,l,k,d]=ppbrk(f);
  16.    coefs=coefs./(ones(d*l,1)*[k:-1:1]);
  17.    if (l<2),
  18.       intgrf=ppmak(breaks,[coefs zeros(d,1)],d);
  19.    else
  20.      % evaluate each integrated polynomial at the right endpoint of its 
  21.      % interval (adapted from  ppual )
  22.       xs=diff(breaks(1:l));index=[1:l-1];
  23.       if (d>1),
  24.          ones(d,1)*xs; xs=ans(:)';
  25.          1+d*ones(d,1)*index+[-d:-1]'*ones(1,l-1); index=ans(:);
  26.       end
  27.       vv=xs.*coefs(index,1)';
  28.       for i=2:k,
  29.          vv = xs.*(vv + coefs(index,i)');
  30.       end
  31.       if (d>1),
  32.          junk=zeros(d,l-1);junk(:)=vv;last=(cumsum([zeros(1,d);junk']))';
  33.       else,
  34.          last=cumsum([0,vv]);
  35.       end
  36.  
  37.       intgrf=ppmak(breaks,[coefs(:,1:k) last(:)],d);
  38.    end
  39.  
  40. elseif (f(1)==11), % the function is in B-repr.
  41.    % Set it up so that it would be correct on the interval [t(1),t(n+k)].
  42.    % There is no way to make it correct everywhere since the integral of
  43.    % a spline over the interval [t(1),t(n+k)] need not be zero.
  44.    [t,a,n,k,d]=spbrk(f); 
  45.  
  46.    index = find(diff(t)>0); % increase multiplicity of last knot to  k
  47.    needed = index(length(index)) - n; % =  k+1 - (n+k - index(length(index));
  48.    if (needed > 0), 
  49.       t = [t t(n+k)*ones(1,needed)];a = [a zeros(1,needed)]; n = n+needed;
  50.    end
  51.  
  52.    intgrf=spmak([t t(n+k)],cumsum((a.*(ones(d,1)*((t(k+[1:n])-t(1:n))/k)))')');
  53. else,  % unknown representation
  54.    error('unknown function type encountered')
  55. end
  56.