Ряды Фурье. Численные методы расчета коэффициентов

Ряды Фурье. Численные методы расчета коэффициентов

МИНИСТЕРСТВО ОБРАЗОВАНИЯ И НАУКИ

РЕСПУБЛИКИ КАЗАХСТАН

Иновационный евразийский университет

КАФЕДРА МАТЕМАТИКИ

Курсовая работа.

Тема: Ряды Фурье. Численные методы расчета коэффициентов.

Курсовая работа студента гр. МТ-31

Нургалиев А. З.

Научный руководитель               

Бокаева М. С.

Дата сдачи курсовой работы _________

Дата защиты _________

                                                                                                             Оценка _________

Павлодар 2006 год.

 Содержание.

1.   Введение.

2.   Разложение периодической функции.

2.1. Периодические величины.

2.2. Гармонический анализ.

3.   Ряды Фурье.

3.1. Методы определения коэффициентов.

3.1.1.   Схема Рунге.

3.1.1.1.   Схема для двенадцати ординат.

3.1.1.2.   Примеры.

3.1.1.3.   Схема для двадцати четырёх ординат.

3.1.2.   Быстрое преобразование Фурье.

3.1.2.1.   Тригонометрическая интерполяция.

3.1.2.2.   Дискретное преобразование Фурье.

3.1.2.3.   Быстрое преобразование Фурье.

3.1.3.   Расчёт коэффициентов на ЭВМ.

3.1.3.1.   Схема Рунге.

3.1.3.2.   Метод трапеций.

3.1.3.3.   Сравнение методов.

4.   Заключение.

5.   Список использованной литературы.

6.   Приложение.

1. Введение.

В курсовой работе рассмотрено различные методы определения коэффициентов рядов Фурье. При разработки данного вопроса было рассмотрено тригонометрическая интерполяция теории и дискретное преобразование рядов Фурье. Также была разработана программа для расчетов коэффициентов на ЭВМ.

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

2. Разложение периодической функции.

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

Таковы например, сила и напряжение переменного тока или – пример паровой машины – путь, скорость и ускорение крейцкопфа, давление пара, касательное усилие в пальце кривошипа и т. д.

            Простейшей из периодических функций (если не считать постоянной) является синусоидальная величина: , где  есть «частота», связанная с периодом Т соотношением:

                                                                                                                               (1)

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

, , , , …,     (2)

которые , если не считать постоянной, имеют частоты

, 2, 3, …,

кратные наименьшей из них, , и периоды

Т, , , …,

то получится периодическая функция (с периодом Т), но уже существенно отличная от величин типа (2).

            Для примера мы воспроизводим здесь сложение трех синусоидальных величин:

;

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

            Теперь естественно поставить обратный вопрос: можно ли данную периодическую функцию  периода Т представить в виде суммы конечного или хотя бы бесконечного множества синусоидальных величин вида (2)?  Как увидим ниже, по отношению к довольно широкому классу функций на этот вопрос можно дать утвердительный ответ, но только если привлечь именно всю бесконечную последовательность величин (2). Для функций этого класса имеет место разложение в «тригонометрический ряд»:

,

причем  суть постоянные, имеющие особые значения для каждой такой функции, а частота  дается формулой (1).

            Геометрически это означает, что график периодической функции получается путем наложения ряда синусоид. Если же истолковать каждую синусоидальную величину механически как представляющую гармоническое колебательное движение, то можно также сказать, что здесь сложное колебание, характеризуемое функцией   или просто ее гармониками (первой, второй и т. д.). Самый же процесс разложения периодической функции на гармоники носит название гармонического анализа.

            Если за независимую переменную выбрать

,

то получится функция от x:

,

тоже периодическая, но со стандартным периодом . Разложение же (3) примет вид

.    (4)

Развернув      члены этого ряда по формуле для синуса суммы и положив

 (n=1,2,3,…),

мы придем к окончательной форме тригонометрического разложения:

,

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

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

3.1.1. Схема Рунге.

Разложение функции в ряд Фурье, или гармонический анализ, оказывается нужным во многих чисто практических вопросах машиноведения, электротехники и пр. Но в этих случаях очень редко приходится непосредственно пользоваться формулами Эйлера-Фурье:

, , ,   (10)

