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