Mercurial > hg > octave-lyh

comparison scripts/polynomial/ppint.m @ 12608:59e2460acae1

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
make piecewise polynomial (pp) functions more compatible
author Kai Habel <kai.habel@gmx.de>
date Wed, 23 Feb 2011 08:11:40 +0100
parents d0b799dafede
children e81ddf9cacd5
comparison
equal deleted inserted replaced
12606:3047363c376d 12608:59e2460acae1
26 26
27 function ppi = ppint (pp, c) 27 function ppi = ppint (pp, c)
28 if (nargin < 1 || nargin > 2) 28 if (nargin < 1 || nargin > 2)
29 print_usage (); 29 print_usage ();
30 endif 30 endif
31 if (! isstruct (pp)) 31 if (! isstruct (pp) && strcmp (pp.form, "pp"))
32 error ("ppint: PP must be a structure"); 32 error ("ppint: PP must be a structure");
33 endif 33 endif
34 34
35 [x, p, n, k, d] = unmkpp (pp); 35 [x, p, n, k, d] = unmkpp (pp);
36 p = reshape (p, [], k); 36 p = reshape (p, [], k);
37 37
38 ## Get piecewise antiderivatives 38 ## Get piecewise antiderivatives
39 pi = p / diag (k:-1:1); 39 pi = p / diag (k:-1:1);
40 k += 1; 40 k += 1;
41 if (nargin == 1) 41 if (nargin == 1)
42 pi(:,k) = 0; 42 pi(:, k) = 0;
43 else 43 else
44 pi(:,k) = repmat (c(:), n, 1); 44 pi(:, k) = repmat (c(:), n, 1);
45 endif 45 endif
46 46
47 ppi = mkpp (x, pi, d); 47 ppi = mkpp (x, pi, d);
48 48
49 ## Adjust constants so the the result is continuous. 49 tmp = -cumsum (ppjumps (ppi), length (d) + 1);
50 50 ppi.coefs(prod(d)+1:end, k) = tmp(:);
51 jumps = reshape (ppjumps (ppi), prod (d), n-1); 51
52 ppi.P(:,2:n,k) -= cumsum (jumps, 2);
53
54 endfunction 52 endfunction
55 53
54 %!shared x,y,pp,ppi
55 %! x=0:8;y=[ones(size(x));x+1];pp=spline(x,y);
56 %! ppi=ppint(pp);
57 %!assert(ppval(ppi,x),[x;0.5*x.^2+x],1e-14)
58 %!assert(ppi.order,5)