Решение задачи Коши для системы ОДУ методом Рунге-Кутты 4 порядка
Последовательный алгоритм
Последовательная сложность	4mn
Объём входных данных	m + 3
Объём выходных данных	(m+1)n
Параллельный алгоритм
Высота ярусно-параллельной формы	O(n)
Ширина ярусно-параллельной формы	O(m)

Основные авторы описания: А.А. Лесничий (1.1,1.2,1.3,1.4,1.5,1.6,1.7,1.8,1.9), Д.А. Алимов (1.1,1.7,1.10,2.4,2.7)

1 Свойства и структура алгоритма

1.1 Общее описание алгоритма

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

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

y(x) = y_{0} + \int\limits_{x_{0}}^{x}f(t,y)dt,

то, применяя различные численные формулы для расчёта интеграла в правой части уравнения, можно получить методы Рунге - Кутты различных порядков.

Общий вид формул методов Рунге - Кутты с шагом сетки h_{n} :

y = f(t,y),\,\,\, y(t_{0}) = y_{0}, y_{n+1} = y_{n} + h_{n+1}\sum\limits_{i=1}^{s}b_{i}K_{n}^{i}, K_{n}^{i} = f(t_{n} + c_{i}h_{n+1}, y_{n} + h_{n+1}\sum\limits_{j=1}^{i-1}a_{ij}K_{n}^{j}).

Конкретный метод Рунге - Кутты определяется набором коэффициентов b_{i}, c_{j}, a_{ij} , которые должны удовлетворять определённым соотношениям.

1.2 Математическое описание алгоритма

1.2.1 Метод Рунге-Кутты 4-го порядка для задачи Коши для ДУ первого порядка

Рассмотрим задачу Коши, где правая часть удовлетворяет условиям теорем существования и единственности решения.

y" = f(x,y),\ a \leq x \leq b;\ y(a) = y^0

Зададим равномерную сетку

x_i = a + ih,\ i = 1,\dots, n,\ h = \frac{b-a}{n}

Введём обозначения y(x_i) = y_i . Получим вычислительную формулу:

\begin{cases} k_1 = hf(x_i,y_i)\\ k_2 = hf(x_i + h/2,y_i + k_1/2)\\ k_3 = hf(x_i + h/2,y_i + k_2/2)\\ k_4 = hf(x_i + h,y_i + k_3)\\ y_{i+1} = y_i + [ k_1 + 2k_2 + 2k_3 + k_4 ]/6 \\ \end{cases}

1.2.2 Метод Рунге-Кутты 4-го порядка для задачи Коши для системы ДУ первого порядка

Численное решение задачи Коши для систем ОДУ 1-го порядка методами Рунге-Кутты ищется по тем же формулам, что и для ОДУ первого порядка. Например, решение методом Рунге-Кутты 4-го порядка можно найти, если положить:

y_i \rightarrow \bar y_i f(x_i,y_i) \rightarrow \bar f(x_i,\bar y_i) k_l \rightarrow \bar k_l \bar k_l = \begin{pmatrix} k^i_{l,1}\\ \vdots\\ k^i_{l,m}\\ \end{pmatrix}

где m – размерность системы, l = 1, \dots, 4 . В результате получим

\begin{cases} \bar k_1 = h\bar f(x_i,\bar y_i)\\ \bar k_2 = h\bar f(x_i + h/2,\bar y_i + \bar k_1/2)\\ \bar k_3 = h\bar f(x_i + h/2,\bar y_i + \bar k_2/2)\\ \bar k_4 = h\bar f(x_i + h,\bar y_i + \bar k_3)\\ \bar y_{i+1} = \bar y_i + [ \bar k_1 + 2\bar k_2 + 2\bar k_3 + \bar k_4 ]/6 \\ \end{cases}

1.3 Вычислительное ядро алгоритма

1.4 Макроструктура алгоритма

Макроструктуру алгоритма составляют вычисления коэффициентов k_{j},\,j=1,..,4 при нахождении значений искомой функции в каждом узле.

1.5 Схема реализации последовательного алгоритма

Приведем возможную реализацию последовательного алгоритма на языке Matlab (правая часть ОДУ задана в виде функции func(x,y))

1 function = RungeKutta4 (y0, a, b, n) 2 h = (b - a ) / n ; 3 x = a : h : b 4 5 y (:, 1 ) = y0 ; 6 for i = 2 : n 7 k1 = h * func (x (i ), y (:, i )) 8 k2 = h * func (x (i ) + h / 2 , y (:, i ) + k1 / 2 ) 9 k3 = h * func (x (i ) + h / 2 , y (:, i ) + k2 / 2 ) 10 k4 = h * func (x (i ) + h , y (:, i ) + k3 ) 11 y (:, i + 1 ) = y (:, i ) + (k1 + 2 * k2 + 2 * k3 + k4 ) / 6 12 end 13 end

1.6 Последовательная сложность алгоритма

Каждый шаг цикла алгоритма состоит из 4 обращений к функции \bar f , 11 умножений и 10 сложений. Так как обращение к функции \bar f является наиболее сложной операцией, то сложность линейного выполнения алгоритма можно определить как 4mn .

1.7 Информационный граф

Опишем граф алгоритма аналитически и графически.
Пусть задана система из m одномерных дифференциальных уравнений. Граф имеет трёхмерную структуру, которая представляет собой связанные между собой ограниченные плоскости (уровни). На каждом i-ом уровне вычисляется i-ое приближение вектора решения системы. Всего таких уровней – n, где n – количество узлов, в которых вычисляется приближённое решение системы дифференциальных уравнений. Выходные данные i-го уровня являются входными данными для i+1-го уровня.
Вершины информационного графа делятся на два типа:

• Первому типу вершин соответствует вычисление одной координаты функции, стоящей в правой части дифференциального уравнения, в точках, подаваемых на вход вершине. Результатом выполнения операции является одна координата вектора k_{l},\,\,\,l = 1,2,3,4 . На каждом уровне таких вершин будет 4m , а их общее количество – 4mn .
• Второму типу вершин соответствует операция a + bc . Эти вершины также можно разделить на две группы. В вершинах первой группы вычисляются выражения, результат которых будет подаваться в качестве аргументов для вершин первого типа, при этом входными данными для этих вершин будут результаты срабатывания m вершин первого типа, расположенных на одном уровне. На каждом уровне вершин первой группы будет 3m , а их общее количество – 3mn . Входными данными для вершин второй группы будет результат срабатывания одной вершины первого типа и одной такой же вершины второго типа. Результатом срабатывания данной вершины на i-ом уровне будут m координат вектора i-го приближения решения системы. На каждом уровне таких вершин будет 4m , а их общее количество – 4mn .

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

