doc/info/interpol.texi

   1 @menu
   2 * Introduction to interpol::
   3 * Functions and Variables for interpol::
   4 @end menu
   5
   6 @node Introduction to interpol, Functions and Variables for interpol, interpol-pkg, interpol-pkg
   7 @section Introduction to interpol
   8
   9 Package @code{interpol} defines the Lagrangian, the linear and the cubic
  10 splines methods for polynomial interpolation.
  11
  12 For comments, bugs or suggestions, please contact me at @var{'mario AT edu DOT xunta DOT es'}.
  13
  14 @opencatbox{Categories:}
  15 @category{Numerical methods}
  16 @category{Share packages}
  17 @category{Package interpol}
  18 @closecatbox
  19
  20 @node Functions and Variables for interpol,  , Introduction to interpol, interpol-pkg
  21 @section Functions and Variables for interpol
  22
  23
  24 @anchor{lagrange}
  25 @deffn {Function} lagrange @
  26 @fname{lagrange} (@var{points}) @
  27 @fname{lagrange} (@var{points}, @var{option})
  28
  29 Computes the polynomial interpolation by the Lagrangian method. Argument @var{points} must be either:
  30
  31 @itemize @bullet
  32 @item
  33 a two column matrix, @code{p:matrix([2,4],[5,6],[9,3])},
  34 @item
  35 a list of pairs, @code{p: [[2,4],[5,6],[9,3]]},
  36 @item
  37 a list of numbers, @code{p: [4,6,3]}, in which case the abscissas will be assigned automatically to 1, 2, 3, etc.
  38 @end itemize
  39
  40 In the first two cases the pairs are ordered with respect to the first coordinate before making computations.
  41
  42 With the @var{option} argument it is possible to select the name for the independent variable, which is @code{'x} by default; to define another one, write something like @code{varname='z}.
  43
  44 Note that when working with high degree polynomials, floating point evaluations are unstable.
  45
  46 See also @mrefcomma{linearinterpol}  @mrefcomma{cspline}  and @mrefdot{ratinterpol}
  47
  48 Examples:
  49
  50 @example
  51 (%i1) load("interpol")$
  52 (%i2) p:[[7,2],[8,2],[1,5],[3,2],[6,7]]$
  53 (%i3) lagrange(p);
  54        (x - 7) (x - 6) (x - 3) (x - 1)
  55 (%o3)  -------------------------------
  56                      35
  57    (x - 8) (x - 6) (x - 3) (x - 1)
  58  - -------------------------------
  59                  12
  60    7 (x - 8) (x - 7) (x - 3) (x - 1)
  61  + ---------------------------------
  62                   30
  63    (x - 8) (x - 7) (x - 6) (x - 1)
  64  - -------------------------------
  65                  60
  66    (x - 8) (x - 7) (x - 6) (x - 3)
  67  + -------------------------------
  68                  84
  69 (%i4) f(x):=''%;
  70                (x - 7) (x - 6) (x - 3) (x - 1)
  71 (%o4)  f(x) := -------------------------------
  72                              35
  73    (x - 8) (x - 6) (x - 3) (x - 1)
  74  - -------------------------------
  75                  12
  76    7 (x - 8) (x - 7) (x - 3) (x - 1)
  77  + ---------------------------------
  78                   30
  79    (x - 8) (x - 7) (x - 6) (x - 1)
  80  - -------------------------------
  81                  60
  82    (x - 8) (x - 7) (x - 6) (x - 3)
  83  + -------------------------------
  84                  84
  85 (%i5) /* Evaluate the polynomial at some points */
  86       expand(map(f,[2.3,5/7,%pi]));
  87                                   4          3           2
  88                     919062  73 %pi    701 %pi    8957 %pi
  89 (%o5)  [- 1.567535, ------, ------- - -------- + ---------
  90                     84035     420       210         420
  91                                              5288 %pi   186
  92                                            - -------- + ---]
  93                                                105       5
  94 (%i6) %,numer;
  95 (%o6) [- 1.567535, 10.9366573451538, 2.89319655125692]
  96 (%i7) load("draw")$  /* load draw package */
  97 (%i8) /* Plot the polynomial together with points */
  98       draw2d(
  99         color      = red,
 100         key        = "Lagrange polynomial",
 101         explicit(f(x),x,0,10),
 102         point_size = 3,
 103         color      = blue,
 104         key        = "Sample points",
 105         points(p))$
 106 (%i9) /* Change variable name */
 107       lagrange(p, varname=w);
 108        (w - 7) (w - 6) (w - 3) (w - 1)
 109 (%o9)  -------------------------------
 110                      35
 111    (w - 8) (w - 6) (w - 3) (w - 1)
 112  - -------------------------------
 113                  12
 114    7 (w - 8) (w - 7) (w - 3) (w - 1)
 115  + ---------------------------------
 116                   30
 117    (w - 8) (w - 7) (w - 6) (w - 1)
 118  - -------------------------------
 119                  60
 120    (w - 8) (w - 7) (w - 6) (w - 3)
 121  + -------------------------------
 122                  84
 123 @end example
 124
 125 @opencatbox{Categories:}
 126 @category{Package interpol}
 127 @closecatbox
 128
 129 @end deffn
 130
 131
 132 @anchor{charfun2}
 133 @deffn {Function} charfun2 (@var{x}, @var{a}, @var{b})
 134
 135 The characteristic or indicator function on the half-open interval @math{[@var{a}, @var{b})},
 136 that is, including @var{a} and excluding @var{b}.
 137
 138 When @math{@var{x} >= @var{a}} and @math{@var{x} < @var{b}} evaluates to @code{true} or @code{false},
 139 @code{charfun2} returns 1 or 0, respectively.
 140
 141 Otherwise, @code{charfun2} returns a partially-evaluated result in terms of @mref{charfun}.
 142
 143 Package @code{interpol} loads this function.
 144
 145 See also @mrefdot{charfun}
 146
 147 Examples:
 148
 149 When @math{@var{x} >= @var{a}} and @math{@var{x} < @var{b}} evaluates to @code{true} or @code{false},
 150 @code{charfun2} returns 1 or 0, respectively.
 151
 152 @c ===beg===
 153 @c load ("interpol") $
 154 @c charfun2 (5, 0, 100);
 155 @c charfun2 (-5, 0, 100);
 156 @c ===end===
 157 @example
 158 (%i1) load ("interpol") $
 159 (%i2) charfun2 (5, 0, 100);
 160 (%o2)                           1
 161 (%i3) charfun2 (-5, 0, 100);
 162 (%o3)                           0
 163 @end example
 164
 165 Otherwise, @code{charfun2} returns a partially-evaluated result in terms of @mref{charfun}.
 166
 167 @c ===beg===
 168 @c load ("interpol") $
 169 @c charfun2 (t, 0, 100);
 170 @c charfun2 (5, u, v);
 171 @c assume (v > u, u > 5);
 172 @c charfun2 (5, u, v);
 173 @c ===end===
 174 @example
 175 (%i1) load ("interpol") $
 176 (%i2) charfun2 (t, 0, 100);
 177 (%o2)            charfun((0 <= t) and (t < 100))
 178 (%i3) charfun2 (5, u, v);
 179 (%o3)             charfun((u <= 5) and (5 < v))
 180 (%i4) assume (v > u, u > 5);
 181 (%o4)                    [v > u, u > 5]
 182 (%i5) charfun2 (5, u, v);
 183 (%o5)                           0
 184 @end example
 185
 186 @opencatbox{Categories:}
 187 @category{Package interpol}
 188 @closecatbox
 189
 190 @end deffn
 191
 192
 193 @anchor{linearinterpol}
 194 @deffn {Function} linearinterpol @
 195 @fname{linearinterpol} (@var{points}) @
 196 @fname{linearinterpol} (@var{points}, @var{option})
 197
 198 Computes the polynomial interpolation by the linear method. Argument @var{points} must be either:
 199
 200 @itemize @bullet
 201 @item
 202 a two column matrix, @code{p:matrix([2,4],[5,6],[9,3])},
 203 @item
 204 a list of pairs, @code{p: [[2,4],[5,6],[9,3]]},
 205 @item
 206 a list of numbers, @code{p: [4,6,3]}, in which case the abscissas will be assigned automatically to 1, 2, 3, etc.
 207 @end itemize
 208
 209 In the first two cases the pairs are ordered with respect to the first coordinate before making computations.
 210
 211 With the @var{option} argument it is possible to select the name for the independent variable, which is @code{'x} by default; to define another one, write something like @code{varname='z}.
 212
 213 See also @mrefcomma{lagrange}  @mrefcomma{cspline}  and @mrefdot{ratinterpol}
 214
 215 Examples:
 216
 217 @example
 218 (%i1) load("interpol")$
 219 (%i2) p: matrix([7,2],[8,3],[1,5],[3,2],[6,7])$
 220 (%i3) linearinterpol(p);
 221         13   3 x
 222 (%o3)  (-- - ---) charfun2(x, minf, 3)
 223         2     2
 224  + (x - 5) charfun2(x, 7, inf) + (37 - 5 x) charfun2(x, 6, 7)
 225     5 x
 226  + (--- - 3) charfun2(x, 3, 6)
 227      3
 228
 229 (%i4) f(x):=''%;
 230                 13   3 x
 231 (%o4)  f(x) := (-- - ---) charfun2(x, minf, 3)
 232                 2     2
 233  + (x - 5) charfun2(x, 7, inf) + (37 - 5 x) charfun2(x, 6, 7)
 234     5 x
 235  + (--- - 3) charfun2(x, 3, 6)
 236      3
 237 (%i5)  /* Evaluate the polynomial at some points */
 238        map(f,[7.3,25/7,%pi]);
 239                             62  5 %pi
 240 (%o5)                 [2.3, --, ----- - 3]
 241                             21    3
 242 (%i6) %,numer;
 243 (%o6)  [2.3, 2.952380952380953, 2.235987755982989]
 244 (%i7) load("draw")$  /* load draw package */
 245 (%i8)  /* Plot the polynomial together with points */
 246        draw2d(
 247          color      = red,
 248          key        = "Linear interpolator",
 249          explicit(f(x),x,-5,20),
 250          point_size = 3,
 251          color      = blue,
 252          key        = "Sample points",
 253          points(args(p)))$
 254 (%i9)  /* Change variable name */
 255        linearinterpol(p, varname='s);
 256        13   3 s
 257 (%o9) (-- - ---) charfun2(s, minf, 3)
 258        2     2
 259  + (s - 5) charfun2(s, 7, inf) + (37 - 5 s) charfun2(s, 6, 7)
 260     5 s
 261  + (--- - 3) charfun2(s, 3, 6)
 262      3
 263 @end example
 264
 265 @opencatbox{Categories:}
 266 @category{Package interpol}
 267 @closecatbox
 268
 269 @end deffn
 270
 271
 272
 273 @anchor{cspline}
 274 @deffn {Function} cspline @
 275 @fname{cspline} (@var{points}) @
 276 @fname{cspline} (@var{points}, @var{option1}, @var{option2}, ...)
 277
 278 Computes the polynomial interpolation by the cubic splines method. Argument @var{points} must be either:
 279
 280 @itemize @bullet
 281 @item
 282 a two column matrix, @code{p:matrix([2,4],[5,6],[9,3])},
 283 @item
 284 a list of pairs, @code{p: [[2,4],[5,6],[9,3]]},
 285 @item
 286 a list of numbers, @code{p: [4,6,3]}, in which case the abscissas will be assigned automatically to 1, 2, 3, etc.
 287 @end itemize
 288
 289 In the first two cases the pairs are ordered with respect to the first coordinate before making computations.
 290
 291 There are three options to fit specific needs:
 292 @itemize @bullet
 293 @item
 294 @code{'d1}, default @code{'unknown}, is the first derivative at @math{x_1}; if it is @code{'unknown}, the second derivative at @math{x_1} is made equal to 0 (natural cubic spline); if it is equal to a number, the second derivative is calculated based on this number.
 295
 296 @item
 297 @code{'dn}, default @code{'unknown}, is the first derivative at @math{x_n}; if it is @code{'unknown}, the second derivative at @math{x_n} is made equal to 0 (natural cubic spline); if it is equal to a number, the second derivative is calculated based on this number.
 298
 299 @item
 300 @code{'varname}, default @code{'x}, is the name of the independent variable.
 301 @end itemize
 302
 303 See also @mrefcomma{lagrange}  @mrefcomma{linearinterpol}  and @mrefdot{ratinterpol}
 304
 305 Examples:
 306 @example
 307 (%i1) load("interpol")$
 308 (%i2) p:[[7,2],[8,2],[1,5],[3,2],[6,7]]$
 309 (%i3) /* Unknown first derivatives at the extremes
 310          is equivalent to natural cubic splines */
 311       cspline(p);
 312               3         2
 313         1159 x    1159 x    6091 x   8283
 314 (%o3)  (------- - ------- - ------ + ----) charfun2(x, minf, 3)
 315          3288      1096      3288    1096
 316             3         2
 317       2587 x    5174 x    494117 x   108928
 318  + (- ------- + ------- - -------- + ------) charfun2(x, 7, inf)
 319        1644       137       1644      137
 320           3          2
 321     4715 x    15209 x    579277 x   199575
 322  + (------- - -------- + -------- - ------) charfun2(x, 6, 7)
 323      1644       274        1644      274
 324             3         2
 325       3287 x    2223 x    48275 x   9609
 326  + (- ------- + ------- - ------- + ----) charfun2(x, 3, 6)
 327        4932       274      1644     274
 328
 329 (%i4) f(x):=''%$
 330 (%i5) /* Some evaluations */
 331       map(f,[2.3,5/7,%pi]), numer;
 332 (%o5) [1.991460766423356, 5.823200187269903, 2.227405312429507]
 333 (%i6) load("draw")$  /* load draw package */
 334 (%i7) /* Plotting interpolating function */
 335       draw2d(
 336         color      = red,
 337         key        = "Cubic splines",
 338         explicit(f(x),x,0,10),
 339         point_size = 3,
 340         color      = blue,
 341         key        = "Sample points",
 342         points(p))$
 343 (%i8) /* New call, but giving values at the derivatives */
 344       cspline(p,d1=0,dn=0);
 345               3          2
 346         1949 x    11437 x    17027 x   1247
 347 (%o8)  (------- - -------- + ------- + ----) charfun2(x, minf, 3)
 348          2256       2256      2256     752
 349             3          2
 350       1547 x    35581 x    68068 x   173546
 351  + (- ------- + -------- - ------- + ------) charfun2(x, 7, inf)
 352         564       564        141      141
 353          3          2
 354     607 x    35147 x    55706 x   38420
 355  + (------ - -------- + ------- - -----) charfun2(x, 6, 7)
 356      188       564        141      47
 357             3         2
 358       3895 x    1807 x    5146 x   2148
 359  + (- ------- + ------- - ------ + ----) charfun2(x, 3, 6)
 360        5076       188      141      47
 361 (%i8) /* Defining new interpolating function */
 362       g(x):=''%$
 363 (%i9) /* Plotting both functions together */
 364       draw2d(
 365         color      = black,
 366         key        = "Cubic splines (default)",
 367         explicit(f(x),x,0,10),
 368         color      = red,
 369         key        = "Cubic splines (d1=0,dn=0)",
 370         explicit(g(x),x,0,10),
 371         point_size = 3,
 372         color      = blue,
 373         key        = "Sample points",
 374         points(p))$
 375 @end example
 376
 377 @opencatbox{Categories:}
 378 @category{Package interpol}
 379 @closecatbox
 380
 381 @end deffn
 382
 383 @anchor{ratinterpol}
 384 @deffn {Function} ratinterpol @
 385 @fname{ratinterpol} (@var{points}, @var{numdeg}) @
 386 @fname{ratinterpol} (@var{points}, @var{numdeg}, @var{option1})
 387
 388 Generates a rational interpolator for data given by @var{points} and the degree of the numerator
 389 being equal to @var{numdeg}; the degree of the denominator is calculated
 390 automatically. Argument @var{points} must be either:
 391
 392 @itemize @bullet
 393 @item
 394 a two column matrix, @code{p:matrix([2,4],[5,6],[9,3])},
 395 @item
 396 a list of pairs, @code{p: [[2,4],[5,6],[9,3]]},
 397 @item
 398 a list of numbers, @code{p: [4,6,3]}, in which case the abscissas will be assigned automatically to 1, 2, 3, etc.
 399 @end itemize
 400
 401 In the first two cases the pairs are ordered with respect to the first coordinate before making computations.
 402
 403 There is one option to fit specific needs:
 404 @itemize @bullet
 405 @item
 406 @code{'varname}, default @code{'x}, is the name of the independent variable.
 407 @end itemize
 408
 409 See also @mrefcomma{lagrange}  @mrefcomma{linearinterpol}  @mrefcomma{cspline}  @mrefcomma{minpack_lsquares}  and @ref{lbfgs-pkg}
 410
 411 Examples:
 412
 413 @example
 414 (%i1) load("interpol")$
 415 (%i2) load("draw")$
 416 (%i3) p:[[7.2,2.5],[8.5,2.1],[1.6,5.1],[3.4,2.4],[6.7,7.9]]$
 417 (%i4) for k:0 thru length(p)-1 do
 418         draw2d(
 419           explicit(ratinterpol(p,k),x,0,9),
 420           point_size = 3,
 421           points(p),
 422           title = concat("Degree of numerator = ",k),
 423           yrange=[0,10])$
 424 @end example
 425
 426 @opencatbox{Categories:}
 427 @category{Package interpol}
 428 @closecatbox
 429
 430 @end deffn