Нахождение локальных экстремумов функций одна из наиболее распространённых задач вычислительной математики, которые приходится решать программистам. Нахождение максимумов и минимумов на заданном интервале лежит в основе решения многих задач по оптимизации.
Один из наиболее предпочтительных способов решения этой математической задачи, использование метода Ньютона (или по-другому метода касательных), так как он обладает наиболее быстрой сходимостью.
Основы метода Ньютона
Метод Ньютона основан на геометрическом смысле производной.
Напомним, что геометрический смысл производной – тангенс угла наклона касательной к графику функции в данной точке.
Таким образом, задав некоторую точку x0 (начальное приближение), значение функции в ней (y0) и значение производной (tgα), можно вычислить координаты точки x1, в которой эта функция обращается в ноль.
Данное обстоятельство и позволяет использовать метод Ньютона для решения математических задач, В частности, нахождения корней уравнений.
На чертеже показан сильно упрощённый вариант, так как обычно для нахождения целевой точки однократного выполнения расчёта недостаточно и требуется его повторение до тех пор, пока погрешность не станет меньше некоторой критической величины. На второй итерации расчёта в качестве x0 используется вычисленное значение x1 и т.д. по формуле.
x_{n+1}=x_{0}-\frac{f'(x)}{f»(x)}
Поэтому метод Ньютона относится к числу итерационных методов.
Важно отметить, что вычисления с помощью итерационных методов всегда имеют некоторую погрешность. Однако, очень часто такие методы единственный способ решения задачи.
Поэтому на каждой итерации вычисляется погрешность.
\Delta x=|x_{n+1}-x_{n}|
Если погрешность меньше максимально допустимой (∈), полученный результат вычислений принимаем в качестве ответа (корня уравнения).
При всех своих достоинствах метод Ньютона имеет ряд недостатков, которые ограничивают его применение. Основные из них:
- Функция должна быть непрерывна и дифференцируема на всём исследуемом отрезке;
- Производная функции должна существовать и быть непрерывна на всём исследуемом интервале:
- Как все итерационные методы, метод Ньютона требует единственности решения на всём исследуемом интервале;
- Если начальное приближение выбрано неправильно (слишком далеко от предполагаемого корня или меньше его), метод может не сойтись.
Условия экстремумов
Нахождение экстремумов основано на том, что в них производная функции (производная первого порядка) равна 0. При этом в точке максимума производная от производной функции (производная второго порядка) меньше нуля, а в точке минимума больше.
Алгоритм метода Ньютона при поиске экстремумов
- Выбираем (для последующих итераций, задаём на основании расчётов) начальное значение (x0). Если вычисляем для него значение производной первого порядка;
- На основании начального значения и значения производной первого порядка вычисляем следующее значение (x1).
- Вычисляем для x1 производные первого и второго порядка
- Проверяем, выполняется ли условие экстремума и условие остановки;
- Если нет, принимаем x1 в качестве x0 и повторяем пункты 2-4.
Если да, принимаем данное значение в качестве решения задачи и определяем, является ли экстремум минимумом или максимумом.
Сама формула расчёта при поиске экстремумов будет иметь вид:
x_{n+1}=x_{n}-\frac{f'(x_n)}{f»(x_n)}
Пример задачи
Составим подпрограмму для поиска максимума и минимума на некотором интервале для следующей функции:
f(x)=cos(x)+sin(x)
Вначале продифференцируем эту функцию.
Производная первого порядка этой функции имеет вид:
f'(x)=-sin(x)+cos(x)
Производная второго порядка имеет вид:
f»(x)=-\left ( cos(x)+sin(x) \right )
Теперь, остаётся только реализовать все необходимые вычисления на языке программирования. Вычисление производных при этом лучше вынести в отдельные функции, чтобы упростить реализацию самого метода.
Ниже приведены примеры исходного кода поиска максимума. Если экстремум не является максимумом, будет возбуждено исключение.
На Delphi
Производная первого порядка:
unction Dfdx(x: double): double; begin result := -1 * sin(x) + cos(x); end;
Производная второго порядка:
function D2fdx2(x: double): double; begin result := -1 * (cos(x) + sin(x)); end;
Сам метод Ньютона:
function Newton(start: double; eps: double): double; var x0, x1, dx: double; begin x0 := start; repeat x1 := x0 - dfdx(x0) / d2fdx2(x0); dx := abs(x1 - x0); x0 := x1; until dx < eps; if (d2fdx2(x1) < 0) then result := x1 else raise Exception.Create('Экстремум на этом интервале не является максимумом!'); end;
На C#
Производная первого порядка:
private double Dfdx(double x) { return -1 * Math.Sin(x) + Math.Cos(x); }
Производная второго порядка:
private double D2fdx2(double x) { return -1 * (Math.Cos(x) + Math.Sin(x)); }
Сам метод Ньютона:
private double Newton(double start, double eps, ref bool isMax) { double x1, dx; double x0 = start; do { x1 = x0 - dfdx(x0) / d2fdx2(x0); dx = Math.Abs(x1 - x0); x0 = x1; } while (dx < eps); if (d2fdx2(x1) < 0) { retun x1; } else { throw new Exception("Экстремум на этом интервале не является максимумом!"); } }
На Java
Реализация на Java практически идентична C#. Все различия обусловлены только особенностями синтаксиса.
Производная первого порядка:
private double dfdx(double x) { return -1 * Math.sin(x) + Math.cos(x); }
Производная второго порядка:
private double d2fdx2(double x) { return -1 * (Math.cos(x) + Math.sin(x)); }
Сам метод Ньютона:
private static double newton(double start, double eps) throws Exception { double x1, dx; double x0 = start; do { x1 = x0 - dfdx(x0) / d2fdx2(x0); dx = Math.abs(x1 - x0); x0 = x1; } while (dx < eps); if (d2fdx2(x1) < 0) { return x1; } else { throw new Exception("Экстремум на этом интервале не является максимумом!"); } }
Нахождение минимума осуществляется аналогичным способом. Только условие определения экстремума изменяется на противоположное.
D2fdx2(x1) > 0