Методы решения нелинейных уравнений в MathCAD
Реализуем для некоторого уравнения 4 наиболее популярных численных метода для решения нелинейных уравнений. При этом мы стремимся именно запрограммировать методы, а не воспользоваться встроенным инструментом Given...Find
или функциями root
, polyroot
. Об этих способах решения почитайте, например, здесь.
Определим функцию уравнения f(x)=0
как функцию пользователя, интервал поиска решения зададим переменными a
и b
. Найти этот интервал можно, например, табличным или графическим методом:

Уравнение и интервал поиска решения
Начальный интервал[a,b]
должен быть таким, чтобы значенияf(a)
иf(b)
имели противоположные знаки. Если искомый корень уравнения окажется единственным на интервале, то совсем хорошо :)
Логика метода дихотомии (возможно, более правильные названия - метод бисекции, метод половинного деления) довольно проста: если на концах выбранного интервала [a,b]
знаки функции совпадают (произведение f(a)*f(b)>0
), то вернуть результат "недопустимый интервал" (вернём в этом случае ответ "бесконечность"), в противном случае до тех пор, пока длина интервала не станет меньше заданной погрешности ε
, будем находить середину текущего интервала c=(a+b)/2
, считать в ней значение функции и проверять, какую из половин отрезка [a,c]
или [c,b]
нужно отбросить для выполнения следующего шага - а именно, ту, в которой знак f(c)
совпадает со знаком функции на левой или правой границе интервала (в листинге - проверка f(a)*f(c)>0
). Для большей точности вернём середину "последнего" интервала [a,b]
, меньшего ε
:

Метод дихотомии в MathCAD
В методе простой итерации исходное уравнение f(x)=0
представляется в эквивалентном виде φ(x)=x
(что, вообще говоря, можно сделать бесконечным числом способов), а затем шаг метода выполняется по формуле xk+1 = φ(xk)
, пока не будет достигнута заданная точность |xk+1-xk|<ε
. Если выбрать φ(x)=x-c*f(x)
, то константу c
целесообразнее всего искать методом релаксации, для которого c=2/(M+m)
, где M
- максимальное из значений первой производной на концах отрезка или в находящихся на нём точках перегиба функции (точках, где f''(x)=0
), а m
- минимальное из таких значений. Вот соответствующий расчёт в MathCAD:

Метод простой итерации (релаксации) в MathCAD
Если заданной сходимости нет в течение 10000 шагов, в подпрограмме предусмотрен аварийный выход.
Численный метод Ньютона решения нелинейного уравнения основан на формуле вида xk+1 = xk-f(xk)/f'(xk)
, обеспечивающей наилучшую сходимость, но требующей дополнительного вычисления производной на каждом шаге. Так как производные для MathCAD - не проблема, можно всё сделать "в лоб":

Метод Ньютона для решения нелинейного уравнения в MathCAD
Видно, что сходимость метода - на 2 порядка выше (погрешность найденного решения ~ 10-8).
Наконец, существует метод хорд, в котором кривая f(x)
заменяется прямой линией (хордой), стягивающей точки (a, f(a))
и (b, f(b))
. Формула этого метода зависит от знака выражения f(a)*f''(a)
, то есть, имеет два варианта:
Если f(a)*f''(a)>0
, то x0=b
, xk+1=a-(f(a)(xk-a))/(f(xk)-f(a))
Если f(a)*f''(a)<0
, то x0=a
, xk+1=xk-(f(xk)*(b-xk))/(f(b)-f(xk))
Вот примерная реализация на MathCAD, как и в предыдущих двух случаях, контролируется максимальное число итераций, равное 10000:

Метод хорд решения нелинейного уравнения в MathCAD
Видно, что сходимость метода оказалась в нашем случае не столь высока.
Подсчитать, сколько шагов какому методу потребовалось, можете сами, немного поменяв выдачу подпрограмм.
Скачать расчёты из статьи в архиве .zip с файлом .xmcd Mathcad 14/15 (34 Кб)
05.09.2013, 15:07 [116892 просмотра]