share/numeric/interpol.mac

   1 /*               COPYRIGHT NOTICE
   2
   3 Copyright (C) 2005-2012 Mario Rodriguez Riotorto
   4
   5 This program is free software; you can redistribute
   6 it and/or modify it under the terms of the
   7 GNU General Public License as published by
   8 the Free Software Foundation; either version 2
   9 of the License, or (at your option) any later version.
  10
  11 This program is distributed in the hope that it
  12 will be useful, but WITHOUT ANY WARRANTY;
  13 without even the implied warranty of MERCHANTABILITY
  14 or FITNESS FOR A PARTICULAR PURPOSE. See the
  15 GNU General Public License for more details at
  16 http://www.gnu.org/copyleft/gpl.html
  17 */
  18
  19
  20 /*             INTRODUCTION
  21
  22 This package defines some interpolation techniques.
  23
  24 For questions, suggestions, bugs and the like, feel free
  25 to contact me at
  26
  27 mario @@@ edu DOT xunta DOT es
  28
  29 */
  30
  31
  32
  33 /* Returns de input in the form of a list of pairs, ordered wrt the first */
  34 /* coordinate. The argument must be either:                               */
  35 /*  a) a two column matrix, p:matrix([2,4],[5,6],[9,3])                   */
  36 /*  b) a list of pairs, p: [[2,4],[5,6],[9,3]]                            */
  37 /*  c) a list of numbers, p: [4,6,3], in which case the abscissas will be */
  38 /*     assigned automatically to 1, 2, 3, etc.                            */
  39 /* This is an uxiliary function for the 'interpol' package.               */
  40 interpol_check_input(data,funame):=
  41  block([n,out],
  42    if not listp(data) and not matrixp(data)
  43       then error("Argument to '",funame,"' must be a list or matrix"),
  44    n: length(data),
  45    if n<2
  46       then error("Argument to '",funame,"' has too few sample points")
  47    elseif listp(data) and every('identity,map(lambda([x], listp(x) and length(x)=2),data))
  48       then out: sort(data)
  49    elseif matrixp(data) and length(data[1]) = 2
  50       then out: sort(args(data))
  51    elseif listp(data) and every('identity,map(lambda([x], not listp(x)),data))
  52       then out: makelist([i,data[i]],i,1,n)
  53       else error("Error in arguments to '",funame,"' function"),
  54    /* controlling duplicated x's */
  55    for i:2 thru n do
  56       if out[i-1][1] = out[i][1]
  57          then error("Duplicated abscissas are not allowed"),
  58    out )$
  59
  60
  61
  62 /* Lagrangian interpolation. The argument must be either:                      */
  63 /*  a) a two column matrix, p:matrix([2,4],[5,6],[9,3])                        */
  64 /*  b) a list of pairs, p: [[2,4],[5,6],[9,3]]                                 */
  65 /*  c) a list of numbers, p: [4,6,3], in which case the abscissas will be      */
  66 /*     assigned automatically to 1, 2, 3, etc.                                 */
  67 /* In cases a) and b) the pairs are ordered wrt the 1st. coordinate before any */
  68 /* computation is made. Option:                                                */
  69 /*   'varname='x: the name of the independent variable                         */
  70 /* Sample session:                                                             */
  71 /* load(interpol);                                                             */
  72 /* p:[[7,2],[8,2],[1,5],[3,2],[6,7]]$                                          */
  73 /* lagrange(p);                                                                */
  74 /* f(x):=''%;                                                                  */
  75 /* map(f,[2.3,5/7,%pi]);                                                       */
  76 /* load(draw)$;                                                                */
  77 /* draw2d(                                                                     */
  78 /*    explicit(f(x),x,0,9),                                                    */
  79 /*    point_size = 3,                                                          */
  80 /*    points(p)) $                                                             */
  81 lagrange(tab,[select]) := block([n,sum:0,prod,options,defaults,ratprint:false,tab2],
  82    tab2: interpol_check_input(tab,"lagrange"),
  83    options:  ['varname],
  84    defaults: ['x],
  85    for i in select do(
  86       aux: ?position(lhs(i),options),
  87       if numberp(aux) and aux <= length(options) and aux >= 1
  88         then defaults[aux]: rhs(i)),
  89    if not symbolp(defaults[1])
  90       then error("Option 'varname' is not correct"),
  91
  92    /* constructing the interpolating polynomial */
  93    n: length(tab2),
  94    for i:1 thru n do(
  95       prod: 1,
  96       for k:1 thru n do
  97          if k#i then prod: prod * (defaults[1]-tab2[k][1]) / (tab2[i][1]-tab2[k][1]),
  98       sum: sum + prod * tab2[i][2] ),
  99    sum )$
 100
 101
 102
 103 /* Characteristic function for intervals. Returns true iif  */
 104 /* z belongs to [l1, l2). This is an auxiliary function to  */
 105 /* be called from linearinterpol and cspline                */
 106 charfun2(z,l1,l2):= charfun(l1 <= z and z < l2)$
 107
 108
 109
 110 /* Linear interpolation. The argument must be either:                          */
 111 /*  a) a two column matrix, p:matrix([2,4],[5,6],[9,3])                        */
 112 /*  b) a list of pairs, p: [[2,4],[5,6],[9,3]]                                 */
 113 /*  c) a list of numbers, p: [4,6,3], in which case the abscissas will be      */
 114 /*     assigned automatically to 1, 2, 3, etc.                                 */
 115 /* In cases a) and b) the pairs are ordered wrt the 1st. coordinate before any */
 116 /* computation is made. Option:                                                */
 117 /*   'varname='x: the name of the independent variable                         */
 118 /* Sample session:                                                             */
 119 /* load(interpol);                                                             */
 120 /* p:[[7,2],[8,2],[1,5],[3,2],[6,7]]$                                          */
 121 /* linearinterpol(p);                                                          */
 122 /* f(x):=''%;                                                                  */
 123 /* map(f,[2.3,5/7,%pi]);                                                       */
 124 /* load(draw)$;                                                                */
 125 /* draw2d(                                                                     */
 126 /*    explicit(f(x),x,0,9),                                                    */
 127 /*    point_size = 3,                                                          */
 128 /*    points(p)) $                                                             */
 129 linearinterpol(tab,[select]) := block([n,s:0,a,b,options, defaults,ratprint:false,tab2],
 130    tab2: interpol_check_input(tab,"linearinterpol"),
 131    options:  ['varname],
 132    defaults: ['x],
 133    for i in select do(
 134       aux: ?position(lhs(i),options),
 135       if numberp(aux) and aux <= length(options) and aux >= 1
 136         then defaults[aux]: rhs(i)),
 137    if not symbolp(defaults[1])
 138       then error("Option 'varname' is not correct"),
 139
 140    /* constructing the interpolating polynomial */
 141    n: length(tab2),
 142    if n=2 /* case of two points */
 143       then s: tab2[2][2] + (tab2[2][2]-tab2[1][2]) *
 144                           (defaults[1]-tab2[2][1]) /
 145                           (tab2[2][1]-tab2[1][1])
 146       else for i:2 thru n do(
 147                if i=2
 148                   then (a: 'minf, b: tab2[i][1])
 149                   else if i=n
 150                           then (a: tab2[i-1][1], b: 'inf)
 151                           else (a: tab2[i-1][1], b: tab2[i][1]),
 152                s: s + apply('charfun2,[defaults[1], a, b]) *
 153                       expand( tab2[i][2] + (tab2[i][2]-tab2[i-1][2]) *
 154                                     (defaults[1]-tab2[i][1]) /
 155                                     (tab2[i][1]-tab2[i-1][1]) )   ),
 156    s )$
 157
 158
 159
 160 /* Cubic splines interpolation. The argument must be either:                           */
 161 /*  a) a two column matrix, p:matrix([2,4],[5,6],[9,3])                                */
 162 /*  b) a list of pairs, p: [[2,4],[5,6],[9,3]]                                         */
 163 /*  c) a list of numbers, p: [4,6,3], in which case the abscissas will be              */
 164 /*     assigned automatically to 1, 2, 3, etc.                                         */
 165 /* In cases a) and b) the pairs are ordered wrt the 1st. coordinate before any         */
 166 /* computation is made. Options:                                                       */
 167 /*   'd1='unknown: 1st derivative at x_1; if it is 'unknown, the second derivative     */
 168 /*         at x_1 is made equal to 0 (natural cubic spline); if it is equal to a       */
 169 /*         number, the second derivative is estimated based on this number             */
 170 /*   'd2='unknown: 1st derivative at x_n; if it is 'unknown, the second derivative     */
 171 /*         at x_n is made equal to 0 (natural cubic spline); if it is equal to a       */
 172 /*         number, the second derivative is estimated based on this number             */
 173 /*   'varname='x: the name of the independent variable                                 */
 174 /* Reference: this algorithm is based on 'Numerical Recipes in C', section 3.3         */
 175 /* Sample session:                                                                     */
 176 /* load(interpol);                                                                     */
 177 /* p:[[7,2],[8,2],[1,5],[3,2],[6,7]]$                                                  */
 178 /* cspline(p); ==> natural cubic spline (second derivatives are zero in both extremes) */
 179 /* f(x):=''%;                                                                          */
 180 /* map(f,[2.3,5/7,%pi]);                                                               */
 181 /* load(draw)$;                                                                        */
 182 /* draw2d(                                                                             */
 183 /*    explicit(f(x),x,0,9),                                                            */
 184 /*    point_size = 3,                                                                  */
 185 /*    points(p)) $                                                                     */
 186 /* cspline(p,d1=0,dn=0);                                                               */
 187 /* g(x):=''%;                                                                          */
 188 /* draw2d(                                                                             */
 189 /*    explicit(g(x),x,0,9),                                                            */
 190 /*    point_size = 3,                                                                  */
 191 /*    points(p)) $                                                                     */
 192 cspline(tab,[select]):= block([options, defaults, n, aux, y2, u, sig, p,
 193                                qn, un, a, b, s:0, aj, bj, cj, dj, ratprint:false,tab2],
 194    tab2: interpol_check_input(tab,"cspline"),
 195    options:  ['d1, 'dn, 'varname],
 196    defaults: ['unknown, 'unknown, 'x],
 197    for i in select do(
 198       aux: ?position(lhs(i),options),
 199       if numberp(aux) and aux <= length(options) and aux >= 1
 200         then defaults[aux]: rhs(i)),
 201    if not numberp(defaults[1]) and  defaults[1] # 'unknown
 202       then error("Option 'd1' is not correct"),
 203    if not numberp(defaults[2]) and  defaults[2] # 'unknown
 204       then error("Option 'dn' is not correct"),
 205    if not symbolp(defaults[3])
 206       then error("Option 'varname' is not correct"),
 207
 208    /* if tab2 contains only two points, linear interpolation */
 209    n: length(tab2),
 210    if n=2 /* case of two points */
 211       then return(ratsimp( tab2[2][2] + (tab2[2][2]-tab2[1][2]) *
 212                                        (defaults[3]-tab2[2][1]) /
 213                                        (tab2[2][1]-tab2[1][1]))),
 214
 215
 216    /* constructing the interpolating polynomial */
 217    y2: makelist(0,i,1,n),
 218    u: makelist(0,i,1,n-1),
 219
 220    /* controlling the lower boundary condition */
 221    if /*d1*/ defaults[1] = 'unknown
 222       then (y2[1]: 0,
 223             u[1]: 0)
 224       else (y2[1]: -1/2,
 225             u[1]: 3 / (tab2[2][1]-tab2[1][1]) *
 226                       ((tab2[2][2] - tab2[1][2])/(tab2[2][1] - tab2[1][1]) - defaults[1]) ),
 227
 228    /* decomposition loop of the triangular algorithm */
 229    for i:2 thru n-1 do (
 230       sig: (tab2[i][1] - tab2[i-1][1]) / (tab2[i+1][1] - tab2[i-1][1]),
 231       p: sig * y2[i-1] + 2,
 232       y2[i]: (sig - 1) / p,
 233       u[i]: (tab2[i+1][2] - tab2[i][2]) /(tab2[i+1][1] - tab2[i][1]) -
 234             (tab2[i][2] - tab2[i-1][2]) /(tab2[i][1] - tab2[i-1][1]),
 235       u[i]: (6 * u[i] / (tab2[i+1][1] - tab2[i-1][1]) - sig * u[i-1]) / p ) ,
 236
 237    /* controlling the upper boundary condition */
 238    if /*dn*/ defaults[2] = 'unknown
 239       then (qn: 0,
 240             un: 0)
 241       else (qn: 1/2,
 242             un: 3 / (tab2[n][1] - tab2[n-1][1]) *
 243                 (defaults[2] - (tab2[n][2] - tab2[n-1][2]) / (tab2[n][1] - tab2[n-1][1]))),
 244    y2[n]: (un - qn * u[n-1]) / (qn * y2[n-1] + 1),
 245
 246    /* backsubstitution loop of the tridiagonal algorithm */
 247    for k: n-1 thru 1 step -1 do
 248       y2[k]: y2[k] * y2[k+1] + u[k],
 249
 250    /* constructing the cubic splines */
 251    for j:2 thru n do (
 252       if j=2
 253           then (a: 'minf, b: tab2[j][1] )
 254           else if j=n
 255                   then (a: tab2[j-1][1], b: 'inf)
 256                   else (a: tab2[j-1][1], b: tab2[j][1]),
 257       /* in the following sentences, defaults[3] is variable's name */
 258       aux: (tab2[j][1] - tab2[j-1][1]),
 259       aj: (tab2[j][1] - defaults[3]) / aux,
 260       bj: (defaults[3] - tab2[j-1][1]) / aux,
 261       aux: aux * aux /6,
 262       cj: (aj^3 - aj) * aux,
 263       dj: (bj^3 - bj) * aux,
 264
 265       s: s + charfun2(defaults[3], a, b) *
 266              expand(aj * tab2[j-1][2] + bj * tab2[j][2] + cj * y2[j-1] + dj * y2[j])  ),
 267    s )$
 268
 269
 270
 271 /* Rational interpolation, with interpolating function of the form             */
 272 /*               r                                                             */
 273 /*           p  x  + ... + p  x + p                                            */
 274 /*            r             1      0                                           */
 275 /*    R(x) = ------------------------                                          */
 276 /*               m                                                             */
 277 /*           q  x  + ... + q  x + q                                            */
 278 /*            m             1      0                                           */
 279 /* The 2nd. argument r is the degree of the numerator (r < sample size). The   */
 280 /* degree of the denominator is calculated as m: sample_size - r - 1.          */
 281 /* The 1st. argument must be either:                                           */
 282 /*  a) a two column matrix, p:matrix([2,4],[5,6],[9,3])                        */
 283 /*  b) a list of pairs, p: [[2,4],[5,6],[9,3]]                                 */
 284 /*  c) a list of numbers, p: [4,6,3], in which case the abscissas will be      */
 285 /*     assigned automatically to 1, 2, 3, etc.                                 */
 286 /* In cases a) and b) the pairs are ordered wrt the 1st. coordinate before any */
 287 /* computation is made. Option:                                                */
 288 /*   'varname='x: the name of the independent variable                         */
 289 /* Sample session:                                                             */
 290 /* load(interpol)$                                                             */
 291 /* load(draw)$                                                                 */
 292 /* p: [[7.2, 2.5], [8.5, 2.1], [1.6, 5.1], [3.4, 2.4], [6.7, 7.9]]$            */
 293 /* for k:0 thru length(p)-1 do                                                 */
 294 /*   draw2d(explicit(ratinterpol(p,k),x,0,9),                                  */
 295 /*          point_size = 3,                                                    */
 296 /*          points(p),                                                         */
 297 /*          title = concat("Degree of numerator = ",k),                        */
 298 /*          yrange=[0,10])$                                                    */
 299 ratinterpol(tab,r,[select]) :=
 300   block([n,m,coef,unk,sol,lovtab,lov,options,defaults,ratprint:false,tab2],
 301     tab2: interpol_check_input(tab,"ratinterpol"),
 302     options:  ['varname],
 303     defaults: ['x],
 304     for i in select do(
 305        aux: ?position(lhs(i),options),
 306        if numberp(aux) and aux <= length(options) and aux >= 1
 307          then defaults[aux]: rhs(i)),
 308     if not symbolp(defaults[1])
 309        then error("Option 'varname' is not correct"),
 310
 311     /* constructing the interpolating rational function */
 312     n: length(tab2),
 313     if not integerp(r) or r > n-1 or r < 0
 314       then error("Degree of numerator must be a positive integer less than sample size"),
 315     m: n - r - 1, /* degree of denominator */
 316
 317     /* coef is the matrix of an homogeneous linear system */
 318     coef: apply(matrix,
 319                 makelist(block([x,y],
 320                            [x,y]: p,
 321                            append([1],
 322                                   makelist(x^k,k,1,r),
 323                                   [-y],
 324                                   makelist(-y*x^k,k,1,m))),
 325                          p, tab2)),
 326     unk: makelist(gensym(), k, r+m+2),
 327     sol: map(last, linsolve(flatten(args(coef.unk)),unk)),
 328     lovtab : listofvars(tab2),
 329     lov: listofvars(sol),
 330     lov: listify(setdifference(setify(lov),setify(lovtab))),
 331     sol: subst(map(lambda([z1,z2], z1=z2), lov, makelist(1, k, length(lov))), sol),
 332     makelist(sol[k],k,1,r+1) . makelist(defaults[1]^k,k,0,r) /
 333         makelist(sol[k],k,r+2,r+2+m) . makelist(defaults[1]^k,k,0,m) )$
 334