Рисунок 1. Информационный граф для системы из двух дифференциальных уравнений.

1.8 Ресурс параллелизма алгоритма

Поскольку в описанной выше вычислительной схеме наиболее трудоемкой является операция расчета правых частей ОДУ при вычислении k_i (i = 1, \dots, 4) , то основное внимание следует уделить распараллеливанию этой операции. Здесь будет применяться подход декомпозиции уравнений системы ОДУ на подсистемы. Поэтому для инициализации рассмотрим следующую схему декомпозиции данных по имеющимся процессорным элементам с локальной памятью: на каждый \mu - ПЭ (процессорный элемент) (\mu = 0, \dots, p-1 ) распределяется n/p дифференциальных уравнений и вектор \bar y_0 . Далее расчеты производятся по следующей схеме:

1. на каждом ПЭ одновременно вычисляются n/p соответствующих компонент вектора \bar k_1 по формуле [ \bar k_1 ]_{\mu} = h[ \bar f(x_i, \bar y_i) ]_{\mu}
2. для обеспечения второго расчетного этапа необходимо провести сборку вектора \bar k_1 целиком на каждом ПЭ. Затем независимо выполняется вычисление компонент вектора \bar k_2 по формуле [ \bar k_2 ]_{\mu} = h[ \bar f(x_i + h/2,\bar y_i + \bar k_1/2)]_{\mu} ;
3. проводится сборка вектора \bar k_2 на каждом ПЭ, вычисляются компоненты вектора \bar k_3:\ [ \bar k_3 ]_{\mu} = h [\bar f(x_i + h/2,\bar y_i + \bar k_2/2)]_{\mu} ;
4. проводится сборка вектора \bar k_3 на каждом ПЭ, вычисляются компоненты вектора \bar k_4:\ [ \bar k_4 ]_{\mu} = h [\bar f(x_i + h,\bar y_i + \bar k_3)]_{\mu} ;
5. рассчитываются с идеальным параллелизмом компоненты вектора \bar y_{i+1}:\ [\bar y_{i+1}]_{\mu} = [\bar y_{i}]_{\mu} + ([ \bar k_1 ]_{\mu} + 2[ \bar k_2 ]_{\mu} + 2[ \bar k_3 ]_{\mu} + [ \bar k_4 ]_{\mu})/6\ и производится сборка вектора \bar y_{i+1} на каждом ПЭ. Если необходимо продолжить вычислительный процесс, то полагается i = i + 1 и осуществляется переход на п. 1

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

На каждом ярусе в данном алгоритме каждый ПЭ производит четыре операции вычисления правых частей ОДУ, шестнадцать операций сложения векторов и умножения вектора на число, а так же р-1 пересылку данных между другими ПЭ, что довольно сильно замедляет выполнение алгоритма и является накладными расходами при распараллеливании алгоритма. При этом количество итераций цикла равно длине вектора x , то есть n . Из вышеуказанного следует, что при классификации по высоте ЯПФ алгоритм имеет сложность O(n) , а при классификации по ширине ЯПФ – O(m) (при условии, что p = m ).

1.9 Входные и выходные данные алгоритма

Входными данными алгоритма являются:

1. вектор y^0 размерности m ;
2. границы временного интервала a и b ;
3. частота дискретизации n ;

Всего размер входных данных: m + 3

Выходными данными являются

1. n векторов \bar y_i размерности m ;
2. вектор \bar x размерности n

Всего размер выходных данных: (m+1)n

1.10 Свойства алгоритма

Перечислим основные свойства алгоритма:

2 Программная реализация алгоритма

2.1 Особенности реализации последовательного алгоритма

2.2 Локальность данных и вычислений

2.2.1 Локальность реализации алгоритма

2.2.1.1 Структура обращений в память и качественная оценка локальности

2.2.1.2 Количественная оценка локальности

2.2.1.3 Анализ на основе теста Apex-Map

2.3 Возможные способы и особенности параллельной реализации алгоритма

2.4 Масштабируемость алгоритма и его реализации

2.4.1 Масштабируемость алгоритма

2.4.2 Масштабируемость реализации алгоритма

Проведём исследование масштабируемости параллельной реализации метода Рунге - Кутты согласно методике . Исследование проводилось на суперкомпьютере "Ломоносов" Суперкомпьютерного комплекса Московского университета . Исследовалась собственная параллельная реализация алгоритма, написанная на языке C++ с использованием стандарта MPI. Сборка программы производилась при помощи компилятора компании Intel версии 15.0 (без указания ключей компиляции) с использованием библиотеки OpenMPI версии 1.8.4 . Расчеты проводились на узлах из группы T-Blade2 с процессорами Intel Xeon 5570 Nehalem 2.93GHz, при этом для каждого MPI - процесса выделялось одно ядро.

Набор и границы значений изменяемых параметров запуска реализации алгоритма:

• число процессоров с шагом по степеням двойки до 64, далее с шагом 10;
• размерность системы .

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

• минимальная эффективность реализации 0.0000032%;
• максимальная эффективность реализации 0.0024%.

Перечислим некоторые особенности тестируемой параллельной реализации:

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

Рисунок 3. Эффективность работы параллельного алгоритма в зависимости от размерности системы и числа процессоров.

Построим оценки масштабируемости протестированной параллельной реализации метода Рунге - Кутты 4-го порядка:

• По числу процессов -0.000000768. При увеличении числа процессов эффективность на рассмотренной области изменений параметров запуска уменьшается, однако уменьшение достаточно медленное, что связано с крайне низкой максимальной эффективностью работы параллельной реализации алгоритма. Уменьшение эффективности при увеличении числа процессов для систем с небольшой размерностью объясняется тем, что при количестве процессов превышающем размерность системы часть из них не будет участвовать в вычислениях. Для систем с большой размерностью увеличение процессов негативно сказывается на эффективности, так как увеличивается количество пересылок между ними, что существенно затормаживает работу.
• По размеру задачи 0.000000763. При увеличении размерности системы при фиксированном количестве процессов эффективность на рассмотренной области изменений параметров запуска увеличивается, причём увеличение наблюдается на почти всей рассматриваемой области. Это объясняется тем, что при фиксированном количестве процессов увеличение размерности системы приводит к увеличению вычислительной нагрузки на каждый процессор, вследствие чего увеличивается время вычислительной работы относительно времени простоя. С увеличением числа процессов скорость возрастания эффективности при увеличении размерности задачи уменьшается, так как количество пересылок становится достаточно большим.
• По двум направлениям 0.00000000174. При одновременном увеличении количества процессов и размерности системы эффективность увеличивается. Скорость возрастания эффективности крайне низкая, что объясняется высокими накладными расходами.

