Турчак Л.И. Основы численных методов. Под ред. В.В.Щенникова (1987) (1095857), страница 17
Текст из файла (страница 17)
Таким образом, погрешности метода Симпсона и формулы (ЗЛЗ), использующей методы прямоугольников и трапеций, имеют одни порядок. Отметим разницу между формулой С нмпсона (3.34) и комбинацией методов прямоугольников и трапеции (3.33) . Эти формулы имеют одинаковую погрешность, однако формула (3,33) требует двукратного вычисления интеграла разиымп методами. Кроме того, для метода Симпсона нужно почти вдвое меньше табличных значений фунт;цип, поскольку для метода прямоугольников нужны дополнительные данные в Блок-схема одного пз простейших лепия определенного интеграла по представлена на рис. 14, Б качеств задаются границы отрезка интегриро Рнс. 14. Ьлок-схема метода Симпсона юлуцелых точках.
алгоритмов вычисматоду Симпсона е исходных данных ванпя а, Ь, погреш- 1С~2 Гл. 3. лпФФНРенцпРОВлнпГ и пнтГГРНРОвлния г = 1, 2. . . и. х;,<х<х;, Выражение для интеграла представим в виде ь и '+1 и х1+1 х=~!И~ =Х ~!(*)ь* Х 1 и(*)н* О 1=1 х1 1=1 х.
1-1 Используя выражение (3,35), в результате вычисления интегралов находим и (а;Ь;+ -В;Й, + -с;В~ ~- -Н~й~~). (3.36) 1=1 Способ вычисления коэффициентов а;, Ь;, е;, И~ описан в гл. 2. Здесь лишь отметим, что а; = у;,. Для практических расчетов формулу (3.36) можно представить в виде и и 1 ~ —, ~ й; (У; 1 + У,) — —, ~ 1г~ (с; 1 + сь), (3,37) 1=1 ";=1 ность е, а также формула для вычисления значений подынтегральной функции у = 1(х). Первоначально отрезок [а, Ь1 разбивается на четыре части с шагом й= =(Ь вЂ” а)/4.
Вычисляется значение интеграла 1~. Потом число шагов удваивается, вычисляется значение 1, с шагом 6/2. Условие окончания счета принимается в виде !1, — Х~! (е. Если это условие не выполнено, происходит новое деление шага пополам и т, д. Отметим, что представленный на рис. 14 алгоритм не является оптимальным. В частности, при вычислении каждого последующего приближения 1, не используются значения функции ~(х), уже найденные на предыдущем этапе.
Более экономичные алгоритмы будут рассмотрены в п. 5. 4. Использование спла11нов. Одним из методов численного интегрирования, особенно эффективным при строго ограниченном числе узлов, является л1етод силайнов, пспользую1цпй интерполяцию сплайнами (см. гл. 2, Я 3, п.
2). Разобьем отрезок интегрирования ~а, Ь] на и частей точками х;. Пусть х; — х;, = Ь; (1=1, 2, ..., и). На каждом элементарном отрезке пнтерполпруем подынтегральную функци1о ~(х) с помощью кубического сплайна: гр;(х) = а; + Ь;(х — х,,) + с,(х — х;,)' + сг;(х — х,,)', (3.35) й 2. ЧИСЛЕННОЕ ИНТЕГРИРОВАНИЕ Анализ этой формулы показывает, что первый член в правоп части совпадает с правой частью формулы (3.30) для метода трапеций. Следовательно, второй член характеризует поправку к методу трапеций, которую дает использование сплайпов. Как следует пэ формулы (3.35), коэффициенты с; выражаются через вторые производные д;(х): сг = ~ я~1 (хг-1) ~~ ~ у1 — 1ю Это позволяет оценить второй член правой части формулы (3.37): Ьз 1 3 '1 , (С; 1 + С,) У~, где у* — вторая производная в некоторой внутренней точке.
Полученная оценка показывает, что добавка к формуле трапеций, которую дает использование сплайнов, компенсирует погрешность самой формулы трапеций. Отметим, что во всех предыдущих методах (см. и. 2, 3) формулы численного интегрирования можно условно записать в виде линейной комбинации табличных значений функции: ь и / (х) Их = —,'~~ сг д .- а 1=0 При использованпи сплайнов такое представление невозмо>кно, поскольку сами коэффициенты сг~ зависят от всех значений у;. Рассмотрев разные методы численного интегрирования, трудно сравнивать их достоинства и недостатки.
Любая попытка такого сравнения непременно поставит перед нами альтернативный вопрос: что больше, Ь'у" или Ь'у"'? Все зависит от самой функции у = ~ (х) и поведения ее производных. Уточнение результатов численного интегрирования можно проводить по-разному. В частности, в представленном на рпс. 14 алгоритме с использованием метода Симпсона проводится сравнение двух значений интеграла 1, и 1„полученных при разбиениях отрезка ~а, 6~ соответственно с шагами Ь и Ь!2, Аналогичный алгоритм можно построить и для других методов.
Здесь мы упомянем другую схему уточнения значения интеграла — процесс Эйткена. Он дает возможность 104 гл. 3. диФФеРннциРОВлние 11 инткгРироВАние оценить погрешность метода О (Й") и указывает алгоритм уточнения результатов. 1'асчет проводится последовательно трн раза прн различных шагах разбиения Ь„ /г„7г„ причем их отношения постоянны: Ь,/Й, = А;/Ь, = а (например, при делении шага пополам о = 0.5).
Пусть в результате численного пптегрпроваппя получены значения шгтеграла 1„1~, 7з. Тогда уточненное значение интеграла вычисляется по формуле а порядок точности используемого метода численного интегрирования определяется соотношением Уточнение значения интеграла можно также проводить методом Рунге — Ромберга (см. ~ 1, и.
5). 5. Адаптивные алгоритмы. Из анализа погрешностей методов численного интегрирования следует, что точность получаемых результатов зависит как от характера изменения подынтегральной функции, так и от шага интегрирования. Будем считать, что величину шага мы задаем. При этом ясно, что для достижения сравцпмой точности при интегрировании слабо меняющейся функции шаг можно выбирать большим, чем при интегрировании резко меняющихся функций. 11а практике нередко встречаются случаи, когда подынтегральная функция меняется по-разному на отдельных участках отрезка интегрирования. Это обстоятельство требует такой организации экономичных численных алгорптмов, прп которой они автоматически приспосабливались бы к характеру изменения функции. Такие алгоритмы называготся адаптивныжи (приспосабливагоигизгися) . Онп позволяют вводить разные значения шага интегрирования на отдельных участках отрезка интегрирования.
Это дает возможность уменьшить машинное время без потери точности результатов расчета. 11одчеркнем, что этот подход используется обычно при задании подынтегральпой функции у=1(х) в виде формулы, а пе в табличном виде. Программа, реализующая адаптивный алгоритм численного интегрирования, входит обычно в виде стандарт ~ 2 чпслвнное пнтГгР11РОВАнин 105 пой подпрограммы в математическое обеспечение Э))М.
Пользователь готовой программы задает границы отрезка интегрирования а, Ь, допустимую абсолютную погрешность е и составляет блок программы для вычисления значения подынтегральной функции ~(х) . Программа вычисляет значение пптеграла 1 с заданной погрешностью е, т. е. 1 — ) ~(х) сЬ (3.38) (е, Разумеется, не для всякой функции можно получить результат с заданной погрешностью. Поэтому в программе может быть предусмотрено сообщсние пользователю о недостижимости заданной погрешности. Интеграл при этом вычисляется с максимально возможной точностью, и программа выдает эту реальную точность.
Рассмотрим принцип работы адаптивного алгоритма. Первоначально отрезок ~а, Ь1 разбиваем на и частей. В дальпейшем каждый такой элемептарпып отрезок делим последовательно пополам. Окончательное число шагов, их расположение и размеры зависят от подынтегральиой функции и допустимой погрешности е. К каждому элементарному отрезку ~х; „х,] применяем формулы численного интегрирования при двух различных его разбиениях. Получаем приолпжения 1;, Х; (1) (2) для интеграла по этому отрезку: 1; = ) УЯЫх.
(3.39) Полученные значения сравниваем и проводим оценку их погрешности. Если погрешность находится в допустимых границах, то одно из этих приближений принимается за значение интеграла по атому элементарному отрезку. В противном случае происходит дальнейшее деление отрезка и вычисление новых приближений. С целью эко-' номии машинного времени точки деления располагаются таким образом, чтобы использовались вычисленные значения функции в точках предыдущего разбиения. Папример, при вычислении интеграла (3.39) по формуле Симпсона отрезок ~х;,, х,1 сначала разбиваем на две части с шагом Ь;!2 и вычисляем значение 1~ ° (1) ~щ гл. 3. ДиФФеРенц11РОВАние и пнтегРИРовхние (а) Потом вычисляем 1, с шагом Ь;/4. Получим выражения 1';" = — ' /(х(,) + 4У Х1 — 1 + —.
+ У(хз) уР) - '(((х.;,) .~- 4( (х,, ~- — ') + + 2/ ..;, + —,.' + Ч; + —,* + 1(;), (3,4() (3.40) Формулу (3.41) можно также получить двукратным применением формулы (3.40) для отрезков [х; „х~,+ Ь/21 и 1х; + г1;/2, х;1. Процесс деления отрезка пополам и вычисления уточненных значений 1; и 1; продолжается до тех пор, (() (а) пока пх разность станет не больше некоторой заданной величины б;, зависящей от е и Ь: 11'," — 1",-'! ~ б,. (3.42) Аналогичная процедура проводится для всех п элеи ментарных отрезков. Велцчина1=,~ 1, принимается в 1=1 качестве искомого значения интеграла.