share/numeric/rtest_interpol.mac

   1 (kill (all),
   2  load (interpol),
   3  float_approx_equal_tolerance : 1e-8,
   4  dat1: matrix([6.5,3],[1.8,-8],[9,2],[-2,4],[1/8,1],[9/12,3],[10,%pi]),
   5  dat2: [[2,4],[7,1],[-7,2/7],[5,9],[8,4/9],[9,4]],
   6  dat3: [8,5.9,4.1,8.6,9.7,3,99.8],
   7  0);
   8 0;
   9
  10
  11 /* Lagrange interpolation */
  12
  13 lagrange(dat1);
  14 3.178762054958412E-5*%pi*(x-9)*(x-6.5)*(x-1.8)*(x-3/4)*(x-1/8)*(x+2)
  15 -1.3795660402574614E-4*(x-10)*(x-6.5)*(x-1.8)*(x-3/4)*(x-1/8)*(x+2)
  16 +2.3412532015037063E-4*(x-10)*(x-9)*(x-1.8)*(x-3/4)*(x-1/8)*(x+2)
  17 +.004313780800478495*(x-10)*(x-9)*(x-6.5)*(x-3/4)*(x-1/8)*(x+2)
  18 +.003788399045316041*(x-10)*(x-9)*(x-6.5)*(x-1.8)*(x-1/8)*(x+2)
  19 -8.045639731099019E-4*(x-10)*(x-9)*(x-6.5)*(x-1.8)*(x-3/4)*(x+2)
  20 +1.605431979101289E-4*(x-10)*(x-9)*(x-6.5)*(x-1.8)*(x-3/4)*(x-1/8);
  21
  22 lagrange(dat2);
  23 ((x-8)*(x-7)*(x-5)*(x-2)*(x+7))/224
  24 -2*(x-9)*(x-7)*(x-5)*(x-2)*(x+7)/1215
  25 +(x-9)*(x-8)*(x-5)*(x-2)*(x+7)/280
  26 +(x-9)*(x-8)*(x-7)*(x-2)*(x+7)/-96
  27 +2*(x-9)*(x-8)*(x-7)*(x-5)*(x+7)/2835
  28 +(x-9)*(x-8)*(x-7)*(x-5)*(x-2)/-1270080;
  29
  30 lagrange(dat3, varname=z);
  31 .1386111111111111*(z-6)*(z-5)*(z-4)*(z-3)*(z-2)*(z-1)
  32 -((z-7)*(z-5)*(z-4)*(z-3)*(z-2)*(z-1))/40
  33 +.2020833333333333*(z-7)*(z-6)*(z-4)*(z-3)*(z-2)*(z-1)
  34 -.2388888888888889*(z-7)*(z-6)*(z-5)*(z-3)*(z-2)*(z-1)
  35 +.08541666666666665*(z-7)*(z-6)*(z-5)*(z-4)*(z-2)*(z-1)
  36 -.04916666666666667*(z-7)*(z-6)*(z-5)*(z-4)*(z-3)*(z-1)
  37 -((2-z)*(z-7)*(z-6)*(z-5)*(z-4)*(z-3))/90;
  38
  39 /* Linear interpolation */
  40
  41 linearinterpol(dat1);
  42 (20/17-24*x/17)*charfun(x<1/8)
  43 +(%pi*x-2*x-9*%pi+20)*charfun(9<=x)+(5.6-0.4*x)*charfun(6.5<=x and x<9)
  44 +(2.340425531914894*x-12.21276595744681)*charfun(1.8<=x and x<6.5)
  45 +(10.85714285714286-10.47619047619048*x)*charfun(3/4<=x and x<1.8)
  46 +(16*x/5+3/5)*charfun(1/8<=x and x<3/4);
  47
  48 linearinterpol(dat2);
  49 (26*x/63+200/63)*charfun(x<2)+((32*x)/9-28)*charfun(8<=x)
  50 +(44/9-5*x/9)*charfun(7<=x and x<8)
  51 +(29-4*x)*charfun(5<=x and x<7)
  52 +(5*x/3+2/3)*charfun(2<=x and x<5);
  53
  54 linearinterpol(dat3,varname=w);
  55 (10.1-2.1*w)*charfun(w<2)+(96.8*w-577.8)
  56 *charfun(6<=w)+(43.2-6.7*w)*charfun(5<=w and w<6)
  57 +(1.1*w+4.2)*charfun(4<=w and w<5)
  58 +(4.5*w-9.4)*charfun(3<=w and w<4)+(9.5-1.8*w)*
  59 charfun(2<=w and w<3);
  60
  61 /* Cubic splines interpolation */
  62
  63 cspline(dat1);
  64 (2.233410484248884E-4*%pi*x^3+0.658739927100864*x^3
  65 +0.00134004629054933*%pi*x^2+3.952439562605184*x^2
  66 +.001671568159305024*%pi*x+3.518491936013175*x
  67 -2.3032045618816613E-4*%pi+.4971450384125281)
  68 *charfun(x<1/8)
  69 +(-0.153921193200841*%pi*x^3-.06508981912631977*x^3
  70 +4.61763579602523*%pi*x^2+1.952694573789593*x^2
  71 -45.02243676705146*%pi*x-21.46185591876961*x+143.3819812688326*%pi
  72 +84.43892093505657)*charfun(9<=x)
  73 +(.07396021366527837*%pi*x^3+.3549365250720062*x^3
  74 -1.535162189359993*%pi*x^2-9.388016719565206*x^2
  75 +10.35274510141554*%pi*x+80.60454572142358*x-22.74356433656842*%pi
  76 -221.760283985523)*charfun(6.5<=x and x<9)
  77 +(-0.00936633903933342*%pi*x^3-.5795316182501081*x^3
  78 +.08970558837993714*%pi*x^2+8.834112075216023*x^2
  79 -.2088954538939998*%pi*x-37.83929144465441*x
  80 +.1399901999355959*%pi+34.86802987431265)*charfun(1.8<=x and x<6.5)
  81 +(0.0163989810743406*%pi*x^3+6.51539068532972*x^3
  82 -.04942714023390259*%pi*x^2-29.47846836411505*x^2
  83 +.04154345761091165*%pi*x+31.12335334614151*x
  84 -.01027314696735098*%pi-6.509557000164901)*charfun(3/4<=x and x<1.8)
  85 +(-.007441723733517281*%pi*x^3-10.14309672394197*x^3
  86 +.004214445583777643*%pi*x^2+8.003128306746245*x^2
  87 +.001312268247651485*%pi*x+3.012155842995544*x
  88 -2.153496265359353E-4*%pi+0.518242375621596)*charfun(1/8<=x and x<3/4);
  89
  90 cspline(dat2,d1=2,varname=r);
  91 (635237*r^3/13492332+436937*r^2/1124361+333821*r/642492+6946987/6746166)*charfun(r<2)
  92 +(-71261*r^3/83286+213783*r^2/9254-941613*r/4627+2720884/4627)
  93 *charfun(8<=r)
  94 +(13907*r^3/83286-39995*r^2/27762-14209*r/1983+383564/5949)
  95 *charfun(7<=r and r<8)
  96 +(81073*r^3/111048-491039*r^2/37016+1199465*r/15864-680383/5288)
  97 *charfun(5<=r and r<7)
  98 +(-165773*r^3/499716+110833*r^2/41643-670655*r/166572+1014497/249858)
  99 *charfun(2<=r and r<5);
 100
 101 cspline(dat3,d1=8/5,dn=1/3,varname=k);
 102 (2.410085470085469*k^3-13.34034188034188*k^2+21.05042735042734*k
 103 -2.120170940170938)*charfun(k<2)
 104 +(-119.4588034188034*k^3+2292.709401709402*k^2-14537.15418803419*k
 105 +30491.4882051282)*charfun(6<=k)
 106 +(61.94307692307693*k^3-972.5244444444445*k^2+5054.24888888889*k
 107 -8691.31794871795)*charfun(5<=k and k<6)
 108 +(-17.01350427350427*k^3+211.8242735042735*k^2-867.4947008547009*k
 109 +1178.254700854701)*charfun(4<=k and k<5)
 110 +(1.710940170940171*k^3-12.86905982905983*k^2+31.27863247863247*k
 111 -20.10974358974358)*charfun(3<=k and k<4)
 112 +(.4697435897435902*k^3-1.698290598290602*k^2-2.233675213675204*k
 113 +13.4025641025641)*charfun(2<=k and k<3);
 114
 115 /* Rational interpolation */
 116
 117 float(ratinterpol(dat1,4));
 118 (.1620328314185472*x^4-2.534071721796399*x^3+11.2288985461703*x^2
 119 -7.58855393365858*x+5.924974869319725)/(x^2-6.159704297522396*x+5.901285404551714);
 120
 121 ratinterpol(dat2,3,varname=t);
 122 (6519*t^3/19843-45257*t^2/19843-399142*t/19843+2751280/19843)
 123 /(t^2-274651*t/19843+925960/19843);
 124
 125 ratinterpol(dat3,5,varname=t);
 126 (-38471*t^5/35760+692581*t^4/35760-1525837*t^3/11920+13508819*t^2/35760
 127 -8456297*t/17880+216979/1490)/(t-1253/149);
 128
 129 (kill (h),
 130  h:[[1,(-(13+2*b)*1)/12],
 131     [2,(-(15+2*b)*1)/7],
 132     [3,(-(17+2*b)*3)/16],
 133     [4,(-(19+2*b)*2)/9],
 134     [5,(-(21+2*b)*1)/4],
 135     [6,(-(23+2*b)*3)/11]],
 136  0);
 137 0;
 138
 139 is (equal (ratinterpol(h,1),
 140            ((32*b^5+1552*b^4+30064*b^3+291128*b^2+1412514*b+2766033)*x
 141              /(32*b+16)+144)
 142              /(x^4-(2*b+53)*x^3/2+(4*b^2+128*b+1283)*x^2/4
 143                -(8*b^3+300*b^2+3974*b+19993)*x/8
 144                -(80*b^4+3440*b^3+56240*b^2+418500*b+1260063)/(16*b+8))));
 145 true;
 146
 147 is (equal (ratinterpol(h,2), (-x^2-(2*b+11)*x/2)/(x+5)));
 148 true;
 149
 150 is (equal (ratinterpol (subst (b=5/7, h), 2), subst (b=5/7, ratinterpol (h, 2))));
 151 true;
 152
 153 (reset(float_approx_equal_tolerance),0);
 154 0$