Гурский Д., Турбина Е. - Вычисления в MathCad 12 (1077322), страница 66
Текст из файла (страница 66)
Чтобы данная ошибка не вызвала сбой в работе программы, нужно сделать так,чтобы функция полинома Р никогда не была равна 0. Для этого необходимо отслеживать значение Р. Если Р окажется равным 0, в качестве значения функции следует возвратить очень малую, однако отличную от 0 величину. Логичным шагом будет увязатьданную величину и TOL (к примеру, как 0.1-T0L). Технически описанный механизмпроще всего реализовать, используя функцию if.
Данная функция по своему назначению аналогична конструкции if-otherwise и, соответственно, она принимает три параметра. Первый параметр — это некоторое условие. Второй параметр — величина, которая должна быть возвращена, если условие окажется истинным. Третий параметр —значение, которое следует возвратить, если условие даст ложь.'last(V) ,last(V)S ЫЬ IP(V,x)J=0s,0.1-TOLj=0,2lp(V,x),,H(V,x)( , ) <- G(V,x)( , )2 P(V,x)'G(V,x)dx2P(V,x)Далее нужно создать цикл, на каждом обороте которого будет вычисляться один корень.
Соответственно, количество итераций, совершаемых данным циклом, должнобыть равно степени полинома.for i e 0.. last(V) - 1Перейдя к нахождению очередного корня, нужно задать начальное приближение а.Чаще всего не имеет значения, что это будет за приближение. Однако в случае комплексных корней оно может влиять на то, сойдется алгоритм или нет. Чтобы исключитьпредвзятость, будем задавать начальное приближение случайным образом в интервале от -100 до 100, используя генератор равномерно распределенных случайных чиселrnd. Задав а, необходимо объявить еще две переменные. Переменная а | ы будет хранитьзначение предыдущего приближения.
Изначально ей можно присвоить любое значение, заметно отличающееся от а. Переменная к — это счетчик итераций, совершенныхпри поиске корня. Первоначально она должна быть равна 0.а<100+ md(200),a l a s t <-2-a,k<-08.1. Решение уравнений• 285Теперь следует создать цикл, который будет реализовывать итеративный механизмпостепенного приближения к корню по формуле Лаггера. Так как количество необходимых итераций заведомо неизвестно, нужно использовать цикл while. В качестве условий его остановки необходимо задать то, что значение полинома в точке приближенияпо модулю должно быть меньше TOL, а также разность между данным и предыдущимприближениями по абсолютной величине должна не превышать TOL.while -,(|P(V,a)|< ТОЬл|а - alastКак уже указывалось выше, в большинстве случаев не имеет значения, какое числовзять в качестве начального приближения.
Но изредка встречаются и исключения (некоторые комплексные корни). Поэтому метод Лаггера может и не сойтись. При этомнужно сменить начальное приближение и осуществить попытку определения корнязаново. Судить о том, что последовательность приближений к корню не сходится, можно по величине счетчика итераций к. Так как скорость сходимости метода Лаггераочень высока, то даже при самом малом значении TOL и самом неудачном выборе начального приближения на то, чтобы найти корень, не должно уйти больше, чем несколько десятков итераций. Если же значение к превысит установленный лимит, тонужно сгенерировать новое приближение, обнулить переменные к и alast, после чегоперейти на новую итерацию посредством оператора continue.
Если же значение к невелико, следует просто увеличить его на 1, зафиксировав очередную итерацию.if k > 100а <- -100+ md(200),alast2-a,k <- 0continuek + 1 otherwiseДалее нужно написать код, который будет, используя формулу Лаггера, на основаниитекущего приближения находить приближение следующее. Величина же старого приближения будет присваиваться переменной a|ast.<- last(V),G a <- G(V,a),H a <- H(V,a)ldSlif(max(m) * O,max(m) ,0.1-TOL)В приведенном фрагменте кода есть два неочевидных момента. Во-первых, нужно ответить на вопрос, зачем создается массив т . Это необходимо, чтобы минимальным количеством кода определить, при каком знаке перед корнем выражение знаменателядроби из формулы Лаггера примет максимальное значение. Занося оба варианта в массив т , мы можем определить, какой из них больше посредством функции max.
Во-вторых, нужно предусмотреть то, что наибольший элемент массива m может быть и нулем. Чтобы избежать ошибки деления на нуль, действуем точно так же, как в случаелокальной функции Р.286 •Глава 8. Решение уравнений и систем уравненийКогда критерии сходимости к корню выполнятся и цикл while прекратит работу, значение последнего приближения должно быть занесено в вектор результатов res. Затемполином должен быть разделен на соответствующий найденному корню линейныймножитель.
Для этого следует задействовать написанную ранее функцию pol_dev. После выполнения этих действий алгоритм будет готов перейти к поиску очередного корня.res. <— aV <- pol_dev(V, a)Когда все корни будут найдены, вектор res должен быть возвращен как результат работы функции. Однако перед этим корни стоит отсортировать в порядке возрастания,воспользовавшись функцией sort.sort (res)Функция LaGuerra готова. Приведем ее полный код:'last(V) ,last(V)Л•TOLLaGuerra(V,TOI):= P(V,x)<-ifj=o2G(V,x)P(V,x),H(V,x) <- G(V,x)2 - dxP(V,x)for i E 0..
last(V) - 1a <- -100+ rnd(200),a l a s t <- 2 a , k < - 0while-,(|P(V,a)| < T O L A | a - a l a s t< TOL,)if k > 100a <- -100 + rnd(200),a l a s t «- 2-a,k <- 0continuek <- k + 1 otherwisen <- last(V),Ga <- G(V,a),H a <- H(V,a)Gaif(max(m) * 0,max(m),0.1TOL)res. <- aV<-pol_dev(V,a)sort (res)8.1. Решение уравнений • 2 8 7Проверим эффективность написанной программы, применив ее по отношению к полиному, корни которого заведомо известны. Сравним выданный ею результат с результатом работы функции polyroots.LaGuerralZ,-144016122Z :=(х+ 1)-(х- 9) -(х- 7)/х + 4)coeffs,x-456186-24.
1 .-0.99999999999998685-1.0000000230252142"-l-9999999999999915i-2-0000000149387862i1.999999999999984912.0000000153338786ipolyroots(Z) =6.99999999999999737.00000003880974398.99999981334185158.9998382448970879V 9.00000018665814499 0001617261872831 .Как видите, наша программа если и уступает функции polyroots, так только в скоростиработы.
Точность же определения ею корней полинома может даже (при малых TOL)превосходить точность, показываемую функцией polyroots. Она способна находить всевиды корней: действительные, комплексные, кратные.Зная суть алгоритма Лаггера, несложно понять, отчего точность определения корнейс повышением степени полинома резко падает. Это связано с тем, что нахождению очередного корня предшествует операция деления полинома на соответствующий определенному ранее корню линейный множитель. Так как приближение к корню имеетограниченную точность, коэффициенты нового полинома будут найдены с ошибкой.Соответственно, ошибка при вычислении нового корня будет складываться из общейпогрешности метода и погрешности, связанной с неточностью коэффициентов.
Поэтому данный корень будет менее точен по сравнению с предыдущим. Таким образом помере вычисления корней точность будет неуклонно падать. Это можно заметить дажев случае полиномов не очень высокой степени. Так, если определить корни полиномастепени 8, то часть из них будет точна до 13-14-го знака мантиссы. Но два-три корнябудут иметь ошибку уже в 4-5-м знаке.Учитывая описанную выше особенность метода Лаггера, его не стоит использовать дляпоиска корней полиномов высокой степени. Для этого лучше применять функцию root,находя каждый корень индивидуально. Или же можно обратиться к методу сопровождающей матрицы, который менее устойчив, однако не склонен к такому выраженномунакоплению ошибки.Метод сопровождающей матрицы мы не будем описывать так подробно, как методЛаггера, однако изложим его основные идеи.Из курса линейной алгебры известно, что любой полином может быть представленв матричном виде как следующий определитель:Р(х)=288 •Глава 8.
Решение уравнений и систем уравненийЗдесь А — так называемая сопровождающая матрица (companion matrix), x — переменная, I — единичная матрица, соразмерная А.Сопровождающая матрица, по сути, хранит коэффициенты полинома. В случае полинома степени m сопровождающая матрица — это матрица размерности mxm, имеющаяследующий вид:lVаmаmаm0000100Оо1оV-iaаА=m1m-2aМожно показать, что корням полинома соответствуют собственные значения сопровождающей матрицы. Следовательно, задача определения корней сводится к задаченахождения собственных значений. Существуют специальные численные методы, которые позволяют решать эту задачу с высокой точностью. Один из этих методов реализован в Mathcad в форме функции eigenvals.
Мы воспользуемся данной функцией,чтобы написать программу, реализующую метод сопровождающей матрицы. Код этойпрограммы будет предельно прост:comp_matr(vector) :=for i e last(vector) - 1.. ОOO,last(vector)-l-i-(vector. + vector,last( vectorfor j e 0.. last(vector) - 2I. . .4-1J+1.Jeigenvals(T)Несложно показать, что точность программы compmatr совпадает с точностью встроенного алгоритма сопровождающей матрицы. Это означает, что они абсолютно идентичны.Z:=(24 -50 35 -10 I ) 1D000338 ^comp_matr(Z) =2.99999999999995652.00000000000001640.99999999999999734/^0.99999999999999734^polyroots(Z) =2.00000000000001642.99999999999995654.0000000000000338,)8.1.4. Графическое решение уравненийГрафическое решение уравнения заключается в определении по графику функции, соответствующей левой части уравнения, при какой величине аргумента данная функция принимает значение, равное правой части уравнения.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
