Этот алгоритм вычисления вещественного корня произвольной функции был взят из книги "Машинные методы математических вычислений" (Дж.Форсайт, М.Малькольм, К.Моулер) и переписан с небольшими изменениями с FORTRANа на С++. class MathFunc1 { public: virtual double operator () ( double ) const = 0; }; inline int sign ( double d ) { if ( d > 0 ) return 1; if ( d < 0 ) return -1; return 0; } bool zeroin ( double ax, double bx, const MathFunc1 & func, double tol, double & res ) { double fa = func(ax); double fb = func(bx); if ( sign(fa) * sign(fb) > 0 ) return false; if ( tol < 0 ) tol = 0.; // Присвоение начальных значений double a = ax; double b = bx; for(;;) { double c = a; double fc = fa; double d = b - a; double e = d; m30: if ( fabs(fc) < fabs(fb) ) { a = b; b = c; c = a; fa = fb; fb = fc; fc = fa; } // Проверка сходимости if ( fb == 0 ) break; const double tol1 = 2. * macheps * fabs(b) + 0.5 * tol; const double xm = 0.5 * ( c - b ); if ( fabs(xm) <= tol1 ) break; // Необходима ли бисекция? if ( fabs(e) < tol1 || fabs(fa) <= fabs(fb) ) goto m70; { double p, q; const double s = fb / fa; if ( a == c ) { // Линейная интерполяция p = 2. * xm * s; q = 1. - s; } else { // Обратная квадратичная интерполяция q = fa / fc; const double r = fb / fc; p = s * ( 2. * xm * q * ( q - r ) - ( b - a ) * ( r - 1. ) ); q = ( q - 1. ) * ( r - 1. ) * ( s - 1. ); } if ( p > 0 ) q = -q; p = fabs ( p ); // Приемлема ли интерполяция if ( p + p >= 3.*xm*q - fabs(tol1*q) ) goto m70; if ( p >= fabs(0.5*e*q) ) goto m70; e = d; d = p / q; goto m80; } m70: e = d = xm; // бисекция // Завершить шаг m80: a = b; fa = fb; if ( fabs(d) > tol1 ) b += d; else { if ( xm > 0 ) b += tol1; else b -= tol1; } fb = func(b); if ( sign(fb) * sign(fc) <= 0 ) goto m30; } res = b; return true; }
Параметры ax и bx функции zeroin задают интервал в котором ищется нуль и в них значения функции func
должны иметь разные знаки. Параметр tol задаёт допустимую погрешность результата res.
Функция zeroin выполняет итерационный процесс, в котором на каждом шаге присутствуют три абсциссы a, b, c.
b - последнее и наилучшее приближение к нулю, a - предыдущее приближение, c - предыдущее или ещё более
раннее приближение, такое, что func(b) и func(c) имеют разные знаки. Во всех случаях b и c ограничивают нуль.
Кроме того, |func(b)| <= |func(c)|. Если длина интервала |b - c| уменьшилась настолько, что выполняется
условие |b - c| <= tol + 4.*macheps*|b|, то итерации прерываются. Кроме tol в проверке участвует macheps, чтобы
подстраховаться на случай, когда заданное значение tol слишком мало. В частности, чтобы заставить zeroin
найти наименьший возможный интервал, нужно задать tol равным нулю.
Исходники находятся в файле mathem.cpp. Наверх |