src/nr/romberg.cpp

Go to the documentation of this file.
00001 //
00002 // C++ Interface: romberg
00003 //
00004 // Description: 
00005 //
00006 //
00007 // Author: Roman Geus <geus@maxwell>, (C) 2004
00008 //
00009 // Copyright: See COPYING file that comes with this distribution
00010 //
00011 //
00012 
00013 #include "romberg.h"
00014 
00015 namespace NR {
00016 
00017 using namespace std;
00018 
00019 void polint(const vector<double>& xa, const vector<double>& ya, const double x, double &y, double &dy) {
00020     int i, m, ns=0;
00021     double den, dif, dift, ho, hp, w;
00022 
00023     int n = xa.size();
00024     vector<double> c(n), d(n);
00025     dif = fabs(x - xa[0]);
00026     for (i = 0; i < n; ++ i) {
00027         if ((dift = fabs(x - xa[i])) < dif) {
00028             ns = i;
00029             dif = dift;
00030         }
00031         c[i] = ya[i];
00032         d[i] = ya[i];
00033     }
00034     y = ya[ns--];
00035     for (m = 1; m < n; ++ m) {
00036         for (i = 0; i < n-m; i ++) {
00037             ho = xa[i] - x;
00038             hp = xa[i+m] - x;
00039             w = c[i+1] - d[i];
00040             den = ho - hp;
00041             assert(den != 0.0);
00042             den = w/den;
00043             d[i] = hp*den;
00044             c[i] = ho*den;
00045         }
00046         y += (dy = (2*(ns+1) < (n-m) ? c[ns+1] : d[ns--]));
00047     }
00048 }
00049 
00050 } // namespace NR

Generated on Fri Oct 26 13:35:12 2007 for FEMAXX (Finite Element Maxwell Eigensolver) by  doxygen 1.4.7