для вычисления коэффициентов разложения. Дело в том, что функции, которые нужно подвергнуть гармоническому анализу, обыкновенно задаются таблицей своих значений или графиком. Таким образом, аналитического выражения функции в нашем распоряжении нет; иногда к самому гармоническому анализу прибегают именно для того, чтобы таким путем получить хотя бы приближенное аналитическое выражение для функции. В этих условиях для вычисления коэффициентов Фурье нужно обратится к приближенным методам. Разумеется, на практике приходится пользоваться лишь немногими первыми членами тригонометрического разложения. Коэффициенты ряда Фурье в большинстве случаев убывают, а с ними быстро падает и влияние далеких гармоник.

            Обычно дается (или снимается с графика) ряд равноотстоящих ординат, т.е. ряд значений функции , отвечающих равноотстоящим значениям аргумента . По этим ординатам величины (10) можно приближенно вычислить, пользуясь методами изложенными выше. Но вычисления здесь оказываются довольно громоздкими, и для того чтобы упростить и, так сказать, автоматизировать их, придумано много различных приемов, один из которых мы и изложим.

3.1.1.1. Схема для двенадцати ординат.

Пусть, скажем, промежуток от 0 до  разделен на k равных частей и пусть известны ординаты

,

отвечающие точкам деления

.

тогда по формуле трапеции имеем (конечно, лишь приближенно!):

.

            Ввиду периодичности нашей функции , значение  можно написать и так:

.

            Аналогично, применяя формулу трапеции к другим интегралам (10), найдем:

или

,

а также

.

            Положим сначала k=12 и будем исходить из двенадцати ординат

,

отвечающих двенадцати ординатам равноотстоящих значениям аргумента:

,

            Все множители, на которые придется умножить эти ординаты, по формулам приведения сведутся к следующим:

.

Именно легко проверить, что

        (11)

Например,

что совпадает с написанным выше выражением.

            Для того чтобы свести выкладки (особенно - умножение) к минимуму, их производят по определенной схеме, предложенной Рунге.

            Сначала выписываются в указанном ниже порядке ординаты и над каждой парой подписанных одна под другой ординат производят сложение и вычитание:

Затем аналогично выписывают эти суммы и разности и снова подвергают их сложению и вычитанию:

  Теперь, получив после всех этих сложений и вычитаний ряд величин , мы можем следующим образом выразить через них искомые коэффициенты:

         (12)

Нетрудно убедится, что эти формулы в точности соответствуют формулам (11).

3.1.1.2. Примеры.

1) Дана некоторая диаграмма касательных усилий (на пальце кривошипа) для некоторой паровой машины. В связи с вопросом о крутильных колебаниях вала представляет интерес выделить гармонические составляющие касательного усилия Т как функции от угла  поворота кривошипа. Сняв с графика двенадцать равноотстоящих ординат, произведем гармонический анализ по указанной схеме:

Теперь по формулам (12):

Таким образом,

Объединим члены, содержащие косинус и синус одного и того же угла:

            Мы видим, что наиболее сильное влияние здесь оказывает вторая гармоника.

2) Для того чтобы дать себе отчет в том, с какой примерно точностью получаются коэффициенты Фурье функции по двенадцати ординатам ее графика, мы приложим изложенный метод к некоторым аналитически заданным функциям и сравним приближенные результаты с точными.

            Сначала рассмотрим функцию , которую в промежутке  задается формулой

,

А для остальных значений x определяется по закону периодичности

.

            Вычислим табличку:

При этом можно использовать легко проверяемое тождество:

По схеме Рунге по этим значениям yнайдем:

b₁=0.608; b₂=0.076; b₃=0.022;

все числа , а с ними и все коэффициенты  оказываются нулями.

            В то же время формулы (10) непосредственно дают (с помощью трехкратного интегрирования по частям):

,

Так что

; ; .

Совпадение превосходное!

3) Однако далеко не всегда получается столь точный результат. В виде второго примера мы возьмем функцию с периодом , которая в промежутке  определяется так:

.

Пользуясь тождеством:

,

составим таблицу:

Тогда по схеме Рунге

числа же  и коэффициенты  - на этот раз нули. Точные значения коэффициентов будут:

в частности,

; ; .

Таким образом, если для первых двух коэффициентов относительная погрешность не превосходит 1,5-2 %, то для последующих она достигнет10% и даже 20%! Ясно, что для повышения этой точности нужно брать больше ординат.

3.1.1.3. Схема для двадцати четырех ординат.

Положим теперь, что даны или сняты с графика двадцать четыре ординаты:

,

отвечающие значениям аргумента:

,

или

.

На этот раз все множители, на которые при приближенном вычислении коэффициентов Фурье приходится умножать ординаты, сведутся к таким:

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

Отметим, что с величинами q и r нет надобности проделывать сложения и вычитания.

