Лабораторная работ №2 (Численное дифференцирование. Численное интегрирование. Быстрое преобразование Фурье) (834487), страница 2
Текст из файла (страница 2)
Получен ответ на вопрос о том, существует лиоптимальный шаг интегрирования для данной формулы,минизимирующий достижимую погрешность.3. С помощью теоремы о корнях многочленов Лежандра, доказаннойв лекциях, выведена квадратура Гаусса, имеющая степень точности 5.Определено количество узлов необходимых для использования такойквадратуры.4. Написана функция gauss_quad5(f ) численного интегрирования функции fс помощью квадратуры Гаусса пятой степени точности.5. Доказано, что квадратура Гаусса имеет степень точности 5.8Быстрое преобразование Фурье (БПФ):1. Написана функция fft_coeff(y_nodes), которая вычисляет и возвращаеткомплексные коэффициенты тригонометрического полинома,интерполирующего узлы y_nodes, равномерно распределенные на отрезке[− π ; π ] .2. Протестирована корректность результатов работы функцииfft_coeff(y_nodes) с помощью БПФ для функции f1(x). С помощьювыкладок из лекций объяснено, как связаны возвращаемые комплексныекоэффициенты (и их индексы) с исходной функцией.3.
Написана функция trigonometric_interpolant(x, coeffs), которая вычисляетзначение тригонометрического полинома с коэффициентами coeffs вточке x.4. Произведена тригонометрическая интерполяция функции f2(x) для N = 2n ,где n ∈ 1, . . . , 8 и выведены результаты в виде графиков.Проанализирована непрерывность функции f2(x) и исходя из графиковсделан вывод о сходимости подобного приближения.5. Проделаны аналогичные шаги для функции f3(x).91. Численное дифференцированиеВыведем общую центральную формулу численного дифференцирования4-го порядка, аппроксимирующую первую производную по 5 узлам, вместе состаточным членом. Для этого допустим, что нам известны значения функцииf(x) в точках x1 - 2h, x1 - h, x1 , x1 + h, x1 + 2h, и разложим ее в ряд Тейлора в точкеx1 :(2 )(3)(4)(5)f ( x1 )2 f ( x1 )3 f ( x 1)4 f ( ξ)5f (x)=f ( x 1)+ f '(x 1)( x−x1 )+( x−x1 ) +( x −x1 ) +( x−x 1) +( x−x1 ) .2624120где ξ ∈ (x1 ; x).
Тогда значения ряда в точках x1 - h и x1 + h будут равны:2f ( x1 −h)=f ( x1 )−hf ' ( x1 )+2f ( x1 +h)= f ( x 1 )+ hf ' ( x 1)+345h (2)h (3)h (4)h ( 5)f ( x 1 )− f ( x 1)+ f ( x 1 )−f (ξ 1 ) ,2624120345h (2)hhh (5 )f ( x 1 )+ f (3) ( x 1 )+ f ( 4) ( x 1 )+f (ξ 1 ) .2624120(1)(2)где ξ1 ∈ (x1 - h; x ) и ξ2 ∈ (x ; x1 + h). Вычтем из равенства (1) равенство (2):3f ( x1 −h)−f ( x 1 + h)=−2 hf ' ( x 1 )−5h ( 3)hf ( x 1 )− ( f (5 ) (ξ 1 )+f ( 5) ( ξ 2 )).360(3)Предположим, что f(x) ∈ C 5 [x1 - h; x1 + h].
Тогда по теореме опромежуточном значении существует такое ξ3 ∈ (x1 - h; x1 + h), что1f 5 ( ξ 3 )= [ f ( 5) ( ξ 1 )+ f (5 ) (ξ 2 )] ,2(4)тогда равенство (3) примет вид:f ( x1 −h)−f ( x 1+ h)=−2 hf ' ( x 1 )−h3 (3)h5f ( x 1 )− ( f (5 ) (ξ 3 )).330(5)Аналогично получим значения ряда в точках x1 - 2h и x1 + 2h:f ( x1 −2h)=f ( x 1 )−2 hf ' ( x 1)+2h2 (2)4 h3 (3 )2h 4 (4 )4 h5 (5)f ( x 1)−f ( x 1)+f ( x 1)−f ( ξ4 ) ,23315(6)f ( x1 +2 h)= f ( x 1)+2 hf ' ( x 1 )+ 2h 2 (2)4 h 3 (3 )2 h 4 (4 )4 h 5 (5)f ( x 1 )+f ( x1 )+f ( x 1)+f (ξ 5) .23315(7)10Аналогично предыдущему выводу вычтем из равенства (6) равенство (7)и воспользуемся теоремой о промежуточном значении.
Получим равенство:816f (x1 −2h)−f (x 1 +2 h)=−4 hf ' (x 1)− h3 f (3 ) (x1 )− h5 f (5) (ξ6 ).315(8)Сложим равенство (5), умноженное на 8, с равенством (8), умноженнымна (-1) и в результате выполнения арифметических операций получим искомуюобщую центральную формулу численного дифференцирования 4-го порядка:4 5 (5 )8 f (x1 −1)−8 f ( x 1 +h)−f (x 1−2 h)+ f (x 1 +2 h)− h f (ξ)5f ' ( x1 )=−12 h(9)f ' ( x 1 )=⇒12211 4 ( 5)f ( x 1 − 2 h)−f ( x 1 − h)+ 0⋅f ( x 1 )+f ( x 1 + h )−f ( x 1 + 2 h )+h f (ξ)12 h3h3h12 h15Формула действительно имеет 4-й порядок точности, так шагдифференцирования h имеет 4-ую степень в выражении остаточного члена.Разработана программа, которая позволяет получить log-log графикзависимости абсолютной погрешности численного дифференцирования от шагадифференцирования h. Для этого была написана функция diff2(x_0, h, f), котораявозвращает значение первой производной функции f на основе центральнойформулы численного дифференцирования 2-го порядка в точке x0 для шагадифференцирования h, а так же функция diff4(x_0, h, f), которая возвращаетзначение первой производной функции f на основе центральной формулычисленного дифференцирования 4-го порядка в точке x0 для шагадифференцирования h.Посчитаны значения производной функцииxf ( x )= x⋅eв точке x0 = 2 длямножества значений h ∈ [10 −16 ; 1] сначала с помощью функции diff2, а затем спомощью функции diff4 .
Для обоих случаев построены log-log графикизависимости абсолютной погрешности численного дифференцирования от шагадифференцирования (см. рисунок 1).11Рис. 1: log-log график зависимости абсолютной погрешности численногодифференцирования от шага дифференцирования h для формулы численногодифференцирования 2-го порядка точности (синие точки), для формулычисленного дифференцирования 4-го порядка точности (оранжевые точки),зелёная прямая — функция E = h2, красная прямая — функция E = h4.На рисунке 1 помимо зависимости абсолютной погрешности численногодифференцирования от шага дифференцирования h прямыми линиямивыведены функции вида E = h2 и E = h4.
Как можно заметить, они параллельнысоответствующим зависимостям. Приходим к выводу, что порядок выведенныхформул дифференцирования на log-log графике совпадает с действительнымпорядком.Найдём оптимальный шаг дифференцирования для двух формул, прикотором абсолютная погрешность минимальна.
Рассмотрим формулудифференцирования второго порядка:12f ( x 1 + 2h)− f ( x 1 ) h 2 (3 )f ' ( x1 +h )=− f (ξ )2h6.(10)Из лекций [1] возьмём формулу (3.27) для оптимального шагадифференцирования для данной формулы:hopt =(3⋅ϵ ( 13 ))M,(11)−16где ϵ=2.2⋅10 - машинная точность, M — максимальное значение f (3).x (3)xM =max ( x⋅e ) =max (e (3+ x )) , x = x = 2, тогда M =e 2⋅5 .0Считаем оптимальный шаг:hopt =(3⋅2.2⋅10−6 1/ 3) ≈2.6⋅10−6 ,25⋅e(12)что практически совпадает со значением на графике на рисунке 1.Для формулы дифференцирования четвертого порядка (9) запишемвыражение полной погрешности:~~~~f (x 1−2 h)−8 f (x1 −h)+ 8 f (x 1 +h)− f (x 1 +2 h) ϵ 1 4 (5)f ' (x1 )−= + h f (ξ)12hh 15,(13)e( x 1 −2 h)−8 e( x 1 −h)+8 e ( x 1 + h)−e ( x 1 +2 h)12где вычислительная погрешностьограниченаϵ . Пустьf (5) ограничена M. Тогда верным является следующеенеравенство:~~~~f (x 1−2 h)−8 f (x 1−h)+8 f (x 1+h)− f (x 1+2 h) ϵ 1 4f ' (x 1)−⩽ + h M .
(14)12 hh 15Отсюда,||дифференцируя правую часть, получаем оптимальный шаг:,(15)hopt =(15⋅ϵ 1 /5 15⋅2.2⋅10−16 1 /5) =() ≈ 4.3⋅10−424⋅M4⋅7⋅eчто также практически совпадает со значением на графике на рисунке 1.13Существование данного минимума связано с тем, что при стремлениишага дифференцирования к нулю, решающей становится вычислительнаяпогрешность, связанная с машинным эпсилон.2. Численное интегрированиеРазработаны две программы.В первой программе написана функция composite_simpson(a, b, n, f)численного интегрирования функции f на интервале [a; b] по n узлам спомощью составной формулы Симпсона. С её помощью рассчитан интегралπ∫ g(x )dxдля множества значений n ∈ [3; 9999].
Построен log-log графикзависимости абсолютной погрешности численного интегрирования от шагаинтегрирования (см. рисунок 2).0На рисунке 2 помимо зависимости абсолютной погрешности численногоинтегрирования от шага дифференцирования h прямой линией выведенафункция вида E = h4. Как можно заметить, она параллельна соответствующейзависимости. Делаем вывод, что полученное значение, полученное с помощьюграфика, совпадает с аналитическим порядком точности составной формулыСимпсона.Для данной формулы существует оптимальный шаг интегрирования,минимизирующий достижимую погрешность, так как помимо погрешностиметода, существует вычислительная погрешность, которая появляется привычислении значений функций в узлах.
Будет существовать такой шагинтегрирования, в котором будет минимальна и погрешность метода, ивычислительная погрешность.Согласно формуле (3.66) из лекции [1] верхняя грань вычислительнойпогрешности не зависит от количества отрезков n или их длины h, что говоритнам о том, что увеличение подотрезков не приводит к дестабилизацииабсолютной погрешности вычислений, погрешность не начинает расти, как этобыло на рисунке 1.14Рис. 2: log-log график зависимости абсолютной погрешности численногоинтегрирования от шага интегрирования h для формулы численногоинтегрирования 4-го порядка точности (синие точки), оранжевая прямая —функция E = h4.Во второй программе написана функция gauss_quad5(f) численногоинтегрирования функции f с помощью квадратуры Гаусса пятой степениточности.По теореме (3.2.8) из лекции[1], которая говорит нам о том, чтоквадратуры Гаусса могут быть построены, если в качестве узлов выбраны корнисоответствующего многочлена Лежандра, выведем квадратуру Гаусса пятойстепени точности.
Пятая степень точности квадратуры означает, что она точноинтегрирует полиномы до степени m = 5 включительно. Полином степени mимеет m + 1 оптимизируемых параметров, а квадратура Гаусса имеет их вколичестве 2n, где n — число слагаемых в квадратуре Гаусса. Исходя из этихрассуждений, запишем уравнение, из которого получим число слагаемых вквадратуре Гаусса:15(16)(m+1)=2 nn= 3,⇒значит для такой квадратуры необходимо 3 узла:1∫ f ( x ) dx≈c 1 ( x 1)+ c2 ( x 2 )+ c3 ( x3 )(17).−1Полином Лежандра третьей степени имеет вид:ϕ 3 (x)=x3 −.(18)Он имеет три корня: x 1=0 ,x 2=−√3535x 3=,√35.Вычислим коэффициенты по формуле:11c i =∫ l i ( x) dx=∫Отсюда∏−1 j= 1, i≠ j−1c 1=nx− x j .x i− x j(19)855c =c =9, 2 9, 3 9.Подставив полученные значения узлов и коэффициентов в формулу (17),получим квадратуру Гаусса степени точности 5:1∫ f ( x ) dx≈c 1 f ( x 1 )+ c 2 f ( x 2)+c 3 f ( x 3 )= 89 f (0 )+ 59 f (−−1√3 53)+ f ( ) .5 95√(20)Докажем,что квадратура Гаусса имеет степень точности 5 с помощью следующеговычислительного эксперимента:Построим последовательность полиномов, имеющих степени 0, 1, 2, 3, 4,5 и 6, используя случайно сгенерированные значения коэффициентовполиномов.f 0 ( x)=−333 ,f 1 ( x )=−333+72 x ,f 2 ( x )=− 333 + 72 x + 526 x2,f 3 ( x )=− 333 + 72 x + 526 x 2 − 641 x 323,f 4 ( x )=−333 + 72 x + 526 x − 641 x − 721 x4,16f 5 ( x )=− 333 + 72 x + 526 x 2 − 641 x 3 − 721 x 4 − 732 x 52345,f 6 ( x )=− 333 + 72 x + 526 x − 641 x − 721 x − 732 x + 995 x6.Проинтегрируем их на интервале [0; 2] аналитически и с помощьювыведенной квадратуры Гаусса.Посчитаем абсолютную погрешность.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
