К вопросу о численном решении алгебраических уравнений

Vladimir_S · 09.04.2014, 12:51

Здесь я предлагаю рассмотреть некоторые методы численного решения нелинейных алгебраических уравнений вида 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

(Продолжение следует)

Vladimir_S · 09.04.2014, 14:36

(продолжение)

3. Метод касательных (Ньютона)

Начинаем с того, что выбираем стартовую точку x0. Этой точкой может быть один из концов интервала (a или b), а может быть просто точка вблизи предполагаемого корня (чем ближе, тем лучше).
(Лирическое отступление. Некоторые идиоты-"преподаватели" требуют от студентов, чтобы этой точкой непременно был один из концов интервала, причем не любой, а именно тот, где знаки функции и ее производной совпадают. Заявляю ответственно: подобные "украшательства" есть полная чушь, свидетельствующая лишь о некомпетентности этих самых...).
Далее, через точку (x0, F(x0)) проводим касательную к кривой y=F(x) до пересечения этой касательной оси абсцисс в точке х1. Таким образом, имеем прямоугольный треугольник с длинами катетов вертикального – F(x0), и горизонтального – расстояние между х0 и х1. Гипотенузой является отрезок касательной.
Если угол между касательной и осью абсцисс, отсчитанный от положительного направления последней, есть α, то из геометрии следует, что (x0-x1)=F(x0)/tg(α). Но тангенс угла наклона касательной есть ни что иное, как производная в точке х0, а потому, окончательно
x1 = x0 – F(x0)/F'(x0)
Далее берем функцию и производную в точке х1 и по той же формуле получаем х2, потом х3 и т.д. Точки х0, х1, х2, х3, ... и будут последовательными приближениями искомого корня.
К сожалению, метод не всегда устойчив. Например, если в стартовой точке производная близка к нулю, то точка х1 рискует оказаться далеко за пределами интервала и попасть в область, где функция либо не определена вовсе, либо успела пройти один или даже несколько экстремумов, и в таком случае процедура будет не приближать, а отдалять получаемые значения х1, х2 и т.д. от корня. Тем важнее становится этап предварительного исследования, о котором писалось выше.

Код:

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;

Function G(z:real):Real;
begin
 G:=(Sin(z)+Cos(z))*2*(Cos(z)-Sin(z))/Exp(Ln(33.5)*2/3)-1;
end;

Procedure EQRoot(e:real; var R:real);
var
 x1,x2,d:real;
begin
 x1:=0.0;
 Repeat
  x2:=x1-F(x1)/G(x1);
  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.8476303072
0.8463956635
0.8463953790
0.8463953790
0.8463953790
0.8463953790
0.8463953790
0.8463953790
0.8463953790
0.8463953790
-------------
0.8463953790 - accurate value

4. Метод простой итерации

Тут хочу отметить сразу: метод, в общем случае, неустойчив, и выдаст ли программа желаемый результат или пойдет "вразнос" – это уж как повезет. Есть некие критерии сходимости, но они сложны и, в основном, избыточны, а потому рассматривать мы их не будем, а остановимся на самом алгоритме.
Первое, что нужно сделать, это преобразовать исходное уравнение F(x)=0 к виду х=Ф(х). Для уравнения, взятого в качестве примера, такое преобразование очевидно (функция Ф(х) там представлена). Дальше, нужно постараться найти число х0, близкое к искомому корню, и подставить его в ПРАВУЮ часть уравнения х=Ф(х). Тогда то, что получится слева, будет следующим приближением х1, каковое мы снова подставляем в правую часть и получаем х2, ну и т.д.

Код:

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);
end;

Procedure EQRoot(e:real; var R:real);
var
 x1,x2,d:real;
begin
 x1:=0.0;
 Repeat
  x2:=F(x1);
  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.8468813641
0.8463839519
0.8463839519
0.8463956466
0.8463953728
0.8463953728
0.8463953792
0.8463953792
0.8463953790
0.8463953790
-------------
0.8463953790 - accurate value

5. Метод "десятичного деления"