Теперь через полученные указанным путем величины k, l, m, n, q и r коэффициентов Фурье выразятся следующим образом:

Дальнейшие коэффициенты по двадцати четырем ординатам получаются с все меньшей точностью.

Нужно обратить внимание на одну подробность. Для получения коэффициентов  и  нужно отдельно вычислить те выражения, которые поставлены в квадратные скобки, а затем сложить их (для нахождения ) и вычесть (для нахождения ). Аналогичное замечание – относительно вычисления коэффициентов  и  .

3.1.2. Быстрое преобразование Фурье.

Тригонометрическая интерполяция. Дискретное преобразование.

            Дискретное преобразование Фурье применяется при решении многих прикладных задач. К ним относится тригонометрическая интерполяция, вычисление сверстки функций, распознавание образов и многие другие. Дискретное преобразование Фурье стало особенно эффективным методом решения прикладных задач после создания быстрого преобразования Фурье.

            Пусть f(x) – периодическая функция с периодом 1 – разложена в ряд Фурье

,                 (1)

причем

.                            (2)

Здесь i – мнимая единица.

            Рассмотрим значение этой функции на сетке из точек , где l, N целые, N фиксировано, и обозначим . Если , где k целое, то , где kl целое. Следовательно,

        (3)

в узлах сетки. Поэтому если функция f(x) рассматривается в узлах сетки , то в соотношении (1) можно привести подобные члены

,     (4)

где

.             (5)

Лемма. При , определяемых (5), соотношение (4) остается в силе, если пределы суммирования [0, N-1] заменить на [m,N-1+m], где m – любое целое.

            Если с самого начала была задана функция, определенная только на сетке, то на этой сетке ее можно также представить в форме (1). Действительно, такую функцию можно продолжить на всю прямую, доопределив ее между узлами сетки путем линейной интерполяции. Для непрерывной кусочно-дифференцируемой функции выполняется (2), поэтому в точках сетки после приведения подобных членов получим (4).

            Определим скалярное произведение для функции на сетке следующим образом:

.

(Множитель  введен для согласованности получаемых соотношений с непрерывным случаем: если f(x) и g(x) – непрерывные функции на отрезке [0,1], то вследствие интегрируемости f(x)g(x) по Риману

при ). Функции  при  образуют ортогональную систему относительно введенного таким образом скалярного произведения. Действительно,

.

При , суммируя геометрическую прогрессию, имеем

(при  знаменатель отличен от 0). Поскольку , то в итоге имеем

 при .             (6)

Умножая (4) скалярно на , получим равенство

             (7)

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

,

поэтому

при  и фиксированном j.

            Покажем, что соотношение

        (8)

в общем случае не имеет места. Пусть . Из (4) получаем , остальные . Таким образом, правая часть (8) есть . Она совпадает с f(x) в точках , но, как правило, далека от нее вне этих точек.

Воспользовавшись утверждением леммы, перепишем (4) в виде

.       (9)

Если f(x) – достаточно гладкая функция, то величины  с ростом j убывают быстро, поэтому  при малых q. Кроме того, при гладкой f(x) величины  и  малы при больших q.

            Напомним, что это приближенное равенство обращается в точное равенство в точках сетки. Способ аппроксимации

Носит название тригонометрической интерполяции. Соотношение (9) называют конечным или дискретным рядом Фурье, а коэффициенты  - дискретными коэффициентами Фурье.

            Игнорирование установленного нами факта о равенстве функций  и  в узлах сетки при  часто являются источником получения неверных соотношений.

            Существует соответствие между задачей приближения функций линейными комбинациями Чебышева и тригонометрическим многочленами. Пусть на отрезке [-1,1] функция f(x) приближается линейными комбинациями . Замена переменных x=cost сводит исходную задачу к задаче приближения функции f(cost) линейной комбинацией .

            Справедливо равенство

.

            Следовательно, задача наилучшего приближения f(x) в норме, соответствующей скалярному произведению , эквивалентна задаче приближения  в норме, соответствующей скалярному произведению . Точно так же существует соответствие в случае задач интерполяции и наилучшего приближения в равномерной метрике. Задача интерполирования функции многочленом по узлам  - нулям многочлена Чебышева - после такой замены сводится к задаче интерполирования функции f(cost) при помощи тригонометрического многочлена  по узлам , образующим равномерную сетку.

3.1.2.3. Быстрое преобразование Фурье.

Осуществление прямого и обратного дискретных преобразований Фурье

