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

ПЕНЗЕНСКАЯ ГОСУДАРСТВЕННАЯ ТЕХНОЛОГИЧЕСКАЯ

АКАДЕМИЯ

Кафедра «Прикладная математика и исследование операций в экономике»

«УТВЕРЖДАЮ»

Зав. кафедрой ______________

«_____»______________20__ г.

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

по дисциплине «Численные методы»

на тему « Численное решение систем дифференциальных уравнений »

Вариант № ___

Автор работы:

Специальность: 080116 («Математические методы в

экономике»)

Обозначение курсовой работы: ПГТА 080116-08КР000.17 ПЗ

Группа: 09ММ

Руководитель: Козлов А.Ю.

Работа защищена «__»_____20__г. Оценка _______________

20___ г.

Задание

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

Для решения систем обыкновенных дифференциальных уравнений в системе MATLAB также имеются необходимые «решатели». Это функции ode23, ode45, odell3, odel5s, ode23s, ode23t и ode23tb.

Функции с суффиксом s предназначены для решения так называемых систем жестких дифференциальных уравнений. Для всех остальных систем дифференциальных уравнений наиболее употребительной является функция ode45, реализующая алгоритм Рунге - Кутты 4 - 5-го порядка (разные порядки точности используются для контроля шага интегрирования).

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

yl' = Fl(х, yl, у2,..., уn);

у2' = F2(х, yl, у2,..., уn);

уn' = Fn(х, yl, у2,..., уn);

то есть состоит из n уравнений, разрешенных относительно первых производных функций yl, у2,..., уn (это все функции от х). Введем вектор-столбцы Y и F, состоящие из yl, у2,..., уn и Fl, F2,..., Fn соответственно. Тогда система дифференциальных уравнений примет следующий векторный вид:

Y' = F(х, Y);

Чтобы применить «решатель» ode45, нужно оформить в виде собственной функции правую часть системы уравнений F (х, Y).

Пусть, к примеру, требуется решить следующую систему уравнений:

yl' = у2 + К * х * х;

у2' = -yl;

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

Для данного примера неизвестная вектор-функция Y состоит из двух элементов:

Y = [yl, у2]

Так как функции yl и у2 вычисляются во многих точках в процессе нахождения решения, то реально yl и у2 являются вектор-столбцами. Вектор F правых частей системы уравнений для К = 0.01 вычисляем с помощью собственной функции MyDifEql:

function F = MyDifEql (x, у)

F = [0.01 * x * x + y(2); -y(l)];

текст которой записываем в файл MyDifEql.m Эта функция вызывается каждый раз, когда требуется вычислить правые части уравнений в конкретной точке х, так что здесь х является скаляром, а вектор у состоит всего из двух элементов. Другая функция – MyDifEq0:

function F = MyDifEq0(x, у) F = [у(2); -y(l)];

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

[X, Y] = ode45('MyDifEq0', [0, 20], [0, 1]);

находящую решение нашей системы дифференциальных уравнений с начальными условиями [0, 1] на отрезке [0, 20]. Поскольку здесь использована функция MyDifEq0 для правых частей уравнений, то мы нашли решение для линейного случая гармонических колебаний. Это решение (как yl = Y(:, l), так и у2 = Y (:, 2)) с помощью команд

plot(X, Y(:, l))

hold on

plot(X, Y(:, 2))

отображается на следующем графике (рис. 1).

Рис. 1

Затем вызовем «решатель» ode45 уже с функцией MyDifEql, соответствующей нелинейному случаю, когда К = 0.01:

[X, Y] = ode45('MyDifEql', [0, 20], [0, 1]);

и отобразим на графике соответствующие ему решения (рис. 2).

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

Рис. 2

Теперь перейдем к так называемым жестким системам дифференциальных уравнений. Не вдаваясь в точные определения (оставим их для специальных математических пособий), отметим лишь, что это такие системы дифференциальных уравнений, решения которых на разных отрезках изменения независимых переменных ведут себя абсолютно по-разному: в одних местах наблюдается чрезвычайно быстрое изменение зависимых переменных, в то время как в других местах имеет место их сверхмедленная эволюция. Это слишком сложный характер поведения для обычных алгоритмов численного решения дифференциальных уравнений. Поэтому для надежного решения жестких систем уравнений применяют специальные методы. Как мы упоминали выше, в системе MATLAB для этих целей предусмотрены «решатели» с суффиксом s (stiff - жесткий) в их именах.

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

Уравнение Ван-дер-Поля имеет следующий вид:

yl' = у2;

у2' = -yl + К * (1 - yl *.yl) * у2;

с начальными условиями yl (0) = 2, у2(0) = 0. Здесь К = 1000-большой коэффициент нелинейности задачи. Составим функцию MyVanDerPol, которая описывает правые части представленных дифференциальных уравнений:

function F = MyVanDerPol(х, у)

F = [ у(2); -у(1) + 1000* (1 - у(1)^л2) * у(2) ];

Для решения дифференциальных уравнений Ван-дер-Поля на отрезке изменения независимой переменной [0, 3000] используем «решатель» odel5s:

[X, Y] = odel5s('MyVanDerPol', [0, 3000], [2, 0]);

после чего визуализируем решение, а именно первый столбец матрицы Y:

plot(X, Y(:, l));

Вот как выглядит это решение на рис. 3:

Рис. 3

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

Задание.

1. Составить программу для решения системдифференциальных уравнений методом Рунге - Кутты 4-5-го порядка в Matlab.

2. Решить систему дифференциальных уравнений методом Рунге - Кутты 4-5-го порядка при заданных начальных и конечных условиях (постоянный h =0.1; 0.01) в Matlab (с помощью составленной программы и решателя Matlab).

3. Оценить точность результатов решения.

4. Построить графики , , , а также график f(x, y, z) в декартовой системе координат.

5. Создать видеофайл решения задачи.

Варианты заданий.

Данные взять из таблицы 1.

Таблица 1

№ п/п	СДУ	Начальные условия	Конечные условия
1.		1.0 0.0 0.0	20.0
2.		0.0 1.0 0.0	15.0
3.		0.0 0.0 1.0	25.0
4.		1.0 0.0 0.0	20.0
5.		0.0 1.0 0.0	15.0
6.		0.0 0.0 1.0	25.0
7.		1.0 0.0 0.0	20.0
8.		0.0 0.5 0.0	15.0
9.		0.5 0.0 0.0	20.0
10.		0.0 0.5 0.0	15.0
11.		0.0 0.0 0.5	25.0
12.		0.0 0.0 0.5	15.0
13.		0.5 0.0 0.0	20.0
14.		0.0 0.0 0.5	20.0
15.		0.5 0.0 0.0	15.0
16.		0.5 0.5 0.0	25.0
17.		0.5 0.5 0.5	20.0
18.		0.5 0.0 0.5	20.0
19.		0.0 0.5 0.5	15.0
20.		-0.5 0.0 0.5	25.0
21.		0.5 0.5 -0.5	20.0
22.		-0.5 0.0 -0.5	25.0