В основном код взят с сайта prografix.narod.ru.
ОГЛАВЛЕНИЕ
Основной элемент - класс MathFunc1, как раз и описывающий функцию одной действительной переменной:
class MathFunc1
{
public:
virtual double operator () (double) = 0;
};
Нельзя создавать объекты этого класса, надо обязательно создавать классы-продолжения, например:
class Polynom5 : public MathFunc1
{
public:
double a, b, c, d, e;
double operator () (double t)
{
return e + t * (d + t * (c + t * (b + t * (a + t))));
}
};
Polynom5 f; // объявляем переменную f типа Polynom5
f.f0 = 1;
f.f1 = 1;
f.f2 = f.f3 = f.f4 = 0;
double value = f(0.3) + f(1.4);
или
class f01 : public MathFunc1
{
public:
double operator () (double t)
{
return sqrt(t);
}
};
f01 g; // объявляем переменную g типа f01
double value = g(0.3) + g(1.4);
Из prografix:
Этот алгоритм вычисления вещественного корня произвольной функции был взят из книги "Машинные методы математических вычислений" (Дж.Форсайт, М.Малькольм, К.Моулер) и переписан с небольшими изменениями с FORTRANа на С++.
// Поиск нуля функции на заданном интервале // Значения функции в ax и bx должны иметь разные знаки // tol - заданная точность, // res - найденный корень bool zeroin(double ax, double bx, MathFunc1 & func, double tol, double & res);
Параметры 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 равным нулю.
На каждом шаге zeroin выбирает очередное приближение из двух кандидатов - один получен алгоритмом бисекции, а другой - алгоритмом интерполяции. Если a, b, c различны, используется обратная квадратичная интерполяция; в противном случае - линейная интерполяция ( метод секущих ). Если точка, полученная интерполяцией, "разумна", то выбирается она; иначе выбирается точка бисекции. Определение "разумности" довольно техническое, но по существу оно означает, что точка находится внутри текущего интервала и не слишком близко к его концам. Следовательно, длина интервала гарантированно убывает на каждом шаге и убывает быстро, если функция хорошо ведёт себя.
Из prografix:
Для того, чтобы найти минимум или максимум функции на заданном интервале [a,b] с точностью eps можно воспользоваться методом "золотого сечения" или комбинированным методом:
// Поиск минимума или максимума функции на заданном интервале double fmin(double a, double b, MathFunc1 & func, double eps); double fmax(double a, double b, MathFunc1 & func, double eps);
Алгоритм fmin я взял из книги "Машинные методы математических вычислений" (Дж.Форсайт, М.Малькольм, К.Моулер) и переписал с небольшими изменениями с FORTRANа на С++. Он представляет собой комбинацию метода "золотого сечения" и параболической интерполяции. Аналогично был сделан алгоритм fmax. В моих тестах этот алгоритм в среднем делал в 2 раза меньше запросов вычисления значений функции, чем метод "золотого сечения". Для разных функций этот показатель у меня колебался от 10 до 90%.
(В моей версии остался только комбинированный метод. С.Х.)
(Этого не было на сайте "ProGrafix". С.Х.)
// Численное интегрирование методом Симпсона с удвоением количества узлов на каждом шаге double Simpson(double a, double b, MathFunc1 & f, double eps);
Отрезок интегрирования разбивается на 2 части, потом на 4, потом на 8 и т.д., пока не получим требуемую точность.
Не задавайте избыточную точность!
Например, для фунции f(x)=x4
для вычисления интеграла на отрезке [0,1] с максимальной точностью потребуется разбиение на 4096 участков,
для фунции f(x)=sqrt(x) на том же отрезке – на 8 млн. участков,
Пример использования:
class f03 : public MathFunc1
{
public:
double operator () (double t)
{
return t*t*t-t+1;
}
};
...
f03 f;
double x1;
if (!zeroin(-2, 2, f, 1e-12, x1)) printf("Error\n");
double max = fmax(-1, 1, f, 1e-12);
double min = fmin(-1, 1, f, 1e-12);
double ing = Simpson(-1, 1, f, 1e-12);
printf(" Root = %10.7f, max = %10.7f, min = %10.7f, ing = %10.7f\n", x1, max, min, ing);
Root = -1.3247180, max = -0.5773503, min = 0.5773503, ing = 2.0000000