WWW.DISSERS.RU

БЕСПЛАТНАЯ ЭЛЕКТРОННАЯ БИБЛИОТЕКА

   Добро пожаловать!


Pages:     || 2 |
Министерство образования Российской Федерации РОСТОВСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ С.С.Михалкович, А.В.Олифер, А.М.Столяр ЧИСЛЕННЫЕ МЕТОДЫ Выпуск IV Интегралы от быстро осциллирующих функций.

Многомерные интегралы Методические указания к выполнению индивидуальных заданий на ЭВМ для студентов 2 курса физического факультета Ростов-на-Дону 2000 Введение В данных указаниях рассматриваются методы приближенного вычисления интегралов от быстро осциллирующих функций и многомерных интегралов. Умение вычислять такие интегралы оказывается принципиально важным для решения некоторых задач радиофизики и ядерной физики соответственно, а также и в других областях.

Вычисление интегралов рассматриваемых типов не всегда удается осуществить с помощью готовых процедур, предлагаемых вычислительными пакетами типа Мaple, Мathcad и т.п. Когда же это оказывается возможным, то всегда существует некоторый риск того, что полученный ответ неверен.

В такой ситуации большую роль играет глубокое понимание сути численных методов, используемых в стандартных процедурах, и способность, при необходимости, самостоятельно запрограммировать нужный метод.

На развитие такого уровня понимания как раз и нацелены описанные ниже практические задания и рекомендации по их решению Общие замечания по поводу составления и тестирования программ для выполнения заданий смотри во Введении к [2].

1. Приближенное интегрирование быстро осциллирующих функций. Формула Филона 1.1. Задание 1) Написать процедуру, реализующую составную квадратурную формулу Филона приближенного интегрирования быстро осциллирующих функций. С помощью этой процедуры найти значение интеграла с точностью не менее 0.01% cогласно правилу Рунге. Зафиксировать N - число подынтервалов, потребовавшихся для достижения заданной точности.

4 2) Применить для вычисления того же самого интеграла составную формулу Симпсона с N подынтервалами.

3) Найти приближенное значение интеграла с помощью пакета Maple и определить, насколько от него отличаются значения, найденные в пунктах 1 и 2.

Тестовый пример:

x ei10x dx = cos 10 25 + 49sin 10 250 -01401909989.

.

( ) ( ) -Рекомендуемая литература: [3], гл.3, §7.

1.2. О методе Квадратурная формула Филона применяется для приближенного вычисления интегралов вида b f x exp ix dx, (1) ( ) ( ) a где выделены относительно медленно меняющийся на интервале [a,b] множитель f x и быстро осциллирующий на [a,b] множитель exp ix.

( ) ( ) Последнее условие формально записывается следующим образом:

b - a >>1. Функция f x предполагается гладкой. Заметим, что часто ( ) ( ) интеграл от быстро осциллирующей функции можно преобразовать к виду (1) с помощью той или иной замены переменной интегрирования.

Вывод формулы Филона аналогичен выводу квадратурных формул Ньютона-Котеса. Поскольку функция f x гладкая, то на [a,b] она прибли( ) жается с известной погрешностью интерполяционным многочленом. Пусть для определенности это интерполяционный многочлен в форме Лагранжа L3 x = P1 x f x1 + P2 x f x2 + P3 x f x3, ( ) ( ) ( ) ( ) ( ) ( ) ( ) построенный по узлам x1 = a, x2 = a + b 2, x3 = b. Как обычно, Pi x, ( ) ( ) i =12,3, суть многочлены второй степени не зависящие от функции f,, Pi x =1, если i = j, и Pi x = 0 иначе. В силу того, что f x L3 x ( ) ( ) ( ) ( ) j j b b f x exp ix dx L3 x exp ix dx = ( ) ( ) ( ) ( ) a a b- a a + b b- a = expi D f x. (2) ( ) j j 22 j=Интеграл от произведения многочлена L3 x на функцию exp ix вычис( ) ( ) лен в явном виде. Он представляет собой линейную комбинацию значений функции f x в точках x1, x, x3 с коэффициентами, не зависящими от ( ) функции f x. Это и есть формула Филона. Аналогично выводятся форму( ) лы более высокого порядка.

Оценка погрешности формулы Филона R f полностью определяется ( ) погрешностью интерполяции функции f x многочленом L3 x и совпада( ) ( ) ет с оценкой погрешности квадратурной формулы Симпсона для интеграла лишь от одной функции f x, без быстро осциллирующего множителя:

