Гурский Д., Турбина Е. - Вычисления в MathCad 12 (1077322), страница 84
Текст из файла (страница 84)
В первом случае нужно проводитьинтегрирование, во втором — применять методы численной оптимизации. И то, и другое само по себе является серьезной задачей. Но можно ли как-то вычислить интегралс заданной точностью, не применяя формул для ошибки?Решение поставленной проблемы можно легко найти, если вспомнить, что с ростомколичества разбиений результаты вычислений любым численным методом образуютравномерно сходящуюся (в случае отсутствия погрешности округления) последовательность.
Данная последовательность более или менее быстро, в зависимости отиспользуемого алгоритма, сходится к истинному значению интеграла. В принципе,абсолютная точность приближения наступает лишь при достижении ширины интервала разбиения бесконечно малого значения. Однако на практике нас всегда интересует определенная, обычно не очень высокая точность результата. Очевидно, что чем ближе находятся два соседних значения интеграла (вычисленные соответственно при пи 2-п разбиениях) к истинному значению, тем на меньшую величину они отличаются.А что если использовать их разность в качестве критерия достижения требуемой точности?Описанная идея реализуется в программе Autolnt, которая вычисляет интеграл по методу трапеций. Данная программа генерирует последовательность приближений к интегралу, удваивая на каждом обороте количество разбиений промежутка.
Техническиэто реализуется с помощью бесконечного цикла while, который прекращает работув том случае, если два последних приближения к интегралу отличаются на величину,меньшую TOL. Сделать же цикл бесконечным очень просто: для этого выражение в егоправом маркере должно всегда давать истину (то есть 1). В принципе, можно простоввести в правый маркер 1.3 6 8 • Глава 10. Вычисление интеграловAutoInt(f,a,b,TOL):=INT_PREV<- (a - b)while In <- n-2.h<b-anINT<-0for i e l..nf[aINT <r- INT + h-(i - 1)] + f(a + h-i)break if (| INT - INT_PREV | < TOL)INTPREV*- INT otherwiseINTПроверим, насколько эффективно работает написанная программа:f(x) := sin(x)Autolnt(f,0,7r, 10" V - 2 = -1.5687316193x10"6AutoIntV.O.n, 10" 8 ) - 2 = -1.5319761015x10"9Autobit(f,0,7i, 10" 1 2 ) - 2 щ -1.2079226508x10" 1 3Программа работает неплохо.
Конечно, она не так эффективна, как встроенные в Mathcadалгоритмы, однако ее можно улучшить. Главный ее недостаток заключается в том, чтопри каждом удвоении количества разбиений значения функции в узловых точках пересчитываются «с нуля», хотя из 2-п значений половина уже была вычислена на предыдущей итерации цикла while. Также, если п>1, то лучше использовать не формулутрапеций, а формулу Симпсона или Буля. Исправить описанные недостатки, сохранивглавную идею программы Autolnt можно, применяя алгоритм Ромберга.Метод Ромберга довольно сложен для понимания, однако мы попробуем максимальнопросто описать его идею.
Для этого лучше всего будет изложить последовательностьдействий при решении задач этим методом.1. В качестве самого первого приближения вычисляется значение площади трапеции,основания которой проведены через границы промежутка интегрирования. Полученное значение заносится в матрицу приближений как:Ro о *- ^ Т ^ № ) + ^ Ь ) )2. Затем запускается бесконечный цикл. На каждом его круге шаг уменьшается вдвое.На первом обороте цикла вычисляется приближение по формуле трапеций из тогоусловия, что интервалов интегрирования уже 2.
Полученное, более точное, приближение интеграла заносится в первую строку нулевого столбца. На втором круге шаг10.6. Численные методы интегрирования• 369будет уменьшен в четыре раза относительно первоначального, и соответственнотрапеций также будет просуммировано 4. Новое приближение будет занесено вовторую строку. И так далее до тех пор, пока не будет выполнено условие остановкицикла.
На языке программирования описанную последовательность действий можно задать так:while 1J+ 1h<-(Ь - а)2R>-J7-UJ,O *- ~2 г ^ +h- ZЧа+H2p"1)]P =1Изучая данный фрагмент программы, вы не найдете уже хорошо нам известнуюформулу метода трапеций. Здесь она заменена на более простую с точки зрения вычислительной работы рекуррентную формулу. Формула позволяет, использовавзначение интеграла, полученное на шаг раньше при вдвое большей ширине интервала разбиения, найти новое приближение.
Достоинство этой формулы в том, чтоона приводит к результату при вдвое меньшем количестве вычислительных операций, чем простая формула трапеций.3. Как было показано выше, точность формул Симпсоиа и Буля значительно выше, чемточность метода трапеций. Поэтому вполне логичным было бы использовать их, а неждать, пока в ходе работы цикла количество разбиений промежутка возрастет до техдесятков тысяч, при которых алгоритм трапеций может дать достаточную точность.В алгоритме Ромберга действительно производится переход к более точным интерполирующим полиномам высокой степени. Причем происходит он таким образом,что никаких приближений интеграла, кроме приближения по методу трапеций,вычислять не приходится.Все дело в том, что в математике доказана связь между приближениями интеграла,рассчитанными различными численными методами.
Так, зная два приближения,полученные по методу трапеций, можно вычислением совершенно простого арифметического выражения получить приближение, которое по точности соответствует формуле Симпсона. Обязательным условием (на основе этого строился вывод)для используемых двух приближений более низкой степени точности является то,что шаг, при котором они были получены, должен отличаться в два раза:Не стоит пугаться этой общей формулы, на самом деле она предельно проста.
Оналишь показывает, что для того чтобы получить приближение но Симпсону с шагомh, следует умножить приближение, полученное гири этом шаге по методу трапеций,на 4, отнять от него элемент того же столбца, соответствующий вычислению интеграла методом трапеций при шаге 2h, и поделить затем результат на 3.
То есть вычисление приближения более высокой степени в том случае, если есть два значения интеграла, рассчитанные по более простому методу, сводится к элементарнойарифметической операции из трех действий.3 7 0 • Глава 10. Вычисление интеграловВычислив приближения по Симпсону, можно получить приближения еще на двапорядка более точные, если применить рекуррентную формулу Буля:B(f,h) =16S(f,h)-S(f,2h)15В общем виде формула уточнения имеет вид:к4-1Проще всего понять ее смысл можно, представив себе соответствующую таблицуприближений. Первый столбец из j элементов будет полностью заполнен приближениями, полученными с помощью рекуррентной формулы трапеций.
Приближений по Симпсону будет на одно меньше, поскольку первое приближение алгоритматрапеций не имеет верхней, вычисленной при вдвое большем шаге, пары. Соответственно в столбце приближений, вычисленных по формуле Буля, не будет двухверхних элементов, в следующем столбце — трех и т. д. В результате будет получена нижняя треугольная матрица, элемент нижнего правого угла которой и является наиболее точным приближением (рис. 10.6).R(J,O)R(J, 1)R(J,2)ФормулатрапецийФормулаСимпсонаФормулаБуляОД 4)ОД, 3)ТретьеприближениеЧетвертоеприближениео1234R(4,0)^~~-~-» R(4,1)^ ^ ~ - > R(4,2) - — = = - R(4,3)^ " " " ^ R(4,4)Рис. 10.6. Схема интегрирования по Ромбергудля пяти приближенийЕсли отнять полученную при вычислении интеграла по Ромбергу матрицу от соразмерной, все ненулевые элементы которой — истинные значения подсчитываемогоинтеграла, поделить на величину истинного значения и умножить на 100, то в результате получится характеристическая матрица процентных ошибок приблизительно следующего вида:38.78940000021.46018-4.7197600005.19406-0.227990.071460001.28842-0.013460.32148 -8.2955х10~4^0.08033 -5.167x108.43455хЮ~4 -2.775х\0~ 41.2375x10"ч1.9x10-75-8.15х 10" 7,-9-5х 102.7хЮ" 7 010.6.
Численные методы интегрирования* 371Из этой матрицы видно, сколь быстро убывают ошибки последовательных приближений, вычисляемых с помощью описанных рекуррентных формул. Самих же вычислений при этом нужно сделать очень немного. Поэтому метод Ромберга не только очень точен, но и весьма быстр.4. В программной реализации алгоритма Ромберга удобно осуществлять одновременное нахождение всех возможных для данного шага приближений (то есть на каждом круге цикла while заполняется одна строка приведенной выше матрицы, начиная со второй):for k e 1..J4kRRI.k-l•U4k-l5. В качестве условия остановки цикла в методе Ромберга обычно используют критерий разности двух последних приближений, лежащих на диагонали матрицы ( R ( ;и R, ,, ,).
В том случае, если разность по модулю окажется меньше TOL, то цикл будет остановлен и в виде ответа будет выдано последнее приближение:break ifF,JRJ-1,J-1<TOLКак видите, реализация метода Ромберга, в общем, несложна. Однако все равно стоитпривести соответствующую программу в «собранном» виде:Romberg(f,a,b,TOL) := J < - 0while 1J < - J +1Rj .02J"1+ h- Vf[a + h-(2p - 1)]p=lfor k e 1..J4 RR' J,k-l~ J-l,k-lRJ,k^break ifk4 -l..<TOL3 7 2 • Глава 10. Вычисление интеграловПроверим, насколько качественно работает написанная программа.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
