X Международная студенческая научная конференция Студенческий научный форум - 2018

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

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

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

Задачи курсовой работы:

изучить теоретические методы решения ОДУ;

разработать параллельный алгоритм метода Рунге-Кутта;

выбор инструментальных средств разработки;

программно реализовать и сделать анализ эффективности метода

Пояснительная записка состоит из введения, трёх глав, заключения, списка использованных источников и приложения.

ТЕОРЕТИЧЕСКАЯ ЧАСТЬ

Разработка и реализация алгоритма решения ОДУ с использованием технологии OpenMP требует тщательного анализа алгоритма и выбора средств разработки. Именно поэтому мы сначала рассмотрим некоторые методы решения ОДУ. Затем рассмотрим стандарт OpenMP для распараллеливания алгоритма. А после более тщательно разберем алгоритм решения ОДУ методом Рунге-Кутта.

1.  
   МЕТОДЫ РЕШЕНИЯ ОДУ

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

Методы их решения подразделяются на три класса:

Аналитические методы.

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

Приближенные методы.

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

Численные методы.

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

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

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

1.  
   СУТЬ МЕТОДА РУНГЕ-КУТТА

Разбор и рассмотрение методов применяемых на практике для решения дифференциальных уравнений мы начнем с их широкой категории известной под общим названием методов Рунге-Кутта

Методы Рунге-Кутта обладают следующими свойствами:

Эти методы являются одноступенчатыми: чтобы найти у_m+1 нужна информация о предыдущей точке (x_my_m)

Они согласуются с рядом Тейлора вплоть до членов порядка hp где степень р различна для различных методов и называется порядковым номером или порядком метода

Они не требуют вычисления производных от f (xy) а требуют вычисления самой функции

Рассмотрим сначала геометрическое построение и выведем некоторые формулы на основе геометрических аналогий. После этого мы подтвердим полученные результаты аналитически

Предположим нам известна точка (x_my_m) на искомой кривой Тогда мы можем провести прямую линию с тангенсом угла наклона у_m=f(x_my_m) которая пройдет через точку (x_my_m) Это построение показано на рисунке 1.2.1 где кривая представляет собой точное но конечно неизвестное решение уравнения а прямая линия L1 построена так как это только что описано

Рис.1.2.1. Геометрическое построение метода Рунге-Кутта

Тогда следующей точкой решения можно считать ту где прямая L₁ пересечет координату проведенную через точку x=x_m+1=x_m+h

Уравнение прямой L₁ выглядит так: y=y_m+y_m(x-x_m), так как y=f(x_my_m) и кроме того x_m+1=x_m+h тогда уравнение примет вид

y_m+1=y_m+h*f(x_my_m) (1)

Ошибка при x=x_m+1 показана в виде отрезка е Очевидно найденное таким образом приближенное значение согласуется с разложением в ряд Тейлора вплоть до членов порядка h так что ошибка ограничения равна

e_t=Кh²

Заметим что хотя точка на рисунке 1.2.1 была показана на кривой в действительности y_m является приближенным значением и не лежит точно на кривой

Формула (1) описывает метод Эйлера один из самых старых и широко известных методов численного интегрирования дифференциальных уравнений. Отметим, что метод Эйлера является одним из методов Рунге-Кутта первого порядка

Рассмотрим исправленный метод Эйлера и модификационный метод Эйлера. В исправленном методе Эйлера мы находим средний тангенс угла наклона касательной для двух точек: (x_my_m) и (x_m+hy_m+hy_m) Последняя точка есть та самая которая в методе Эйлера обозначалась (x_m+1y_m+1) Геометрический процесс нахождения точки (x_m+1y_m+1) можно проследить по рис1.2.2. С помощью метода Эйлера находится точка (x_m+hy_m+hy_m) лежащая на прямой L₁. В этой точке снова вычисляется тангенс дает прямую. Наконец через точку (x_my_m) мы проводим прямую L параллельную L Точка в которой прямая L пересечется с ординатой восстановленной из x=x_m+1=x_m+h и будет искомой точкой (x_m+1y_m+1)

Тангенс угла наклона прямой L и прямой L равен:

Ф(x_my_mh)=½[f(x_my_m)+f(x_m+hy_m+y_mh)] (2)

где y_m=f(x_my_m) (3)

Уравнение линии L при этом записывается в виде:

y=y_m+(x-x_m)Ф(x_my_mh) так что

y_m+1=y_m+hФ(x_my_mh) (4)

Рис.1.2.2.

Соотношения (2) (3) и (4) описывают исправленный метод Эйлера

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

f(xy)=f(x_my_m)+(x-x_m)f/x+(y-y_m)f/x+ (5)

где частные производные вычисляются при x=x_m и y=y_m

