1 #ifndef lint 2 static char *sccsid = "@(#)spline.c 4.5 (Berkeley) 12/02/87"; 3 #endif 4 5 #include <stdio.h> 6 #include <math.h> 7 8 #define NP 1000 9 #define INF HUGE 10 11 struct proj { int lbf,ubf; float a,b,lb,ub,quant,mult,val[NP]; } x,y; 12 float *diag, *r; 13 float dx = 1.; 14 float ni = 100.; 15 int n; 16 int auta; 17 int periodic; 18 float konst = 0.0; 19 float zero = 0.; 20 21 /* Spline fit technique 22 let x,y be vectors of abscissas and ordinates 23 h be vector of differences hi=xi-xi-1 24 y" be vector of 2nd derivs of approx function 25 If the points are numbered 0,1,2,...,n+1 then y" satisfies 26 (R W Hamming, Numerical Methods for Engineers and Scientists, 27 2nd Ed, p349ff) 28 hiy"i-1+2(hi+hi+1)y"i+hi+1y"i+1 29 30 = 6[(yi+1-yi)/hi+1-(yi-yi-1)/hi] i=1,2,...,n 31 32 where y"0 = y"n+1 = 0 33 This is a symmetric tridiagonal system of the form 34 35 | a1 h2 | |y"1| |b1| 36 | h2 a2 h3 | |y"2| |b2| 37 | h3 a3 h4 | |y"3| = |b3| 38 | . | | .| | .| 39 | . | | .| | .| 40 It can be triangularized into 41 | d1 h2 | |y"1| |r1| 42 | d2 h3 | |y"2| |r2| 43 | d3 h4 | |y"3| = |r3| 44 | . | | .| | .| 45 | . | | .| | .| 46 where 47 d1 = a1 48 49 r0 = 0 50 51 di = ai - hi2/di-1 1<i<_n 52 53 ri = bi - hiri-1/di-1i 1<_i<_n 54 55 the back solution is 56 y"n = rn/dn 57 58 y"i = (ri-hi+1y"i+1)/di 1<_i<n 59 60 superficially, di and ri don't have to be stored for they can be 61 recalculated backward by the formulas 62 63 di-1 = hi2/(ai-di) 1<i<_n 64 65 ri-1 = (bi-ri)di-1/hi 1<i<_n 66 67 unhappily it turns out that the recursion forward for d 68 is quite strongly geometrically convergent--and is wildly 69 unstable going backward. 70 There's similar trouble with r, so the intermediate 71 results must be kept. 72 73 Note that n-1 in the program below plays the role of n+1 in the theory 74 75 Other boundary conditions_________________________ 76 77 The boundary conditions are easily generalized to handle 78 79 y0" = ky1", yn+1" = kyn" 80 81 for some constant k. The above analysis was for k = 0; 82 k = 1 fits parabolas perfectly as well as stright lines; 83 k = 1/2 has been recommended as somehow pleasant. 84 85 All that is necessary is to add h1 to a1 and hn+1 to an. 86 87 88 Periodic case_____________ 89 90 To do this, add 1 more row and column thus 91 92 | a1 h2 h1 | |y1"| |b1| 93 | h2 a2 h3 | |y2"| |b2| 94 | h3 a4 h4 | |y3"| |b3| 95 | | | .| = | .| 96 | . | | .| | .| 97 | h1 h0 a0 | | .| | .| 98 99 where h0=_ hn+1 100 101 The same diagonalization procedure works, except for 102 the effect of the 2 corner elements. Let si be the part 103 of the last element in the ith "diagonalized" row that 104 arises from the extra top corner element. 105 106 s1 = h1 107 108 si = -si-1hi/di-1 2<_i<_n+1 109 110 After "diagonalizing", the lower corner element remains. 111 Call ti the bottom element that appears in the ith colomn 112 as the bottom element to its left is eliminated 113 114 t1 = h1 115 116 ti = -ti-1hi/di-1 117 118 Evidently ti = si. 119 Elimination along the bottom row 120 introduces further corrections to the bottom right element 121 and to the last element of the right hand side. 122 Call these corrections u and v. 123 124 u1 = v1 = 0 125 126 ui = ui-1-si-1*ti-1/di-1 127 128 vi = vi-1-ri-1*ti-1/di-1 2<_i<_n+1 129 130 The back solution is now obtained as follows 131 132 y"n+1 = (rn+1+vn+1)/(dn+1+sn+1+tn+1+un+1) 133 134 y"i = (ri-hi+1*yi+1-si*yn+1)/di 1<_i<_n 135 136 Interpolation in the interval xi<_x<_xi+1 is by the formula 137 138 y = yix+ + yi+1x- -(h2i+1/6)[y"i(x+-x+3)+y"i+1(x--x-3)] 139 where 140 x+ = xi+1-x 141 142 x- = x-xi 143 */ 144 145 float 146 rhs(i){ 147 int i_; 148 double zz; 149 i_ = i==n-1?0:i; 150 zz = (y.val[i]-y.val[i-1])/(x.val[i]-x.val[i-1]); 151 return(6*((y.val[i_+1]-y.val[i_])/(x.val[i+1]-x.val[i]) - zz)); 152 } 153 154 spline(){ 155 float d,s,u,v,hi,hi1; 156 float h; 157 float D2yi,D2yi1,D2yn1,x0,x1,yy,a; 158 int end; 159 float corr; 160 int i,j,m; 161 if(n<3) return(0); 162 if(periodic) konst = 0; 163 d = 1; 164 r[0] = 0; 165 s = periodic?-1:0; 166 for(i=0;++i<n-!periodic;){ /* triangularize */ 167 hi = x.val[i]-x.val[i-1]; 168 hi1 = i==n-1?x.val[1]-x.val[0]: 169 x.val[i+1]-x.val[i]; 170 if(hi1*hi<=0) return(0); 171 u = i==1?zero:u-s*s/d; 172 v = i==1?zero:v-s*r[i-1]/d; 173 r[i] = rhs(i)-hi*r[i-1]/d; 174 s = -hi*s/d; 175 a = 2*(hi+hi1); 176 if(i==1) a += konst*hi; 177 if(i==n-2) a += konst*hi1; 178 diag[i] = d = i==1? a: 179 a - hi*hi/d; 180 } 181 D2yi = D2yn1 = 0; 182 for(i=n-!periodic;--i>=0;){ /* back substitute */ 183 end = i==n-1; 184 hi1 = end?x.val[1]-x.val[0]: 185 x.val[i+1]-x.val[i]; 186 D2yi1 = D2yi; 187 if(i>0){ 188 hi = x.val[i]-x.val[i-1]; 189 corr = end?2*s+u:zero; 190 D2yi = (end*v+r[i]-hi1*D2yi1-s*D2yn1)/ 191 (diag[i]+corr); 192 if(end) D2yn1 = D2yi; 193 if(i>1){ 194 a = 2*(hi+hi1); 195 if(i==1) a += konst*hi; 196 if(i==n-2) a += konst*hi1; 197 d = diag[i-1]; 198 s = -s*d/hi; 199 }} 200 else D2yi = D2yn1; 201 if(!periodic) { 202 if(i==0) D2yi = konst*D2yi1; 203 if(i==n-2) D2yi1 = konst*D2yi; 204 } 205 if(end) continue; 206 m = hi1>0?ni:-ni; 207 m = 1.001*m*hi1/(x.ub-x.lb); 208 if(m<=0) m = 1; 209 h = hi1/m; 210 for(j=m;j>0||i==0&&j==0;j--){ /* interpolate */ 211 x0 = (m-j)*h/hi1; 212 x1 = j*h/hi1; 213 yy = D2yi*(x0-x0*x0*x0)+D2yi1*(x1-x1*x1*x1); 214 yy = y.val[i]*x0+y.val[i+1]*x1 -hi1*hi1*yy/6; 215 printf("%f ",x.val[i]+j*h); 216 printf("%f\n",yy); 217 } 218 } 219 return(1); 220 } 221 readin() { 222 for(n=0;n<NP;n++){ 223 if(auta) x.val[n] = n*dx+x.lb; 224 else if(!getfloat(&x.val[n])) break; 225 if(!getfloat(&y.val[n])) break; } } 226 227 getfloat(p) 228 float *p;{ 229 char buf[30]; 230 register c; 231 int i; 232 extern double atof(); 233 for(;;){ 234 c = getchar(); 235 if (c==EOF) { 236 *buf = '\0'; 237 return(0); 238 } 239 *buf = c; 240 switch(*buf){ 241 case ' ': 242 case '\t': 243 case '\n': 244 continue;} 245 break;} 246 for(i=1;i<30;i++){ 247 c = getchar(); 248 if (c==EOF) { 249 buf[i] = '\0'; 250 break; 251 } 252 buf[i] = c; 253 if('0'<=c && c<='9') continue; 254 switch(c) { 255 case '.': 256 case '+': 257 case '-': 258 case 'E': 259 case 'e': 260 continue;} 261 break; } 262 buf[i] = ' '; 263 *p = atof(buf); 264 return(1); } 265 266 getlim(p) 267 struct proj *p; { 268 int i; 269 for(i=0;i<n;i++) { 270 if(!p->lbf && p->lb>(p->val[i])) p->lb = p->val[i]; 271 if(!p->ubf && p->ub<(p->val[i])) p->ub = p->val[i]; } 272 } 273 274 275 main(argc,argv) 276 char *argv[];{ 277 extern char *malloc(); 278 int i; 279 x.lbf = x.ubf = y.lbf = y.ubf = 0; 280 x.lb = INF; 281 x.ub = -INF; 282 y.lb = INF; 283 y.ub = -INF; 284 while(--argc > 0) { 285 argv++; 286 again: switch(argv[0][0]) { 287 case '-': 288 argv[0]++; 289 goto again; 290 case 'a': 291 auta = 1; 292 numb(&dx,&argc,&argv); 293 break; 294 case 'k': 295 numb(&konst,&argc,&argv); 296 break; 297 case 'n': 298 numb(&ni,&argc,&argv); 299 break; 300 case 'p': 301 periodic = 1; 302 break; 303 case 'x': 304 if(!numb(&x.lb,&argc,&argv)) break; 305 x.lbf = 1; 306 if(!numb(&x.ub,&argc,&argv)) break; 307 x.ubf = 1; 308 break; 309 default: 310 fprintf(stderr, "Bad agrument\n"); 311 exit(1); 312 } 313 } 314 if(auta&&!x.lbf) x.lb = 0; 315 readin(); 316 getlim(&x); 317 getlim(&y); 318 i = (n+1)*sizeof(dx); 319 diag = (float *)malloc((unsigned)i); 320 r = (float *)malloc((unsigned)i); 321 if(r==NULL||!spline()) for(i=0;i<n;i++){ 322 printf("%f ",x.val[i]); 323 printf("%f\n",y.val[i]); } 324 exit(0); 325 } 326 numb(np,argcp,argvp) 327 int *argcp; 328 float *np; 329 char ***argvp;{ 330 double atof(); 331 char c; 332 if(*argcp<=1) return(0); 333 c = (*argvp)[1][0]; 334 if(!('0'<=c&&c<='9' || c=='-' || c== '.' )) return(0); 335 *np = atof((*argvp)[1]); 336 (*argcp)--; 337 (*argvp)++; 338 return(1); } 339