Метод дихотомии для решения нелинейных уравнений

Рассмотрим простейший метод уточнения значения корня с заданной точностью - метод деления отрезка пополам (дихотомии или бисекций). Если определен интервал нахождения корня [a,b], то этот алгоритм состоит из:

1. Задания значений и вычисления значений функции на концах отрезка u=f(a), v=f(b).

2. Организации цикла, в котором последовательно выбранный отрезок делится пополам и осуществляется выбор того из двух отрезков, на котором функция меняет знак.

Выбор нужного отрезка можно реализовать так. Определяется середина отрезка и значение функции в этой точке w=f(x). Если произведение функций u*w < 0, то интервал [a,b] сужается справа заменой b=x, v=w, иначе - слева заменой a=x, u=w. Изобразите данные ситуации на графике и разберитесь с предлагаемым способом выбора требуемого отрезка.

Реализацию метода дихотомии можно провести в Excel. Рассмотрим методику на примере уравнения . Начальный интервал неопределенности отрезок [1, 3], заданная точность =0,001. Для нахождения приближенного корня уравнения понадобилось выполнить 19 шагов.

Таблица 1. Вспомогательная таблица для вычисления корней нелинейного уравнения методом дихотомии

Решение уравнения можно произвести в пакете Mathcad. Ниже приведена функция для вычисления корней методом дихотомии в данном пакете.

Корень уравнения с использованием данной функции будет следующим 2.094551 и достигнут за 34 шага.

На языке Pascal для решения уравнения методом дихотомии необходимо использовать несколько подпрограмм-функций

function X_Dich(a,b,eps:real):real;

var u,v,wx:real;

begin

u:=f(a); v:=f(b);

repeat

x:=0.5*(a+b); w:=f(x);

if (u*w<0.0)

then begin b:=x; v:=w end

else begin a:=x; u:=w end;

until abs(w) < eps;

X_Dich:=x

end;

Здесь цикл заканчивается при условии, что функция f(x) попадает «в полосу шума» .

Выбор нужного отрезка можно осуществить и с помощью знаковой функции Sign(x): равной единице, если x>0, минус единице, если x<0 и нулю, если x=0.

function Sign(x:real):integer;

begin

Sign:=0;

If x<0.0 then Sign:=-1;

If x>0.0 then Sign:= 1;

end;

При уменьшении интервала [a,b] точка a всегда остается слева от искомого корня и знак функции f(a) не изменяется. Поэтому в процессе половинного деления можно сравнивать его со знаком функции f(x). При их совподении точку a передвигаем в середину отрезка, в противном случае передвигаем точку b.

function X_Dich(a,b,eps:real):real;

var s:integer; w:real;

begin s:=Sign(f(a));

repeat

x:=0.5*(a+b); w:= f(x);

if Sign(w)=s

then a:=x

else b:=x;

until abs(w) < eps;

X_Dich:=x

end;

Использование знаковой функции наглядно показывает, что в методе половинного деления поведение функции f(x) употребляется пассивно, оно влияет только на выбор нужного интервала.

Изучите приведенные алгоритмы, проведите отладку программы и сделайте её тестирование. В качестве теста используйте любое уравнение с известным решением. В частности, можно рекомендовать уравнение Валлиса , имеющего один вещественный корень 2.09455.

Подпрограмму можно улучшить, вводя новый тип

Type fun=function(x:real):real; {$F+}

и изменив заголовок функции

function X_Dich(a,b,eps:real;f:fun):real;

Теперь можно использовать названия других функций левой части урвнения (3.1), не делая изменений в подпрограмме X_Dich.

Рекомендуется обезопасить подпрограмму от возможного «зацикливания», установив счетчик и соответствующие условия, например так:

if k>kmax

then begin writeln(‘Error01’); Exit end;

Введение счетчика полезно так же и при проведении вычислительных экспериментов по определению числа итераций, требуемых для достижения различной точности, наример, , изучении условий остановки цикла (3.4,3.5), сравнении различных методов уточнения корня и др.

Полезно в своей библиотеке иметь не только подпрограммы-функции, но и процедуры, составленные на основе приведенных алгоритмов. Напишите их, выполните отладку и тестирование.

Приведем рекурсивный вариант метода дихотомии:

function X_Dich(a,b,eps:real):real;

var x:real;

begin

if keypressed then Halt;

if f(a)*f(b) > 0.0

then begin writeln(Error02); Exit end

else begin x:=0.5*(a+b);

if abs(f(x)) > eps

then if f(a)*f(x) <0.0

then X_Dich:=X_Dich(a,x,eps)

else X_Dich:=X_Dich(x,b,eps)

else X_Dich:=x

end

end;

Изучите его, проведите отладку и тестирование программы, сделайте сравнение.

Выбор очередной точки в середине отрезка не является единственным вариантом. Можно в качестве такой точки выбрать случайное число, заменив оператор x:=0.5*(a+b) на

x:=a+(b-a)*random, предварительно инициализируя датчик случайных чисел Randomise. Проведите соответствующие расчеты и сравните требуемое число итераций для достижения заданной точности.

Более совершенный метод выбора точки деления отрезка [a,b] – метод хорд, в котором в качестве x выбирается точка пересечения с осью абсцис прямой y=Ax+B (хорды), проведенной через концы интервала u=f(a) и v=f(b).

a b x

u

Рис. 3.2. Графическая иллюстрация метода хорд

Из рисунка видно, что

, где . (3.9)

Описанные методы являются линейно сходящимися или, как говорят, сходящимися со скоростью геометрической прогрессии. В самом деле, абсолютные погрешности связаны соотношением , где знаменатель С=0.5. Метод половинного деления имеет среднюю скорость сходимости равную ln2, в то время как метод хорд в зависимости от свойств функции может иметь как меньшую, так и большую среднюю скорость сходимости. Рекомендуется исследовать этот вопрос экспериментально.

Не нашли, что искали? Воспользуйтесь поиском: