Технический форум
Вернуться   Технический форум > Программирование > Форум программистов > Помощь студентам


Ответ
 
Опции темы Опции просмотра
Старый 09.04.2014, 12:51   #1 (permalink)
Vladimir_S
Специалист
 
Аватар для Vladimir_S
 
Регистрация: 27.08.2008
Адрес: Санкт-Петербург
Сообщений: 23,721
Сказал(а) спасибо: 122
Поблагодарили 268 раз(а) в 90 сообщениях
Репутация: 60099
По умолчанию К вопросу о численном решении алгебраических уравнений

Здесь я предлагаю рассмотреть некоторые методы численного решения нелинейных алгебраических уравнений вида F(x)=0, т.е. задачу, с которой часто сталкиваются наши будущие светила науки и техники, а ныне студенты. Будут рассмотрены:

1. Метод половинного деления.
2. Метод хорд (секущих).
3. Метод касательных (Ньютона).
4. Метод простой итерации.
5. Метод "десятичного деления" (условное название).

Программы написаны на языке Паскаль. Для всех методов структура программ идентична: функция F(х) (а для метода Ньютона также и функция G(x)=F'(x)) задается в виде отдельной подпрограммы, что позволяет легко заменить ее требуемой, сам процесс нахождения корня оформлен в виде процедуры, а тело программы сводится к циклу по разным точностям определения корня в направлении повышения итоговой точности.

iter.jpg

ВНИМАНИЕ – ВАЖНО!!!
Прежде, чем приступить к написанию программы вычисления корня ЛЮБЫМ методом, надлежит провести небольшое исследование, сводящееся к тому, чтобы:
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
(Продолжение следует)
__________________
With Mozilla Firefox - straight to communism!
Vladimir_S на форуме   Ответить с цитированием

Старый 09.04.2014, 12:51
Helpmaster
Member
 
Аватар для Helpmaster
 
Регистрация: 08.03.2016
Сообщений: 0

Пожалуйста, посмотрите еще несколько тем по вашей проблеме

Решить систему дифференциальных уравнений
Решить систему линейных уравнений методом Зейделя
Решение уравнений методом Ньютона
Система для решения уравнений методом хорд на C/C++
Решение системы уравнений в Экселе
Подскажите, по вопросу определения расстояния

Старый 09.04.2014, 14:36   #2 (permalink)
Vladimir_S
Специалист
 
Аватар для Vladimir_S
 
Регистрация: 27.08.2008
Адрес: Санкт-Петербург
Сообщений: 23,721
Сказал(а) спасибо: 122
Поблагодарили 268 раз(а) в 90 сообщениях
Репутация: 60099
По умолчанию

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

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
__________________
With Mozilla Firefox - straight to communism!
Vladimir_S на форуме   Ответить с цитированием
Ads

Яндекс

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

Опции темы
Опции просмотра

Ваши права в разделе
Вы не можете создавать новые темы
Вы не можете отвечать в темах
Вы не можете прикреплять вложения
Вы не можете редактировать свои сообщения

BB коды Вкл.
Смайлы Вкл.
[IMG] код Выкл.
HTML код Выкл.
Trackbacks are Вкл.
Pingbacks are Вкл.
Refbacks are Выкл.




Часовой пояс GMT +4, время: 12:56.


Powered by vBulletin® Version 6.2.5.
Copyright ©2000 - 2014, Jelsoft Enterprises Ltd.