Гурский Д., Турбина Е. - Вычисления в MathCad 12 (1077322), страница 83
Текст из файла (страница 83)
Иллюстрация метода средних прямоугольников показана на рис. 10.5, а.3 6 2 • Глава 10. Вычисление интеграловаОвРис. 10.5. Иллюстрация основных методов численного интегрированияСовсем не трудно вывести общую формулу данного метода:а + h-i н—2Здесь п — количество разбиений интервала интегрирования, h — ширина среднего прямоугольника. Для того же, чтобы реализовать в системе Mathcad такой простой алгоритм, которым является метод средних прямоугольников, даже не придется прибегатьк использованию языка программирования:Проверим, правильно ли работает заданная функция:func(x) := sin(x)-xfunc(x)dx float, 10 -> .265488009Pij func,0,-,1000] =0.2654879007Проверка показывает, что при больших п метод средних прямоугольников можетбыть весьма точен.
Так, в рассмотренном примере 1000 разбиений дали результат, верный до пяти знаков. Для большинства практических задач это более нежели достаточная точность. А вычислить функцию 1000 раз и просуммировать полученные значения — это для современного компьютера дело, требующее долей секунды. Поэтомуможно было бы остановиться на методе средних прямоугольников как на вполне приемлемом для подсчета интегралов на небольшом интервале. Однако если с помощьюрассматриваемого алгоритма понадобилось бы вычислить с той же точностью тройнойинтеграл, то количество разбиений возросло бы до 10 000 000. А подсчитать столькораз значение сложной функции — это уже дело нескольких часов (а лет 20 назад этобыла бы вообще малоразрешимая задача).
Приведенные цифры иллюстрируют очевидную необходимость найти какой-то другой метод численного интегрирования, ко-10.6. Численные методы интегрирования* 363торый, при сохранении точности, требовал бы меньшего количества узловых значенийфункции.Уменьшить количество разбиений можно прежде всего, увеличивая точность приближения функции легко интегрируемой аналитически кривой. Наиболее простая из таких кривых — это прямая.Метод, основанный на приближении фрагментов кривой функции секущими, проходящими через границы интервалов разбиения, называется методом трапеций (см.
рис. 10.5, б).Точность данного метода близка к точности метода средних прямоугольников, и величина ошибки пропорциональна квадрату ширины интервала разбиения. Общая формула метода трапеций будет вполне очевидна, если вспомнить, каким образом в геометрии вычисляется площадь этой фигуры:1(i-l)]]Задать в Mathcad функцию, вычисляющую интеграл по методу трапеций, столь жепросто, как и в случае метода средних прямоугольников:bjt(f,*,b,n):-V^•i+al+Jb-a-b-afunc(x) := sin(x)-xfimc(x)dx float, 10 -> .265488009Int ftmc,0,-,1000 =0.2654882245I3Точность результата вычисления интеграла по методу трапеций, полностью в соответствии с теоретическими выводами, оказалась крайне близкой к точности метода средних прямоугольников.
Соответственно точность данного метода нас вряд ли можетудовлетворить.Невысокая точность метода трапеций легко объясняется тем, что все-таки приближение кривой функции ломаной из отрезков секущих является довольно грубым дажедля значительного количества разбиений. Чтобы уменьшить ошибку интерполяции,попробуем использовать для приближения функции полином второй степени — параболу.Идея вычисления интеграла с помощью парабол довольно проста: промежуток разбивается на неширокие интервалы и на каждом из них вычисляется площадь, ограниченнаяинтерполирующей функцию параболой (см. рис. 10.5, в). Чтобы задать приближающую параболу, нужно знать коэффициенты ее уравнения.
Конечно, можно просто разложить в средней точке интервала функцию в ряд Тейлора — но это, пожалуй, слишком сложно и не очень надежно. Гораздо проще можно справиться с данной задачей,использовав возможности Mathcad по решению систем уравнений. Так как уравнениепараболы имеет три коэффициента, то, исходя из известного правила алгебры, чтобыих найти, нужно задать систему соответственно из трех уравнений.
Сделать это можно,3 6 4 • Глава 10. Вычисление интеграловвычислив значение приближаемой функции в трех точках интервала. Наиболее просто выбрать в качестве таковых границы и середину интервала разбиения. В результате получим следующую систему уравнений:АхО2 + ВхО + С=ГОА-(хО+ h ) 2 + В(хО + h) + C= flА(хО+ 2-h)2 + В-(хО+ 2-h) + C = f2Здесь хО — левая граница интервала; h — шаг, соответствующий расстоянию междуузловыми точками; fO, fl, f2 — значения функции на левой границе, в середине, на правой границе интервала; А, В, С — коэффициенты уравнения параболы.Система является линейной относительно коэффициентов, поэтому у нее будет только одно решение.
Это означает, что через три точки можно провести одну единственную параболу (справедливо и обобщение — через п точек проходит только один полином степени п-1). Решить заданную систему очень просто или посредством оператораsolve, или матричным способом. Ввиду объемности выражений для коэффициентов,создадим переменные А, В и С, которые будут на них ссылаться:ВхО+С=ЮА-?COF:=А-(хО+ h) + B-(xO+ h) + C = flsolve, А, В, С_А-(хО + 2-h)2 + B-(xO+ 2-h) + C = f2_В := (COF 1 ),А :=С :=После того как коэффициенты приближающей параболы будут определены, элементарно можно найти площадь, которую ограничивает ее фрагмент на интервале интегрирования:xOf2hА-х + В-х+ Сdxsimplify -» --h-(fO+ 4-flхОПолученная общая формула называется формулой Симпсона, и с ее помощью можнонаходить значение интеграла, и не решая систему уравнений для каждого интерваларазбиения.
Достаточно вычислить три узловых значения интегрируемой функции.Реализация метода Симпсона в Mathcad будет выглядеть следующим образом:пъ*3-пfunc(x) :=sin(x)-x10.6. Численные методы интегрированияfunc(x)dxfloat,15 -> .26548800861814* 365Simd func,0,-,100 = 0.265488008580604v•'озуПри количестве разбиений, в 10 раз меньшем, чем было использовано для расчета интеграла методами средних прямоугольников и трапеций, точность результата оказалась на 5 (!) порядков выше. А это означает, что метод Симпсона можно уже использовать и для вычислений кратных интегралов.Приближать интегрируемую функцию можно полиномами и более высокой степени.Обычно для этого используются полиномы 3-й и 4-й степени (описываются соответственно формулой Симпсона 3/8 и формулой Буля). Получить формулы для них можно точно так же, как мы получили формулу Симпсона.
В общем случае, чем выше степеньинтерполирующего полинома, тем точнее приближение. Так, ошибка при вычисленииинтеграла по формуле Буля убывает пропорционально 6-й степени величины интервала разбиений (формула Симпсона на два порядка менее точная). Из-за ряда технических сложностей полиномы выше 4-й степени используются редко.Чтобы более объективно сравнить эффективность работы рассмотренных методов,определим подбором то количество разбиений, которое требуется для каждого из них,чтобы приблизить интеграл с заведомо известным аналитическим решением с точностью до 10-го знака мантиссы:Г f(x) dx -> 2f(x) := sin(x)Pr(f,0,rc,100000)-2 = 8.2235551701614895x10"ПInt (f,O,Ji, 130000) - 2 = -9.7301278145778269x10"ПSimp (f,0,rt,170)-2 = 8.099x ю " 1 1Полученные цифры не могут не впечатлять: при использовании метода параболы разбиений требуется почти в 1000 раз меньше, чем для методов средних прямоугольников и трапеций. Для достижения же стандартной точности (0.001) будет достаточнымразбить промежуток всего на 4 (!) интервала:Simp (f,0,Ji,4) = 2.00026917Благодаря такой высокой эффективности формулы Симпсона до появления компьютерной техники ее использовали для подсчета интегралов, которые нельзя было вычислить аналитически, на бумаге.
Кстати, если использовать, например, формулу Буля,то разбиений потребуется еще меньше.В теории численных методов имеются формулы, используя которые можно предсказать,какой будет порядок у ошибки при вычислении интеграла от функции f(x) на промежутке от а до b с шагом h. Соответствующая формула для метода трапеций имеет вид:2(Ъ - а)123 6 6 • Глава 10. Вычисление интеграловФормула для ошибки метода Симпсона будет схожей:(Ь - а)dxOVh 7 =180Данные формулы показывают, что ошибка зависит прежде всего от величины шага.Также на нее будет влиять ширина интервала интегрирования (чем больше шагов придется проделать алгоритму, тем сильнее будут сказываться накапливающиеся погрешности) и скорость изменения функции (чем она выше, тем больше будет погрешность).Количественно скорость изменения функции можно узнать, вычислив производную.Однако производная дает мгновенную скорость, а нам нужно знать среднюю скоростьизменения функции на промежутке интегрирования.
Поэтому необходимо найти среднее значение производной на промежутке. Именно эта величина фигурирует в приведенных выше формулах.Формулы ошибки методов численного интегрирования важны в связи с тем, что онипозволяют определить, насколько малым должен быть шаг, чтобы была достигнутанеобходимая точность. Например, попробуем узнать, какой величины нужно сделатьшаг, чтобы вычислить интеграл от f(x)=sin(x) на промежутке от 0 до я методом Симпсона с точностью до 5-го знака.Выразив из формулы ошибки метода Симпсона h, получим следующее неравенство:180TOLh<(Ъ - а)-±-Лх)dxВсе входящие в данное выражение величины, кроме среднего значения производной,нам известны. Чтобы найти недостающую величину, рассчитаем четвертую производную аналитически, после чего найдем среднее значение полученной функции на промежутке, используя следующее замечательное свойство интеграла:1b-af(x)dxРасчет:—-jsin(x) ~^ sin(x)dxsin(x) dx —>• —0Зная среднее значение производной на промежутке, находим величину шага и количество шагов::sin(x) —> sin(x)dxi Г S1sin(x).
dx —> —n10.6. Численные методы интегрирования* 367Проводим интегрирование по Симпсону, опираясь на полученный выше результат.Ответ проверяем посредством аналитического интегрирования:f(x) := sin(x)rnSimp(f,0,71,10]) = 2.00000000065sin(x)dx-> 2J0На практике при определении необходимой величины шага среднее значение производной обычно не вычисляется ввиду того, что для этого необходимо провести интегрирование. Вместо него в формулу для ошибки подставляется наибольшее значение(по модулю) производной на промежутке. Неравенство для h от этого не нарушается,только величина шага получается слегка заниженной.
Так, с учетом, что на промежутке от 0 до л; максимальное значение синуса — 1, мы получим следующий результат:18010=0.028о.О28Как видите, величина шага при точном нахождении среднего значения производнойи замене его наибольшим ее значением различается мало — всего на 10 %.Проверка показала, что формулы для определения шага интегрирования вполне эффективны. Однако их применение при реализации алгоритмов численного интегрирования невозможно из-за сложностей, с которыми связано вычисление среднего значения производной или даже ее максимальной величины.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