Является составной частью решения многих задач решения многих задач. Непосредственное осуществление этих преобразований по формулам (4), (7) требует  арифметических операций. Рассмотрим вопрос о возможности сокращение этого числа. Для определенности речь пойдет о вычислении коэффициентов  по заданным значениям функции. Идея построения алгоритмов быстрого преобразования Фурье опирается то, что при составном N в слагаемых правой части (7) можно выделить группы, которые входят в выражения различных коэффициентов . Вычисляя каждую группу только один раз, можно значительно сократить число операций.

            Рассмотрим сначала случай . Представим q, j, лежащие в пределах , в виде , где . Имеем цепочку соотношений

.

Из равенства

и предыдущего соотношения получим

,

где

.

Непосредственное вычисление всех  требует  арифметических операций, а последующее вычисление  - еще  операций. Поэтому при  общее число операций составит . Точно так же при  строится алгоритм вычисления совокупности значений , для которого общее число операций не превосходит , здесь С – постоянная, не зависящая от N. Выпишем соответствующие расчетные формулы для наиболее употребительного случая . Представим числа q, l в виде

,

где . Величину  представим в виде

,

где s - целое, равное сумме всех слагаемых вида , которых . Очевидно, что , поэтому

После перегруппировки слагаемых имеем

Это соотношение можно записать в виде последовательности рекуррентных соотношений

где

Переход от каждой совокупности к совокупности  требует O(N) арифметических и логических операций; всего таких шагов r, поэтому общее число операций имеет порядок .

Вычисление при помощи совокупностей  дает меньшее накопление вычислительной погрешности по сравнению с формулами (3.7). Определенные удобства имеются также при вычислении экспонент, входящих в расчетные формулы. При вычислении величин  используются значения , . В частности, при m=1 величина  принимает значения +1 или -1. Для вычисления значений  потребуются еще значения  при нечетных j, удовлетворяющих неравенству . Их можно вычислить через уже вычисленные до этого величины, в частности, при помощи соотношений

где, в свою очередь,

при .

В ряде случаев удается еще уменьшить число операций. Один из таких случаев упоминался выше: дана вещественная функция , известная в точках ; требуется найти коэффициенты интерполяционного многочлена

.

Другой случай: при четном N заданы значения функции

в точках ; нужно определить коэффициенты .

3.1.3. Расчет коэффициентов на ЭВМ.

Было запрограммировано два метода расчета коэффициентов на языке Паскаль:

1)   по схеме Рунге;

2)   метод трапеций.

3.1.3.1. Схема Рунге.

Расчет ведется для двенадцати орт. Для большего количества ординат алгоритм остается аналогичным с небольшими корректировками в основной части программы (необходимо заменить вычислительные формулы для коэффициентов). См. приложение 1.

3.1.3.2. Метод трапеций.

Метод трапеций был выбран по причине того, что схема Рунге основана на вычислении коэффициентов Фурье методом трапеций и является лишь результатом удачной манипуляции.. См. приложение 2.

3.1.3.3. Сравнение методов.

Если сравнивать две программы то необходимо заметить, что причиной того что мы отказались во второй программе от непосредственного применения схемы Рунге заключается в том, что она является довольно громоздкой и, несмотря на то, что схема  Рунге требует меньшее количество вычислительных операций (+2n), чем  прямой метод трапеций (), в то же время при вычислении на ЭВМ затрачивается большой объем памяти для хранения промежуточных данных (u,v,p,…). 

Метод Рунге скорее удобен для вычисления вручную, но менее актуален в программировании.

            Если говорить о нахождении более оптимального метода расчета коэффициентов Фурье на ЭВМ, то таким является вышеописанное быстрое преобразование Фурье. Он позволяет сократить количество операций до . В сравнении с вышеописанными методами он является более приемлемым. Способы его алгоритмизации были разработаны и подробно описаны в работе «Numerical recipes in C: The art of scientific computing»-Cambridge unv.,1992.

            Сам алгоритм лишь упоминается в курсовой работе, потому что количество операций Б.П.Ф. сопоставимо со С.Р., только Б.П.Ф. является более гибким (в С.Р необходимо вводить n кратное 12-ти значений функции, а чтобы уменьшить погрешность необходимо вносить изменения в основную программу для увлечения количества исходных данных). 

В дальнейшем я надеюсь продолжить изучение и разработку методов определения коэффициентов Фурье.  

4. Заключение.

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

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

5. Список литературы.

1.   Фихтенгольц Г. М. «Курс дифференциального и интегрального исчисления»(III том) – Москва, 1970г.

2.   БахваловН.С.  «Численые методы» - Москва, 2002г.