Подставляя в формулу (5) x=x_m+h и y=y_m+hy_m и, используя выражение (3) для y_m, получаем

f(x_m+hy_m+hy_m)=f+hf_x+hff_y+O(h²)

где снова функция f и ее производные вычисляются в точке (x_my_m) Подставляя результат в (4) и производя необходимые преобразования, получаем:

Ф(x_my_mh)=f+h/2(f_x+ff_y)+O(h²)

Подставим полученное выражение в (4) и сравним с рядом Тейлора

y_m+1=y_m+hf+h²/2(f_x+ff_y)+O(h³)

Как видим исправленный метод Эйлера согласуется с разложением в ряд Тейлора вплоть до членов степени h² таким образом являясь методом Рунге-Кутты второго порядка

1.  
   ЧИСЛЕННОЕ РЕШЕНИЕ МЕТОДА РУНГЕ-КУТТА

Методы Рунге-Кутта находят широкое применение при решении ДУ. Наибольшее применение нашел метод 4-го порядка.

(3)

_ (4)

_ (5)

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

Общепринятый метод 4-го порядка:

_ (6)

(7)

(8)

(9)

(10)

Ошибка формулы (10) пропорциональна h⁵.

Этот метод намного более точен, чем методы Эйлера, но требует и большего объема вычислений: положение точки (x_i+1, y_i+1) определяется в результате 4-кратного вычисления значения функции f (x,y). С появлением ЭВМ этот недостаток перестал быть существенным и метод Рунге-Кутта 4-го порядка применяется на практике чрезвычайно широко.

Число микроотрезков [x_i; x_i+1], на которые разбивается исходный отрезок [x₀;xn], определяется требуемой точностью вычислений. Для достижения нужной точности задача решается несколько раз при последовательно удваиваемом числе микроотрезков n. Точность считается достигнутой, если при начальном и удвоенном числе n значения y_i и y_2i (в совпадающих точках x) отличаются не более чем на заданную величину:

, i =0, ..,n, (11)

где p – порядок точности метода.

Метод Ругне-Кутта обладает следующими свойствами:

Метод является одноступенчатым (чтобы найти _, нужна информация о предыдущей точке_, _)

Не требует вычисления производных от f(x,y), а требует вычисления самой функции

Имеет небольшую погрешность

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

1.  
   ОБЗОР ТЕХНОЛОГИИ OPENMP

OpenMP (Open Multi-Processing) - это набор директив компилятора, библиотечных процедур и переменных окружения, которые предназначены для программирования многопоточных приложений на многопроцессорных системах с общей памятью (SMP-системах).

Первый стандарт OpenMP был разработан в 1997 г. как API, ориентированный на написание легко переносимых многопоточных приложений. Сначала он был основан на языке Fortran, но позднее включил в себя и языки Си и Си++.

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

Разработку спецификации OpenMP ведут несколько крупных производителей вычислительной техники и программного обеспечения, чья работа регулируется некоммерческой организацией "OpenMP Architecture Review Board" (ARB) [1].

В OpenMP используется модель параллельного выполнения "ветвление-слияние". Программа OpenMP начинается как единственный поток выполнения, называемый начальным потоком. Когда поток встречает параллельную конструкцию, он создает новую группу потоков, состоящую из себя и некоторого числа дополнительных потоков, и становится главным в новой группе. Все члены новой группы (включая главный) выполняют код внутри параллельной конструкции. В конце параллельной конструкции имеется неявный барьер. После параллельной конструкции выполнение пользовательского кода продолжает только главный поток. В параллельный регион могут быть вложены другие параллельные регионы, в которых каждый поток первоначального региона становится основным для своей группы потоков. Вложенные регионы могут в свою очередь включать регионы более глубокого уровня вложенности.

Число потоков в группе, выполняющихся параллельно, можно контролировать несколькими способами. Такими как:

Использование переменной окружения OMP_NUM_THREADS;

Вызов процедуры omp_set_num_threads();

Использование выражения num_threads в сочетании с директивой parallel.

1.  
   OPENMP И ДРУГИЕ ТЕХНОЛОГИИ ПАРАЛЛЕЛЬНОГО ПРОГРАММИРОВАНИЯ

На данный момент считается, что наиболее гибким, переносимым и общепринятым интерфейсом параллельного программирования является MPI (интерфейс передачи сообщений). Однако модель передачи сообщений:

недостаточно эффективна на SMP-системах;

относительно сложна в освоении, так как требует мышления в "невычислительных" терминах.

POSIX-интерфейс для организации нитей (Pthreads) поддерживается широко (практически на всех UNIX-системах), однако по многим причинам не подходит для практического параллельного программирования:

нет поддержки Fortran;

слишком низкий уровень;

