![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to spline.c CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [Research Unix] / researchv10no / cmd|
Return to spline.c CVS log | Up to [Research Unix] / researchv10no / cmd |
1.1 root 1: #include <stdio.h>
2: #include <ctype.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-18
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: {
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: {
153: float d,s,u,v,hi,hi1;
154: float h;
155: float D2yi,D2yi1,D2yn1,x0,x1,yy,a;
156: int end;
157: float corr;
158: int i,j,m;
159: if(n<3) return(0);
160: if(periodic) konst = 0;
161: d = 1;
162: r[0] = 0;
163: s = periodic?-1:0;
164: for(i=0;++i<n-!periodic;){ /* triangularize */
165: hi = x.val[i]-x.val[i-1];
166: hi1 = i==n-1?x.val[1]-x.val[0]:
167: x.val[i+1]-x.val[i];
168: if(hi1*hi<=0) return(0);
169: u = i==1?zero:u-s*s/d;
170: v = i==1?zero:v-s*r[i-1]/d;
171: r[i] = rhs(i)-hi*r[i-1]/d;
172: s = -hi*s/d;
173: a = 2*(hi+hi1);
174: if(i==1) a += konst*hi;
175: if(i==n-2) a += konst*hi1;
176: diag[i] = d = i==1? a:
177: a - hi*hi/d;
178: }
179: D2yi = D2yn1 = 0;
180: for(i=n-!periodic;--i>=0;){ /* back substitute */
181: end = i==n-1;
182: hi1 = end?x.val[1]-x.val[0]:
183: x.val[i+1]-x.val[i];
184: D2yi1 = D2yi;
185: if(i>0){
186: hi = x.val[i]-x.val[i-1];
187: corr = end?2*s+u:zero;
188: D2yi = (r[i]-hi1*D2yi1-s*D2yn1+end*v)/
189: (diag[i]+corr);
190: if(end) D2yn1 = D2yi;
191: if(i>1){
192: a = 2*(hi+hi1);
193: if(i==1) a += konst*hi;
194: if(i==n-2) a += konst*hi1;
195: d = diag[i-1];
196: s = -s*d/hi;
197: }
198: }
199: else D2yi = D2yn1;
200: if(!periodic) {
201: if(i==0) D2yi = konst*D2yi1;
202: if(i==n-2) D2yi1 = konst*D2yi;
203: }
204: if(end) continue;
205: m = hi1>0?ni:-ni;
206: m = 1.001*m*hi1/(x.ub-x.lb);
207: if(m<=0) m = 1;
208: h = hi1/m;
209: for(j=m;j>0||i==0&&j==0;j--){ /* interpolate */
210: x0 = (m-j)*h/hi1;
211: x1 = j*h/hi1;
212: yy = D2yi*(x0-x0*x0*x0)+D2yi1*(x1-x1*x1*x1);
213: yy = y.val[i]*x0+y.val[i+1]*x1 -hi1*hi1*yy/6;
214: printf("%f ",x.val[i]+j*h);
215: printf("%f\n",yy);
216: }
217: }
218: return(1);
219: }
220: readin()
221: {
222: for(n=0;n<NP;n++){
223: if(auta) x.val[n] = n*dx+x.lb;
224: else if(scanf("%f",&x.val[n])!=1)
225: break;
226: if(scanf("%f",&y.val[n])!=1)
227: break;
228: }
229: }
230:
231: getlim(p)
232: struct proj *p;
233: {
234: int i;
235: for(i=0;i<n;i++) {
236: if(!p->lbf && p->lb>(p->val[i])) p->lb = p->val[i];
237: if(!p->ubf && p->ub<(p->val[i])) p->ub = p->val[i];
238: }
239: }
240:
241:
242: main(argc,argv)
243: char *argv[];
244: {
245: extern char *malloc();
246: int i;
247: x.lbf = x.ubf = y.lbf = y.ubf = 0;
248: x.lb = INF;
249: x.ub = -INF;
250: y.lb = INF;
251: y.ub = -INF;
252: while(--argc > 0) {
253: argv++;
254: again: switch(argv[0][0]) {
255: case '-':
256: argv[0]++;
257: goto again;
258: case 'a':
259: auta = 1;
260: numb(&dx,&argc,&argv);
261: break;
262: case 'k':
263: numb(&konst,&argc,&argv);
264: break;
265: case 'n':
266: numb(&ni,&argc,&argv);
267: break;
268: case 'p':
269: periodic = 1;
270: break;
271: case 'x':
272: if(!numb(&x.lb,&argc,&argv)) break;
273: x.lbf = 1;
274: if(!numb(&x.ub,&argc,&argv)) break;
275: x.ubf = 1;
276: break;
277: default:
278: fprintf(stderr, "spline: bad argument\n");
279: exit(1);
280: }
281: }
282: if(auta&&!x.lbf) x.lb = 0;
283: readin();
284: getlim(&x);
285: getlim(&y);
286: i = (n+1)*sizeof(dx);
287: diag = (float *)malloc((unsigned)i);
288: r = (float *)malloc((unsigned)i);
289: if(r==NULL||!spline()) for(i=0;i<n;i++){
290: printf("%f ",x.val[i]);
291: printf("%f\n",y.val[i]);
292: }
293: }
294:
295: numb(np, argcp, argvp)
296: int *argcp;
297: float *np;
298: register char ***argvp;
299: {
300: register char c;
301: extern double atof();
302:
303: if(*argcp <= 1)
304: return(0);
305: while((c=(*argvp)[1][0]) == '+')
306: (*argvp)[1]++;
307: if(!(isdigit(c) || c=='-'&&(*argvp)[1][1]<'A' || c=='.'))
308: return(0);
309: *np = atof((*argvp)[1]);
310: (*argcp)--;
311: (*argvp)++;
312: return(1);
313: }
314:
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.