Курс лекций для студентов специальности 140104 «Промышленная теплоэнергетика» москва 2010

Вид материала

Содержание

9.7. Метод Рунге - Кутта 4-го порядка
9.8. Решение дифференциальных уравнений высоких порядков
9.9. Решение дифференциальных уравнений второго порядка
9.10. Решение дифференциальных уравнений m-го порядка методом Рунге-Кутта (4-го порядка)

Подобный материал:

1 ... 14 15 16 17 18 19 20 21 22

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

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

В этом методе на каждом шаге интегрирования дифференциальных уравнений искомая функция y(x) аппроксимируется рядом Тейлора (9.4), содержащим члены ряда с h⁴:

В результате ошибка на каждом шаге имеет порядок h⁵.

Для сохранения членов ряда, содержащих h₂,h₃,h₄ необходимо определить вторую y", третью y"' и четвертую y⁽⁴⁾ производные функции y(x). Эти производные аппроксимируем разделенными разностями второго, третьего и четвертого порядков соответственно.

В результате для получения значения функции y_i+1 по методу Рунге-Кутта выполняется следующая последовательность вычислительных операций:

Вывод формулы не приведен. Предоставляется возможность вывод формул выполнить самостоятельно.

Алгоритм метода Рунге-Кутта (4-го порядка) можно построить в виде двух программных модулей: основной программы и подпрограммы Rk4, реализующей метод (рис 9.14).

Рис. 9.14. Схема алгоритма метода Рунге-Кутта 4-го порядка.

Здесь

(x,y)-при вводе начальная точка, далее текущие значения табличной функции,

h-шаг интегрирования дифференциального уравнения,

b-конец интервала интегрирования.

^

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

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

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

Любое дифференциальное уравнение m-го порядка

(9.8)

можно свести к системе, состоящей из m уравнений первого порядка при помощи замен.

Заменим:

В результате дифференциальное уравнение m-го порядка (9.8) сводится к системе, состоящей из m дифференциальных уравнений первого порядка:

(9.9)

Решением системы (9.2), а значит и дифференциального уравнения m-го порядка (9.1) является m табличных функций

.

^

9.9. Решение дифференциальных уравнений второго порядка

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

Общий вид дифференциальных уравнений второго порядка:

(9.10)

Нормальная форма дифференциальных уравнений второго порядка:

(9.11)

Пример

Уравнение в общем виде

Его нормальная форма

(9.12)

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

Заменим y₁=y',

Тогда y'₁=y".

В результате уравнение (9.11) сводится к системе, состоящей из двух дифференциальных уравнений первого порядка:

Для примера (12.12) эта система имеет вид:

(9.13)

Решением этой системы являются две функции y(x) и y₁(x),

где

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

Дана система

(9.14)

Даны два начальных условия:

Необходимо проинтегрировать систему на участке [a, b] с шагом h.

В численных методах задача Коши для системы (9.14) ставится следующим образом:

Найти табличные функции

т.е. найти таблицу

i	x	y	y₁
0	x₀	y₀	(y₁)₀
1	x₁	y₁	(y₁)₁
2	x₂	y₂	(y₁)₂
3	x₃	y₃	(y₁)₃
...	...	...	...
n	x_n	y_n	(y₁)_n

Здесь

h - шаг интегрирования дифференциального уравнения,

a=x₀ - начало участка интегрирования уравнения,

b=x_n - конец участка,

n=(b-a)/h - число шагов интегрирования уравнения.

На графике решением задачи Коши для системы, состоящей из двух дифференциальных уравнений первого порядка, является совокупность узловых точек (рис. 9.15).

При этом на каждом шаге, т.е. для каждого значения x_i решением являются две узловые точки с координатами (x_i, y_i), (x_i, (y₁)_i).

Рис. 9.15.

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

Для вычисления правых частей уравнений системы (9.14) необходимо сформировать подпрограмму PRAV.

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

Схема алгоритма решения системы (9.13) представлена на рис 9.16.

Рис. 9.16. Схема алгоритма решения системы (12.6)

Здесь

h - шаг интегрирования дифференциального уравнения,

b - конец участка,

n - число шагов интегрирования уравнения,

