Здесь я предлагаю рассмотреть некоторые методы численного решения нелинейных алгебраических уравнений вида F(x)=0, т.е. задачу, с которой часто сталкиваются наши будущие светила науки и техники, а ныне студенты. Будут рассмотрены:
1. Метод половинного деления.
2. Метод хорд (секущих).
3. Метод касательных (Ньютона).
4. Метод простой итерации.
5. Метод "десятичного деления" (условное название).
Программы написаны на языке Паскаль. Для всех методов структура программ идентична: функция F(х) (а для метода Ньютона также и функция G(x)=F'(x)) задается в виде отдельной подпрограммы, что позволяет легко заменить ее требуемой, сам процесс нахождения корня оформлен в виде процедуры, а тело программы сводится к циклу по разным точностям определения корня в направлении повышения итоговой точности.
ВНИМАНИЕ – ВАЖНО!!!
Прежде, чем приступить к написанию программы вычисления корня ЛЮБЫМ методом, надлежит провести небольшое исследование, сводящееся к тому, чтобы:
1. Найти некое предварительное грубо-приближенное значение корня.
2. Выбрать интервал значений аргумента такой, чтобы он:
а) содержал внутри себя ЕДИНСТВЕННЫЙ корень;
б) (для методов хорд и касательных) функция внутри него не имела экстремумов (производная не должна менять знак), то есть, иными словами, чтобы внутри выбранного интервала функция была либо монотонно-возрастающей, либо монотонно-убывающей.
Проще всего это исследование осуществить путем построения графика функции F(x).
Договоримся обозначать границы интервала буквами
a и
b. Для функции, взятой в качестве примера, положим
a = 0.0,
b = 2.0.
Во всех рассмотренных примерах точность постепенно нарастает от 0.1 до 0.0000000001, чтобы можно было наглядно увидеть, как результат вычислений сходится к точному значению корня.
Итак, "начнем, благословясь".
1. Метод половинного деления
Его еще иногда не совсем корректно называют "методом дихотомии". На мой взгляд, метод неудобный, устаревший но, увы, страстно любимый преподавателями информатики. Суть: исходный интервал разбиваем пополам, далее определяем, в какой из половинок находится наш корень (а это проще всего узнать, рассмотрев произведение значений функции на концах половинок: в нужной оно окажется отрицательным), затем делим выбранную половинку пополам и т.д., пока длина интервала не окажется меньше установленной наперед точности. Метод устойчив.
Код:
Var
Eps,Res:Real;
i:Integer;
Function F(x:real):Real;
begin
F:=Sqr(Sin(x)+Cos(x))/Exp(Ln(33.5)*2/3)+Sqrt(3/7)-x;
end;
Procedure EQRoot(e:real; var R:real);
var
a,b,c:real;
begin
a:=0.0;
b:=2.0;
Repeat
c:=(a+b)/2;
if F(a)*F(c)<0 then b:=c else a:=c;
Until b-a<e;
R:=c;
end;
Begin
Eps:=0.1;
for i:=1 to 10 do
begin
EQRoot(Eps,Res);
Writeln(Res:12:10);
Eps:=Eps/10;
end;
Readln
End.
Результат:
Код:
0.8125000000
0.8515625000
0.8466796875
0.8463745117
0.8463973999
0.8463945389
0.8463954329
0.8463953808
0.8463953799
0.8463953790
-------------
0.8463953790 - accurate value
2. Метод хорд (секущих)
Тут так: одну из точек начального интервала, например, точку
a, "закрепляем". Проводим прямую (хорду) между точками с координатами (
a, F(a)) и (
b, F(b)). Естественно, поскольку эти две точки лежат по разные стороны оси абсцисс, наша хорда пересечет эту ось в точке
x1, и координата точки пересечения будет первым приближением решения. Далее ею мы заменим точку
b и проведем хорду между точками (
a, F(a)) и (
x1, F(x1)). Получим новую точку пересечения
x2. Это следующее приближение. И так повторяем до тех пор, пока абсолютная величина разности между предыдущим и последующим приближениями не станет меньше установленной точности. При выполнении оговоренных выше условий выбора интервала, метод устойчив.
Код:
Var
Eps,Res:Real;
i:Integer;
Function F(z:real):Real;
begin
F:=Sqr(Sin(z)+Cos(z))/Exp(Ln(33.5)*2/3)+Sqrt(3/7)-z;
end;
Procedure EQRoot(e:real; var R:real);
var
x1,x2,a,d:real;
begin
a:=0.0;
x1:=2.0;
Repeat
x2:=x1-F(x1)*(x1-a)/(F(x1)-F(a));
d:=abs(x1-x2);
x1:=x2;
Until d<e;
R:=x2;
end;
Begin
Eps:=0.1;
for i:=1 to 10 do
begin
EQRoot(Eps,Res);
Writeln(Res:12:10);
Eps:=Eps/10;
end;
Readln
End.
Результат:
Код:
0.8435500972
0.8468323673
0.8463282561
0.8464056892
0.8463956223
0.8463953417
0.8463953848
0.8463953781
0.8463953790
0.8463953790
-------------
0.8463953790 - accurate value
Полной устойчивости метода хорд (секущих) можно добиться, если точки, через которые проводится хорда, выбирать так, чтобы они всегда располагались по разные стороны оси абсцисс, т.е. заменять найденным текущим значением корня либо нижний, либо верхний предел интервала.
Код:
Var
Eps,Res:Real;
i:Integer;
Function F(z:real):Real;
begin
F:=Sqr(Sin(z)+Cos(z))/Exp(Ln(33.5)*2/3)+Sqrt(3/7)-z;
end;
Procedure EQRoot(e:real; var R:real);
var
x1,x2,a,b,d:real;
begin
a:=0.0;
b:=2.0;
x1:=b;
Repeat
x2:=b-F(b)*(b-a)/(F(b)-F(a));
d:=abs(x1-x2);
if F(a)*F(x2)<0 then b:=x2 else a:=x2;
x1:=x2;
Until d<e;
R:=x2;
end;
Begin
Eps:=0.1;
for i:=1 to 10 do
begin
EQRoot(Eps,Res);
Writeln(Res:12:10);
Eps:=Eps/10;
end;
Readln
End.
Результат:
Код:
0.8448594303
0.8462309819
0.8463778020
0.8463934999
0.8463951781
0.8463953576
0.8463953767
0.8463953788
0.8463953790
0.8463953790
-------------
0.8463953790 - accurate value
(Продолжение следует)