Этот метод в отличие от ранее рассмотренных не требует предварительно указывать интервал, в котором располагается корень уравнения. Для начала работы требуется задать лишь одну начальную точку x0, расположенную вблизи от предполагаемого корня. Направление поиска определяется из этой точки с помощью линейной экстраполяции f(x). Таким образом, при начале расчета из заданной точки x0 определяется точка x1, затем из точки x1 рассчитывается x2 и так далее. Продолжение этого процесса далее дает последовательность чисел x0, x1, x2, x3, …, xi, … последовательно приближающихся к корню уравнения. Для получения итерационной формулы метода Ньютона воспользуемся разложением функции f(x) в окрестности точки xi в ряд Тейлора:
где f'(x),if"(x) иif'"(x) –i первая, вторая и третья производные от функ-
ции f (x) по x.
Сократим (10), отбросив слагаемые, содержащие ∆x во второй и более
высоких степенях. Тогда
Полагая далее, что в окрестностях xi имеется точка xi+1 = xi + ∆x, в которой
функция f (+1x)=i(f+x∆i) xравна нулю, получим линейное уравнение
из которого найдем xi+1:
Это соотношение является итерационной формулой метода Ньютона. Алгоритм метода Ньютона представлен на рис. 7. Получаемые методом Ньютона точки xi образуют ряд чисел x0, x1, x2, x3, …, который сходится к точному решению, то есть к корню уравнения. Из (11) следует, что каждый шаг метода Ньютона требует большего объема вычислений чем, например, метод половинного деления, так как приходится находить значение не только функции f(x), но и ее производной. Несмотря на это метод Ньютона и его модификации широко используются на практике. Это обусловлено, во-первых, тем, что он не требует задания отрезка [a, b], содержащего корень, а может стартовать от одной начальной точки. Во-вторых, он имеет более высокую скорость сходимости, чем ранее рассмотренные методы. Теоретически можно показать, что метод Ньютона позволяет получить квадратичную сходимость. Это означает, что на каждой итерации погрешность (отклонение очередного приближения xi от точного решения) уменьшается по квадратичному закону, то есть количество верных значащих цифр решения удваивается.
Если на очередном шаге достигнута погрешность не более 0,5 то за пять-шесть итераций она уменьшится до величины порядка 2–64, что сопоставимо с погрешностью вычислений на ЭВМ. В методе половинного деления для достижения такой же погрешности количество итераций потребовалось бы увеличить более чем на порядок. На рис. 8, а представлен ход решения методом Ньютона в графическом виде.
При использовании метода Ньютона следует учитывать ряд его особенностей. Одна из них состоит в необходимости правильного выбора начального приближения. Чтобы понять, как влияет выбор начальной точки на работу метода, попробуйте графически найти решение для рис. 8, начав его из точки x0 = a. Метод Ньютона обладает локальной сходимостью, то есть способен найти корень, если начальное приближение задано в некоторой малой его окрестности. Если же начальное приближение взято неудачно и функция немонотонна, метод может дать расходящуюся последовательность xi (см. п. 1.5). Другая проблема заключается в том, что производная f'(x) в (11) находится в знаменателе. Это означает, что f'(x) не должна обращаться в ноль, так как в противном случае итерационная формула перестает работать. Трудности могут возникнуть и в том случае, если f'(x)не равна нулю, но достаточно мала, вследствие чего результат деления f (x)f'(x) может оказаться неприемлемо большим. Во многих математических пакетах, например, в MathCAD и MATLAB эти проблемы решаются применением комбинированных алгоритмов, сочетающих достоинства различных методов, например, метода половинного деления и метода Ньютона. Первый обеспечивает устойчивую сходимость и используется на начальном этапе решения, а после некоторого числа итераций включается второй, быстрее приближающийся к корню уравнения.
A=0.01; R=10;
f=@(L)0.14.*(0.418+0.086.*(L./R).^2)./(R.*(L.^0.5).*((1-0.086.*(L./R).^2)).^0.5)-A;
L=21:0.1:33.9;
y=f(L);
figure(1);
plot(L,y); grid on;
xlabel('L'); ylabel('y');
L=31;%начальное приближение
dL=1e-1;%малое приращение аргумента
TOL=1e-1; %точность поиска решений
for i=1:1000
fL=f(L);
df=(f(L+dL)-f(L))/dL;%производная
Li=L-f(L)/df;
if abs(Li-L)<TOL*10.^-3
break;
end
L=Li;
end
fprintf('otvet: L=%g, f(L)=%g (%i iteraciy) \n', Li,fL,i)
%стандартное решение
L1=fzero(f,32);
fprintf('стандартное решение: L=%g, f(L)=%g (%i iteraciy) /n', Li,fL,i)
otvet: L=32.2641, f(L)=7.67933e-008 (6 iteraciy)
стандартное решение: L=32.2641, f(L)=7.67933e-008 (6 iteraciy) /n>>