Гурский Д., Турбина Е. - Вычисления в MathCad 12 (1077322), страница 73
Текст из файла (страница 73)
Аналогично тому, как функция считается равной нулю, если ее величина по модулю меньше TOL,так и ограничивающее неравенство будет считаться соблюденным, если отклонение отнего не превышает определенного порога. Порог этот задается специальной системнойконстантой CTOL (Constraint Tolerance — Точность ограничения). К примеру, приСТ01.= 10~3неравенство х>5 окажется невыполненным, если х=5,099. По умолчанию CTOLравняется 10~3, ее предельная величина — 10 17. Задать ее можно точно так же, как TOL:переопределением на рабочем поле или открыв вкладку Build-in Variables (Встроенныепеременные) окна Worksheet Options (Параметры документа).
Менять CTOL в сторонууменьшения нужно тогда, когда результат должен быть получен с максимальной точностью. Однако нужно учитывать, что при этом может возрасти время расчета вплотьдо того, что система посчитает, что корня не существует. Поэтому стоит воздержатьсяот «профилактического» присваивания CTOL предельного значения.Довольно очевиден следующий вывод: простые уравнения и системы лучше решатьсимвольно, более сложные — численно. Вообще, численные методы куда более универсальны и надежны.
Однако у них есть один весьма существенный недостаток: корниопределяются не точно, а с заданной погрешностью вычислений. Предельно высоким|516уровнем точности работы функции find является 10 -10" . Для очень многих практически важных задач такой точности явно недостаточно.
Кроме того, для всех итерационных численных методов характерно накопление ошибки, которое может довольно316 •Глава 8. Решение уравнений и систем уравненийсерьезно исказить результат. Так, например, при решении систем уравнений аналитической химии для определения рН какого-то вещества или смеси веществ, искомоезначение концентрации ионов водорода обычно весьма мало.
Поэтому, при использовании численных методов, оно зачастую просто округляется до нуля. А так как рН —+это отрицательный десятичный логарифм от концентрации Н , то вполне понятно, чтопри таком округлении его значение правильно вычислено быть не может. Если же длярешения этой задачи использовать символьный процессор, то ответ будет найден безкаких-либо трудностей. Это связано с тем, что символьный процессор при приблизительном решении систем алгебраических уравнений использует арифметику длинныхмантисс и предельный уровень точности увеличивается до 10~20. Однако, увы, такойподход применим только к простейшим, алгебраическим системам.В следующем примере предложено решение подобной задачи.
Изучив его, вы поймете преимущество используемой символьным процессором арифметики длинныхмантисс, а также важность верного выбора численного метода. Если вы не владеете элементами химии, то просто рассмотрите решение данной системы с точки зрения трудностей, которые могут возникнуть при малой величине параметров или корней уравнений.Пример 8.33. Определение рН миллимолярного раствора гидрофосфатанатрияНачальные условия (к, — константа диссоциации гидрофосфата, к2 — дигидрофосфата, ки.
— константа диссоциации воды, С — концентрация гидрофосфата):K 3 :=4.8-1013.-14к 2 :=6.2-10С:=10-3Приближения для численного метода (z — концентрация фосфата, у — гидрофосфата, х — дигидрофосфата, о — ОН', h — Н + ):z:= 0о := 0h := 0х := 0у := ОРешение системы уравнений, описывающей происходящие процессы диссоциации и гидролиза:Giveno-h = k-=wx+y + z = Cx+2y+3z-o-2C=0Ответ численного метода, использующегося по умолчанию (неверный в рамках рассматриваемойзадачи):,-4З.ЗЗЗхЮ.-43.519x10find(x,y,z,o,h) =.-43.148x10-1.852x1v0-5\8.2.
Решение систем уравнений* 317Ответ метода Левенберга (точность в рамках данной задачи приемлемая):.-44.842x10.-61.389x10find(x,y,z,o,h) =—45.144x10г-53.013x10-101.778x10Решение символьного процессора (наиболее точное):5.3152187228378056243 -10-1.3905082583061470611fnd(x,y,z,o,h)-410~ 64.6986863597452558463 • 10~ 4-6.165323630925497779110~5-1.6219748708469446330 10"'°9.967676533828310467610-44.718737001435861О878 10-1.6184142415843580576 • 10~ 71.3871883846497920911 -10~ 63.394188O41327389O456 \0~65.2673911147176409913 - Ю " 4-9.9337346534150365771 • 10~ 45.4865411328177990344 -10~5-1.00667О7385386202768 - 1 0 " "1.8226419447008066607 - 1 0 " 1 0Определение рН:рН := log (l .8226419447008066607-10~10)pH = 9.739Правильный выбор используемого метода играет первостепенную роль для верногои легкого решения уравнения или системы уравнений.
Однако решаете ли вы задачучисленно или используете возможности символьного процессора, очень внимательноотноситесь к результату вычислений и при малейших сомнениях обязательно делайтепроверку. Этот принцип проверки очень важен, поэтому мы еще не раз его повторимв этой книге. Постарайтесь ему следовать, и вы не допустите ошибку или оплошностьпочти наверняка.Все употребляемые в Mathcad алгоритмы решения систем нелинейных уравнений относятся к так называемым градиентным, то есть связаны с приближением функцийс помощью производных. Однако они могут отличаться целым рядом технических тонкостей (влияющих на эффективность алгоритмов): способом вычисления производнойили решения возникающей на определенном этапе системы линейных уравнений, условием принятия точки в качестве корня или приемом разрешения проблемы сингулярности якобиана.
Но все алгоритмы решения систем уравнений в Mathcad можноотнести к модификациям метода Ньютона.Чтобы разобраться, в чем заключается данный метод и с чем могут быть связаны проблемы при его использовании, сначала рассмотрим простейший случай поиска корнейодного уравнения с одним неизвестным. В этом случае принцип работы метода Ньютона довольно схож с уже рассмотренным ранее методом секущих и заключается в следующем: через определенную пользователем начальную точку проводится касательнаяк кривой функции, нуль которой нужно найти.
В точке, в которой касательная пересекает ось X, вычисляется значение функции. Если оно меньше (по модулю) установленного минимума точности TOL, то данная точка определяется как корень. В противном случае через точку функции, которая соответствует полученному приближению,3 1 8 •:• Глава 8. Решение уравнений и систем уравненийпроводится вторая касательная. Точка ее пересечения с осью X определяется как второе приближение. Если значение f(x) в ней больше по модулю, чем TOL, проводится третьякасательная.
И так далее до тех пор, пока условие корня не будет выполнено (рис. 8.13).Рис. 8.13. Иллюстрация метода Ньютона. Достаточное приближение к корню достигаетсяна пятой итерацииТеперь, зная идею метода Ньютона, попробуем самостоятельно написать программупоиска корней нелинейных уравнений. Для этого сначала следует решить проблемычисто технического плана: во-первых, нужно придумать, каким образом мы будем проводить касательные через точки приближений, и, во-вторых, как затем находить точкиих пересечения с осью X.В общем случае условие пересечения прямой линии (каковой и является касательная)с осью X можно определить как:k-x+b = 0Из этого уравнения требуется найти х, однако чтобы это сделать, нужно знать величины констант к и Ь.
Определить одну из них предельно просто: ведь, как известно, к численно равна тангенсу угла наклона прямой к оси X. В случае же касательной тангенсэтот равен значению производной функции в точке касания. То есть:dxЗная к, можно очень легко найти и Ь, исходя из того факта, что в точке касания значения функции касательной и исследуемой функции совпадают. То есть:^ f (vx n /) . x n + b = f ( x n )'dxВыражаем b:Подставляем полученные выражения для констант в первую формулу:if(xn).(x-xn)+f(xn)=08.2.
Решение систем уравнений • 3 1 9Из этого простого уравнения с одним неизвестным находим требуемую точку пересечения касательной с осью X:х=хп-<»•)dxДанная формула носит название формулы Ньютона. Зная ее, реализовать соответствующий численный метод предельно просто:Newton(f,a,TOL) :=(п <- 0х <- а |df(x)<--fl:x)dxwhile 1х <— х + TOL on error 1 + dfnc |nn\njx - ffx V df(xn+1return (xn n j if f(xn)|<TOLA \ffo-f ^< TOLreturn error("Dont con\ergeto a solution!" ) if n > 100В коде программы Newton есть несколько тонкостей, которые стоит прокомментировать.О Производную от f(x) необходимо задать как локальную функцию df(x).
Это связано с тем, что если использовать оператор дифференцирования в выражении, соответствующем формуле Ньютона, напрямую, то возникнет ошибка (так как системавоспримет f(xn) как постоянную величину, а не функцию от х).•Если очередное приближение совпадет с экстремумом, то возникнет ошибка деления на нуль.
Обойти эту проблему просто. Для этого нужно, при возникновенииошибки, прибавить к текущему приближению небольшую величину (например,TOL). Точка нового приближения сместится в сторону от точки экстремума, и поэтому значение производной в ней уже не будет равно нулю.• Метод Ньютона является локально сходящимся. Это означает, что то, сойдется онк корню или нет, зависит от правильности выбора начального приближения. Важнопредусмотреть ситуацию, в которой неудачное приближение породит расходящуюся последовательность приближений.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