Вообще-то официально такого не существует, более того, это метод выдумал я сам для своих целей. Вряд ли, конечно, это "открытие века", скорее всего он имеется в природе под каким-то другим названием, просто я не нашел. Мне он очень нравится, поэтому расскажу о нем здесь: вдруг кого-то заинтересует.
Возьмем в качестве стартовой точки нижнюю границу интервала a (верхняя b нам не понадобится) и еще зададим шаг h. Вот тут важно: этот шаг непременно должен быть степенью десяти (положительной, отрицательной либо нулевой), т.е. 10, 1000, 0.1, 1, 0.000001 и т.п. – в зависимости от задачи. Теперь, стартовав с x=a, начинаем увеличивать x с шагом h до тех пор, пока знак функции F(x) не изменится. Как только это произойдет, "отступим" на один шаг, уменьшим h ровно в 10 раз и повторим процедуру. Потом еще раз. И так будем повторять, пока не добьемся требуемой точности. В чем "изюминка" данного подхода? А в том, что каждое прохождение добавляет значащую цифру следующего разряда. То есть сколько нам надо их иметь, столько раз и следует прокручивать цикл. По-моему, очень даже удобно.
Метод абсолютно устойчив.

Код:

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
 x,a,h:real;
begin
 x:=0.0;
 a:=0.0;
 h:=1.0;
 Repeat
  repeat
   x:=x+h;
  until F(a)*F(x)<0;
  x:=x-h;
  h:=h/10;
 Until h<e;
 R:=x;
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.8000000000
0.8400000000
0.8460000000
0.8463000000
0.8463900000
0.8463950000
0.8463953000
0.8463953700
0.8463953790
0.8463953790
-------------
0.8463953790 - accurate value

09.04.2014, 12:51	#1 (permalink)
Vladimir_S Специалист Регистрация: 27.08.2008 Адрес: Санкт-Петербург Сообщений: 27,807 Сказал(а) спасибо: 340 Поблагодарили 583 раз(а) в 208 сообщениях Репутация: 113184	К вопросу о численном решении алгебраических уравнений Здесь я предлагаю рассмотреть некоторые методы численного решения нелинейных алгебраических уравнений вида 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* (Продолжение следует)