x, y, y₁ - при вводе начальные значения, далее - текущие значения табличной функции.
^

9.10. Решение дифференциальных уравнений m-го порядка методом Рунге-Кутта (4-го порядка)

Как уже было сказано, любое дифференциальное уравнение m-го порядка

(9.15)

сводится к системе, состоящей из m дифференциальных уравнений 1-го порядка

(9.6)

Численным решением системы (9.9), а значит и дифференциального уравнения m-го порядка (9.8) является m табличных функций

т.е. функция y(x) и все ее производные, включая производную (m-1)-го порядка.

При этом каждая из табличных функций определяется на промежутке [a, b] с шагом h и включает n узловых точек. Таким образом, численным решением уравнения (9.8) или системы (9.9) является матрица порядка

, (табл. 9.8)

i	x	y	y₁=y',	y₁=y''₁,	:	y_m-1=y^(m-1)
0	x₀	y₀	(y₁)₀	(y₂)₀	:	(y_m-1)₀
1	x₁	y₁	(y₁)₁	(y₂)₁	:	(y_m-1)₁
2	x₂	y₂	(y₁)₂	(y₂)₂	:	(y_m-1)₂
3	x₃	y₃	(y₁)₃	(y₂)₃	:	(y_m-1)₃
:	:	:	:	:	:	:
n	x_n	y_n	(y₁)_n	(y₂)_n	:	(y_m-1)_n

где

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

n = (b-a)/h - количество шагов интегрирования, равно количеству строк матрицы.

Каждый j-й столбец матрицы - это массив решений одной j-й табличной функции по всем n шагам интегрирования.

Каждая i-ая строка матрицы - это массив решений m табличных функций на одном i-ом шаге интегрирования.

На графике решением дифференциального уравнения m-го порядка (9.8) является совокупность

узловых точек. При этом каждому шагу интегрирования, т.е. каждому значению x_i,

, соответствуют m узловых точек скоординатами

При построении алгоритма задачи будем как и ранее расчет вести по шагам интегрирования, т.е. в цикле по

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

(9.17)

где

(9.18)

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

Рассмотрим пример (9.19).

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

(9.19)

С учетом обозначений (9.18) имеем:

Тогда дифференциальное уравнение (9.19) сводится к системе:

(9.20)

Вычисление правых частей уравнений системы (9.20) будем выполнять в подпрограмме PRAV. При этом подпрограмма PRAV будет иметь вид

Обращение к подпрограмме

PRAV(m, x, Y, F),

где

m-порядок системы,

x- значение x на i-м шаге интегрирования,

Y=[y(1), y(2), y(3), : , y(m)] - входной массив длиной m,

F=[f(1), f(2),:, f(m)] - результат, массив значений правых частей уравнений системы (9.10) длиной m.

Программу решения системы дифференциальных уравнений реализуем в виде 3-х программных модулей:

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

;

- подпрограммы RGK, в которой на каждом i-м шаге реализуется метод Рунге Кутта (4-го порядка) для системы дифференциальных уравнений m-го порядка. Здесь на каждом шаге интегрирования вычисляется массив решений Y длиной m. Для этого внутри подпрограммы RGK организуем циклы по j, где

;

- подпрограммы PRAV, обращение к которой осуществляется из подпрограммы RGK для вычисления массива F длиной m - значений правых частей уравнений системы.

Схема алгоритма основной программы представлена на рис.9.17.

Рис. 9.17. Схема алгоритма основной программы

Здесь

m-порядок системы,

h-шаг интегрирования,

n-количество шагов интегрирования,

x-начальное и далее - текущее значение x,

Y-массив длинной m, куда заносим начальные и далее - текущие значения решений системы на одном шаге интегрирования.

В подпрограмме RGK для вычисления элементов массива Y, используем те же формулы, что и для решения одного дифференциального уравнения 1-го порядка методом Рунге-Кутта (4-го порядка), но с учетом поправки на массивы.

Тогда

Здесь

Y- массив решений длиной m,

Y1 - рабочий массив длиной m,

T - рабочая матрица порядка (

).

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

Схема алгоритма подпрограммы RGK представлена на рис.9.18

Рис. 9.18. Схема алгоритма подпрограммы RGK

Вопросы для самопроверки