( ) b b R f = f x - L3 x exp ix dx f x - L3 x dx ( ) ( ) ( ) ( ( ) ( )) ( ) a a b )( ) max f x.

( )( - a ab, [ ] Составная квадратурная формула Филона строится так же, как и составная квадратурная формула Симпсона: исходный интервал разбивается на подынтервалы, на каждом из которых применяется формула (2). Итоговая расчетная формула, однако, несколько сложнее, чем в случае метода Симпсона. Это связано, в частности, с тем, что множитель exp i a + b / 2) для ( ) () каждого подынтервала свой.

1.3. Расчетные формулы В формуле (2) коэффициенты D1 p = p-3 2pcos p - sin p 2 - p2 + i p2 cos p - psin p, ( ) ( ) ( ) ( ) ( ) ( ) ( ) [] D2 p = p-3 4 sin p - 4 pcos p, ( ) ( ) ( ) [] D3 p = p-3 2pcos p + sin p p2 - 2 + i psin p - p2 cos p.

( ) ( ) ( ) ( ) ( ) ( ) ( ) [] Для приближенного интегрирования функций вида f x cos x ( ) ( ) ( f x sin x ) с вещественной функцией f x, очевидно, следует исполь( ) ( ) ( ) зовать вещественную (мнимую) части формулы (2). Например, b b - a cos a + b ReD j b - a f x j f x cos x dx = ( ) ( ) ( ) 22 j= a a + b b - a - sin ImD f x, ( ) j j j= где Re ( Im ) обозначают вещественную и мнимую части соответственно.



2. Приближенное интегрирование методом Монте-Карло 2.1. Задание Написать программу для вычисления интеграла методом Монте-Карло. Организовать вывод приближенных значений интеграла для различного числа испытаний в файл. Построить по этим данным график в среде Maple, отметив на графике точное значение интеграла (см. Рис.2). Сравнить результаты, получающиеся при использовании трех последовательностей случайных чисел, сгенерированных стандартным генератором псевдослучайных чисел, и генератором, описанным в приложении А, для “простейшего” и “геометрического” вариантов метода.

Тестовый пример:

2 2 1- x1 1- x2 1- x6 dx = 4 / 3, (3) ( )( ) ( ) ( )G где область G - шестимерный куб со стороной равной двум с центром в начале координат:

G = = x1,x2,...,x6 R | |xi|1, i = 1,..,6.

( ) {x } Рекомендуемая литература: [3], гл.5, §8; [4], гл. 3, §2.

2.2. О методе Метод Монте-Карло, или метод статистических испытаний, применяется главным образом для приближенного вычисления многомерных интегралов. Дело в том, что с ростом размерности интеграла d количество узлов для его приближенного вычисления по формулам типа Ньютона-Котеса растет пропорционально N d, где N число узлов, достаточное для приближенного вычисления одномерного “сечения” данного интеграла с требуемой точностью. В то же время для метода Монте-Карло число узлов, необходимое для вычисления подынтегральной функции, не зависит от порядка интеграла.

Поэтому, для достаточно больших значений d метод Монте-Карло оказывается единственным реально применимым вычислительным методом. Другими “показаниями” к применению данного метода являются сложная граница области интегрирования, отсутствие областей с резкими изменениями подынтегральной функции и умеренные требования к точности вычислений. Последнее обстоятельство является весьма существенным, так как ме-1/тод Монте-Карло имеет невысокую скорость сходимости, порядка N (см. (5)).

В любом случае, прежде чем принимать решение о применении того или иного метода необходимо попытаться понизить размерность исходного интеграла, например, за счет тех или иных симметрий области интегрирования и подынтегральной функции, и исследовать характер поведения подынтегральной функции, при необходимости оценив возможность выделения особенностей (ср. с [2]).

Характерной особенностью метода Монте-Карло является вероятностный характер получаемого с его помощью ответа: с ростом N повышается вероятность того, что вычисленное значение близк к истинному значению интеграла. Это утверждение справедливо лишь для “истинно случайныx” последовательностей чисел, используемых для теоретического обоснования метода. Поэтому для получения правильных результатов при практическом использовании метода Монте-Карло совершенно необходимо высокое качество генератора псевдослучайных чисел.

