Гурский Д., Турбина Е. - Вычисления в MathCad 12 (1077322), страница 74
Текст из файла (страница 74)
В программе Newton отслеживается количество проделанных итераций. Если оно превышает 100, то возвращается сообщение об ошибке.Проверим, насколько эффективна программа Newton, решив уравнение, корень которого можно точно найти с использованием оператора solve:320 •Глава 8. Решение уравнений и систем уравненийf(x) :=sin(x)хsin(x) - - solve ,x -> 2.772604708265991234x/_3ч (2.7726047083006953^ х тЛ . 1Л-15Л (2.77260470826599АNewton^, 4,10 J =Newton(f,4,10)=\Проверка показывает, что метод Ньютона сходится очень быстро. Каждая итерацияувеличивает точность в среднем на два знака мантиссы. Это больше, чем в случае таких методов, как метод секущих или Больцано. Но этот вывод вовсе не означает, чтодля решения всех уравнений лучше использовать вычислительный блок Given-Find (реализующий для нахождения корней одного уравнения с одним неизвестным схожийс описанным алгоритм). Все дело в том, что для метода Ньютона характерны те же недостатки, что и для метода секущих (трудности в поиске решений при наличии экстремумов, перегибов, точек разрыва и особенностей поведения функции на бесконечности).
Однако, помимо этого, к описанным недостаткам добавляются еще и многочисленныепроблемы, возникающие в связи с использованием производной. Поэтому в случае сложных уравнений лучше не рисковать и применять для решения более надежный алгоритм секущих (а еще лучше глобально сходящиеся методы интервалов локализациикорня). Вообще, на практике метод Ньютона редко используется (ввиду описанныхнедостатков) для решения одиночного уравнения. Однако очень важно его расширение для систем уравнений. Конечно, в нем не исчезают присущие этому алгоритму слабости и недостатки, но для поиска решений систем уравнений просто пока не созданоничего лучшего.Попробуем самостоятельно расширить метод Ньютона на случай систем уравнений.Для начала разберемся с наиболее простым случаем и придумаем идею алгоритма решения системы из двух нелинейных уравнений. Для этого вспомним, что любое такоеуравнение — это частный случай (точка нуля) некоторой функции z=f(x, у), которуюможно представить в пространстве как поверхность.
Очевидно, что решением системыуравнений будет являться точка па линии пересечения (или просто точка касания)поверхностей, соответствующих каждому из уравнений. С другой стороны, значенияобеих функций в такой точке должны быть равны нулю.Итак, с тем, какой геометрический смысл имеет точка корня системы из двух уравнений,мы разобрались. Но как можно найти такую точку? Чтобы ответить на этот вопрос, вспомним, как мы справились с этой проблемой в двумерном случае. Тогда удалось найтикорни благодаря тому, что сложные кривые были заменены на простую и техничнуюпрямую.
Первая мысль, которая возникает, — это просто повторить этот проверенныйи такой эффективный шаг. Однако проблема заключается в том, что поверхность нельзяаппроксимировать прямой! Это то же самое, что заменить кривую точкой: глупо и бессмысленно. Поверхность должна быть заменена каким-то, по возможности более простым и легко исследуемым, но объемным объектом. Очевидно, что таким объектом является плоскость. Не менее очевидно, что две плоскости, приближающие в даннойокрестности некоторые поверхности, почти наверняка пересекутся.
Результатом ихпересечения будет прямая. А, точка, в которой данная прямая пересечет плоскость X0Y,и будет либо корнем (в лучшем случае), либо приближением.Нарисовать или тем более представить работу алгоритма, который бы реализовывал описанную выше идею, очень сложно. Поэтому просто поверим собственной интуиции.8.2.
Решение систем уравнений•> 3 2 1Первый вопрос, на который мы должны ответить, это как можно провести касательнуюплоскость к поверхности в некоторой точке. Сделать это очень просто: для этого записываем разложение соответствующей функции в ряд Тейлора по двум первым членам.Математической тонкостью этого действия в случае функции двух переменных является то, что разложение должно быть проведено по обеим переменным.f(x,y) = f ( x n , y n ) + ^ ( х п , У п ) ( х - х д ) + д-^п,Уп)(у- Уп)Аналогичное выражение можно записать и для второго уравнения системы:Р(х,у) = р(х„,У п )+ | Р ( Х П . У П ) ( Х - х п ) + ^ уТак как значение обеих функций в точке корня равно 0, то имеем следующую системулинейных уравнений:- Х п ) + | Р ( Х п ' У п ) ( У - Уп-f(x ,y )J-f(x ,х-хU8-{(х ,- 3уИз этой системы нам нужно найти значения х и у, которые являются координатамиточки следующего приближения (или корня).
Однако, раскрывая скобки, мы лишьусложним вид системы. Поэтому на данном этапе сделаем замены:X = х - хпПроизведя замену, перепишем полученную систему из двух линейных уравнений в матричном виде:Левая матрица в произведении — это так называемый якобиан. В данном случае ониграет роль матрицы коэффициентов.Решить систему можно или через обратную матрицу, или (техничнее) используя функцию LsoLve. Когда корни системы будут найдены, выражаем из них координаты точкинового приближения.
Если это приближение удовлетворит критериям нуля функцийf(x,y) и р(х,у), то возвращаем его как решение. В противном случае переходим к следующей итерации.Важно учесть то, что якобиан может быть вырожден. Это означает, что в выбраннойточке частные производные двух функций равны, а следовательно, аппроксимирующие плоскости будут параллельны и никакой точки приближения получено не будет.Также возможна ситуация, что одна из плоскостей будет параллельна плоскости X0Y.В этом случае одна из строк якобиана будет нулевой, и он не будет нести необходимойинформации. Как же таких проблем избежать? Наиболее простой вариант — не попадать3 2 2 •:• Глава 8.
Решение уравнений и систем уравненийв точки, в которых якобиан вырожден. А сделать это можно очень просто: в том случае,если такое попадание все-таки происходит, следует просто изменить рабочие значенияприближений, прибавляя к ним некоторую (небольшую) величину, например TOLРассмотрим теперь, как можно решить систему из N уравнений. В общем-то, это можно сделать точно так же, как в случае системы из двух уравнений: решая систему из Nлинейных уравнений, в качестве матрицы коэффициентов которой выступает якобиан, вектора неизвестных — вектор разностей новых и старых приближений, вектораправых частей — вектор значений соответствующих уравнениям функций в точке текущего приближения, взятых с противоположным знаком.Язык программирования Mathcad слишком примитивен, чтобы на нем можно былореализовать функцию, способную решать системы из произвольного количества уравнений.
Поэтому мы ограничимся созданием функции, позволяющей решать системыиз двух нелинейных уравнений.sys(f,p,xO,yO,TOL):= (п <- 0 з^ <- хО у 0 <- уО-z(x,y) dF v (z,x,y) <- - z ( x , y )УЭхдуwhile>TOLv>TOLreturn error("Dont converge to a solutions" ) if n > 100dFx(P'Xn'yn)d(x n + 1 <-x n +TOL y n + ] <-y n +TOL n < - n + l continue) if | j | = 0Z <— 1 solveLH-OJJ<- Z_ - x Jу , < - Z , - yJ0n n+11 n(Xnynn < - n + A)nПопробуем решить с помощью созданного алгоритма следующую систему уравнений:22х + у =16у = sin(x)Для данной системы довольно легко подобрать начальные приближения. Это связанос тем, что в обоих уравнениях можно выразить одну переменную через другую, по причине чего система может быть визуализирована двумя кривыми. Корням же будут соответствовать пересечения этих кривых (рис.
8.14).8.2. Решение систем уравнений* 323-5- 1Рис. 8.14. У системы будет два решения. Первое будет расположено вблизи точки (4, -1),второе — (-4, 1)С определением корней рассматриваемой системы, ввиду того, что мы знаем приблизительное расположение решений, программа sys легко справится:f(x,y) : = 1 6 - x2-у2p(x,y):=y-sin(x)f 3.9358754812sys(f,pA-l,1013)=-0.7133611965 sys (f,p, -4,l,104-3.93587548120.71336119654JВ основе функции find лежат куда более сложные алгоритмы, чем созданный нами.Однако базовая концепция у них та же — метод Ньютона. Поэтому мы можем легкопредсказать, какие у них будут слабые места. Чтобы решение было найдено, начальныеприближения должны быть близки к корням. Если между приближением и корнемокажется экстремум, разрыв или область медленного изменения функции, то численный метод, скорее всего, не сойдется.
Главная сложность заключается в том, чтобы такие приближения найти. Задача эта трудная и, если в системе больше двух уравнений,практически нерешаемая. По этой причине использовать численные методы решениясистем уравнений крайне сложно. Увы, но универсальных и надежных алгоритмовпока не существует.8.2.4. Приближенное решение систем уравненийОчень часто приходится сталкиваться с системами, не имеющими решения. Однакоиногда необходимо выяснить, при каких значениях переменных система не согласована наименее, иначе говоря, имеет минимальную невязку.Для решения задач такого рода в Mathcad имеется специальная функция minerr(xl,...).По особенностям своего применения она абсолютно идентична функции find, то естьтакже требует ввода ключевого слова Given и определения начальных приближений.При своей работе функция minerr использует те же численные алгоритмы, что и функция find.