нет поддержки параллелизма по данным;

механизм нитей изначально разрабатывался не для целей организации вычислительного параллелизма.

OpenMP можно рассматривать как высокоуровневую надстройку над Pthreads (или аналогичными библиотеками нитей). Перечислим преимущества OpenMP.

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

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

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

Одним из достоинств OpenMP его разработчики считают поддержку так называемых "orphan" (оторванных) директив, то есть директивы синхронизации и распределения работы могут не входить непосредственно в лексический контекст параллельной области.

2. ПРАКТИЧЕСКАЯ ЧАСТЬ

1.  
   РАЗРАБОТКА АЛГОРИТМА

Алгоритм Рунге-Кутты 4-го порядка является итерационным алгоритмом. В нём каждый шаг цикла зависит от предыдущего. В самом шаге, каждая последующая операция также зависит от результата предыдущей операции. Таким образом, упразднить цикл до независимых и соответственно параллельных операций не является возможным. Однако, так как на результат функции f влияет только i-ый параметр вектора y, то для каждого i можно производить независимое вычисление. В итоге, если вычислять результаты выполнения функции f одновременно для всех i-ых параметров вектора y, получим алгоритмическую сложность равную 4nC.

Опишем алгоритм метода Рунге-Кутта по пунктам:

Вводим исходные значения x₀,y₀,xn.

Задаем значение n, например, 10.

Рассчитываем значение h=(xn-x₀)/n.

Последовательно для i=1,2,…,n определяем:

k₁=F(x_i,y_i),k₂=F(x_i+,y_i+),k₃=F(x_i+,y_i+),k₄=F(x_i+h,y_i+k₃h),

y_i+1=y_i+(k₁+2k₂+2k₃+k₄).

Получаем приближенные значения.

1.  
   РАЗРАБОТКА ПРОГРАММЫ

1.  
   1.  
      ВЫБОР ЯЗЫКА ПРОГРАММИРОВАНИЯ

Разработка программы велась в среде Unix-подобной операционной системы Ubuntu 16.04 LTS для процессоров архитектуры Intel i386. Языком программирования для написания программы был выбран C++. Выбор языка обусловлен тем, что библиотеки MPI в стандартной поставке идут только для двух языков — C, С++ и Fortran. При этом C++ является более понятным и знакомым языком программирования.

1.  
   1.  
      НАПИСАНИЕ ПРОГРАММЫ

Рассмотрим код программы, реализующий решение однородных дифференциальных уравнений методом Рунге-Кутты с использованием технологии OpenMP.

Сначала происходит инициализация переменных.

Листинг 1 «Инициализация переменных»

double func( int tip,double x,double y,double a,double b,double c, double d, double e, double f )

{

double s=0;

switch (tip) {

case 1: {s = a+b*(y*c*sin(d*x))-(e*y*f*y); break;}

case 2: {s =a*cos(b*x+c*y)+d*(e*x-f*y); break;}

case 3: {s=((a*cos(b*x)/(x+c))-(d*y*e*y)*f); break;}

case 4: {s = a*(b*x+c*y)/(e*f)*d; break;}

default:{ s =0; }

}

Обозначим ключевые переменные.

tip – представленные для выбора функции.

x, y – переменные функций.

a, b, c, d, e, f – коэффициенты функций.

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

Листинг 2 «Инициализация переменных»

int myrank, size, kol, *tips= new int[100], buffer=1;

int i=0, n=0, k, num;

double h, ot1, ot2, k1=0.0, k2=0.0, k3=0.0, k4=0.0, *a= new double[100], *b=new double[100], *c=new double[100], *d=new double[100], *e=new double[100], *f=new double[100];

Обозначим ключевые переменные.

myrank, size, kol – переменные для расчета размерности на этапе выполнения программы.

*tips= newint[100] – в динамической памяти создается объект tips.

ot1, ot2 – начало и конец отрезка функции.

k1, k2, k3, k4 – значение полученные при вычислении функции.

*a= new double[100], *b=new double[100], *c=new double[100], *d=new double[100], *e=new double[100], *f=new double[100] – в динамической памяти запоминается новое значение объектов a, b, c, d, e, f.

Представим главную директиву библиотеки OpenMP parallel, которая создает параллельную область в представленном ниже блоке. Параметр private директивы parallel – ограничивает доступ к переменным только конкретным потоком. Параметр shared директивы parallel – разрешает доступ к переменным. Количеством потоков можно легко управлять с помощью аргумента num_threads.

Листинг 3 «Директива parallel»

#pragma omp parallel for private(k,i,x,y,j,k1,k2,k3,k4) shared(a,b,c,d,e,f,h,y2) num_threads(num)

for (k=0; k

Код для цитирования: Скопировать