Итак, пусть имеется d - мерный интеграл d I= f x dx, G R. (4) ( ) G Предположим, что элементы последовательности i, i=1..N, - это наугад { } выбранные точки области G. Тогда, согласно “простейшему” варианту метода f - f I mes G f ± mes G, (5) ( ) ( ) N где mes( G) - мера области G, N N 1 f f i, f f i. (6) ( ) ( ) N N i=1 i=Второй член в формуле (5), как обычно при записи приближенного числа, характеризует (вероятную) ошибку приближения.

Суть формул (5), (6) проста. Точное значение интеграла I, разделенного на mes G, - это математическое ожидание случайной величины f x, по( ) ( ) рожденной случайной величиной x, равномерно распределенной в области G. Формула N 1 f x dx f i ( ) ( ) mes G N ( ) i=G (см. (5)) - это обычное приближение математического ожидания выборочным средним значением.

Приведем также “геометрический” вариант метода. Пусть для простоты интеграл (4) - одномерный: G = a,b, и 0 f x M при x G (Рис. 1).

( ) [ ] Выберем наугад N точек в прямоугольнике a,b M.

[ ] y M y = f x ( ) a b x Рис.1. Отношение площади криволинейной трапеции, ограниченной графиком функции y = f x, к площади прямоугольника b - a M при( ) ( ) ближенно равно доле точек, выбранных случайным образом в прямоугольнике, оказавшихся ниже графика функции. Чем больше взято точек, тем приближение точнее.

Пусть Nf из них оказалось ниже графика функции f x. Тогда значение ( ) интеграла, равное площади криволинейной трапеции, пропорционально отношению N N. Коэффициент пропорциональности равен величине f b ( - a M - площади прямоугольника.

) В многомерном случае поступают аналогичным образом. Пусть, 0 f x M, x G (Если, m f x M, m < 0, x G, то можно пе( ) ( ) рейти к функции g x f x - m ). Берут наугад N точек i, i=1..N, в ( ) ( ) { } области G. Для каждой такой точки выбирается случайное число i из интервала 0,M. Пусть Nf - это количество точек i, для которых [ ] { } i < f i. Тогда ( ) N f I= f x dx mes G M. (7) ( ) ( ) N G Проиллюстрируем оба варианта метода: “простейший” и “геометрический”, на тестовом примере (3). Сначала применялась формула (5). В данном случае мера mes G области G равна объему шестимерного куба со ( ) стороной равной двум, т.е. 64. Таким образом, приближенное значение интеграла вычислялось по формуле N I 64 f( i ). (8) N i=На рисунке 2 приведены приближенные значения тестового интеграла в зависимости от количества точек, выбранных с помощью встроенного и “внешнего” генераторов случайных чисел соответственно. Каждая точка - это шесть последовательных случайных чисел из интервала -11. Обратите, [ ] внимание на меньшие колебания с ростом N приближенных значений интеграла, найденных с помощью внешнего датчика, и на важность удачного выбора значения числа: RandSeed (для встроенного генератора) или idummy (для ran1), в зависимости от которого вычисляется та или иная случайная последовательность.





Рассмотрим теперь применение “геометрического” варианта метода (Рис.3). Так как подынтегральная функция в интеграле (2.1) принимает значения из интервала 01 на G, то, в соответствие с формулой (7),, [ ] N f 2 2 1- x1 1- x2 1- x6 dx 64. (9) ( )( ) ( ) N G Из сравнения результатов расчетов по формулам (8), (9) видно, что “геометрический” вариант метода Монте-Карло вычисления многомерных интегралов менее точен, чем “простейший”.

В любом случае сходимость достаточно медленная. Существуют различные подходы для ее ускорения. Прежде всего необходимо определить области имеются ли области резкого изменения подынтегральной функции.

Для таких областей интегрирование по методу Монте-Карло дает наибольшую ошибку. В простейшем случае интегралы по этим областям вычисляют отдельно. По существу это отражает естественное стремление брать больше случайных точек там, где функция быстро изменяется. Иногда из этих же соображений с самого начала используют подходящее неравномерное распределения случайных точек. Очень полезным оказывается выделение особенностей с помощью явно интегрируемых функций.

А Б Рис.2. Значения интеграла (3), найденные с помощью формулы (8) (“простейший” вариант метода), для различного числа случайно выбранных точек из области G. Каждая кривая отвечает отдельной последовательности случайных чисел, (А) сгенерированных с помощью датчика случайных чисел встроенного в Turbo Pascal, v.7.0 (Borland Intl) ;

RandSeed=1, 48322, -13049; (Б) датчика Парка-Миллера ran1 (приложение А); idummy=-1, -48332, -13049. Горизонтальная линия - точное значение интеграла.

A Б Рис.3. Значения интеграла (3), найденные с помощью формулы (9) (“геометрический” вариант метода Монте-Карло), для различного числа случайно выбранных точек из области G; (А) кривые полученные с помощью встроенного генератора случайных чисел; стандартная переменная RandSeed=1, 48322, -13049; (Б) с помощью генератора ran1;

idummy=-1, -48332, -13049.

Приложение А Существует немало хороших алгоритмов вычисления “случайных” чисел (см., например, [5] ). В основе многих из них лежит линейный конгруэнтный метод, согласно которому последовательность случайных чисел получается из соотношения Xn+1 =(a Xn +c) mod m, n>=0, для некоторых специально подобранных значений a, c, m; значение X0 - достаточно произвольно. Ниже приведен фрагмент кода, реализующий с небольшими модификациями один из вариантов данного метода, - “минимальный” генератор случайных чисел Парка-Миллера (Park, Miller) c a==16807, с=0, и m=231 -1= 2147483647 ([6]).

const IA=16807; IM=2147483647; AM=1.0/IM;

IQ=127773; IR=2836; MASK=123459876;

NTAB = 32;

NDIV = (1+(IM-1) div NTAB);

EPS = 1.2e-7; RNMX = (1.0-EPS);

function ran1 (var idum:longint):real;

{ начальное значение idum должно быть ОТРИЦАТЕЛЬНЫМ; } { между последоватнльными обращениями к ran1 значение idum } { не должно меняться } const iy: longint = 0;

iv: array [0..NTAB-1] of longint =(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0);

var k,j:longint;

temp:double;

begin if ((idum <= 0) or (iy=0)) then begin if (-idum)<1 then idum:=1 else idum:=-idum;

for j:=NTAB+7 downto 0 do begin k:=idum div IQ;

idum:=IA*(idum-k*IQ)-IR*k;

if (idum < 0) then idum:=idum+IM;

if (j

end;

iy:=iv[0];

end;

k:=idum div IQ;

idum:=IA*(idum-k*IQ)-IR*k;

if (idum < 0) then idum:=idum+IM;

j:=iy div NDIV;

iy:=iv[j];

iv[j]:=idum;

temp:=AM*iy;

if (temp>RNMX) then ran1:=RNMX else ran1:=temp;

end;

Варианты заданий по теме “Интегралы от быстро осциллирующих функций” 2 2 1. x e-2x sin 40x dx 2. x e-2x cos 40x dx ( ) ( ) 0 1 3. ln 1+ x sin 60x dx 4. 1+ x cos 60x dx ( ) ( ) ( ) 0 2 2 5. x sin 40x dx 6. 1+ x cos 60x dx ( ) ( ) -1 1 2 7. x e-x sin 45x dx 8. x sin 50x dx ( ) ( ) -2 1 3 9. x e-x cos 40x dx 10. x e-2x sin 40x dx ( ) ( ) -2 2 4 11. 1+ x sin 40x dx 12. x e2x sin 40x dx ( ) ( ) -2 -2 13. x e2x cos 40x dx 14. ln x + 3 cos 50x dx ( ) ( ) ( ) -2 -2 15. x e2x cos 50x dx 16. xe2x cos 45x dx ( ) ( ) -2 -3 17. e-3x cos 45x dx 18. xe-2x cos 30x dx ( ) ( ) -1 2 19. 1+ x cos 45x dx 20. ln 1+ 2x cos 45x dx ( ) ( ) ( ) -2 2 21. xe-2x cos 50x dx 22. 1+ x sin 45x dx ( ) ( ) 0 2 23. ln 4 - x sin 30x dx 24. 1+ x sin 50x dx ( ) ( ) ( ) 1 -2 25. xe-x cos 60x dx 26. xe-2x cos 50x dx ( ) ( ) 0 Варианты заданий по теме “Многомерные интегралы” Во всех заданиях интегралы рассматриваются по области G, представляющей собой четырехмерный куб со стороной s и с центром в точке С=(c1, c2, c3,c4 ).

Pages:     || 2 |










© 2011 www.dissers.ru - «Бесплатная электронная библиотека»

Материалы этого сайта размещены для ознакомления, все права принадлежат их авторам.
Если Вы не согласны с тем, что Ваш материал размещён на этом сайте, пожалуйста, напишите нам, мы в течении 1-2 рабочих дней удалим его.