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