1626435587-55f52a4de97976f3c6215fa7c103f544 (844241), страница 13
Текст из файла (страница 13)
Для комбинации двух переключателей (флагов) используется операция «побитовое ИЛИ» — стандартный приём в языке C. Например, чтобы скомбинировать исчерпывающий способ оптимизации (FFTW_EXHAUSTIVE) с разрешением наперезапись входных значений (FFTW_DESTROY_INPUT), следует передавать в качестве пятого параметра функции fftw_plan_dft_1d значениеFFTW_EXHAUSTIVE | FFTW_DESTROY_INPUT.Заметим, что построение плана может занимать достаточно много времени, особенно при использовании больших массивов данныхсовместно с опциями FFTW_PATIENT или FFTW_EXHAUSTIVE. Например, для таблиц длиной = 220 ≈ 106 создание плана одномерного преобразования Фурье на процессоре Intel(R) Core(TM) i7-47703,40 ГГц занимает более получаса при использовании уровня оптимизации FFTW_EXHAUSTIVE и 7,5 секунд при использовании FFTW_MEASURE.Для экономии времени и повышения быстродействия программы библиотека FFTW позволяет сохранить информацию о найденном оптимальном способе выполнения преобразования Фурье в файл.
Данныймеханизм носит название(мудрость) и подробно описан в разделе «Words of Wisdom — Saving Plans» официальной документацииFFTW [A9].wisdomТаблица 110FFTW_MEASUREFFTW_PATIENTFFTW_EXHAUSTIVE=2План, сек БПФ, мкс0,0280,376,5720=2План, сек БПФ, мс7,53060023186023Необходимость использования механизма wisdom возникает только для сохранения планов между запусками программы. В пределаходного запуска использование «накопленной мудрости» (информациио наиболее эффективном алгоритме выполнения дискретного преобразования Фурье) происходит автоматически, без участия пользователя FFTW.
При необходимости создать в программе план того жетипа и с тем же числом точек , что был создан ранее, повторныйвызов fftw_plan_dft_1d выполняется за пренебрежимо малое время.64(Использование нескольких планов FFTW в одной программе встречается достаточно часто для выполнения прямого и обратного преобразования Фурье либо для преобразования различных величин, сохранённых в разных массивах: поскольку fftw_plan содержит информациюо входном и выходном массиве данных, использование нескольких парвходных и выходных массивов предполагает использование различныхпланов FFTW27 .)Решив проблемы, связанные с установкой библиотеки FFTW, еёподключением на этапе компиляции и компоновки программы, разобравшись с выделением памяти, созданием планов и соответствующими операциями финализации (освобождения ресурсов в конце работы программы), нам осталось приложить последнее усилие собственно для выполнения преобразования Фурье.
Хотя с точки зрения программирования данный шаг тривиален и сводится к вызову функции void fftw_execute(const fftw_plan plan), это место зачастуюстановится источникомошибок в расчётах. Типичная ошибка связана с использованием различных определений дискретного преобразования Фурье, используемых пользователем и разработчиками библиотеки. Действительно, в соответствии с выражением(56) на с. 33, возможны различные определения, отличающиеся знаками в показателе экспонент и значениями коэффициентов + , − .В FFTW преобразование со знаком «минус» в экспоненте считаетсяпрямым (DFT_FORWARD), со знаком «плюс» — обратным (DFT_BACKWARD).Константы + = − = 1, так что суперпозиция прямого и обратногопреобразований не даёт единицу, но эквивалентна (с точностью до погрешности вычислений) домножению таблицы значений на её размер.
Чтобы избежать трудно обнаруживаемых ошибок, перед написанием программы необходимо перепроверить свои обозначения и привестиих в соответствие с использованными в FFTW.В заключение заметим, что возможности библиотеки FFTW значительно превосходят объём данной главы, не претендующей на всеобъемлющее руководство, но призванной лишь познакомить читателя сFFTW и облегчить начало работы с этим замечательным и эффективным инструментом.
Более подробную информацию, включая описаниемногомерных преобразований, различных типов преобразований вещественных данных, использования параллельных вычислений, а такжеответы на целый ряд других вопросов по использованию FFTW можнонайти на официальном сайте библиотеки [A9] в разделе Documentation.математическихОб альтернативном решении см. раздел «New-array Execute Functions» официальной документации FFTW [A9].2765Упражнения1) Выберите параметры временно́й сетки (число узлов, ширина)для моделирования распространения импульсов длительностью10−12 с при требуемом уровне точности не хуже 0,1%. Как следуетизменить параметры, если в процессе распространения ожидаетсяувеличение ширины спектра импульсов в десять раз?2) Постройте спектр мощности гауссовского импульса exp(−2 ) и пары импульсов exp(−2 ) + exp(−( − 5)2 ), объясните вид графиков.3) Постройте спектр мощности гармонического сигнала () = cos 0 при разных 0 , используя прямоугольное окно и окноХанна.4) Какое число членов ряда Фурье позволит аппроксимироватьосцилляции температуры воздуха возле НГУ [A10] за последние3 дня (10 дней) с точностью 5∘ C? 10∘ C? Как убывает амплитудакоэффициентов Фурье с ростом частоты? (Для повышения точности аппроксимации следует использовать сумму линейной функции и конечного ряда Фурье.)5) Используя библиотеку FFTW, выполните дискретное преобразование Фурье от функции () = 1 + + 2 на промежутке0 ≤ < 2 .
Напечатайте и проанализируйте таблицу значенийфункции и её спектра, используя = 8. Выполните обратноепреобразование Фурье, сравните с исходной таблицей значенийфункции.6) Сравните время выполнения 1000 преобразований Фурье с использованием библиотеки FFTW для таблиц длиной = 4093,4095 и 4096. Объсните результат. (Для вычисления времени выполнения программного кода на языке C можно использоватьфункцию clock, объявленную в заголовочном файле <time.h>.)7) Вычисляяконечную разность, найдите производную по време∫︀ни (, ) , где (, ) — численное решение однородногоуравнения теплопроводности − = 0, построенное по явнойсхеме (89).
Сравните график зависимости ( ˜(, ))/|˜ (, )| от саналитическим решением при разных значениях отношения /ℎ2 .8) Решите аналогичную задачу, используя неявную схему Кранка–Николсона (100). Какие спектральные компоненты затухают медленнее всего? Как зависит ответ от величины шага ?663. Уравнение теплопроводностиВ данной главе мы познакомимся с основами применения разностных схем для решения уравнений в частных производных на примерепараболического уравнения второго порядка, описывающего теплопроводность и диффузию.
Начнём с наиболее простого случая, когда имеется всего одна пространственная координата: = 2 + (, ).(86)(Здесь и далее мы будем использовать обозначения и для дифференциальных операторов и соответственно.) Если неоднородность (, ) имеет достаточно простой вид, уравнение (86) может бытьпроинтегрировано аналитически. Однако при включении в (86) нелинейных членов, зависимости коэффициента теплопроводности от координат и времени, а также при поиске решения в областях сложной формы оказываются востребованными численные методы.
Исключительнов целях упрощения изложения мы будем рассматривать применениечисленных методов на наиболее простых примерах, допускающих также точное аналитическое решение. При этом следует иметь в виду, чтотакое упрощение не принципиально, и рассматриваемые методы могутбыть использованы в более сложных случаях, когда получение точногорешения невозможно.Будем искать решение уравнения (86) в области 0 ≤ ≤ 1, > 0.
Начальное распределение температуры (, 0) = 0 (). Граничные условия будут рассмотрены ниже, в п. 3.1. Численное решение уравнения(86) будем строить на равномерной сетке с шагом ℎ = 1/ по пространственной координате и шагом по времени, рис. 11 (а). Условимсяиспользовать (, ) для обозначения точного решения (86), а (, ) —для численного решения.
При этом для краткости будем обозначатьзначения функций в узлах сетки соответствующей буквой с двумя индексами: нижним (пространственным) и верхним (временны́м). Например, ≡ ( , ), ≡ ( , ) и т. п. Пространственный индекс пробегает значения от 0 до , что соответствует 0 ≤ ≤ 1.3.1. Граничные условияУсловия на границе будем задавать в одном из двух видов, известных как задача Дирихле ((0, ) = 0 (), (1, ) = 1 ()) и задача Неймана ( (0, ) = 0 (), (1, ) = 1 (), где ≡ ). Легко увидеть,что заменой неизвестной функции можно свести каждую из этих задач(либо их линейную суперпозицию) к аналогичной задаче с нулевыми67(а)t(б)tm+1tm0(в)xj-1 xjxj+1x(г)(а) Равномерная сетка для построения численного решения; шаблоны численных схем: (б) явной, (в) неявной, (г) Кранка — НиколсонаРис.
11.условиями на границе. Так, в случае задачи Дирихле можно сделатьподстановку:(, ) = (1 − )0 () + 1 () + (, ),где функция удовлетворяет уравнению (, ) = 2 (, ) + (, ) − (1 − )0′ () − 1′ ()с нулевыми граничными условиями (0, ) = (1, ) = 0. Аналогично,для задачи Неймана можно использовать подстановку:2+ 0 () + (, ), (0, ) = (1, ) = 0.2Использование нулевых граничных условий позволяет упростить процедуру численного решения.Чтобы учесть граничные условия при решении задачи Дирихле,очевидно, следует положить 0 = 0, = 0 для всех = 1, 2, .
. . Какприменить условие = 0 при решении задачи Неймана? Наиболеепростой способ состоит в замене производной на уже известную намконечную разность: ( = 0) ≈ (1 − 0 )/ℎ. Однако использованиеданного соотношения приводит к появлению погрешности (ℎ). Какмы увидим дальше, сравнительно легко можно построить разностнуюсхему второго порядка точности по шагу сетки ℎ, поэтому применение формулы первого порядка для производной нежелательно: этоснизит порядок точности численного решения. В этой связи выведемразностную формулу, аппроксимирующую производную на краю сеткисо вторым порядком точности. Для этого запишем невязку междузначением производной ′ (0) и разностной формулой с неопределёнными коэффициентами:(, ) = (1 () − 0 ()) = ′ (0) −0 + 1 + 2.ℎ68Наличие ℎ в знаменателе формулы позволяет сделать коэффициенты, и безразмерными.
Разложив решение () в ряд Тейлора, выберемкоэффициенты так, чтобы занулить невязку в первых трёх порядкахпо шагу сетки ℎ: = −0)︀++ℎ (︀+ ′ (0)(1 − − 2) − ′′ (0) + 4 + (ℎ2 ).ℎ2Разрешая систему линейных уравнений относительно коэффициентов, и , получаем искомую разностную формулу: ′ (0) =−30 + 41 − 2+ (ℎ2 ).2ℎ(87)Действуя аналогичным образом либо записывая полученное выражение (87) для производной (1 − ) при = 0, несложно получитьвыражение для производной на правом краю сетки: ′ ( ) =3 − 4 −1 + −2+ (ℎ2 ).2ℎ(88)Выражения (87) и (88) понадобятся нам для построения численныхсхем второго порядка точности для решения задачи Неймана.3.2.
Явная схемаЗаменяя в уравнении (86) производные и 2 в точке ( , ) конечными разностями, несложно выписать наиболее простую схему длячисленного решения уравнения теплопроводности:+1 − +1− 2 + −1=+ .ℎ2(89)Если значения изестны на -м слое по времени, используя соотношение (89) и граничные условия, можно вычислить значения численного решения +1 на следующем, ( + 1)-м, слое по времени поформулам, ввиду чего (89) относят к классучисленныхсхем. Например, для задачи Дирихле имеем:явнымявных+1= 0,(90))︀ (︀ +1 = + 2 +1− 2 + −1+ , 1 ≤ ≤ − 1.(91)ℎПри решении задачи Неймана вначале нужно применить формулу (91)для вычисления +1 во внутренних точках 1 ≤ ≤ − 1, а затем использовать граничные условия, вычислив значения искомого решения690+1 = 0,+1на границах 0+1 и с использованием полученных выше разностных соотношений (87) и (88).Соотношение (91) связывает одно неизвестное значение сеточнойфункции +1 с тремя известными значениями −1, и +1.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
