Mercurial > hg > octave-lyh

annotate scripts/polynomial/spline.m @ 13141:e81ddf9cacd5

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
maint: untabify and remove trailing whitespace from source files * bicg.m, gmres.m, pkg.m: Untabify and remove trailing whitespace. * libcruft/Makefile.am, libcruft/blas-xtra/cdotc3.f, libcruft/blas-xtra/cmatm3.f, libcruft/blas-xtra/ddot3.f, libcruft/blas-xtra/dmatm3.f, libcruft/blas-xtra/sdot3.f, libcruft/blas-xtra/smatm3.f, libcruft/blas-xtra/zdotc3.f, libcruft/blas-xtra/zmatm3.f, libcruft/lapack-xtra/crsf2csf.f, libcruft/lapack-xtra/zrsf2csf.f, liboctave/Array.cc, liboctave/DASPK-opts.in, liboctave/DASRT-opts.in, liboctave/DASSL-opts.in, liboctave/LSODE-opts.in, liboctave/Makefile.a,mliboctave/Quad-opts.in, liboctave/Sparse-perm-op-defs.h, scripts/Makefile.a,mscripts/deprecated/glpkmex.m, scripts/general/blkdiag.m, scripts/general/interp1.m, scripts/general/profshow.m, scripts/general/quadl.m, scripts/general/triplequad.m, scripts/help/__makeinfo__.m, scripts/io/strread.m, scripts/io/textread.m, scripts/io/textscan.m, scripts/linear-algebra/rank.m, scripts/miscellaneous/gzip.m, scripts/miscellaneous/private/__xzip__.m, scripts/miscellaneous/tempdir.m, scripts/miscellaneous/unpack.m, scripts/pkg/pkg.m, scripts/plot/allchild.m, scripts/plot/ancestor.m, scripts/plot/cla.m, scripts/plot/clf.m, scripts/plot/findall.m, scripts/plot/findobj.m, scripts/plot/gca.m, scripts/plot/gcf.m, scripts/plot/hggroup.m, scripts/plot/isfigure.m, scripts/plot/ishghandle.m, scripts/plot/legend.m, scripts/plot/line.m, scripts/plot/loglog.m, scripts/plot/patch.m, scripts/plot/print.m, scripts/plot/private/__quiver__.m, scripts/plot/private/__scatter__.m, scripts/plot/rectangle.m, scripts/plot/semilogx.m, scripts/plot/semilogy.m, scripts/plot/surface.m, scripts/plot/text.m, scripts/plot/title.m, scripts/plot/trisurf.m, scripts/plot/view.m, scripts/plot/whitebg.m, scripts/plot/xlabel.m, scripts/plot/xlim.m, scripts/plot/ylabel.m, scripts/plot/ylim.m, scripts/plot/zlabel.m, scripts/plot/zlim.m, scripts/polynomial/mkpp.m, scripts/polynomial/polygcd.m, scripts/polynomial/ppint.m, scripts/polynomial/ppjumps.m, scripts/polynomial/ppval.m, scripts/set/setxor.m, scripts/sparse/bicgstab.m, scripts/sparse/cgs.m, scripts/sparse/spconvert.m, scripts/specfun/nthroot.m, scripts/strings/strmatch.m, scripts/strings/untabify.m, scripts/testfun/demo.m, scripts/testfun/example.m, src/DLD-FUNCTIONS/filter.cc, src/DLD-FUNCTIONS/mgorth.cc, src/DLD-FUNCTIONS/quadcc.cc, src/DLD-FUNCTIONS/str2double.cc, src/Makefile.a,msrc/gl-render.cc, src/gl2ps-renderer.cc, src/graphics.cc, src/octave-config.cc.in, src/octave-config.in, src/ov-class.h, src/ov-fcn.h, src/profiler.cc, src/profiler.h, src/pt-binop.cc, src/pt-unop.cc, src/symtab.cc, src/txt-eng-ft.cc: Remove trailing whitespace.
author John W. Eaton <jwe@octave.org>
date Thu, 15 Sep 2011 12:51:10 -0400
parents 59e2460acae1
children 614505385171
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
11523
fd0a3ac60b0e update copyright notices
John W. Eaton <jwe@octave.org>
parents: 11196
diff changeset
1 ## Copyright (C) 2000-2011 Kai Habel
5824
448f9982e7fb [project @ 2006-05-19 06:53:31 by jwe]
jwe
parents:
diff changeset
2 ## Copyright (C) 2006 David Bateman
448f9982e7fb [project @ 2006-05-19 06:53:31 by jwe]
jwe
parents:
diff changeset
3 ##
448f9982e7fb [project @ 2006-05-19 06:53:31 by jwe]
jwe
parents:
diff changeset
4 ## This file is part of Octave.
448f9982e7fb [project @ 2006-05-19 06:53:31 by jwe]
jwe
parents:
diff changeset
5 ##
448f9982e7fb [project @ 2006-05-19 06:53:31 by jwe]
jwe
parents:
diff changeset
6 ## Octave is free software; you can redistribute it and/or modify it
448f9982e7fb [project @ 2006-05-19 06:53:31 by jwe]
jwe
parents:
diff changeset
7 ## under the terms of the GNU General Public License as published by
7016
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe]
jwe
parents: 7001
diff changeset
8 ## the Free Software Foundation; either version 3 of the License, or (at
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe]
jwe
parents: 7001
diff changeset
9 ## your option) any later version.
5824
448f9982e7fb [project @ 2006-05-19 06:53:31 by jwe]
jwe
parents:
diff changeset
10 ##
448f9982e7fb [project @ 2006-05-19 06:53:31 by jwe]
jwe
parents:
diff changeset
11 ## Octave is distributed in the hope that it will be useful, but
448f9982e7fb [project @ 2006-05-19 06:53:31 by jwe]
jwe
parents:
diff changeset
12 ## WITHOUT ANY WARRANTY; without even the implied warranty of
448f9982e7fb [project @ 2006-05-19 06:53:31 by jwe]
jwe
parents:
diff changeset
13 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
448f9982e7fb [project @ 2006-05-19 06:53:31 by jwe]
jwe
parents:
diff changeset
14 ## General Public License for more details.
448f9982e7fb [project @ 2006-05-19 06:53:31 by jwe]
jwe
parents:
diff changeset
15 ##
448f9982e7fb [project @ 2006-05-19 06:53:31 by jwe]
jwe
parents:
diff changeset
16 ## You should have received a copy of the GNU General Public License
7016
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe]
jwe
parents: 7001
diff changeset
17 ## along with Octave; see the file COPYING. If not, see
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe]
jwe
parents: 7001
diff changeset
18 ## <http://www.gnu.org/licenses/>.
5824
448f9982e7fb [project @ 2006-05-19 06:53:31 by jwe]
jwe
parents:
diff changeset
19
448f9982e7fb [project @ 2006-05-19 06:53:31 by jwe]
jwe
parents:
diff changeset
20 ## -*- texinfo -*-
10793
be55736a0783 Grammarcheck the documentation from m-files.
Rik <octave@nomad.inbox5.com>
parents: 10658
diff changeset
21 ## @deftypefn {Function File} {@var{pp} =} spline (@var{x}, @var{y})
7650
eb7bdde776f2 Texinfo fixes
John W. Eaton <jwe@octave.org>
parents: 7017
diff changeset
22 ## @deftypefnx {Function File} {@var{yi} =} spline (@var{x}, @var{y}, @var{xi})
5824
448f9982e7fb [project @ 2006-05-19 06:53:31 by jwe]
jwe
parents:
diff changeset
23 ##
11587
c792872f8942 all script files: untabify and strip trailing whitespace
John W. Eaton <jwe@octave.org>
parents: 11536
diff changeset
24 ## Return the cubic spline interpolant of @var{y} at points @var{x}.
11536
702dbd0c53f5 Add undocumented ppder, ppint, ppjumps functions to documentation.
Rik <octave@nomad.inbox5.com>
parents: 11523
diff changeset
25 ## If called with two arguments, @code{spline} returns the piecewise
8602
827f0285a201 spline.m doc fix
carlo@guglielmo.local
parents: 8517
diff changeset
26 ## polynomial @var{pp} that may later be used with @code{ppval} to
827f0285a201 spline.m doc fix
carlo@guglielmo.local
parents: 8517
diff changeset
27 ## evaluate the polynomial at specific points.
11587
c792872f8942 all script files: untabify and strip trailing whitespace
John W. Eaton <jwe@octave.org>
parents: 11536
diff changeset
28 ## If called with a third input argument, @code{spline} evaluates the
9051
1bf0ce0930be Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents: 8920
diff changeset
29 ## spline at the points @var{xi}. There is an equivalence
8602
827f0285a201 spline.m doc fix
carlo@guglielmo.local
parents: 8517
diff changeset
30 ## between @code{ppval (spline (@var{x}, @var{y}), @var{xi})} and
827f0285a201 spline.m doc fix
carlo@guglielmo.local
parents: 8517
diff changeset
31 ## @code{spline (@var{x}, @var{y}, @var{xi})}.
5824
448f9982e7fb [project @ 2006-05-19 06:53:31 by jwe]
jwe
parents:
diff changeset
32 ##
448f9982e7fb [project @ 2006-05-19 06:53:31 by jwe]
jwe
parents:
diff changeset
33 ## The variable @var{x} must be a vector of length @var{n}, and @var{y}
9051
1bf0ce0930be Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents: 8920
diff changeset
34 ## can be either a vector or array. In the case where @var{y} is a
5824
448f9982e7fb [project @ 2006-05-19 06:53:31 by jwe]
jwe
parents:
diff changeset
35 ## vector, it can have a length of either @var{n} or @code{@var{n} + 2}.
448f9982e7fb [project @ 2006-05-19 06:53:31 by jwe]
jwe
parents:
diff changeset
36 ## If the length of @var{y} is @var{n}, then the 'not-a-knot' end
9051
1bf0ce0930be Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents: 8920
diff changeset
37 ## condition is used. If the length of @var{y} is @code{@var{n} + 2},
8602
827f0285a201 spline.m doc fix
carlo@guglielmo.local
parents: 8517
diff changeset
38 ## then the first and last values of the vector @var{y} are the values
827f0285a201 spline.m doc fix
carlo@guglielmo.local
parents: 8517
diff changeset
39 ## of the first derivative of the cubic spline at the end-points.
5824
448f9982e7fb [project @ 2006-05-19 06:53:31 by jwe]
jwe
parents:
diff changeset
40 ##
448f9982e7fb [project @ 2006-05-19 06:53:31 by jwe]
jwe
parents:
diff changeset
41 ## If @var{y} is an array, then the size of @var{y} must have the form
448f9982e7fb [project @ 2006-05-19 06:53:31 by jwe]
jwe
parents:
diff changeset
42 ## @tex
448f9982e7fb [project @ 2006-05-19 06:53:31 by jwe]
jwe
parents:
diff changeset
43 ## $$[s_1, s_2, \cdots, s_k, n]$$
448f9982e7fb [project @ 2006-05-19 06:53:31 by jwe]
jwe
parents:
diff changeset
44 ## @end tex
8517
81d6ab3ac93c Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents: 8310
diff changeset
45 ## @ifnottex
5824
448f9982e7fb [project @ 2006-05-19 06:53:31 by jwe]
jwe
parents:
diff changeset
46 ## @code{[@var{s1}, @var{s2}, @dots{}, @var{sk}, @var{n}]}
8517
81d6ab3ac93c Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents: 8310
diff changeset
47 ## @end ifnottex
5824
448f9982e7fb [project @ 2006-05-19 06:53:31 by jwe]
jwe
parents:
diff changeset
48 ## or
448f9982e7fb [project @ 2006-05-19 06:53:31 by jwe]
jwe
parents:
diff changeset
49 ## @tex
8828
8463d1a2e544 Doc fixes.
Brian Gough <bjg@network-theory.co.uk>
parents: 8602
diff changeset
50 ## $$[s_1, s_2, \cdots, s_k, n + 2].$$
5824
448f9982e7fb [project @ 2006-05-19 06:53:31 by jwe]
jwe
parents:
diff changeset
51 ## @end tex
8517
81d6ab3ac93c Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents: 8310
diff changeset
52 ## @ifnottex
5824
448f9982e7fb [project @ 2006-05-19 06:53:31 by jwe]
jwe
parents:
diff changeset
53 ## @code{[@var{s1}, @var{s2}, @dots{}, @var{sk}, @var{n} + 2]}.
8517
81d6ab3ac93c Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents: 8310
diff changeset
54 ## @end ifnottex
8602
827f0285a201 spline.m doc fix
carlo@guglielmo.local
parents: 8517
diff changeset
55 ## The array is then reshaped internally to a matrix where the leading
11587
c792872f8942 all script files: untabify and strip trailing whitespace
John W. Eaton <jwe@octave.org>
parents: 11536
diff changeset
56 ## dimension is given by
5824
448f9982e7fb [project @ 2006-05-19 06:53:31 by jwe]
jwe
parents:
diff changeset
57 ## @tex
448f9982e7fb [project @ 2006-05-19 06:53:31 by jwe]
jwe
parents:
diff changeset
58 ## $$s_1 s_2 \cdots s_k$$
448f9982e7fb [project @ 2006-05-19 06:53:31 by jwe]
jwe
parents:
diff changeset
59 ## @end tex
8517
81d6ab3ac93c Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents: 8310
diff changeset
60 ## @ifnottex
5824
448f9982e7fb [project @ 2006-05-19 06:53:31 by jwe]
jwe
parents:
diff changeset
61 ## @code{@var{s1} * @var{s2} * @dots{} * @var{sk}}
8517
81d6ab3ac93c Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents: 8310
diff changeset
62 ## @end ifnottex
9051
1bf0ce0930be Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents: 8920
diff changeset
63 ## and each row of this matrix is then treated separately. Note that this
5824
448f9982e7fb [project @ 2006-05-19 06:53:31 by jwe]
jwe
parents:
diff changeset
64 ## is exactly the opposite treatment than @code{interp1} and is done
7001
8b0cfeb06365 [project @ 2007-10-10 18:02:59 by jwe]
jwe
parents: 6702
diff changeset
65 ## for compatibility.
5824
448f9982e7fb [project @ 2006-05-19 06:53:31 by jwe]
jwe
parents:
diff changeset
66 ## @seealso{ppval, mkpp, unmkpp}
448f9982e7fb [project @ 2006-05-19 06:53:31 by jwe]
jwe
parents:
diff changeset
67 ## @end deftypefn
448f9982e7fb [project @ 2006-05-19 06:53:31 by jwe]
jwe
parents:
diff changeset
68
448f9982e7fb [project @ 2006-05-19 06:53:31 by jwe]
jwe
parents:
diff changeset
69 ## This code is based on csape.m from octave-forge, but has been
448f9982e7fb [project @ 2006-05-19 06:53:31 by jwe]
jwe
parents:
diff changeset
70 ## modified to use the sparse solver code in octave that itself allows
448f9982e7fb [project @ 2006-05-19 06:53:31 by jwe]
jwe
parents:
diff changeset
71 ## special casing of tri-diagonal matrices, modified for NDArrays and
448f9982e7fb [project @ 2006-05-19 06:53:31 by jwe]
jwe
parents:
diff changeset
72 ## for the treatment of vectors y 2 elements longer than x as complete
448f9982e7fb [project @ 2006-05-19 06:53:31 by jwe]
jwe
parents:
diff changeset
73 ## splines.
448f9982e7fb [project @ 2006-05-19 06:53:31 by jwe]
jwe
parents:
diff changeset
74
448f9982e7fb [project @ 2006-05-19 06:53:31 by jwe]
jwe
parents:
diff changeset
75 function ret = spline (x, y, xi)
448f9982e7fb [project @ 2006-05-19 06:53:31 by jwe]
jwe
parents:
diff changeset
76
448f9982e7fb [project @ 2006-05-19 06:53:31 by jwe]
jwe
parents:
diff changeset
77 x = x(:);
448f9982e7fb [project @ 2006-05-19 06:53:31 by jwe]
jwe
parents:
diff changeset
78 n = length (x);
12444
07e102029d2a spline.m: Allow length(x) == 2 and unsorted x values.
Ben Abbott <bpabbott@mac.com>
parents: 11587
diff changeset
79 if (n < 2)
07e102029d2a spline.m: Allow length(x) == 2 and unsorted x values.
Ben Abbott <bpabbott@mac.com>
parents: 11587
diff changeset
80 error ("spline: requires at least 2 points");
5824
448f9982e7fb [project @ 2006-05-19 06:53:31 by jwe]
jwe
parents:
diff changeset
81 endif
448f9982e7fb [project @ 2006-05-19 06:53:31 by jwe]
jwe
parents:
diff changeset
82
448f9982e7fb [project @ 2006-05-19 06:53:31 by jwe]
jwe
parents:
diff changeset
83 ## Check the size and shape of y
448f9982e7fb [project @ 2006-05-19 06:53:31 by jwe]
jwe
parents:
diff changeset
84 ndy = ndims (y);
448f9982e7fb [project @ 2006-05-19 06:53:31 by jwe]
jwe
parents:
diff changeset
85 szy = size (y);
12608
59e2460acae1 make piecewise polynomial (pp) functions more compatible
Kai Habel <kai.habel@gmx.de>
parents: 12444
diff changeset
86 if (ndy == 2 && (szy(1) == n || szy(2) == n))
59e2460acae1 make piecewise polynomial (pp) functions more compatible
Kai Habel <kai.habel@gmx.de>
parents: 12444
diff changeset
87 if (szy(2) == n)
6014
437f9086b967 [project @ 2006-09-30 15:27:06 by jwe]
jwe
parents: 5838
diff changeset
88 a = y.';
5824
448f9982e7fb [project @ 2006-05-19 06:53:31 by jwe]
jwe
parents:
diff changeset
89 else
448f9982e7fb [project @ 2006-05-19 06:53:31 by jwe]
jwe
parents:
diff changeset
90 a = y;
448f9982e7fb [project @ 2006-05-19 06:53:31 by jwe]
jwe
parents:
diff changeset
91 szy = fliplr (szy);
448f9982e7fb [project @ 2006-05-19 06:53:31 by jwe]
jwe
parents:
diff changeset
92 endif
448f9982e7fb [project @ 2006-05-19 06:53:31 by jwe]
jwe
parents:
diff changeset
93 else
12608
59e2460acae1 make piecewise polynomial (pp) functions more compatible
Kai Habel <kai.habel@gmx.de>
parents: 12444
diff changeset
94 a = shiftdim (reshape (y, [prod(szy(1:end-1)), szy(end)]), 1);
5824
448f9982e7fb [project @ 2006-05-19 06:53:31 by jwe]
jwe
parents:
diff changeset
95 endif
11587
c792872f8942 all script files: untabify and strip trailing whitespace
John W. Eaton <jwe@octave.org>
parents: 11536
diff changeset
96
c792872f8942 all script files: untabify and strip trailing whitespace
John W. Eaton <jwe@octave.org>
parents: 11536
diff changeset
97 for k = (1:columns (a))(any (isnan (a)))
c792872f8942 all script files: untabify and strip trailing whitespace
John W. Eaton <jwe@octave.org>
parents: 11536
diff changeset
98 ok = ! isnan (a(:,k));
c792872f8942 all script files: untabify and strip trailing whitespace
John W. Eaton <jwe@octave.org>
parents: 11536
diff changeset
99 a(!ok,k) = spline (x(ok), a(ok,k), x(!ok));
c792872f8942 all script files: untabify and strip trailing whitespace
John W. Eaton <jwe@octave.org>
parents: 11536
diff changeset
100 endfor
c792872f8942 all script files: untabify and strip trailing whitespace
John W. Eaton <jwe@octave.org>
parents: 11536
diff changeset
101
5824
448f9982e7fb [project @ 2006-05-19 06:53:31 by jwe]
jwe
parents:
diff changeset
102 complete = false;
448f9982e7fb [project @ 2006-05-19 06:53:31 by jwe]
jwe
parents:
diff changeset
103 if (size (a, 1) == n + 2)
448f9982e7fb [project @ 2006-05-19 06:53:31 by jwe]
jwe
parents:
diff changeset
104 complete = true;
448f9982e7fb [project @ 2006-05-19 06:53:31 by jwe]
jwe
parents:
diff changeset
105 dfs = a(1,:);
448f9982e7fb [project @ 2006-05-19 06:53:31 by jwe]
jwe
parents:
diff changeset
106 dfe = a(end,:);
448f9982e7fb [project @ 2006-05-19 06:53:31 by jwe]
jwe
parents:
diff changeset
107 a = a(2:end-1,:);
448f9982e7fb [project @ 2006-05-19 06:53:31 by jwe]
jwe
parents:
diff changeset
108 endif
448f9982e7fb [project @ 2006-05-19 06:53:31 by jwe]
jwe
parents:
diff changeset
109
12444
07e102029d2a spline.m: Allow length(x) == 2 and unsorted x values.
Ben Abbott <bpabbott@mac.com>
parents: 11587
diff changeset
110 if (~issorted (x))
07e102029d2a spline.m: Allow length(x) == 2 and unsorted x values.
Ben Abbott <bpabbott@mac.com>
parents: 11587
diff changeset
111 [x, idx] = sort(x);
07e102029d2a spline.m: Allow length(x) == 2 and unsorted x values.
Ben Abbott <bpabbott@mac.com>
parents: 11587
diff changeset
112 a = a(idx,:);
07e102029d2a spline.m: Allow length(x) == 2 and unsorted x values.
Ben Abbott <bpabbott@mac.com>
parents: 11587
diff changeset
113 endif
07e102029d2a spline.m: Allow length(x) == 2 and unsorted x values.
Ben Abbott <bpabbott@mac.com>
parents: 11587
diff changeset
114
5824
448f9982e7fb [project @ 2006-05-19 06:53:31 by jwe]
jwe
parents:
diff changeset
115 b = c = zeros (size (a));
448f9982e7fb [project @ 2006-05-19 06:53:31 by jwe]
jwe
parents:
diff changeset
116 h = diff (x);
5838
376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe]
jwe
parents: 5837
diff changeset
117 idx = ones (columns (a), 1);
5824
448f9982e7fb [project @ 2006-05-19 06:53:31 by jwe]
jwe
parents:
diff changeset
118
448f9982e7fb [project @ 2006-05-19 06:53:31 by jwe]
jwe
parents:
diff changeset
119 if (complete)
448f9982e7fb [project @ 2006-05-19 06:53:31 by jwe]
jwe
parents:
diff changeset
120
12444
07e102029d2a spline.m: Allow length(x) == 2 and unsorted x values.
Ben Abbott <bpabbott@mac.com>
parents: 11587
diff changeset
121 if (n == 2)
07e102029d2a spline.m: Allow length(x) == 2 and unsorted x values.
Ben Abbott <bpabbott@mac.com>
parents: 11587
diff changeset
122 d = (dfs + dfe) / (x(2) - x(1)) ^ 2 + ...
12608
59e2460acae1 make piecewise polynomial (pp) functions more compatible
Kai Habel <kai.habel@gmx.de>
parents: 12444
diff changeset
123 2 * (a(1,:) - a(2,:)) / (x(2) - x(1)) ^ 3;
12444
07e102029d2a spline.m: Allow length(x) == 2 and unsorted x values.
Ben Abbott <bpabbott@mac.com>
parents: 11587
diff changeset
124 c = (-2 * dfs - dfe) / (x(2) - x(1)) - ...
12608
59e2460acae1 make piecewise polynomial (pp) functions more compatible
Kai Habel <kai.habel@gmx.de>
parents: 12444
diff changeset
125 3 * (a(1,:) - a(2,:)) / (x(2) - x(1)) ^ 2;
12444
07e102029d2a spline.m: Allow length(x) == 2 and unsorted x values.
Ben Abbott <bpabbott@mac.com>
parents: 11587
diff changeset
126 b = dfs;
07e102029d2a spline.m: Allow length(x) == 2 and unsorted x values.
Ben Abbott <bpabbott@mac.com>
parents: 11587
diff changeset
127 a = a(1,:);
5824
448f9982e7fb [project @ 2006-05-19 06:53:31 by jwe]
jwe
parents:
diff changeset
128
12444
07e102029d2a spline.m: Allow length(x) == 2 and unsorted x values.
Ben Abbott <bpabbott@mac.com>
parents: 11587
diff changeset
129 d = d(1:n-1,:);
07e102029d2a spline.m: Allow length(x) == 2 and unsorted x values.
Ben Abbott <bpabbott@mac.com>
parents: 11587
diff changeset
130 c = c(1:n-1,:);
07e102029d2a spline.m: Allow length(x) == 2 and unsorted x values.
Ben Abbott <bpabbott@mac.com>
parents: 11587
diff changeset
131 b = b(1:n-1,:);
07e102029d2a spline.m: Allow length(x) == 2 and unsorted x values.
Ben Abbott <bpabbott@mac.com>
parents: 11587
diff changeset
132 a = a(1:n-1,:);
07e102029d2a spline.m: Allow length(x) == 2 and unsorted x values.
Ben Abbott <bpabbott@mac.com>
parents: 11587
diff changeset
133 else
07e102029d2a spline.m: Allow length(x) == 2 and unsorted x values.
Ben Abbott <bpabbott@mac.com>
parents: 11587
diff changeset
134 if (n == 3)
12608
59e2460acae1 make piecewise polynomial (pp) functions more compatible
Kai Habel <kai.habel@gmx.de>
parents: 12444
diff changeset
135 dg = 1.5 * h(1) - 0.5 * h(2);
12444
07e102029d2a spline.m: Allow length(x) == 2 and unsorted x values.
Ben Abbott <bpabbott@mac.com>
parents: 11587
diff changeset
136 c(2:n-1,:) = 1/dg(1);
07e102029d2a spline.m: Allow length(x) == 2 and unsorted x values.
Ben Abbott <bpabbott@mac.com>
parents: 11587
diff changeset
137 else
07e102029d2a spline.m: Allow length(x) == 2 and unsorted x values.
Ben Abbott <bpabbott@mac.com>
parents: 11587
diff changeset
138 dg = 2 * (h(1:n-2) .+ h(2:n-1));
07e102029d2a spline.m: Allow length(x) == 2 and unsorted x values.
Ben Abbott <bpabbott@mac.com>
parents: 11587
diff changeset
139 dg(1) = dg(1) - 0.5 * h(1);
07e102029d2a spline.m: Allow length(x) == 2 and unsorted x values.
Ben Abbott <bpabbott@mac.com>
parents: 11587
diff changeset
140 dg(n-2) = dg(n-2) - 0.5 * h(n-1);
5824
448f9982e7fb [project @ 2006-05-19 06:53:31 by jwe]
jwe
parents:
diff changeset
141
12444
07e102029d2a spline.m: Allow length(x) == 2 and unsorted x values.
Ben Abbott <bpabbott@mac.com>
parents: 11587
diff changeset
142 e = h(2:n-2);
07e102029d2a spline.m: Allow length(x) == 2 and unsorted x values.
Ben Abbott <bpabbott@mac.com>
parents: 11587
diff changeset
143
07e102029d2a spline.m: Allow length(x) == 2 and unsorted x values.
Ben Abbott <bpabbott@mac.com>
parents: 11587
diff changeset
144 g = 3 * diff (a(2:n,:)) ./ h(2:n-1,idx) ...
5838
376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe]
jwe
parents: 5837
diff changeset
145 - 3 * diff (a(1:n-1,:)) ./ h(1:n-2,idx);
12444
07e102029d2a spline.m: Allow length(x) == 2 and unsorted x values.
Ben Abbott <bpabbott@mac.com>
parents: 11587
diff changeset
146 g(1,:) = 3 * (a(3,:) - a(2,:)) / h(2) ...
5824
448f9982e7fb [project @ 2006-05-19 06:53:31 by jwe]
jwe
parents:
diff changeset
147 - 3 / 2 * (3 * (a(2,:) - a(1,:)) / h(1) - dfs);
12444
07e102029d2a spline.m: Allow length(x) == 2 and unsorted x values.
Ben Abbott <bpabbott@mac.com>
parents: 11587
diff changeset
148 g(n-2,:) = 3 / 2 * (3 * (a(n,:) - a(n-1,:)) / h(n-1) - dfe) ...
5838
376e02b2ce70 [project @ 2006-06-01 20:23:53 by jwe]
jwe
parents: 5837
diff changeset
149 - 3 * (a(n-1,:) - a(n-2,:)) / h(n-2);
5824
448f9982e7fb [project @ 2006-05-19 06:53:31 by jwe]
jwe
parents:
diff changeset
150
12444
07e102029d2a spline.m: Allow length(x) == 2 and unsorted x values.
Ben Abbott <bpabbott@mac.com>
parents: 11587
diff changeset
151 c(2:n-1,:) = spdiags ([[e(:); 0], dg, [0; e(:)]],
10549
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9768
diff changeset
152 [-1, 0, 1], n-2, n-2) \ g;
12444
07e102029d2a spline.m: Allow length(x) == 2 and unsorted x values.
Ben Abbott <bpabbott@mac.com>
parents: 11587
diff changeset
153 endif
5824
448f9982e7fb [project @ 2006-05-19 06:53:31 by jwe]
jwe
parents:
diff changeset
154
12444
07e102029d2a spline.m: Allow length(x) == 2 and unsorted x values.
Ben Abbott <bpabbott@mac.com>
parents: 11587
diff changeset
155 c(1,:) = (3 / h(1) * (a(2,:) - a(1,:)) - 3 * dfs
12608
59e2460acae1 make piecewise polynomial (pp) functions more compatible
Kai Habel <kai.habel@gmx.de>
parents: 12444
diff changeset
156 - c(2,:) * h(1)) / (2 * h(1));
12444
07e102029d2a spline.m: Allow length(x) == 2 and unsorted x values.
Ben Abbott <bpabbott@mac.com>
parents: 11587
diff changeset
157 c(n,:) = - (3 / h(n-1) * (a(n,:) - a(n-1,:)) - 3 * dfe
12608
59e2460acae1 make piecewise polynomial (pp) functions more compatible
Kai Habel <kai.habel@gmx.de>
parents: 12444
diff changeset
158 + c(n-1,:) * h(n-1)) / (2 * h(n-1));
12444
07e102029d2a spline.m: Allow length(x) == 2 and unsorted x values.
Ben Abbott <bpabbott@mac.com>
parents: 11587
diff changeset
159 b(1:n-1,:) = diff (a) ./ h(1:n-1, idx) ...
10549
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9768
diff changeset
160 - h(1:n-1,idx) / 3 .* (c(2:n,:) + 2 * c(1:n-1,:));
12444
07e102029d2a spline.m: Allow length(x) == 2 and unsorted x values.
Ben Abbott <bpabbott@mac.com>
parents: 11587
diff changeset
161 d = diff (c) ./ (3 * h(1:n-1, idx));
5824
448f9982e7fb [project @ 2006-05-19 06:53:31 by jwe]
jwe
parents:
diff changeset
162
12444
07e102029d2a spline.m: Allow length(x) == 2 and unsorted x values.
Ben Abbott <bpabbott@mac.com>
parents: 11587
diff changeset
163 d = d(1:n-1,:);
07e102029d2a spline.m: Allow length(x) == 2 and unsorted x values.
Ben Abbott <bpabbott@mac.com>
parents: 11587
diff changeset
164 c = c(1:n-1,:);
07e102029d2a spline.m: Allow length(x) == 2 and unsorted x values.
Ben Abbott <bpabbott@mac.com>
parents: 11587
diff changeset
165 b = b(1:n-1,:);
07e102029d2a spline.m: Allow length(x) == 2 and unsorted x values.
Ben Abbott <bpabbott@mac.com>
parents: 11587
diff changeset
166 a = a(1:n-1,:);
07e102029d2a spline.m: Allow length(x) == 2 and unsorted x values.
Ben Abbott <bpabbott@mac.com>
parents: 11587
diff changeset
167 endif
5824
448f9982e7fb [project @ 2006-05-19 06:53:31 by jwe]
jwe
parents:
diff changeset
168 else
448f9982e7fb [project @ 2006-05-19 06:53:31 by jwe]
jwe
parents:
diff changeset
169
12444
07e102029d2a spline.m: Allow length(x) == 2 and unsorted x values.
Ben Abbott <bpabbott@mac.com>
parents: 11587
diff changeset
170 if (n == 2)
07e102029d2a spline.m: Allow length(x) == 2 and unsorted x values.
Ben Abbott <bpabbott@mac.com>
parents: 11587
diff changeset
171 b = (a(2,:) - a(1,:)) / (x(2) - x(1));
07e102029d2a spline.m: Allow length(x) == 2 and unsorted x values.
Ben Abbott <bpabbott@mac.com>
parents: 11587
diff changeset
172 a = a(1,:);
07e102029d2a spline.m: Allow length(x) == 2 and unsorted x values.
Ben Abbott <bpabbott@mac.com>
parents: 11587
diff changeset
173 d = [];
07e102029d2a spline.m: Allow length(x) == 2 and unsorted x values.
Ben Abbott <bpabbott@mac.com>
parents: 11587
diff changeset
174 c = [];
07e102029d2a spline.m: Allow length(x) == 2 and unsorted x values.
Ben Abbott <bpabbott@mac.com>
parents: 11587
diff changeset
175 b = b(1:n-1,:);
07e102029d2a spline.m: Allow length(x) == 2 and unsorted x values.
Ben Abbott <bpabbott@mac.com>
parents: 11587
diff changeset
176 a = a(1:n-1,:);
07e102029d2a spline.m: Allow length(x) == 2 and unsorted x values.
Ben Abbott <bpabbott@mac.com>
parents: 11587
diff changeset
177 elseif (n == 3)
5824
448f9982e7fb [project @ 2006-05-19 06:53:31 by jwe]
jwe
parents:
diff changeset
178
11196
d17cb8a1271d spline.m: Fix algorithm for input with 3 elements (bug #31098)
Marco Caliari <marco.caliari@univr.it>
parents: 10793
diff changeset
179 n = 2;
d17cb8a1271d spline.m: Fix algorithm for input with 3 elements (bug #31098)
Marco Caliari <marco.caliari@univr.it>
parents: 10793
diff changeset
180 c = (a(1,:) - a(3,:)) / ((x(3) - x(1)) * (x(2) - x(3))) ...
d17cb8a1271d spline.m: Fix algorithm for input with 3 elements (bug #31098)
Marco Caliari <marco.caliari@univr.it>
parents: 10793
diff changeset
181 + (a(2,:) - a(1,:)) / ((x(2) - x(1)) * (x(2) - x(3)));
d17cb8a1271d spline.m: Fix algorithm for input with 3 elements (bug #31098)
Marco Caliari <marco.caliari@univr.it>
parents: 10793
diff changeset
182 b = (a(2,:) - a(1,:)) * (x(3) - x(1)) ...
d17cb8a1271d spline.m: Fix algorithm for input with 3 elements (bug #31098)
Marco Caliari <marco.caliari@univr.it>
parents: 10793
diff changeset
183 / ((x(2) - x(1)) * (x(3) - x(2))) ...
d17cb8a1271d spline.m: Fix algorithm for input with 3 elements (bug #31098)
Marco Caliari <marco.caliari@univr.it>
parents: 10793
diff changeset
184 + (a(1,:) - a(3,:)) * (x(2) - x(1)) ...
d17cb8a1271d spline.m: Fix algorithm for input with 3 elements (bug #31098)
Marco Caliari <marco.caliari@univr.it>
parents: 10793
diff changeset
185 / ((x(3) - x(1)) * (x(3) - x(2)));
d17cb8a1271d spline.m: Fix algorithm for input with 3 elements (bug #31098)
Marco Caliari <marco.caliari@univr.it>
parents: 10793
diff changeset
186 a = a(1,:);
12444
07e102029d2a spline.m: Allow length(x) == 2 and unsorted x values.
Ben Abbott <bpabbott@mac.com>
parents: 11587
diff changeset
187 d = [];
11196
d17cb8a1271d spline.m: Fix algorithm for input with 3 elements (bug #31098)
Marco Caliari <marco.caliari@univr.it>
parents: 10793
diff changeset
188 x = [min(x), max(x)];
d17cb8a1271d spline.m: Fix algorithm for input with 3 elements (bug #31098)
Marco Caliari <marco.caliari@univr.it>
parents: 10793
diff changeset
189
12444
07e102029d2a spline.m: Allow length(x) == 2 and unsorted x values.
Ben Abbott <bpabbott@mac.com>
parents: 11587
diff changeset
190 c = c(1:n-1,:);
07e102029d2a spline.m: Allow length(x) == 2 and unsorted x values.
Ben Abbott <bpabbott@mac.com>
parents: 11587
diff changeset
191 b = b(1:n-1,:);
07e102029d2a spline.m: Allow length(x) == 2 and unsorted x values.
Ben Abbott <bpabbott@mac.com>
parents: 11587
diff changeset
192 a = a(1:n-1,:);
11196
d17cb8a1271d spline.m: Fix algorithm for input with 3 elements (bug #31098)
Marco Caliari <marco.caliari@univr.it>
parents: 10793
diff changeset
193 else
5824
448f9982e7fb [project @ 2006-05-19 06:53:31 by jwe]
jwe
parents:
diff changeset
194
11196
d17cb8a1271d spline.m: Fix algorithm for input with 3 elements (bug #31098)
Marco Caliari <marco.caliari@univr.it>
parents: 10793
diff changeset
195 g = zeros (n-2, columns (a));
d17cb8a1271d spline.m: Fix algorithm for input with 3 elements (bug #31098)
Marco Caliari <marco.caliari@univr.it>
parents: 10793
diff changeset
196 g(1,:) = 3 / (h(1) + h(2)) ...
d17cb8a1271d spline.m: Fix algorithm for input with 3 elements (bug #31098)
Marco Caliari <marco.caliari@univr.it>
parents: 10793
diff changeset
197 * (a(3,:) - a(2,:) - h(2) / h(1) * (a(2,:) - a(1,:)));
d17cb8a1271d spline.m: Fix algorithm for input with 3 elements (bug #31098)
Marco Caliari <marco.caliari@univr.it>
parents: 10793
diff changeset
198 g(n-2,:) = 3 / (h(n-1) + h(n-2)) ...
d17cb8a1271d spline.m: Fix algorithm for input with 3 elements (bug #31098)
Marco Caliari <marco.caliari@univr.it>
parents: 10793
diff changeset
199 * (h(n-2) / h(n-1) * (a(n,:) - a(n-1,:)) - (a(n-1,:) - a(n-2,:)));
d17cb8a1271d spline.m: Fix algorithm for input with 3 elements (bug #31098)
Marco Caliari <marco.caliari@univr.it>
parents: 10793
diff changeset
200
d17cb8a1271d spline.m: Fix algorithm for input with 3 elements (bug #31098)
Marco Caliari <marco.caliari@univr.it>
parents: 10793
diff changeset
201 if (n > 4)
5824
448f9982e7fb [project @ 2006-05-19 06:53:31 by jwe]
jwe
parents:
diff changeset
202
11196
d17cb8a1271d spline.m: Fix algorithm for input with 3 elements (bug #31098)
Marco Caliari <marco.caliari@univr.it>
parents: 10793
diff changeset
203 g(2:n - 3,:) = 3 * diff (a(3:n-1,:)) ./ h(3:n-2,idx) ...
d17cb8a1271d spline.m: Fix algorithm for input with 3 elements (bug #31098)
Marco Caliari <marco.caliari@univr.it>
parents: 10793
diff changeset
204 - 3 * diff (a(2:n-2,:)) ./ h(2:n - 3,idx);
5824
448f9982e7fb [project @ 2006-05-19 06:53:31 by jwe]
jwe
parents:
diff changeset
205
11196
d17cb8a1271d spline.m: Fix algorithm for input with 3 elements (bug #31098)
Marco Caliari <marco.caliari@univr.it>
parents: 10793
diff changeset
206 dg = 2 * (h(1:n-2) .+ h(2:n-1));
d17cb8a1271d spline.m: Fix algorithm for input with 3 elements (bug #31098)
Marco Caliari <marco.caliari@univr.it>
parents: 10793
diff changeset
207 dg(1) = dg(1) - h(1);
d17cb8a1271d spline.m: Fix algorithm for input with 3 elements (bug #31098)
Marco Caliari <marco.caliari@univr.it>
parents: 10793
diff changeset
208 dg(n-2) = dg(n-2) - h(n-1);
d17cb8a1271d spline.m: Fix algorithm for input with 3 elements (bug #31098)
Marco Caliari <marco.caliari@univr.it>
parents: 10793
diff changeset
209
d17cb8a1271d spline.m: Fix algorithm for input with 3 elements (bug #31098)
Marco Caliari <marco.caliari@univr.it>
parents: 10793
diff changeset
210 ldg = udg = h(2:n-2);
d17cb8a1271d spline.m: Fix algorithm for input with 3 elements (bug #31098)
Marco Caliari <marco.caliari@univr.it>
parents: 10793
diff changeset
211 udg(1) = udg(1) - h(1);
d17cb8a1271d spline.m: Fix algorithm for input with 3 elements (bug #31098)
Marco Caliari <marco.caliari@univr.it>
parents: 10793
diff changeset
212 ldg(n - 3) = ldg(n-3) - h(n-1);
d17cb8a1271d spline.m: Fix algorithm for input with 3 elements (bug #31098)
Marco Caliari <marco.caliari@univr.it>
parents: 10793
diff changeset
213 c(2:n-1,:) = spdiags ([[ldg(:); 0], dg, [0; udg(:)]],
10549
95c3e38098bf Untabify .m scripts
Rik <code@nomad.inbox5.com>
parents: 9768
diff changeset
214 [-1, 0, 1], n-2, n-2) \ g;
5824
448f9982e7fb [project @ 2006-05-19 06:53:31 by jwe]
jwe
parents:
diff changeset
215
11196
d17cb8a1271d spline.m: Fix algorithm for input with 3 elements (bug #31098)
Marco Caliari <marco.caliari@univr.it>
parents: 10793
diff changeset
216 elseif (n == 4)
d17cb8a1271d spline.m: Fix algorithm for input with 3 elements (bug #31098)
Marco Caliari <marco.caliari@univr.it>
parents: 10793
diff changeset
217
d17cb8a1271d spline.m: Fix algorithm for input with 3 elements (bug #31098)
Marco Caliari <marco.caliari@univr.it>
parents: 10793
diff changeset
218 dg = [h(1) + 2 * h(2); 2 * h(2) + h(3)];
d17cb8a1271d spline.m: Fix algorithm for input with 3 elements (bug #31098)
Marco Caliari <marco.caliari@univr.it>
parents: 10793
diff changeset
219 ldg = h(2) - h(3);
d17cb8a1271d spline.m: Fix algorithm for input with 3 elements (bug #31098)
Marco Caliari <marco.caliari@univr.it>
parents: 10793
diff changeset
220 udg = h(2) - h(1);
d17cb8a1271d spline.m: Fix algorithm for input with 3 elements (bug #31098)
Marco Caliari <marco.caliari@univr.it>
parents: 10793
diff changeset
221 c(2:n-1,:) = spdiags ([[ldg(:);0], dg, [0; udg(:)]],
d17cb8a1271d spline.m: Fix algorithm for input with 3 elements (bug #31098)
Marco Caliari <marco.caliari@univr.it>
parents: 10793
diff changeset
222 [-1, 0, 1], n-2, n-2) \ g;
5824
448f9982e7fb [project @ 2006-05-19 06:53:31 by jwe]
jwe
parents:
diff changeset
223
11196
d17cb8a1271d spline.m: Fix algorithm for input with 3 elements (bug #31098)
Marco Caliari <marco.caliari@univr.it>
parents: 10793
diff changeset
224 endif
d17cb8a1271d spline.m: Fix algorithm for input with 3 elements (bug #31098)
Marco Caliari <marco.caliari@univr.it>
parents: 10793
diff changeset
225
d17cb8a1271d spline.m: Fix algorithm for input with 3 elements (bug #31098)
Marco Caliari <marco.caliari@univr.it>
parents: 10793
diff changeset
226 c(1,:) = c(2,:) + h(1) / h(2) * (c(2,:) - c(3,:));
d17cb8a1271d spline.m: Fix algorithm for input with 3 elements (bug #31098)
Marco Caliari <marco.caliari@univr.it>
parents: 10793
diff changeset
227 c(n,:) = c(n-1,:) + h(n-1) / h(n-2) * (c(n-1,:) - c(n-2,:));
d17cb8a1271d spline.m: Fix algorithm for input with 3 elements (bug #31098)
Marco Caliari <marco.caliari@univr.it>
parents: 10793
diff changeset
228 b = diff (a) ./ h(1:n-1, idx) ...
d17cb8a1271d spline.m: Fix algorithm for input with 3 elements (bug #31098)
Marco Caliari <marco.caliari@univr.it>
parents: 10793
diff changeset
229 - h(1:n-1, idx) / 3 .* (c(2:n,:) + 2 * c(1:n-1,:));
d17cb8a1271d spline.m: Fix algorithm for input with 3 elements (bug #31098)
Marco Caliari <marco.caliari@univr.it>
parents: 10793
diff changeset
230 d = diff (c) ./ (3 * h(1:n-1, idx));
5824
448f9982e7fb [project @ 2006-05-19 06:53:31 by jwe]
jwe
parents:
diff changeset
231
12608
59e2460acae1 make piecewise polynomial (pp) functions more compatible
Kai Habel <kai.habel@gmx.de>
parents: 12444
diff changeset
232 d = d(1:n-1,:);d = d.'(:);
59e2460acae1 make piecewise polynomial (pp) functions more compatible
Kai Habel <kai.habel@gmx.de>
parents: 12444
diff changeset
233 c = c(1:n-1,:);c = c.'(:);
59e2460acae1 make piecewise polynomial (pp) functions more compatible
Kai Habel <kai.habel@gmx.de>
parents: 12444
diff changeset
234 b = b(1:n-1,:);b = b.'(:);
59e2460acae1 make piecewise polynomial (pp) functions more compatible
Kai Habel <kai.habel@gmx.de>
parents: 12444
diff changeset
235 a = a(1:n-1,:);a = a.'(:);
5824
448f9982e7fb [project @ 2006-05-19 06:53:31 by jwe]
jwe
parents:
diff changeset
236 endif
448f9982e7fb [project @ 2006-05-19 06:53:31 by jwe]
jwe
parents:
diff changeset
237
448f9982e7fb [project @ 2006-05-19 06:53:31 by jwe]
jwe
parents:
diff changeset
238 endif
12608
59e2460acae1 make piecewise polynomial (pp) functions more compatible
Kai Habel <kai.habel@gmx.de>
parents: 12444
diff changeset
239 ret = mkpp (x, cat (2, d, c, b, a), szy(1:end-1));
5824
448f9982e7fb [project @ 2006-05-19 06:53:31 by jwe]
jwe
parents:
diff changeset
240
448f9982e7fb [project @ 2006-05-19 06:53:31 by jwe]
jwe
parents:
diff changeset
241 if (nargin == 3)
448f9982e7fb [project @ 2006-05-19 06:53:31 by jwe]
jwe
parents:
diff changeset
242 ret = ppval (ret, xi);
448f9982e7fb [project @ 2006-05-19 06:53:31 by jwe]
jwe
parents:
diff changeset
243 endif
448f9982e7fb [project @ 2006-05-19 06:53:31 by jwe]
jwe
parents:
diff changeset
244
448f9982e7fb [project @ 2006-05-19 06:53:31 by jwe]
jwe
parents:
diff changeset
245 endfunction
448f9982e7fb [project @ 2006-05-19 06:53:31 by jwe]
jwe
parents:
diff changeset
246
448f9982e7fb [project @ 2006-05-19 06:53:31 by jwe]
jwe
parents:
diff changeset
247 %!demo
448f9982e7fb [project @ 2006-05-19 06:53:31 by jwe]
jwe
parents:
diff changeset
248 %! x = 0:10; y = sin(x);
448f9982e7fb [project @ 2006-05-19 06:53:31 by jwe]
jwe
parents:
diff changeset
249 %! xspline = 0:0.1:10; yspline = spline(x,y,xspline);
448f9982e7fb [project @ 2006-05-19 06:53:31 by jwe]
jwe
parents:
diff changeset
250 %! title("spline fit to points from sin(x)");
6702
b2391d403ed2 [project @ 2007-06-12 21:39:26 by dbateman]
dbateman
parents: 6686
diff changeset
251 %! plot(xspline,sin(xspline),"r",xspline,yspline,"g-",x,y,"b+");
b2391d403ed2 [project @ 2007-06-12 21:39:26 by dbateman]
dbateman
parents: 6686
diff changeset
252 %! legend("original","interpolation","interpolation points");
5824
448f9982e7fb [project @ 2006-05-19 06:53:31 by jwe]
jwe
parents:
diff changeset
253 %! %--------------------------------------------------------
448f9982e7fb [project @ 2006-05-19 06:53:31 by jwe]
jwe
parents:
diff changeset
254 %! % confirm that interpolated function matches the original
448f9982e7fb [project @ 2006-05-19 06:53:31 by jwe]
jwe
parents:
diff changeset
255
6686
2aad75fcc93a [project @ 2007-06-03 20:58:28 by dbateman]
dbateman
parents: 6248
diff changeset
256 %!shared x,y,abserr
2aad75fcc93a [project @ 2007-06-03 20:58:28 by dbateman]
dbateman
parents: 6248
diff changeset
257 %! x = [0:10]; y = sin(x); abserr = 1e-14;
2aad75fcc93a [project @ 2007-06-03 20:58:28 by dbateman]
dbateman
parents: 6248
diff changeset
258 %!assert (spline(x,y,x), y, abserr);
2aad75fcc93a [project @ 2007-06-03 20:58:28 by dbateman]
dbateman
parents: 6248
diff changeset
259 %!assert (spline(x,y,x'), y', abserr);
2aad75fcc93a [project @ 2007-06-03 20:58:28 by dbateman]
dbateman
parents: 6248
diff changeset
260 %!assert (spline(x',y',x'), y', abserr);
2aad75fcc93a [project @ 2007-06-03 20:58:28 by dbateman]
dbateman
parents: 6248
diff changeset
261 %!assert (spline(x',y',x), y, abserr);
5824
448f9982e7fb [project @ 2006-05-19 06:53:31 by jwe]
jwe
parents:
diff changeset
262 %!assert (isempty(spline(x',y',[])));
448f9982e7fb [project @ 2006-05-19 06:53:31 by jwe]
jwe
parents:
diff changeset
263 %!assert (isempty(spline(x,y,[])));
6686
2aad75fcc93a [project @ 2007-06-03 20:58:28 by dbateman]
dbateman
parents: 6248
diff changeset
264 %!assert (spline(x,[y;y],x), [spline(x,y,x);spline(x,y,x)],abserr)
12608
59e2460acae1 make piecewise polynomial (pp) functions more compatible
Kai Habel <kai.habel@gmx.de>
parents: 12444
diff changeset
265 %!assert (spline(x,[y;y],x'), [spline(x,y,x);spline(x,y,x)],abserr)
59e2460acae1 make piecewise polynomial (pp) functions more compatible
Kai Habel <kai.habel@gmx.de>
parents: 12444
diff changeset
266 %!assert (spline(x',[y;y],x), [spline(x,y,x);spline(x,y,x)],abserr)
59e2460acae1 make piecewise polynomial (pp) functions more compatible
Kai Habel <kai.habel@gmx.de>
parents: 12444
diff changeset
267 %!assert (spline(x',[y;y],x'), [spline(x,y,x);spline(x,y,x)],abserr)
6014
437f9086b967 [project @ 2006-09-30 15:27:06 by jwe]
jwe
parents: 5838
diff changeset
268 %! y = cos(x) + i*sin(x);
6686
2aad75fcc93a [project @ 2007-06-03 20:58:28 by dbateman]
dbateman
parents: 6248
diff changeset
269 %!assert (spline(x,y,x), y, abserr)
2aad75fcc93a [project @ 2007-06-03 20:58:28 by dbateman]
dbateman
parents: 6248
diff changeset
270 %!assert (real(spline(x,y,x)), real(y), abserr);
2aad75fcc93a [project @ 2007-06-03 20:58:28 by dbateman]
dbateman
parents: 6248
diff changeset
271 %!assert (real(spline(x,y,x.')), real(y).', abserr);
2aad75fcc93a [project @ 2007-06-03 20:58:28 by dbateman]
dbateman
parents: 6248
diff changeset
272 %!assert (real(spline(x.',y.',x.')), real(y).', abserr);
2aad75fcc93a [project @ 2007-06-03 20:58:28 by dbateman]
dbateman
parents: 6248
diff changeset
273 %!assert (real(spline(x.',y,x)), real(y), abserr);
2aad75fcc93a [project @ 2007-06-03 20:58:28 by dbateman]
dbateman
parents: 6248
diff changeset
274 %!assert (imag(spline(x,y,x)), imag(y), abserr);
2aad75fcc93a [project @ 2007-06-03 20:58:28 by dbateman]
dbateman
parents: 6248
diff changeset
275 %!assert (imag(spline(x,y,x.')), imag(y).', abserr);
2aad75fcc93a [project @ 2007-06-03 20:58:28 by dbateman]
dbateman
parents: 6248
diff changeset
276 %!assert (imag(spline(x.',y.',x.')), imag(y).', abserr);
2aad75fcc93a [project @ 2007-06-03 20:58:28 by dbateman]
dbateman
parents: 6248
diff changeset
277 %!assert (imag(spline(x.',y,x)), imag(y), abserr);
10658
c66a4657d764 spline.m: Ignore NaNs within input vectors.
Ben Abbott <bpabbott@mac.com>
parents: 10549
diff changeset
278 %!test
c66a4657d764 spline.m: Ignore NaNs within input vectors.
Ben Abbott <bpabbott@mac.com>
parents: 10549
diff changeset
279 %! xnan = 5;
c66a4657d764 spline.m: Ignore NaNs within input vectors.
Ben Abbott <bpabbott@mac.com>
parents: 10549
diff changeset
280 %! y(x==xnan) = NaN;
c66a4657d764 spline.m: Ignore NaNs within input vectors.
Ben Abbott <bpabbott@mac.com>
parents: 10549
diff changeset
281 %! ok = ! isnan (y);
c66a4657d764 spline.m: Ignore NaNs within input vectors.
Ben Abbott <bpabbott@mac.com>
parents: 10549
diff changeset
282 %! assert (spline (x, y, x(ok)), y(ok), abserr);
c66a4657d764 spline.m: Ignore NaNs within input vectors.
Ben Abbott <bpabbott@mac.com>
parents: 10549
diff changeset
283 %!test
c66a4657d764 spline.m: Ignore NaNs within input vectors.
Ben Abbott <bpabbott@mac.com>
parents: 10549
diff changeset
284 %! ok = ! isnan (y);
c66a4657d764 spline.m: Ignore NaNs within input vectors.
Ben Abbott <bpabbott@mac.com>
parents: 10549
diff changeset
285 %! assert (! isnan (spline (x, y, x(!ok))));
12444
07e102029d2a spline.m: Allow length(x) == 2 and unsorted x values.
Ben Abbott <bpabbott@mac.com>
parents: 11587
diff changeset
286 %!test
07e102029d2a spline.m: Allow length(x) == 2 and unsorted x values.
Ben Abbott <bpabbott@mac.com>
parents: 11587
diff changeset
287 %! x = [1,2];
07e102029d2a spline.m: Allow length(x) == 2 and unsorted x values.
Ben Abbott <bpabbott@mac.com>
parents: 11587
diff changeset
288 %! y = [1,4];
07e102029d2a spline.m: Allow length(x) == 2 and unsorted x values.
Ben Abbott <bpabbott@mac.com>
parents: 11587
diff changeset
289 %! assert (spline (x,y,x), [1,4], abserr);
07e102029d2a spline.m: Allow length(x) == 2 and unsorted x values.
Ben Abbott <bpabbott@mac.com>
parents: 11587
diff changeset
290 %!test
07e102029d2a spline.m: Allow length(x) == 2 and unsorted x values.
Ben Abbott <bpabbott@mac.com>
parents: 11587
diff changeset
291 %! x = [2,1];
07e102029d2a spline.m: Allow length(x) == 2 and unsorted x values.
Ben Abbott <bpabbott@mac.com>
parents: 11587
diff changeset
292 %! y = [1,4];
07e102029d2a spline.m: Allow length(x) == 2 and unsorted x values.
Ben Abbott <bpabbott@mac.com>
parents: 11587
diff changeset
293 %! assert (spline (x,y,x), [1,4], abserr);
07e102029d2a spline.m: Allow length(x) == 2 and unsorted x values.
Ben Abbott <bpabbott@mac.com>
parents: 11587
diff changeset
294 %!test
07e102029d2a spline.m: Allow length(x) == 2 and unsorted x values.
Ben Abbott <bpabbott@mac.com>
parents: 11587
diff changeset
295 %! x = [1,2];
07e102029d2a spline.m: Allow length(x) == 2 and unsorted x values.
Ben Abbott <bpabbott@mac.com>
parents: 11587
diff changeset
296 %! y = [1,2,3,4];
07e102029d2a spline.m: Allow length(x) == 2 and unsorted x values.
Ben Abbott <bpabbott@mac.com>
parents: 11587
diff changeset
297 %! pp = spline (x,y);
07e102029d2a spline.m: Allow length(x) == 2 and unsorted x values.
Ben Abbott <bpabbott@mac.com>
parents: 11587
diff changeset
298 %! [x,P] = unmkpp (pp);
07e102029d2a spline.m: Allow length(x) == 2 and unsorted x values.
Ben Abbott <bpabbott@mac.com>
parents: 11587
diff changeset
299 %! assert (norm (P-[3,-3,1,2]), 0, abserr);
07e102029d2a spline.m: Allow length(x) == 2 and unsorted x values.
Ben Abbott <bpabbott@mac.com>
parents: 11587
diff changeset
300 %!test
07e102029d2a spline.m: Allow length(x) == 2 and unsorted x values.
Ben Abbott <bpabbott@mac.com>
parents: 11587
diff changeset
301 %! x = [2,1];
07e102029d2a spline.m: Allow length(x) == 2 and unsorted x values.
Ben Abbott <bpabbott@mac.com>
parents: 11587
diff changeset
302 %! y = [1,2,3,4];
07e102029d2a spline.m: Allow length(x) == 2 and unsorted x values.
Ben Abbott <bpabbott@mac.com>
parents: 11587
diff changeset
303 %! pp = spline (x,y);
07e102029d2a spline.m: Allow length(x) == 2 and unsorted x values.
Ben Abbott <bpabbott@mac.com>
parents: 11587
diff changeset
304 %! [x,P] = unmkpp (pp);
07e102029d2a spline.m: Allow length(x) == 2 and unsorted x values.
Ben Abbott <bpabbott@mac.com>
parents: 11587
diff changeset
305 %! assert (norm (P-[7,-9,1,3]), 0, abserr);