Методы решения нелинейных уравнений в 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 [113605 просмотров]