Курсовая работа (т): Изучение возможностей массивно-параллельных вычислений в применении к задачам математической физики

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

Сейчас необходимо набрать массив входных данных для преобразования. Как оговаривалось в начале этого пункта, эти данные формально должны быть комплексными, но с нулевой мнимой частью. Для работы с комплексными данными на CUDA существует подходящий тип , определённый как структура из двух элементов, первый из которых несет информацию о действительной части, а второй - о мнимой. Сначала создадим указатель на массив типа  и выделим место:

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

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

Создаём CUFFT план:

Выполняем обратное дискретное преобразование Фурье:

Копируем данные из памяти видеокарты:

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

Абсолютно аналогичные действия приведут нас к нахождению коэффициентов для мембраны. В этом случае предлагается функцию начального положения задавать как произведение двух - для координаты  и координаты , тогда интеграл (3.2.19) распадается на два, и появляется возможность рассчитать коэффициенты для двух осей независимо друг от друга, а потом найти окончательные просто их перемножив. Такой подход был выбран исключительно из соображений экономии времени на преобразовании программы одномерного случая в программу двумерного. Если же по какой-то причине представить функцию  в виде такого произведения не представляется возможным, то программа подлежит доработке, а вместо двух одномерных преобразований Фурье должно быть выполнено одно двумерное.

Вот окончательный набор операций, который приведёт к получению коэффициентов для задачи о мембране:

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

Не стоит забывать, что несмотря на всю эффективность вычисления с помощью FFT, использованная замена интегралов суммами даёт не самые точные значения. Поэтому в следующей главе будет изложен более точный метод вычисления интегралов - с помощью квадратур Лежандра-Гаусса.

5. ВЫЧИСЛЕНИЕ МЕТОДОМ ЛЕЖАНДРА-ГАУССА

5.1 Квадратуры Лежандра-Гаусса

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

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

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

Для удобства будем рассматривать интеграл на отрезке , поиск результатов для других отрезков интегрирования будет осуществляется простым пересчётом. Представляя  последовательно как , и полагая, что  должно равняться нулю, запишем:

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

Представим, что нам известен некоторый полином  степени

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

Так как степень произведения полиномов  меньше, чем  - максимальная степень, для которой применима формула Гаусса, выражение  не несёт в себе противоречий.

Из  следует вывод о том, что полином  ортогонален всякому полиному степени . Но это значит, что с точностью до некоторого множителя полином  совпадает с полиномом Лежандра соответствующей степени. Действительно, полиномы Лежандра, определяемые выражением:

образуют систему, ортогональную на , это их известное свойство /1/. Известно также, что всякая система полиномов, обладающая свойством ортогональности на некотором отрезке, определена однозначно с точностью до множителя. Поэтому можно утверждать, что искомый полином  представим в виде:

а это значит, что необходимые нам узлы  могут быть найдены как корни соответсвующего полинома Лежандра.

Осталось получить веса квадратурной формулы. Для этого воспользуемся системой некоторых полиномов  степени , для каждого из которых выполняется:

Тогда очевидно:

Чтобы получить все  весовых коэффициентов, удобнее всего положить  и представить  в виде:

тогда  будет выполняться. Запишем окончательное выражение для весов:

В связи с тем, что узлы и веса вообще никак не зависят от начального условия - функции произведём их расчёт заблаговременно.

5.2 Расчёт узлов и весовых коэффициентов

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

Прежде всего необходимо подключить пакет ортогональных полиномов, который позволит использовать полиномы Лежандра как уже определённые функции: