Пример 3.3 выполнения задания

Уточнение корней уравнений

Для численного решения алгебраических уравнений разработано множество итерационных методов (методов последовательного приближения к точному значению) уточнения корня. Задача ставится так: при заданном одном или двух (зависит от метода) начальных приближениях корня уравнения F(X)=0 получить приближение корня с заданной точностью ε.

Требуемая точность ε определяет условие завершения итерационного процесса, которое задаётся отношением |X_n-X_n-1| < ε, где X_n и X_n-1 – соседние приближения корня, полученные на (n-1) –м и n–м шагах его уточнения, а начальные (грубые) приближения корней можно найти, например, по результатам табулирования функции F(X).

Метод простых итераций

Для уточнения корня уравнения вида F(X)=0 его следует преобразовать к уравнению X=G(X). Исходными данными для уточнения корня являются требуемая точность ε и начальное приближение X0. Очередное приближение X1 корня вычисляется на основе текущего приближения X0 по формуле X1=G(X0) (на первом шаге уточнения корня Х0 представляет начальное приближение), после чего X0 получает значение X1 и процесс повторяется, пока модуль разности между Х0 и Х1 больше ε.

Применение метода приводит к решению, если |G'(Х)|<1 внутри интервала, содержащем корень уравнения.

Пример. Пусть известно, что при заданном начальном приближении корня X0 метод простых итераций обеспечит получение решения уравнения X=(X-0,1)⁴ +0,1. Тогда для нахождения корня с заданной точностью ε можно использовать следующий фрагмент программы.

uses

SysUtils, Math;

var

X0, X1, Eps, dX: Extended;

begin

ReadLn(X0,Eps);

repeat

X1:=IntPower(X0-0.1, 4)+0.1;

dX:=Abs(X0-X1);

X0:=X1

until dX<Eps;

WriteLn('С точностью ',Eps, ' корень равен ',X0);

............

end.

Метод не всегда обеспечивает нахождение корня. Так, при ε =10^-3 и X0<1,1, где в окрестности корня 0,1 |G'(X)| = |4(X-0,1)³ |<1 будет найден этот корень (как будет проходить уточнение корня, показано стрелками на рис.3.1). Но в окрестности корня 1,1 (1,1 – второй корень уравнения), где |G'(X)| >1 каждый шаг процесса будет приводить к удалению от корня, что, в конечном счете, приведет к переполнению разрядной сетки машины и аварийному останову.

Чтобы найти корень уравнения X=G(X) при условии |G'(X)|>1, его следует преобразовать к виду X=H(X), где H(X) - обратная относительно G(X) функция, и использовать для поиска корня. С учетом сказанного, изменим текст фрагмента программы, чтобы была возможность находить все корни уравнения X=(X-0,1)⁴ +0,1, при любых начальных приближениях, когда |G'(X)|¹0.

uses

SysUtils, Math;

var

X0, X1, Eps, dX, P: Extended;

begin

ReadLn(X0,Eps);

P:=Abs(4*IntPower(X0-0.1,3));

if P=0 then

begin

WriteLn('Нарушено условие применимости метода!');

ReadLn;

Halt;

end

else if P<1 then

//Использование исходной функции

repeat

X1:=IntPower(X0-0.1, 4)+0.1;

dX:=Abs(X0-X1);

X0:=X1

until dX<Eps

else if P>1 then

//Использование обратной функции

repeat

X1:=Power(X0-0.1,0.25)+0.1 ;

dX:=Abs(X0-X1);

X0:=X1

until dX<Eps;

WriteLn('С точностью ',Eps, ' корень равен ',X0);

. . . . . . .

end.