2.5 Динамические характеристики и эффективность реализации алгоритма

2.6 Выводы для классов архитектур

2.7 Существующие реализации алгоритма

Метод 4-го порядка является самым часто используемым из всех схем Рунге - Кутты, поэтому существует множество его программных последовательных реализаций как коммерческих, так и бесплатных. Одной из самых известных программных реализаций является ode45 в среде MATLAB, у которой имеется большое количество возможностей для настройки численного расчёта, что делает её достаточно удобной для использования. Также метод Рунге - Кутты 4-го порядка реализован в параллельной библиотеке PETSc, которая является свободно распространяемой. Несмотря на то, что в этой библиотеки большинство методов реализовано с использованием параллельных алгоритмов, метод Рунге - Кутта в ней реализован последовательным. Довольно трудно найти параллельную реализацию метода в случае рассмотрения общей задачи Коши, но существует множество научных работ, в которых описываются параллельные реализации метода Рунге - Кутты для конкретных классов систем, например, в работе А.В. Старченко описывается параллельная реализация для линейных систем.

Пусть дано дифференциальное уравнение первого порядка

с начальным условием

y(x 0) = y 0 .

Выберем шаг h и введём обозначения:

x i = x 0 + i . h и y i = y(x i) ,

где i = 0, 1, 2, …

x i – узлы сетки,

y i - значение интегральной функции в узлах .

Аналогично описанным выше методам производится решение дифференциального уравнения. Отличие состоит в делении шага на 4 части.

Согласно методу Рунге – Кутта четвёртого порядка, последовательные значения y i искомой функции y определяются по формуле:

, i = 0, 1, 2, …

а числа k 1 (i) , k 2 (i) , k 3 (i) , k 4 (i) на каждом шаге вычисляются по формулам:

Это явный четырёхэтапный метод четвёртого порядка точности.

Методы Рунге – Кутта легко программируются и обладают значительной точностью и устойчивостью для широкого круга задач.

На рисунке 6 приведена блок-схема процедуры RUNGE(X0, XK, Y0, N, Y) для решения задачи Коши описанным выше методом Рунге – Кутта.

i = 0, … , N-1

K1 = h * F(x, Yi)

K2 = h * F(x + h/2, Yi + K1 / 2)

K3 = h * F(x + h/2, Yi + K2 / 2)

K4 = h * F(x + h, Yi + K3)

K = (K1 + 2*K2 + 2*K3 + K4) / 6

Рисунок 6 - Блок-схема процедуры RUNGE

На рисунке 7 приведена блок-схема алгоритма основной программы для решения задачи Коши и получения результатов с фиксированным количеством отрезков разбиения N. В основной программе происходит обращение к процедуре RUNGE(X0, XK, Y0, N, Y), вычисляющей значения искомой функции y j в точках x j методом Рунге – Кутта.

Исходными данными в данной задаче являются:

X0, XK – начальное и конечное значения независимой переменной;

Y0 – значение y 0 из начального условия y(x 0) = y 0 ;

N – количество отрезков разбиения.

Результаты работы программы выводятся в виде двух столбцов:

X – массив значений узлов сетки;

Y – массив значений искомого решения в соответствующих узлах сетки.

Ввод X0, XK, Y0, N

RUNGE(X0,XK,Y0,N,Y)

i=0…N

Вывод X, Y i

Рисунок 7 - Блок-схема алгоритма основной программы для решения задачи Коши с фиксированным количеством отрезков разбиения N

Решение дифференциальных уравнений в среде MathCad

Рисунок 8 - Пример решения дифференциального уравнения методом

Рунге-Кутта 4 порядка в среде MathCad.

Построение графиков функций в среде Visual Basic

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

Задача.

Построить график функции y=sin(x) на отрезке . Шаг табулирования принять равным h.

Решение.

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

Рисунок 9 - Расположение основных компонентов в окне General

Компонент Picture Box () используется в качестве контейнера для построения графика. Он представляет собой матрицу из точек (пикселей), причём имеется возможность управлять цветом каждой отдельной точки. Координаты любой точки определяются парой целых чисел – ее порядковым номером в строке Х и порядковым номером строки внутри объекта Y. Таким образом, координаты левого верхнего угла компонента (0, 0). Число точек в строке и число строк определяются размером компонента.

Рисунок 10 - Координаты объекта PictureBox

На рис. 10 показано расположение осей и координаты угловых точек объекта.

Компонент Line () используется для построения осей и отрезков ломаной графика функции.

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

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

Рассмотрим алгоритм более подробно.

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

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

Ввод исходных данных осуществляется в рассматриваемой программе при нажатии на командную кнопку. Очень часто ввод данных реализуется с помощью компонента TextBox.

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

Количество строк рассчитывается по формуле k=n+2, где k – количество строк, а n – количество отрезков табулирования. Число строк должно быть больше количества отрезков на 2, так как необходимо учесть начальную точку (нулевую) и строку для записи заголовков столбцов страницы.

В самой процедуре табулирования можно совместить два действия – табулирование и расчет экстремумов. Такой вариант решения приведен в листинге программы на рисунке13.

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

Коэффициенты масштабирования графика рассчитываются по следующим формулам:

где kx - коэффициент масштабирования по оси ОХ,

NPX – количество пикселей объекта PictureBox, отводимых для построения графика по горизонтали,

a – начальное значение отрезка аргумента функции,

b – конечное значение отрезка аргумента функции.

,

где Ky - коэффициент масштабирования по оси ОY,

NPY – количество пикселей объекта PictureBox, отводимых для построения графика по вертикали,

min – минимальное значение функции,

max – максимальное значение функции.

Перевод математических координат текущей точки в экранные производится по формулам:

zx = Round(ox + (x(i) - a) * kx),

zy = Round(oy - (y(i) - Min) * ky),

где zx, zy – экранные координаты текущей точки,

ox, oy - координаты точки пересечения осей в компоненте pictureBox,

x(i), y(i) – математические координаты текущей точки,

kx, ky – коэффициенты масштабирования.

В формуле расчета экранной координаты ординаты текущей точки используется знак «минус» для учета противоположного направления осей (на экране и на графике).

Листинг программы построения графика функции приведен на рисунке 13.

Примеры форм с результатами работы программы для различных исходных данных приведены на рисунках 14 и 15.

Рисунок 11 - Пример разработки формы

Рисунок 12 - Алгоритм построения графика функции

