00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013 #include "CubicSpline.h"
00014
00015 CubicSpline::CubicSpline(){;}
00016
00017
00018 void CubicSpline::operator=( const CubicSpline & rhs )
00019 {
00020 yp0=rhs.yp0;
00021 ypend=rhs.ypend;
00022 n=rhs.n;
00023 f=rhs.f;
00024 dfdx=rhs.dfdx;
00025 ddfddx=rhs.ddfddx;
00026
00027 x_list=rhs.x_list;
00028 y_list=rhs.y_list;
00029 yp_list=rhs.yp_list;
00030 a0_list=rhs.a0_list;
00031 a1_list=rhs.a1_list;
00032 a2_list=rhs.a2_list;
00033 a3_list=rhs.a3_list;
00034 }
00035
00036
00037 void CubicSpline::initializeWithFunctionValues(const Array1D<double> &x_input,
00038 const Array1D<double> &y_input,
00039 double y_prime_0, double y_prime_end)
00040 {
00041
00042
00043 n=x_input.dim1();
00044
00045 x_list.allocate(n);
00046 y_list.allocate(n);
00047 yp_list.allocate(n);
00048 a0_list.allocate(n-1);
00049 a1_list.allocate(n-1);
00050 a2_list.allocate(n-1);
00051 a3_list.allocate(n-1);
00052
00053 x_list = x_input;
00054 y_list = y_input;
00055 yp_list = 0.0;
00056 a0_list = 0.0;
00057 a1_list = 0.0;
00058 a2_list = 0.0;
00059 a3_list = 0.0;
00060
00061
00062
00063 yp0 = y_prime_0;
00064 ypend = y_prime_end;
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075 double h,d,h2,d2;
00076 row_right.allocate(n);
00077 row_middle.allocate(n);
00078 row_left.allocate(n);
00079 col_rhs.allocate(n);
00080
00081 for(int i=1;i<n-1;i++)
00082 {
00083 h=x_list(i)-x_list(i-1);
00084 d=x_list(i+1)-x_list(i);
00085
00086 row_right(i) = 1/d;
00087 row_middle(i) = 2/h+2/d;
00088 row_left(i) = 1/h;
00089 }
00090
00091
00092
00093 row_right(0) = 0;
00094 row_middle(0) = 1;
00095 row_left(0) = 0;
00096 row_right(n-1) = 0;
00097 row_middle(n-1) = 1;
00098 row_left(n-1) = 0;
00099
00100
00101
00102 for(int i=1;i<n-1;i++)
00103 {
00104 h=x_list(i)-x_list(i-1);
00105 d=x_list(i+1)-x_list(i);
00106 h2=h*h;
00107 d2=d*d;
00108
00109 col_rhs(i) = 3.0*( y_list(i+1)/d2
00110 + y_list(i) *(1/h2-1/d2)
00111 - y_list(i-1)/h2);
00112 }
00113
00114
00115
00116 col_rhs(0) = yp0;
00117 col_rhs(n-1) = ypend;
00118
00119
00120
00121 solve_tridiagonal_system();
00122
00123
00124
00125 for(int i=0;i<n;i++)
00126 {
00127 yp_list(i)=col_rhs(i);
00128 }
00129
00130
00131
00132
00133
00134
00135
00136
00137
00138
00139
00140 double fj,fj1,dfj,dfj1;
00141
00142 for(int i=0;i<n-1;i++)
00143 {
00144 d = x_list(i+1)-x_list(i);
00145 fj = y_list(i);
00146 fj1 = y_list(i+1);
00147 dfj = yp_list(i);
00148 dfj1 = yp_list(i+1);
00149 a0_list(i)=fj;
00150 a1_list(i)=dfj;
00151 a2_list(i)=3/d/d*(fj1-fj)-1/d*(dfj1+2*dfj);
00152 a3_list(i)=2/d/d/d*(fj-fj1)+1/d/d*(dfj+dfj1);
00153 }
00154 }
00155
00156
00157
00158
00159 void CubicSpline::initializeWithDerivativeValues(const Array1D<double> &x_input,
00160 const Array1D<double> &yp_input,
00161 double y_0)
00162 {
00163
00164
00165 n=x_input.dim1();
00166
00167 x_list.allocate(n);
00168 y_list.allocate(n);
00169 yp_list.allocate(n);
00170 a0_list.allocate(n-1);
00171 a1_list.allocate(n-1);
00172 a2_list.allocate(n-1);
00173 a3_list.allocate(n-1);
00174
00175 x_list = x_input;
00176 y_list = 0.0;
00177 yp_list = yp_input;
00178 a0_list = 0.0;
00179 a1_list = 0.0;
00180 a2_list = 0.0;
00181 a3_list = 0.0;
00182
00183
00184
00185 y0 = y_0;
00186
00187
00188
00189
00190
00191
00192
00193
00194
00195
00196 double h,d,h2,d2;
00197 row_right.allocate(n);
00198 row_middle.allocate(n);
00199 row_left.allocate(n);
00200 col_rhs.allocate(n);
00201
00202 for(int i=1;i<n-1;i++)
00203 {
00204 h=x_list(i)-x_list(i-1);
00205 d=x_list(i+1)-x_list(i);
00206 h2=h*h;
00207 d2=d*d;
00208
00209 row_right(i) = 3/d2;
00210 row_middle(i) = 3/h2-3/d2;
00211 row_left(i) = -3/h2;
00212 }
00213
00214
00215
00216 row_right(0) = 0;
00217 row_middle(0) = 1;
00218 row_left(0) = 0;
00219 row_right(n-1) = 0;
00220 row_middle(n-1) = 1;
00221 row_left(n-1) = -1;
00222
00223
00224
00225 for(int i=1;i<n-1;i++)
00226 {
00227 h=x_list(i)-x_list(i-1);
00228 d=x_list(i+1)-x_list(i);
00229
00230 col_rhs(i) = 3.0*( yp_list(i+1)/d
00231 + yp_list(i) *(2/h+2/d)
00232 + yp_list(i-1)/h);
00233 }
00234
00235
00236
00237 col_rhs(0) = y0;
00238 col_rhs(n-1) = 0.0;
00239
00240
00241
00242 solve_tridiagonal_system2();
00243
00244
00245
00246 for(int i=0;i<n;i++)
00247 {
00248 y_list(i)=col_rhs(i);
00249 }
00250
00251
00252
00253
00254
00255
00256
00257
00258
00259
00260
00261 double fj,fj1,dfj,dfj1;
00262
00263 for(int i=0;i<n-1;i++)
00264 {
00265 d = x_list(i+1)-x_list(i);
00266 fj = y_list(i);
00267 fj1 = y_list(i+1);
00268 dfj = yp_list(i);
00269 dfj1 = yp_list(i+1);
00270 a0_list(i)=fj;
00271 a1_list(i)=dfj;
00272 a2_list(i)=3/d/d*(fj1-fj)-1/d*(dfj1+2*dfj);
00273 a3_list(i)=2/d/d/d*(fj-fj1)+1/d/d*(dfj+dfj1);
00274 }
00275 }
00276
00277
00278 void CubicSpline::solve_tridiagonal_system()
00279 {
00280
00281
00282
00283 double factor;
00284
00285 for(int i=0;i<=(n-2);i++)
00286 {
00287 factor = row_left(i+1)/row_middle(i);
00288 row_middle(i+1) = row_middle(i+1) -
00289 row_right(i)*factor;
00290 col_rhs(i+1) = col_rhs(i+1) - col_rhs(i)*factor;
00291 row_left(i+1) = 0.0;
00292 }
00293
00294
00295
00296
00297 for(int i=n-1;i>=1;i--)
00298 {
00299 factor = row_right(i-1)/row_middle(i);
00300 col_rhs(i-1) = col_rhs(i-1) - col_rhs(i)*factor;
00301 row_right(i-1) = 0.0;
00302 }
00303
00304
00305
00306 for(int i=0;i<n;i++)
00307 {
00308 col_rhs(i) = col_rhs(i)/row_middle(i);
00309 row_middle(i) = 1.0;
00310 }
00311 }
00312
00313
00314
00315
00316 void CubicSpline::solve_tridiagonal_system2()
00317 {
00318
00319
00320
00321 double factor,b,c,d,e,f,g;
00322
00323 for(int i=n-1;i>=1;i--)
00324 {
00325 b=row_middle(i-1);
00326 c=row_right(i-1);
00327 d=col_rhs(i-1);
00328 e=row_left(i);
00329 f=row_middle(i);
00330 g=col_rhs(i);
00331
00332 factor = c/f;
00333 row_middle(i-1) = b - e*factor;
00334 col_rhs(i-1) = d - g*factor;
00335 row_right(i-1) = 0.0;
00336 }
00337
00338
00339
00340
00341 for(int i=0;i<=n-2;i++)
00342 {
00343 b=row_middle(i);
00344 d=col_rhs(i);
00345 e=row_left(i+1);
00346 f=row_middle(i+1);
00347 g=col_rhs(i+1);
00348 factor = e/b;
00349
00350 col_rhs(i+1) = g-d*factor;
00351 row_left(i+1) = 0.0;
00352 }
00353
00354
00355
00356 for(int i=0;i<n;i++)
00357 {
00358 col_rhs(i) = col_rhs(i)/row_middle(i);
00359 row_middle(i) = 1.0;
00360 }
00361 }
00362
00363 void CubicSpline::evaluate(double x)
00364 {
00365
00366 int klo,khi,k;
00367
00368
00369
00370 if(x>=x_list(n-1))
00371 {
00372 x=x_list(n-1);
00373 }
00374 else if(x<=x_list(0))
00375 {
00376 x=x_list(0);
00377 }
00378
00379 klo=0;
00380 khi=n-1;
00381
00382 while (khi-klo > 1)
00383 {
00384 k=(khi+klo) >> 1;
00385 if (x_list(k) > x) khi=k;
00386 else klo=k;
00387 }
00388
00389 evaluate(x,klo);
00390 }
00391
00392 double CubicSpline::integrate(int indexStart, int indexStop)
00393 {
00394 if(indexStart<0 || indexStop>n || indexStart>indexStop)
00395 {
00396 cerr << "ERROR: bad indices to integrate spline\n";
00397 exit(1);
00398 }
00399
00400 double value = 0;
00401 double xr, xl, IL;
00402 double a0, a1, a2, a3;
00403
00404 for(int i=indexStart; i<indexStop-1; i++)
00405 {
00406
00407
00408 xr = 0.0;
00409 xl = x_list(i+1)-x_list(i);
00410 a0=a0_list(i);
00411 a1=a1_list(i);
00412 a2=a2_list(i);
00413 a3=a3_list(i);
00414
00415
00416 IL = a0*xl + a1*xl*xl/2.0 + a2*xl*xl*xl/3.0 + a3*xl*xl*xl*xl/4.0;
00417
00418
00419 value += IL;
00420 }
00421 return value;
00422 }
00423
00424 void CubicSpline::evaluate(double x, int index)
00425 {
00426 int klo = index;
00427 int khi = index + 1;
00428
00429 if(khi>n || khi<1 || klo>(n-1) || klo<0)
00430 {
00431 cerr << "Bad value for index in CubicSpline::evaluate(double x, "
00432 << "int index)" << endl;
00433 cerr << "index= " << index << endl;
00434 cerr << "x= " << x << endl;
00435 exit(1);
00436 }
00437
00438 double r,r2,a0,a1,a2,a3;
00439
00440 r=x-x_list(klo);
00441 r2=r*r;
00442 a0=a0_list(klo);
00443 a1=a1_list(klo);
00444 a2=a2_list(klo);
00445 a3=a3_list(klo);
00446
00447 f = a0+a1*r+a2*r2+a3*r2*r;
00448
00449 dfdx = a1+2.0*a2*r+3.0*a3*r2;
00450
00451 ddfddx = 2.0*a2+6.0*a3*r;
00452 }
00453
00454 double CubicSpline::getFunctionValue()
00455 {
00456 return f;
00457 }
00458
00459 double CubicSpline::getFirstDerivativeValue()
00460 {
00461 return dfdx;
00462 }
00463
00464 double CubicSpline::getSecondDerivativeValue()
00465 {
00466 return ddfddx;
00467 }
00468
00469 void CubicSpline::toXML(ostream& strm)
00470 {
00471 strm << "<CubicSpline>" << endl;
00472 strm << "<x_list>" << endl;
00473 for(int i=0;i<x_list.dim1();i++)
00474 {
00475 strm << x_list(i) << "\t";
00476 }
00477 strm << "</x_list>" << endl;
00478
00479 strm << "<y_list>" << endl;
00480 for(int i=0;i<y_list.dim1();i++)
00481 {
00482 strm << y_list(i) << "\t";
00483 }
00484 strm << "</y_list>" << endl;
00485
00486 strm << "<yp_list>" << endl;
00487 for(int i=0;i<yp_list.dim1();i++)
00488 {
00489 strm << yp_list(i) << "\t";
00490 }
00491 strm << "</yp_list>" << endl;
00492
00493 strm << "<a0_list>" << endl;
00494 for(int i=0;i<a0_list.dim1();i++)
00495 {
00496 strm << a0_list(i) << "\t";
00497 }
00498 strm << "</a0_list>" << endl;
00499
00500 strm << "<a1_list>" << endl;
00501 for(int i=0;i<a1_list.dim1();i++)
00502 {
00503 strm << a1_list(i) << "\t";
00504 }
00505 strm << "</a1_list>" << endl;
00506
00507 strm << "<a2_list>" << endl;
00508 for(int i=0;i<a2_list.dim1();i++)
00509 {
00510 strm << a2_list(i) << "\t";
00511 }
00512 strm << "</a2_list>" << endl;
00513
00514 strm << "<a3_list>" << endl;
00515 for(int i=0;i<a3_list.dim1();i++)
00516 {
00517 strm << a3_list(i) << "\t";
00518 }
00519 strm << "</a3_list>" << endl;
00520
00521 strm << "<f>\t" << f << "\t</f>" << endl;
00522 strm << "<dfdx>\t" << dfdx << "\t</dfdx>" << endl;
00523 strm << "<ddfddx>\t" << ddfddx << "\t</ddfddx>" << endl;
00524 strm << "</CubicSpline>" << endl;
00525 }
00526
00527 ostream& operator<<(ostream & strm, CubicSpline & rhs)
00528 {
00529 rhs.toXML(strm);
00530 return strm;
00531 }
00532
00533