09.04.2014, 14:36	#2 (permalink)
Vladimir_S Специалист Регистрация: 27.08.2008 Адрес: Санкт-Петербург Сообщений: 27,807 Сказал(а) спасибо: 340 Поблагодарили 583 раз(а) в 208 сообщениях Репутация: 113184	(продолжение) 3. Метод касательных (Ньютона) Начинаем с того, что выбираем стартовую точку x0. Этой точкой может быть один из концов интервала (a или b), а может быть просто точка вблизи предполагаемого корня (чем ближе, тем лучше). (Лирическое отступление. Некоторые идиоты-"преподаватели" требуют от студентов, чтобы этой точкой непременно был один из концов интервала, причем не любой, а именно тот, где знаки функции и ее производной совпадают. Заявляю ответственно: подобные "украшательства" есть полная чушь, свидетельствующая лишь о некомпетентности этих самых...). Далее, через точку (x0, F(x0)) проводим касательную к кривой y=F(x) до пересечения этой касательной оси абсцисс в точке х1. Таким образом, имеем прямоугольный треугольник с длинами катетов вертикального – F(x0), и горизонтального – расстояние между х0 и х1. Гипотенузой является отрезок касательной. Если угол между касательной и осью абсцисс, отсчитанный от положительного направления последней, есть α, то из геометрии следует, что (x0-x1)=F(x0)/tg(α). Но тангенс угла наклона касательной есть ни что иное, как производная в точке х0, а потому, окончательно x1 = x0 – F(x0)/F'(x0) Далее берем функцию и производную в точке х1 и по той же формуле получаем х2, потом х3 и т.д. Точки х0, х1, х2, х3, ... и будут последовательными приближениями искомого корня. К сожалению, метод не всегда устойчив. Например, если в стартовой точке производная близка к нулю, то точка х1 рискует оказаться далеко за пределами интервала и попасть в область, где функция либо не определена вовсе, либо успела пройти один или даже несколько экстремумов, и в таком случае процедура будет не приближать, а отдалять получаемые значения х1, х2 и т.д. от корня. Тем важнее становится этап предварительного исследования, о котором писалось выше. Код: 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; Function G(z:real):Real; begin G:=(Sin(z)+Cos(z))2(Cos(z)-Sin(z))/Exp(Ln(33.5)2/3)-1; end; Procedure EQRoot(e:real; var R:real); var x1,x2,d:real; begin x1:=0.0; Repeat x2:=x1-F(x1)/G(x1); 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.8476303072 0.8463956635 0.8463953790 0.8463953790 0.8463953790 0.8463953790 0.8463953790 0.8463953790 0.8463953790 0.8463953790 ------------- 0.8463953790 - accurate value 4. Метод простой итерации Тут хочу отметить сразу: метод, в общем случае, неустойчив, и выдаст ли программа желаемый результат или пойдет "вразнос" – это уж как повезет. Есть некие критерии сходимости, но они сложны и, в основном, избыточны, а потому рассматривать мы их не будем, а остановимся на самом алгоритме. Первое, что нужно сделать, это преобразовать исходное уравнение F(x)=0 к виду х=Ф(х). Для уравнения, взятого в качестве примера, такое преобразование очевидно (функция Ф(х) там представлена). Дальше, нужно постараться найти число х0, близкое к искомому корню, и подставить его в ПРАВУЮ часть уравнения х=Ф(х). Тогда то, что получится слева, будет следующим приближением х1, каковое мы снова подставляем в правую часть и получаем х2, ну и т.д. Код: 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); end; Procedure EQRoot(e:real; var R:real); var x1,x2,d:real; begin x1:=0.0; Repeat x2:=F(x1); 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.8468813641 0.8463839519 0.8463839519 0.8463956466 0.8463953728 0.8463953728 0.8463953792 0.8463953792 0.8463953790 0.8463953790 ------------- 0.8463953790 - accurate value* 5. Метод "десятичного деления" Вообще-то официально такого не существует, более того, это метод выдумал я сам для своих целей. Вряд ли, конечно, это "открытие века", скорее всего он имеется в природе под каким-то другим названием, просто я не нашел. Мне он очень нравится, поэтому расскажу о нем здесь: вдруг кого-то заинтересует. Возьмем в качестве стартовой точки нижнюю границу интервала a (верхняя b нам не понадобится) и еще зададим шаг h. Вот тут важно: этот шаг непременно должен быть степенью десяти (положительной, отрицательной либо нулевой), т.е. 10, 1000, 0.1, 1, 0.000001 и т.п. – в зависимости от задачи. Теперь, стартовав с x=a, начинаем увеличивать x с шагом h до тех пор, пока знак функции F(x) не изменится. Как только это произойдет, "отступим" на один шаг, уменьшим h ровно в 10 раз и повторим процедуру. Потом еще раз. И так будем повторять, пока не добьемся требуемой точности. В чем "изюминка" данного подхода? А в том, что каждое прохождение добавляет значащую цифру следующего разряда. То есть сколько нам надо их иметь, столько раз и следует прокручивать цикл. По-моему, очень даже удобно. Метод абсолютно устойчив. Код: 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 x,a,h:real; begin x:=0.0; a:=0.0; h:=1.0; Repeat repeat x:=x+h; until F(a)F(x)<0; x:=x-h; h:=h/10; Until h<e; R:=x; 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.8000000000 0.8400000000 0.8460000000 0.8463000000 0.8463900000 0.8463950000 0.8463953000 0.8463953700 0.8463953790 0.8463953790 ------------- 0.8463953790 - accurate value

Опции темы
Версия для печати Отправить по электронной почте
Опции просмотра
Линейный вид Комбинированный вид Древовидный вид

Ads
Яндекс Member Регистрация: 31.10.2006 Сообщений: 40200 Записей в дневнике: 0 Сказал(а) спасибо: 0 Поблагодарили 0 раз(а) в 0 сообщениях Репутация: 55070