Введение в вычислительную физику - Федоренко Р.П.
ISBN 5-7417-0002-0
Скачать (прямая ссылка):
о
где f(t) — гладкая функция. Регуляризация состоит в том, что проделывается тождественное преобразование
і і і
J /(<) Г1'2 dt = J [/(0 - <р(0] Г112 dt -I- J «р(0 Г112 dt.
О 0 0
Функция <p(t) выбирается такой, чтобы первый интеграл правой ча-
сти не содержал особенности и при небольшом объеме вычислений достаточно точно определялся хотя бы по формуле Симпсона. Второй интеграл особенность содержит, но вычисляется аналитически.
54
ОСНОВЫ ВЫЧИСЛИТЕЛЬНОЙ МАТЕМАТИКИ
В данном случае цель будет достигнута, если в качестве <р(0 взять
отрезок ряда Тейлора /(О в точке t = 0. Это приводит к вычис-
лению
j dt + /(0) j j= dt + /40)S VTdt.
О 0 0
В примере, с которого мы начали (/(0 = cos 0> приходим к вычислению
і і
I (cos / — I) t~112 dt + J t~112 dt. о о
Второе слагаемое есть 2, первое вычислим по формуле Симпсона: сначала с шагом 0.5, что даст значение 1.807967, затем с шагом 0.25, что даст значение 1.808850. Эти вычисления «стоили» всего четырех вычислений подынтегральной функции. Поучительно сравнить их с тем, сколько вычислений этой функции потребуется при «студенческом» рецепте (для достижения такой же точности).
Вычисление интегралов от быстроосциллирующих функций.
Начнем с простого примера. Пусть требуется вычислить
Jt
J е~г sin kt dt о
при большом значении к, например к = 100. Интегралы типа
J /(0 sin kt dt,
где /(0 — гладкая функция, часто приходится вычислять в некоторых разделах физики. Сложность задачи состоит в том, что подынтегральная функция совершает большое число колебаний. Вычисление интеграла по стандартной формуле Симпсона, конечно, возможно, но требует сетки с очень малым шагом: каждая волна должна быть описана некоторым числом узлов сетки, а волн много.
Дело осложняется еще и тем, что вычисление должно проводиться с высокой точностью, так как результат есть сумма большого числа близких величин с противоположными знаками (интегралов от отдельных волн подынтегральной функции), происходит сильное сокращение знаков и для обеспечения точности остатка (результата) отдельные слагаемые должны вычисляться с существенно более высокой точностью. Для вычисления подобных интегралов используется следующий прием: гладкая функция f(t) аппроксимируется некоторой другой гладкой функцией f(t), такой, чтобы интеграл от /(О sin kt вычислялся аналитически.
§ 41
ВЫЧИСЛЕНИЕ ОПРЕДЕЛЕННЫХ ИНТЕГРАЛОВ
55
Итак, дело сводится к тождественному преобразованию
JC 31 JC
^ /(О sin kt dt = ^ f(t) sin k:t dt + $[/(0 — 7(01 sin kt dt.
ООО
Второе слагаемое является малым и отбрасывается. Правда, если оценить отбрасываемую величину, опираясь только на оценку типа 1/(0 -/(01 < Е> те- в дан-
16
8 4 8 4 8 4 8
16 8 16 8 16 8
8 4 8 4 8 4 8
16 8 16 R 16 3
1
1
1 4 2
Рис. 7
4 2 б
4 2 4
ном случае величинои ЛЕ, ничего хорошего (даже если E — точная оценка погрешности аппроксимации) не получится, так как величина ле может оказаться значительно большей интересующего нас интеграла. На самом деле погрешность существенно меньше. Это ведь интеграл от гладкой функции, не превосходящей е, умноженной на быстроосциллирующую функцию. Естественно ожидать, что погрешность будет во столько раз меньше результата, во сколько раз |/ — /| меньше /. При
/(О = е~г интеграл вычисляется аналитически. Поучительно вычислить его приближенно, заменив функцию е~1 ее интерполяционным полиномом всего лишь второй степени. Интересно, совпадает ли точность результата с ожиданиями? На этой идее построена формула механической квадратуры, аналогичная формуле4 Симпсона. Интервал интегрирования разбивается на четное число шагов длиной х, на каждой паре шагов функция заменяется интерполяционным ПОЛИНОМОМ второй степени L2(O) интегралы от L2(0 sin kt вычисляются аналитически, полученные выражения суммируются.
Вычисление многомерных интегралов. Метод Монте-Карло.
Рассмотрим задачу вычисления интеграла по многомерному кубу -
Iii
J J ...( /(X1, х2, ..., хп) dx{ dx2... dxn.
0 0 о
Нетрудно и здесь построить формулы механических квадратур, аналогичные, например, формуле Симпсона^ Проще всего такие формулы получить, используя процедуру повторного интегрирования, т.е. заменяя многомерный интеграл на равный ему
і 1 і
J dxt J dx2 ... J f(xt, х2, ..., хп) dxn.
0 0 0
56
ОСНОВЫ ВЫЧИСЛИТЕЛЬНОЙ МАТЕМАТИКИ
[Ч.І
В двумерном случае для «элементарной ячейки» имеем коэффициенты (аналогичные коэффициентам 1,4, 1 в одномерной формуле Симпсона), изображенные на рис. 7а. Суммируя интегралы по элементарным ячейкам, получаем коэффициенты, частично изображенные на рис. 16. Эти коэффициенты поставлены около точек сетки, покрывающей двумерный квадрат. Вычисление интеграла по такой квадратурной формуле состоит в вычислении значений / в узлах сетки, умножении на соответствующий коэффициент и суммировании результатов.
Аналогичные формулы можно построить и в кубах бблыыих размерностей, но пользы от этого мало. Дело в том, что при росте размерности объем вычислений катастрофически растет. Операция интегрирования, справедливо считающаяся одной из самых элементарных в анализе, практически является одной из самых сложных, точнее трудоемких. Поэтому описанные выше квадратуры, основанные на регулярных сетках, используются в практических вычислениях лишь для двумерных и трехмерных пространств.