Rem Описание переменных

Dim x() As Single, y() As Single

Private a As Single

Классический метод Рунге - Кутта 4 порядка

Тема 6.3. Метод Рунге-Кутгы

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

Формально, методами Рунге - Кутта являются модифицированный и исправленный метод Эйлера, они представляют собой схемы второго порядка точности. Существуют стандартные схемы третьего порядка, не получившие широкого распространения. Наиболее часто используется и реализована в различных математических пакетах (Maple, MathCAD, Maxima) стандартная схема четвёртого порядка (см. ниже). Иногда при выполнении расчётов с повышенной точностью применяются схемы пятого и шестого порядков.

Метод Рунге - Кутта 4 порядка столь широко распространен, что его часто называют просто метод Рунге - Кутта.

Рассмотрим задачу Коши . Тогда значение в следующей точке вычисляется по формуле:

Величина шага сетки по .

Этот метод имеет 4 порядок, т.е. ошибка на каждом шаге составляет O(h5), а суммарная ошибка на конечном интервале интегрирования O(h4).

Семейство прямых методов Рунге - Кутта является обобщением метода Рунге - Кутта 4 порядка. Оно задается формулами

Конкретный метод определяется числом s и коэффициентами bi,aij и ci. Эти коэффициенты часто упорядочивают в таблицу

Отправить свою хорошую работу в базу знаний просто. Используйте форму, расположенную ниже

хорошую работу на сайт">

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

Размещено на http://www.allbest.ru/

ФЕДЕРАЛЬНОЕ Государственное АВТОНОМНОЕ образовательное учреждение Высшего профессионального образования

«БЕЛГОРОДСКИЙ ГОСУДАРСТВЕННЫЙ НАЦИОНАЛЬНЫЙ

ИССЛЕДОВАТЕЛЬСКИЙ УНИВЕРСИТЕТ»

Институт инженерных технологий и естественных наук

Кафедра математического и программного обеспечения информационных систем

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

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

по дисциплине «Методы вычислений»

студента очного формы обучения

направления подготовки 010500.62

«Математическое обеспечение и администрирование

информационных систем»

3 курса группы 07001302

Данькова Николая Алексеевича

БЕЛГОРОД 2016

Введение

1. Теоретическая часть

1.1 Обыкновенное дифференциальное уравнение первого порядка. Задача Коши

1.2 Суть метода Рунге-Кутта

1.3 Выбор среды разработки

2. Практическая часть

2.1 Программная реализация метода Рунге-Кутта 4-го порядка

3. Тестирование

3.1 Пример

Заключение

Приложение

Введение

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

= f(x, y) ,

где у - искомая функция, х - независимая переменная, f(x,y) - непрерывная функция от х и у. Однако получить аналитическое решение этого уравнения для достаточно произвольной функции f не удается, и только для некоторых частных случаев, с которыми можно ознакомиться в справочной литературе.

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

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

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

Курсовая работа состоит из 3 разделов, содержит 6 рисунков, 3 листинга, 1 приложение и 18 страниц.

1. Теоретическая часть

1.1 Обыкновенное дифференциальное уравнение первого порядка. Задача Коши

Для простоты рассмотрим двумерное пространство переменных х и у и некоторое открытое множество G, принадлежащее ему. Пусть на этом открытом множестве определена непрерывно дифференцируемая функция f(х, у) и задано уравнение

= f(x, y) (1)