3.   Зедгинидзе Г.П., Гогсадзе Р.Ш. «Математические модели в измерительной технике» - Москва, 1970г.

Приложение 1.

{Схема Рунге для 12-ти орт:}

Program MetodRunge;

uses crt;

type ord1=array [0..11] of real;

var Y,U,V,S:ord1;

    A:array [0..3] of real;

    B:array [1..3] of real;

    i,k:integer;

{процедура расчета сумм и разностей значений функции:}

Procedure SummaRaznost(X:ord1;var Sum:ord1;var Raz:ord1;j:integer);

var m:integer;

          begin

               m:=j*2+1;

               for i:=1 to j do

                   begin

                   sum[i]:=X[i]+X[m];

                   raz[i]:=X[i]-X[m];

                   m:=m-1;

                   end;

               sum[j+1]:=X[j+1];

          end;

begin

clrscr;

{Ввод данных:}

writeln('Введите через пробел значения 12-ти, равноотстоящих на pi/6, значений функции');

for i:=0 to 11 do read(Y[i]);

{Расчет значений u и v:}

SummaRaznost(Y,U,V,5);

U[0]:=Y[0];

{Сдвиг всех элементов:u на одну ячейку вправо, чтобы на использовать нулевой элемент матрицы(вообще нулевой элемент будет использоваться только в матрице Y)}

for i:=7 downto 1 do

    U[i]:=U[i-1];

{Расчет значений s и d (значения d заносятся в матрицу Y):}

SummaRaznost(U,S,Y,3);

{Расчет 0-го и 2-го коэффициентов:a, на основе полученных значений s:}

A[0]:=(S[1]+S[2]+S[3]+S[4])/12;

A[2]:=(S[1]-S[4]+0.5*(S[2]-S[3]))/6;

{ Расчет 1-го и 3-го коэффициентов:a, на основе полученных значений d:}

A[1]:=(Y[1]+0.886*Y[2]+0.5*Y[3])/6;

A[3]:=(Y[1]-Y[3])/6;

{Расчет значений sigma и delta(начения sigma заносятся в матрицу Y,   delta в U):}

SummaRaznost(V,Y,U,2);

{Расчет 1 и 3-го коэффициентов:b, ­ на основе полученных значений sigma:}

B[1]:=(0.5*Y[1]+0.886*Y[2]+Y[3])/6;

B[3]:=(Y[1]-Y[3])/6;

{ Расчет 2-го коэффициентов:b, ­ на основе полученных значений delta:}

B[2]:=0.886*(U[1]+U[2])/6;

{Вывод разложения функции в ряд Фурье:}

writeln('Ответ:');

write('T=',A[0]:7:3);

for i:=1 to 3 do begin

                 if A[i]<0 then write(A[i]:7:3)

                    else write('+',A[i]:7:3);

                 write('cos',i,'x');

                 if B[i]<0 then write(B[i]:7:3)

                    else write('+',B[i]:7:3);

                 write('sin',i,'x');

                 end;

end.

Приложение 2.

{Метод трапеций:}

Program MetTrapecyi;

uses crt;

const pi=3.14;

type ord=array [0..5] of real;

var A,B:ord;

    Y:array [0..23] of real;

    h,eps:real;

    m,i,k:integer;

{Функция расчета m-го коэффициента:а}

function af(n:integer;m:integer):real;

var res:real;

    begin

    res:=0;

       for i:=0 to (n-1) do

           res:=res+y[i]*cos(m*i*h);

       af:=res*2/n;

    end;

{Функция расчета m-го коэффициента b :}

function bf(n:integer;m:integer):real;

var res:real;

    begin

    res:=0;

       for i:=0 to (n-1) do

           res:=res+y[i]*(sin(m*i*h));

       bf:=res*2/n;

    end;

begin

clrscr;

writeln('интервал интегрирования: от 0 до 2pi');

{Ввод данных:}

writeln('Введите количество шагов ');

read(k);

writeln(' Введите значения функции с шагом 2pi/',k);

for i:=0 to (k-1) do

    read(Y[i]);

{h-шаг метода}

h:=(2*pi)/k;

            for m:=0 to 5 do begin

                A[m]:=af(k,m);

                B[m]:=bf(k,m);

                end;

{Вывод результата.}

writeln('Отает: ');

writeln('a0=',A[0]/2);

for i:=1 to 5 do

    writeln('a',i,'=',A[i]:5:4);

for i:=1 to 5 do

    writeln('b',i,'=',B[i]:5:4);

end.