Гурский Д., Турбина Е. - Вычисления в MathCad 12 (1077322), страница 62
Текст из файла (страница 62)
Остальные2 6 8 • Глава 8. Решение уравнений и систем уравненийметоды группы основаны на той же основной идее, однако они используют разногорода хитрые приемы, увеличивающие скорость сходимости алгоритма. Поэтому естьсмысл рассмотреть сначала метод Больцано в его классическом варианте, а затем обсудить открытые Риддером и Брентом улучшения.Алгоритм метода Больцано следующий.1. Задается функция, нули которой необходимо найти. Определяется интервал от а0до bQ, на котором локализован корень. На то, что корень находится именно на этоминтервале, будет указывать то, что функция на его границах принимает значенияразных знаков. Этот признак будет работать, если функция непрерывна на промежутке [а0 b j , а также если корень не является корнем типа касания.2.
Определяется координата середины интервала [aQ b j с„.3. Проверяется, в какой из границ интервала функция принимает противоположныйзнак тому знаку, который она имеет в средней точке. Тут возможны три варианта.• Если f(a0) и f(c0) имеют разные знаки, то корень находится на промежутке [а0, с0].Значит, можно сузить интервал локализации, взяв в качестве его правой границы с 0 вместо Ьо.• Если f(b0) и f(c0) имеют разные знаки, то корень находится на интервале [с0, b j .Сужаем интервал, взяв в качестве левой границы с 0 вместо а0.• Если со=О, то с0является искомой точкой корня. Возвращаем ее как результат.4. Повторяем описанные выше действия до тех пор, пока длина интервала [au, b j неуменьшится до 2-TOL. При этом, если взять в качестве приближения к корню среднюю точку интервала, погрешность не превысит TOL. Узнать же, сколько итерацийдолжен совершить алгоритм, чтобы интервал локализации корня сузился до 2-TOLочень просто.
На каждой итерации ширина интервала уменьшается в два раза. Исходя из этого можно записать следующее неравенство (здесь N — необходимое количество итераций):< 2-TOLN2Решив данное элементарное неравенство, получим:N^ln(b-a)-In(TOL)|l1п(2)Таким образом, точность поиска корней методом Больцано не зависит от особенностей функции, а определяется только количеством итераций.Крупным недостатком метода Больцано является его медленная сходимость. Так, к примеру, чтобы найти корень с точностью до 0.00001 при ширине интервала 10, необходима21 итерация.
Это довольно много. Существуют методы, в которых проблема медленной сходимости метода Больцано решается за счет того, что учитываются особенностиповедения функции. Это метод ложного положения, а также используемые в Mathcadметоды Риддера и Брента. Увы, но в случае этих методов уже нельзя однозначно оценить, сколько необходимо итераций для достижения требуемой точности приближения корня. Поэтому в качестве критериев сходимости используются величина функции в точке последнего приближения и (или) разность двух последних приближений.Чтобы вам было проще разобраться в сути метода Больцано, приведем его реализациюна языке программирования Mathcad:8.1.
Решение уравненийBolcano(f,a,b,TOL) :=error("Bad interval!") if* 269v ffa) = 0 v f(b) =. ln(b - a) - lnfTOL)N<-floor —*' +1ln(2)for i e 1.. N + 1с <- (a + b) -i- 2return с if i = N + 1return с if f(c) = 0b < - c if«a)•адa <— с otherwiseПроверка показывает, что функция Bolcano работает ничуть не хуже функции root:f(x) := sin(x)Bolcano(f, 1,4,10" 5 )=l .00000010371853л Bolcano(f, 1,4,10~ 10 )=l .00000000001226тсВ приведенной реализации метода Больцано есть несколько технических моментов,на которые нужно обратить внимание.• Ключевым шагом метода Больцано является определение того, принимает функцияв двух точках значения одного или разных знаков. Чтобы это выяснить, разделим значение функции в точке на модуль этого значения.
В результате мы получим 1 или - 1 ,то есть, по сути, выделим знак. Проделав аналогичную операцию для второй точкии сравнив полученные значения, мы узнаем, имеется ли на промежутке корень.Выделить знак числа также можно, используя встроенную функцию Mathcad sign.•Определяя количество необходимых для достижения требуемой точности итераций по выведенной выше формуле, результат необходимо округлить до ближайшего меньшего целого числа (подумайте, почему). Выполняет эту операцию функцияfloor. Соответственно, округление до ближайшего большего целого числа производит функция ceil.• Перед тем как запустить работу алгоритма Больцано, нужно проверить, принимаетли функция на концах интервала значения разных знаков. Если это не так, останавливаем работу кода и возвращаем сообщение об ошибке посредством функции error.Метод Больцано имеет важное историческое и теоретическое значение, однакона практике для решения уравнений он не используется.
Причина — его медленнаясходимость. То есть, чтобы получить посредством него решение нужной точности, следует проделать значительный объем вычислительной работы. Существует несколькометодов, которые, используя ту же основную идею, что и метод Больцано, обладают лучшей сходимостью. Наиболее простой из них — метод ложного положения (Regula Falsi).Отличие метода ложного положения от метода Больцано заключается в том, как вычисляется очередное приближение. В методе Больцано для этого находится среднееарифметическое границ интервала локализации корня: c n =(a n +b n )/2. В методе ложного положения через точки (an, f(an)) и (bn, f(bn)) проводится секущая (рис.
8.5). Точкой2 7 0 •:• Глава 8. Решение уравнений и систем уравненийсп, в которой секущая пересекает ось X, заменяют одну из границ интервала локализации корня, сужая его. Приближение считается достаточно точным, если I f(c n )| <TOLи (или) I cn-cn.,| <TOL (простой формулы, позволяющей определить, сколько необходимо итераций для достижения требуемой точности, в случае метода ложного положения, в отличие от метода Больцано, нет).Рис. 8.5. Иллюстрация метода ложного положения. Неплохая точность достигается уже навторой итерацииМетод ложного положения сходится быстрее метода Больцано, но есть и более эффективные методы.
По умолчанию функция root применяет метод Риддера. Основноеулучшение в нем, по сравнению с методом Regula Falsi, довольно простое, и основанооно на распространенном в теории численных методов приеме — использовании интерполирующего функцию полинома (точно так же осуществляется, например, переход отметода секущих к методу Мюллера, от метода трапеций — к методу Симпсона, от метода Эйлера — к методу Рунге-Кутта).
Если излагать суть метода Риддера в деталях, тона каждой итерации проделываются следующие шаги. Сначала вычисляется значение функции в средней точке интервала c n =(a n +b n )/2. Затем через имеющиеся три точки (граничные точки интервала и средняя точка) проводится интерполирующая функцию парабола. Так как через К точек проходит только один полином степени К-1(попробуйте это доказать), то интерполирующая парабола будет уникальной. Несложно вывести, что данная парабола будет определяться следующей функцией (подобный вывод приводится в гл. 10 при рассмотрении метода Симпсона):Pn(x)=f(bn)x2-2f(Cn>x+f(aЕсли интервал (an, b n ) достаточно узок, то поведение интерполирующей параболы будет весьма схоже с поведением функции.
Поэтому логично в качестве приближенияк корню использовать точку, в которой парабола пересекает ось X. Найти эту точкунесложно, положив Р п (х)=0 и решив полученное квадратное уравнение. В результате будет получено два корня, однако только один из них будет принадлежать интервалу (an, b n ). Риддер вывел формулу, которая дает возможность сразу вычислить нужныйкорень (при ее выводе считалось, что уравнение задано не относительно переменной х,са относительно экспоненциальной величины е ).
Данная формула имеет вид:8.1. Решение уравнений* 271f(c n ).sign(f(c n )-f(, n ))rn= спПосле нахождения приближения, проверяется, соответствует ли оно используемымкритериям сходимости к корню (в Mathcad это I f(rn)| <TOL). Если нет, то осуществляются стандартные для всех методов интервалов локализации корня шаги. А именно:приближением заменяется та граница интервала (an, b n ), которая имеет с ним одини тот же знак. После этого алгоритм переходит на новую итерацию.Выведенное Риддером уравнение имеет несколько замечательных свойств.
Во-первых,приближение всегда оказывается внутри интервала (an, b n ). Во-вторых, вычисленнаяна основании него последовательность приближений сходится довольно быстро. Каждое последующее приближение будет в среднем на две значимые цифры более точнымпо сравнению с предыдущим. В-третьих, метод Риддера очень устойчив.Реализовать метод Риддера на языке программирования Mathcad несложно. Напишемсоответствующий алгоритм так, чтобы в качестве результата он возвращал не толькоприближение к корню, но и число, показывающее, за сколько итераций было найденорешение. Это необходимо для объективной оценки эффективности программы.Ridder(f,a,b,TOL):= error("Bad Interval!") if sign(a)= sign(b) v f(a) = Ov f(b) = 0(r p r e v <-an«-o)while 1с <- (a + b) т 2rbcitc-a).^*-™return (r n) if (|f(r)| < T O L A | r - r p r e v< TOL)v f(r) = 0b <- г if sign(a) * sign(r)a <— r otherwise0prevРешим с помощью написанной программы уравнение с очевидным корнем при разномзначении TOL:f(x):=x3-l3Ridder(f,0,4,10~ )=(1.000004682614623) Ridder(f,O,4,lCT5)=(1.0OO0OO09556532 4)Ridder(f, 0,4,10~8)=(1-ОО00ОО00ОО398 6)Ridder(f, 0,4,1 СГ10)=(1.000000000000817)Проанализируем полученные результаты.