Согласно теореме существования и единственности для любой точки (x 0 ,y 0) ?G найдется решение у = у(х), определенное на некотором интервале (х 0 -д, х 0 +д), удовлетворяющее условию y(x 0) = y 0, такое, что точки (x,y(x)) ?G и y` x ? f(x, y(x)), причем это решение будет единственным. Задача для уравнения (1) с начальным условием у(х 0) = y 0 (задача Коши) состоит в нахождении функции у(х), обращающей и уравнение (1), и начальное условие в тождество. Допустим, что значения, которые принимает независимое переменное х, принадлежат интервалу (Х 0 , X N) и запишем задачу Коши:

(2)

Разобьём отрезок [Х 0 , X N ] на N частей так, что x n +1 - х n = h n ,

n = 0, … ,N-1. В дальнейшем, не ограничивая общности, рассмотрим случай, когда разбиение равномерное, т.е. все h n = h = = const,

n = 0 ,… ,N-1.

1.2 Суть метода Рунге-Кутта

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

(3)

(4)

(5)

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

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

(6)

(7)

(8)

(9)

(10)

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

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

Число микроотрезков , на которые разбивается исходный отрезок , определяется требуемой точностью вычислений. Для достижения нужной точности задача решается несколько раз при последовательно удваиваемом числе микроотрезков n. Точность считается достигнутой, если при начальном и удвоенном числе n значения y i и y 2i (в совпадающих точках x) отличаются не более чем на заданную величину:

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

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

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

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

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

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

1.3 Выбор среды разработки

C++ Builder-- программный продукт, инструмент быстрой разработки приложений (RAD), интегрированная среда программирования (IDE), система, используемая программистами для разработки программного обеспечения на языке программирования C++. Данный продукт позволяет создавать как консольные приложения, так и приложения с графическим интерфейсом.

Microsoft Visual Studio -- линейка продуктов компании Microsoft, включающих интегрированную среду разработки программного обеспечения и ряд других инструментальных средств. С помощью данного продукта можно разрабатывать консольные приложения, приложения с графическим интерфейсом, а также веб-сайты, веб-приложения, веб-службы как в родном, так и в управляемом кодах для всех платформ, поддерживаемых Windows, Windows Mobile, Windows CE, .NET Framework, Xbox, Windows Phone .NET Compact Framework и Silverlight.

Для выполнения поставленной задачи был выбран программный продукт C++ Builder. Так как является более простым в использовании и соответствует всем необходимым требованиям для создания консольного приложения.

2. Практическая часть

2.1 Программная реализация метода Ркнге-Кутта 4-го порядка

дифференциальный уравнение программирование коши

Разработка программы начинается с описания функций. Для этого мы используем оператор switch.

Листинг 1 «описание функций»

double s=0;

switch (tip){

case 1: {

break; }

case 2: {

break; }

case 3: {

break; }

case 4: {

s = a*(b*x+c*y)/(e*f)*d;

break; }

default:

{ s =0; }

}

return s;

}

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

Когда все нужные данные получены, мы переходим непосредственно к решению ОДУ методом Рунге - Кутта 4-го порядка.

Листинг 2 «программная реализация решения ОДУ методом Рунге - Кутта

4-го порядка»

for (i=0;i<=n;i++) {

x=x[i]+h;

cout<<"y"<

}

3. Тестирование

Запустив программу, мы увидим уравнения, предлагаемые для выбора.

Рис. 1 «выбор уравнения»

Если ввести номер не соответствующий представленным номерам уравнений программа отреагирует на это.

Рис. 2 «ввод неверного параметра»

Выбрав необходимое уравнение, вводим коэффициенты.

Рис. 3 «ввод коэффициентов»

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

Рис. 4 «ввод необходимых параметров»

После окончания вычисления программа выводит решение.

Рис. 5 «вывод результата»

3.1 Пример

Решить задачу Коши:

на отрезке . Найти решение на равномерной сетке с шагом 0.1

Решение. Так как f(x ,y) = х + у, то получаем

= + ,

= ++ ,

= ++ ,

= +h+ ,

= +() ,

= +h ,

для значений i = 1, 2, 3, 4.

Полагая =0, = 1, последовательно находим:

при i = 1

= 0,1(0 +1) = 0,1 ,

= 0,1(0 + 0,05 +1 + 0,05) = 0,11 ,

= 0,1(0 + 0,05 +1 + 0,055) = 0,1105 ,

= 0,1(0 + 0,1 +1 + 0,1105) = 0,121050 ,

= 1 + *(0,1 + 2*0,11+2*0,1105 + 0,12105) = 1,110342 ,

= 0+0,1=0.1 ,

при i = 2

= 0,1*(0 + 1,110342) = 0,121034 ,

= 0,1*(0,1 + 0,05 +1,110342 + 0,0605171) = 0,1320859 ,

= 0,1*(0,1 + 0,05+1,110342 + 0,06604295) = 0,1326385 ,

= 0,1*(0,1 + 0,1 + 1,110342 + 0,11326385)= 0,1442980 ,

= 1,110342 + *(0,121034 + 2*0,1320859+2*0,1326385 + 0,1442980) = 1,242805 ,

= 0,1+0,1 ,

Далее получаем:

при i = 3 = 0,3, =1,399717,

при i = 4 = 0.4, = 1,583648.

Погрешность полученного решения не превышает величины

|y 4 - ц(x 4) | ? 0.000001.

Для наглядности в таблице 1 приведены численные решения одной и той же задачи Коши методами Эйлера, Эйлера-Коши и Рунге-Кутта.

Таблица 1 Численные решения задачи коши разными методами

		Значения, найденного методом	Точное решение ц(x i)=2- x i -1
			Эйлера - Коши	Рунге - Кутта

Теперь сравним полученные результаты с расчетами нашей программы.

Рис. 6 «результат работы программы»

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

Заключение

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

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

Список использованных источников

1. Березин И.С., Жидков Н.П., Методы вычислений: Т.2 - М.: ГИФМЛ, 1960. - 620 с.

2. Бахвалов Н. С., Жидков Н. П., Кобельков Г. М. Численные методы. - М.: Бином, 2001 - с. 363-375.

3. Копченова Н.В., Марон И.А., Вычислительная математика в примерах и задачах - М.: Наука, 1972. - 368 с.

4. https://ru.wikipedia.org/wiki/Microsoft_Visual_Studio

5. https://ru.wikipedia.org/wiki/C%2B%2B_Builder

Приложение 1

#include

#include

#include

#include

using namespace std;

char* rus(const char* text) {

char *buffRus=new char;

CharToOem(text, buffRus);

return buffRus;

}

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; }

}

return s;

}

int main() {

int tip, i=0, n=0;

double h=0.0, ot1=1.0, ot2=0.0, k1=0.0, k2=0.0, k3=0.0, k4=0.0, a=1.0, b=1.0, c=1.0, d=1.0, e=1.0, f=1.0, res=0.0;

bool flag=0;

cout<

cout<<" 1. y"=a+b*(y*c*sin(d*x))-(e*y*f*y)\n 2. y"=a*cos(b*x+c*y)+d*(e*x-f*y)\n 3. y"=((a*cos(b*x)/(x+c))-(d*y*e*y)*f)\n 4. y"=a*(b*x+c*y)/(e*f)*d\n";

while (!flag) {

cout<

cin>>tip;

if((tip == 1) || (tip == 2) || (tip == 3) || (tip == 4)){

flag=1; }

else {

cout <

}

}

cout<

cout <<" a= ";cin>>a;

cout <<" b= ";cin>>b;

cout <<" c= ";cin>>c;

cout <<" d= ";cin>>d;

cout <<" e= ";cin>>e;

cout <<" f= ";cin>>f;

cout<

while (ot1>ot2) {

cout<

return 0;

}

Приложение 2

Размещено на Allbest.ru

Подобные документы

Анализ предметной области объектно-ориентированного программирования. Языки Delphi, Object Pascal - объектно-ориентированная среда программирования. Основные алгоритмические решения. Решение дифференциального уравнения методом Рунге-Кутта в среде Excel.

курсовая работа , добавлен 02.04.2011


Составление программы на алгоритмическом языке Turbo Pascal. Разработка блок-схемы алгоритма её решения. Составление исходной Pascal-программы и реализация вычислений по составленной программе. Применение методов Рунге-Кутта и Рунге-Кутта-Мерсона.

курсовая работа , добавлен 17.09.2009


Суть метода Рунге-Кутта и его свойства. Решение дифференциальных уравнений первого порядка. Вычислительный блок Given/Odesolve. Встроенные функции rkfixed, Rkadapt, Bulstoer. Решения линейных алгебраических уравнений в среде MathCad и Microsoft Excel.

курсовая работа , добавлен 02.06.2014


Обзор методов решения в Excel. Рекурентные формулы метода Эйлера. Метод Рунге-Кутта четвертого порядка для решения уравнения первого порядка. Метод Эйлера с шагом h/2. Решение дифференциальных уравнений с помощью Mathcad. Модифицированный метод Эйлера.

курсовая работа , добавлен 18.01.2011


Математическое описание задачи решения обыкновенного дифференциального уравнения численным явным методом Рунге-Кутта, разработка схемы алгоритма и написание программы в среде программирования Microsoft Visual Studio 2010. Тестирование работы программы.

курсовая работа , добавлен 22.01.2014


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

курсовая работа , добавлен 05.11.2011


Решение дифференциальных уравнений первого порядка. Варианты методов Рунге-Кутта различных порядков. Основные методы численного решения задачи Коши. Повышение точности вычислений и итерационный метод уточнения. Дискретная числовая последовательность.

лабораторная работа , добавлен 14.05.2012


Численные методы решения задачи Коши для обыкновенных дифференциальных уравнений: Эйлера, Рунге-Кутта, Адамса и Рунге. Техники приближенного решения данных уравнений: метод конечных разностей, разностной прогонки, коллокаций; анализ результатов.

курсовая работа , добавлен 14.01.2014


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

курсовая работа , добавлен 11.07.2012


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

и информатики

Уральский технический институт связи и информатики

Кафедра физики, прикладной математики и информатики.

КУРСОВАЯ РАБОТА

по информатике:

Визуализация численных методов.

Решение обыкновенных дифференциальных уравнений.

Выполнил: Арапов

Гр. МЕ-71

Проверил: Минина Е. Е.

Екатеринбург

2008 г.

Введение………………………………………………………………….3

1. Постановка задачи…………………………………………………….4

2. Описание методов решения…………………………………………..5

2. 1. Суть задачи………………………………………………………….5

2. 2. Геометрический смысл задачи…………………………………….5

2. 3. Численные методы решения задачи Коши……………………….6

2. 4. Метод Эйлера ……………………..…………………………….6

2. 5. Метод Рунге-Кутта 4-го порядка………………………………….8

2. 6. Решение поставленной задачи методами Эйлера и

Рунге-Кутта 4-го порядка …………………………………….9

2. 6. 1. Метод Эйлера ………………….. ……….……………………9

2. 6. 2. Метод Рунге-Кутта 4-го порядка ……………………………10

3. Алгоритм решения задачи…………………………………………...11

3. 1. Алгоритмы подпрограмм.………………………………………....11

3. 1. 1. Подпрограмма метода Эйлера……………………..………..11

3. 1. 2. Подпрограмма метода Рунге-Кутта 4-го порядка ………..12

3. 1. 3. Подпрограмма общего решения ………. . …………………1 3

3. 2. Алгоритм функции…………………………………………………1 3

3. 3. Алгоритм программы………………………………………………14

4. Форма программы…………………………………………………….17

5. Листинг программы…………………………………………………..18

6. Решение задачи в MathCad …………………………………………..20

Заключение………………………………………………………………22

Введение

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

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

Для достижения цели я поставил перед собой следующие задачи:

1. Написать программу для решения данного дифференциального уравнения двумя численными методами в программе Visual Basic .
2. Проверить решение с помощью приложения MathCad .
3. Сравнить полученные разными методами результаты с общим решением.

1. Постановка задачи

Решить методами Эйлера и Эйлера модифицированного задачу Коши для дифференциального уравнения 1-го порядка на отрезке [ X 0 ; X k ] с шагом h и начальным условием: Y (X 0 ) = Y 0 .

Ответ должен быть получен в виде таблицы результатов:

	Y (1)	Y (2)
	Y 0 (1)	Y 0 (2)	Y(X 0 )
	Y 1 (1)	Y 1 (2)	Y(X 1 )
	Y k (1)	Y k (2)	Y(X k )

Где Y (1) , Y (2) – решения, полученные различными численными методами, Y T – точное решение дифференциального уравнения.

Возможно представление результатов решения не в виде таблицы, а в виде списков.

Данные таблицы визуализировать на форме в виде графиков.

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

Дифференциальное уравнение					Общее решение
y’ = - x·y/(x+1)	1 ,2		0, 1		y=c· (x+1) ·exp(-x)

2. Описание методов решения

2. 1. Суть задачи

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

Задачу Коши можно сформулировать следующим образом:

Пусть дано дифференциальное уравнение и начальное условие y (x 0 ) = у 0 . Требуется найти функцию у(x ), удовлетворяющую как указанному уравнению, так и начальному условию.

Численное решение задачи Коши сводится к табулированию искомой функции.

График решения дифференциального уравнения называется интегральной кривой.

2. 2. Геометрический смысл задачи

y ’ = f (x , y ) - тангенс угла наклона касательной к графику решения в точке (х, у) к оси 0Х, - угловой коэффициент (рис. 1).

Рисунок 1. Геометрический смысл задачи Коши.

Существование решения:

Если правая часть f (x , y ) непрерывна в некоторой области R , определяемой неравенствами

| x - x 0 | < а; | y - y 0 | < b ,

то существует, по меньшей мере, одно решение у = у(х), определённое в окрестности |х – х 0 | < h , где h - положительное число.

Это решение единственно, если в R выполнено условие Липшица

| f (x , y )- f (x , y )| ≤ N | y - y |(x , y ),

где N - некоторая постоянная (константа Липшица), зависящая, в общем случае, от а и b . Если f (x , у) имеет ограниченную производную

f ’ y (x , y ) в R , то можно положить N = мах | f ’ y (х, у)| при (х, y ) принадлежащим R .

2. 3. Численные методы решения задачи Коши

При использовании численных методов выполняется замена отрезка [х 0 , X ] - области непрерывного изменения аргумента х множеством. состоящего из конечного числа точек х 0 < х 1 < ... < x n = Х - сеткой.

При этом x i называют узлами сетки.

Во многих методах используются равномерные сетки с шагом:

Задача Коши, определённая ранее на непрерывном отрезке [х 0 , X ], заменяется её дискретным аналогом - системой уравнений, решая которую можно последовательно найти значения y 1 , y 2 ,…, y n - приближённые значения функции в узлах сетки.

Численное решение задачи Коши широко применяется в различных областях науки и техники, и число разработанных для него методов достаточно велико. Эти методы могут быть разделены на следующие группы.

Одношаговые методы, в которых для нахождения следующей точки на
кривой у = f (x ) требуется информация лишь об одном предыдущем шаге.
Одношаговыми являются метод Эйлера и методы Рунге - Кутта.

Методы прогноза и коррекции (многошаговые), в которых для отыскания следующей точки кривой у = f (x ) требуется информация более чем об одной из предыдущих точек. Чтобы получить достаточно точное численное значение, часто прибегают к итерации. К числу таких методов относятся методы Милна, Адамса - Башфорта и Хемминга.

Явные методы, в которых функция Ф не зависит от y n +1 .

Неявные методы, в которых функция Ф зависит от y n +1 .

2.4 Метод Эйлера.

Иногда этот метод называют методом Рунге-Кутта первого порядка точности.

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

Пусть дано дифференциальное уравнение первого порядка:

Y’ = f(x, y)

с начальным условием

y (x 0 ) = y 0

Выберем шаг h и введём обозначения:

x i = х 0 + ih и y i = y (x i ), где i = 0, 1, 2, ...,

x i - узлы сетки,

y i - значение интегральной функции в узлах.

Иллюстрации к решению приведены на рисунке 2.

Проведем прямую АВ через точку (x i , y i ) под углом α. При этом tg α = f (x i , y i )

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

Тогда y i +1 = y i + Δy

Приравняем правые части tg α = f (x i , y i ) и. Получим

Отсюда Δ у = h ∙ f (x i , y i ).

Подставим в это выражение формулу y i +1 = y i + Δy , а затем преобразуем его. В результате получаем формулу расчета очередной точки интегральной функции:

Рисунок 2. Метод Эйлера.

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

Блок-схема процедуры решения дифференциального уравнения методом Эйлера приведена на рисунке 3.

F (x , у) - заданная функция – должна

быть описана отдельно.

Входные параметры:

Х0, XK —начальное и конечное

значения независимой переменной;

Y 0 – значение y 0 из начального условия

y (x 0 ) = y 0 ;

Выходные параметры:

У - массив значений искомого решения

в узлах сетки.

Рисунок 3. Блок-схема процедуры решения дифференциального уравнения методом Эйлера.

Метод Эйлера - один из простейших методов численного решения обыкновенных дифференциальных уравнений. Но существенным его недостатком является большая погрешность вычислений. На рисунке 2 погрешность вычислений для i -г o шага обозначена ε. С каждым шагом погрешность вычислений увеличивается.

2.5 Метод Рунге-Кутта 4 порядка

Пусть дано дифференциальное уравнение первого порядка с начальным условием y (x 0 )= y 0. Выберем шаг h и введем обозначения:

x i = x 0 + ih и y i = y (x i ), где i = 0, 1, 2, ... .

Аналогично описанному выше методу производится решение

дифференциального уравнения. Отличие состоит в делении шага на 4 части.

Согласно методу Рунге-Кутта четвертого порядка, последовательные значения y i искомой функции y определяются по формуле:

y i+1 = y i +∆y i где i = 0, 1, 2 ...

∆y=(k1+2*k2+2*k3+k4)/6

a числа k 1, k 2 , k 3, k 4 на каждом шаге вычисляются по формулам:

k 1 = h * f (x i , y i )

k2 = f (x i +h/2, y i +k1 /2)*h

k3 = F(x i +h/2, y i +k2 /2)*h

k4 = F(x i +h, y i +k3) * h

Это явный четырехэтапный метод 4 порядка точности.

Блок-схема процедуры решения дифференциального уравнения методом Рунге-Кутта приведена на рисунке 6.

F (x , у) - заданная функция - должна

быть описана отдельно.

Входные параметры:
Х0, X К - начальное и конечное

значения независимой

переменной;

Y 0 – значение y 0 из начального условия

y (x 0 )= y 0 ;

N - количество отрезков разбиения;

Выходные параметры:

Y - массив значений искомого решения

в узлах сетки.

2. 6. Решение поставленной задачи методами Эйлера и Рунге-Кутта 4 порядка

2. 6. 1. Метод Эйлера

1. Строим оси координат;

2. Отмечаем A (1,2; 1) – первую точку интегральной кривой;

4. Строим касательную l 0 в точке А под углом α 0 ;

5. Находим х 1 по формуле: x i = х 0 + ih , где h – шаг интегрирования

x 1 =1,2 + 1 · 0,1 = 1,3

6. Проводим прямую x = x 1 = 0,1 до пересечения с прямой l 0 , отмечаем точку B (x 1 ; y 1 );

7. Ищем y точки B :

Из прямоугольного треугольника ABC ,

Δ y = y 1 – y 0 ,

Δx = x 1 – x 0 = h ,

f (x 0 ; y 0 ) = (y 1 – y 0 )/ h =>

y 1 = y 0 + h · (f (x 0 ; y 0 )) = 1 + 0,1 · f (1,2;1) = 1-0,545454 = 0,945

Следовательно, точка B имеет координаты (1,3; 0,945).

2. 6 .2. Метод Рунге-Кутта 4 порядка

1. Строим оси координат;

2. Отмечаем А(1,2; 1) – первую точку интегральной кривой;

3. Ищем угол наклона касательной к графику в точке A :

4. Строим касательную l 0 в точке А под углом α 0 ;

5. Находим х 1 по формуле: x i = х 0 + ih

x 1 = 1 ,2 + 1 · 0, 1 = 1, 3 ;

1. Находим по формулам:

k1=0,1·f(1,2;1)=0,1·(-0,545454)= - 0,0545

k2=0,1· f(1,2+0,1/2;1-0,0545/2)= -0,0540

K3=0,1· f(1,2+0,1/2;1-0,0540/2)= - 0,0540

K4=0,1· f(1,2+0,1;1-0,0540)= - 0,0535

∆ y 1 =(- 0,0545+2·(-0,0540)+2·(-0,0540) - 0,0535)/6= - 0,054

∆ y 2 =1- 0,054=0,946

Следовательно, следующая точка графика решения имеет координаты (1,3; 0,946)

3. Алгоритм решения задачи

3. 1. Алгоритмы подпрограмм

3.1.1 Подпрограмма метода Эйлера

3.1.2 Подпрограмма метода Рунге-Кутта 4 порядка

3. 1. 3. Подпрограмма общего решения

3. 2. Алгоритм функции

3. 3. Алгоритм программы

4. Форма программы

1. Листинг программы

Dim j() As Single

Dim x() As Single

Dim y() As Single

Dim o() As Single

Private n, i As Integer

Private xk, x0, kx, ky As Single

Private k, k1, k2, k3, k4 As Single

Private h, max, min, y0 As Single

Private Function f(a, b As Single) As Single

f = -a * b / (a + 1)

End Function

Private Function f1(x As Single) As Single

End Function

Private Sub Eiler()

ReDim x(n)

ReDim j(n)

j(0) = y0

For i = 0 To n

x(i) = x0 + h * i

Next i

For i = 0 To n - 1

Next i

End Sub

Private Sub Runge()

ReDim y(n)

y(0) = y0

For i = 0 To n

x(i) = x0 + h * i

Next i

For i = 0 To n - 1

k1 = h * f(x(i), y(i))

y(i + 1) = y(i) + k

MSFlexGrid1.TextMatrix(1, 2) = Str(y0)

Next i

End Sub

Private Sub Obhee()

ReDim o(n)

For i = 0 To n

o(0) = y0

x(i) = x0 + h * i

o(i) = f1(x(i))

Next i

End Sub

Private Sub Command1_Click()

x0 = Val(Text1.Text)

xk = Val(Text2.Text)

h = Val(Text4.Text)

y0 = Val(Text3.Text)

n = (xk - x0) / h

Label6.Caption = Str(x0)

Label5.Caption = Str(xk)

MSFlexGrid1.Rows = n + 2

MSFlexGrid1.TextMatrix(0, 1) = "Ýéëåð"

MSFlexGrid1.TextMatrix(0, 2) = "Ðóíãå-Êóòò"

MSFlexGrid1.TextMatrix(0, 3) = "Îáùåå ðåøåíèå"

Eiler

Runge

Obhee

max = y0

min = y0

For i = 0 To n

If j(i) > max Then

max = j(i)

End If

If j(i) < min Then

min = j(i)

End If

If y(i) > max Then

max = y(i)

End If

If y(i) < min Then

min = y(i)

End If

If o(i) > max Then

max = o(i)

End If

If o(i) < min Then

min = o(i)

End If

Next i

Label4.Caption = Str(max)

Label7.Caption = Str(min)

Picture1.Cls

For i = 1 To n - 1

X1 = 720 + Round(kx * (x(i - 1) - x0))

X2 = 720 + Round(kx * (x(i) - x0))

X1 = 720 + Round(kx * (x(i - 1) - x0))

X2 = 720 + Round(kx * (x(i) - x0))

Next i

End Sub

Private Sub Command2_Click()

End Sub

1. Решение задачи в MathCad

Эйлер

Рунге-Кутт

Общее

Заключение

Из двух методов (Эйлера и Рунге-Кутта) по полученным результатам точнее (сравненивая с общим решением) оказался метод Рунге-Кутта. Это объясняется тем что, ведь в отличие от метода Эйлера в методе Рунге-Кутта шаг делится не на 4 отрезка, в результате чего погрешность метода становится меньше.

По завершению курсовой работы я выполнил все поставленные задачи: написал программу для решения данного дифференциального уравнения двумя численными методами в программе Visual Basic , проверил решение с помощью приложения MathCad ,сравнил полученные разными методами результаты с общим решением. Я считаю что полностью достигнул поставленную цель данной курсовой работы.

tg(α) = f(x,y)

Y i+1 = Y i + h ∙ F(x, Y i )

x = X0 + i ∙ h

i = 0, …, N - 1

h = (Xk – X0)/N

Конец

MSFlexGrid1.TextMatrix(i + 2, 2) = Str(y(i + 1))

MSFlexGrid1.TextMatrix(1, 2) = Str(y (0))

kx = (6600 - 720) / (xk - x0)

ky = (6600 - 1120) / (max - min)

Label4.Caption = Str(max)

Label7.Caption = Str(min)

o(i)=min

o(i)

MSFlexGrid1.TextMatrix(i + 1, 3) = Str(o(i))

y(i)=max

Eiler(X0, Xk, Y0, N, Y)

x i +1

х i

y i+1

y=y(x)

k2=h*F(x+h/2, Y i +k1/2)

K1=h*F(x,Y i )

y(i) > max

j(i)=min

x = X0 + i ∙ h

i = 0, …, N-1

h = (Xk – X0)/N

Rynge4(X0, Xk, Y0, N, Y)

Command2

o(i)=max

o(i) > max

k=(k1+2*k2+2*k3+k4)/6

k4= h*F(x+h, Y i +k3)

k3= h*F(x+h/2, Y i +k2/2)

Y i+1 = Y i +k

k1 = h * f(x(i), y(i))

k2 = h * f(x(i) + h / 2, y(i) + k1 / 2)

k3 = h * f(x(i) + h / 2, y(i) + k2 / 2)

k4 = h * f(x(i) + h, y(i) + k3)

k = (k1 + 2 * k2 + 2 * k3 + k4) / 6

y(i + 1) = y(i) + k

i = 0, …, n-1

x(i) = x0 + h * i

i = 0, …, n

ReDim y(n)

g(0) = y0

Начало

y0, x0,xk,h

n = (xk - x0) / h

max = y0

min = y0

j(i) > max

i = 0, … n

Eiler

Runge

Obshee

MSFlexGrid1.Rows = n + 2

MSFlexGrid1.TextMatrix(0, 0) = "x"

MSFlexGrid1.TextMatrix(0, 1) = " Эйлер "

MSFlexGrid1.TextMatrix(0, 2) = " Рунге - Кутта "

MSFlexGrid1.TextMatrix(0, 3) = " Общее решение "

Label6.Caption = Str(x0)

Label5.Caption = Str(xk)

j(i)=max

X1 = 720 + Round(kx * (x(i - 1) - x0))

X2 = 720 + Round(kx * (x(i) - x0))

Y1 = 6600 - Round(ky * (o(i - 1) - min))

Y2 = 6600 - Round(ky * (o(i) - min))

Конец

j(i)

X1 = 720 + Round(kx * (x(i - 1) - x0))

X2 = 720 + Round(kx * (x(i) - x0))

Y1 = 6600 - Round(ky * (y(i - 1) - min))

Y2 = 6600 - Round(ky * (y(i) - min))

X1 = 720 + Round(kx * (x(i - 1) - x0))

X2 = 720 + Round(kx * (x(i) - x0))

Y1 = 6600 - Round(ky * (j(i - 1) - min))

Y2 = 6600 - Round(ky * (j(i) - min))

i = 1 , …, n-1

Конец

f1 = y0 / ((x0 + 1) * Exp(-x0)) * (x + 1) * Exp(-x)

f1(x)

MSFlexGrid1.TextMatrix(1, 0) = Str(x0)

MSFlexGrid1.TextMatrix(i + 2, 0) = Str(x(i + 1))

MSFlexGrid1.TextMatrix(i + 2, 1) = Str(j(i + 1))

MSFlexGrid1.TextMatrix(1, 1) = Str(j(0))

j(i + 1) = j(i) + h * f(x(i), j(i))

i = 0, …, n-1

x(i) = x0 + h * i

i = 0, …, n

ReDim x(n)

ReDim j(n)

j(0) = y0

Eiler

Конец

x(i) = x0 + h * i

o (i) = f1 (x(i))

i = 0, …, n

ReDim o (n)

o(0) = y0

Obhee

Конец

MSFlexGrid16

Picture1

f = -a * b / (a + 1)

f(a,b)

Runge

Labe71

Text2

Text1

Labe41

Labe31

Label1

Text3

Labe21

Label6

Text4

Command1

Label4

Labe51

Labe91

Label10

Конец

Picture1.Line (X1, Y1)-(X2, Y2), RGB(0, 200, 0)

Picture1.Line (X1, Y1)-(X2, Y2), RGB(500, 70, 90)

Picture1.Line (X1, Y1)-(X2, Y2), RGB(400, 100, 12)

y(i)

y(i)=min