В некоторых случаях возможно зацикливание – бесконечное выполнение цикла программы (например, для уравнения X =1/X) или медленная сходимость процесса (например, для уравнения X = 1/(X-10^-6).

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

Пример. Составить фрагмент программы для решения следующей задачи. Найти корень уравнения X = 1/(X-10^-6) с заданной точностью ε. Если за заданное число N шагов точность не будет достигнута, то прекратить выполнение программы и вывести с пояснениями модуль разности между двумя последними приближениями, значение ε и значение N, иначе вывести найденное значение корня и число шагов, за которое оно было найдено.

var

X0, X1, Eps:Extended;

I, N:Integer;

F:Boolean;

begin

ReadLn(X1,Eps,N);

F:=False;

for I:=1 to N do

begin

X0:=X1;

X1:=1/(X0-1E-6);

F:=Abs(X0-X1)<Eps;

if F then

//Решение найдено

break;//Выход из цикла

end;

if F then

WriteLn('Корень X0 = ',X0:12,' найден за ',i,' шагов')

else

WriteLn('Заданная точность ',Eps:0,' не достигнута за '

,N,' шагов');

. . . . . . .

end.

Выполнив эту программу при X0=0,5, ε =10^-3, N=5000000 можно убедиться, что после 5000000 итераций заданная точность достигнута не будет.

Метод половинного деления

Исходными данными для уточнения корня уравнения вида F(X)=0 являются требуемая точность ε и два начальных приближения: XL и XR, между которыми должен находиться корень. Поэтому необходимым условием применения метода является истинность отношения F(XL)·F(XR)<0, то есть метод не пригоден в тех случаях, когда график F(X) лишь касается оси абсцисс, не пересекая её, например, в случае уравнения X²=0. Один шаг итерационного процесса уточнения корня состоит в перемещении правой (XR) или левой (XL) границы отрезка (XL,XR) в его середину в соответствии со следующим правилом: если знак F((XR+XL)/2) совпадает со знаком F(XL), то XL получит значение (XR+XL)/2, иначе это значение получит XR (см. рис. 3.2). Процесс повторяется, пока модуль разности между ХR и ХL больше ε.

Пример. Составить фрагмент программы уточнения корня уравнения arctg(X)-X=0 с заданной точностью ε при начальных приближениях корня ХL и ХR методом половинного деления.

ReadLn(XL,XR,Eps);

YL:=ArcTan(XL)-XL;

repeat

X:=(XL+XR)/2;

Y:= ArcTan(X)-X;

if Y*YL>0 then

XL:=X

else

XR:=X;

until Abs(XR-XL)<Eps;

WriteLn('Корень уравнения равен ', X:12);

Метод касательных

Исходными данными для уточнения корня уравнения вида F(X)=0 являются требуемая точность ε и начальное приближение X0. Необходимым условием применения метода является истинность отношения F(X0)·F''(X0)>0. Один шаг итерационного процесса уточнения корня состоит в вычислении очередного приближения по формуле X1= X0-F(X0)/F'(X0), после чего X0 получает значение X1 (см. рис. 3.3). Процесс повторяется, пока модуль разности между Х0 и Х1 больше ε.

Пример. Составить фрагмент программы уточнения корня уравнения
(X-0,1)⁴ –X + 0,1=0
с заданной точностью ε при начальном приближении корня Х0.

ReadLn(X0, Eps);

repeat

dX:=(IntPower(X0-0.1, 4)-X0+0.1)/(4*IntPower(X0-0.1, 3)-1);

X1:=X0-dX;

X0:=X1

until Abs(dX)<Eps;

WriteLn('Корень уравнения равен ', X0:12);

В этом фрагменте использовалась найденное заранее выражение 4(X-0,1)³ –1 первой производной для (X-0,1)⁴ –X + 0,1. С точки зрения объема и точности вычислений такое решение предпочтительнее использования для этих целей разностного отношения, как, например, в следующем операторе

dX1:=(IntPower(X0-0.1, 4) - X0+0.1)

/(( (IntPower(X0+1e-8-0.1, 4) - X0+1e-8+0.1)

-(IntPower(X0-0.1,4) - X0+0.1) )/1e-8);

где для вычисления приближенного хначения производной использовалась формула и ΔX=10^-8.

Пример 3.3 выполнения задания

Программа составлена по условию варианта задания №30 (см. ниже). В реализации метода касательных используются выражения производной 9Х²-10Х+1 и второй производной 18Х-10 выражения, входящего в уравнение 3Х³-5Х²+Х+0,4=0.

program Project1;

{$APPTYPE CONSOLE}

uses

SysUtils, Math;

var

X0,X1,Eps,dX,X,XL,XR,YL,Y:Extended;

i,N:Integer;

begin

//Метод касательных

WriteLn('Метод касательных');

Write('Введите X0, Eps и N: ');

ReadLn(X0,Eps,N);

//Проверка применимости метода

if (3*IntPower(X0,3)-5*Sqr(X0)+X0+0.4)*(18*X0-10) <= 0 then

WriteLn('Не выполнено условие применимости'

,' метода касательных!')

else

begin

i:=0;

repeat

i:=i+1;

dX := (3*IntPower(X0,3)-5*Sqr(X0)+X0+0.4)

/(9* Sqr(X0)-10*X0+1);

X1:=X0-dX;

X0:=X1

until (Abs(dX)<Eps) or (i=N);

if Abs(dX)<Eps then

WriteLn('Корень Х уравнения ',X0:14,' найден за ' ,i

,' шагов, Y(X) = '

,3* IntPower(X0,3)-5* Sqr(X0)+X0+0.4:14)

else

WriteLn('Корень уравнения не найден!');

end;

WriteLn;

//Метод половинного деления

WriteLn('Метод половинного деления');

Write('Введите XL, XR, Eps и N: ');

ReadLn(XL, XR, Eps, N);

//Проверка применимости метода

if (3*IntPower(XL,3)-5*Sqr(XL)+XL+0.4)

*(3*IntPower(XR,3)-5*Sqr(XR)+XR+0.4) > 0 then

WriteLn('Не выполнено условие применимости'

,' метода половинного деления!')

else

begin

YL:=3*IntPower(XL,3)-5*Sqr(XL)+XL+0.4;

i:=0;

repeat

X:=(XL+XR)/2;

Y:= 3*IntPower(X,3)-5*Sqr(X)+X+0.4;

if Y*YL>0 then

XL:=X

else

XR:=X;

inc(i);

until (Abs(XR-XL)<Eps) or (i=N);

if Abs(XR-XL)<Eps then

WriteLn('Корень уравнения ',X:14,' найден за '

,i,' шагов, Y(X) = ',Y:14)

else

WriteLn('Корень уравнения не найден');

end;

ReadLn;

end.