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

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

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

$Описание: y''=f(x,y,y'),\\ y'''=f(x,y,y',y''),\\ \ldots\\ y^{(m)}=f(x,y,y',y''',\ldots,y^{(m-1)}).$

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

$Описание: y^{(m)}=f(x,y,y',y",\ldots,y^{(m-1)})$

( 12.8)

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

Заменим:

$Описание: y_1=y',\\ y_2=y''=y'_1,\\ y_3=y'''=y'_2,\\ \ldots\\ y_m=y^{(m)}=y'_{(m-1)}.$

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

$Описание: \left\{ \begin{array}{l} y'=y_1,\\ y'_1=y_2,\\ y'_2=y_3,\\ \ldots\\ y'_{m-1}=f(x,y,y_1,y_2,\ldots,y_{m-1}). \end{array} \right.$

( 12.9)

Решением системы (12.2), а значит и дифференциального уравнения m -го порядка (12.1) является m табличных функций $Описание: y, y_1=y', y_2=y_1^{''}, \ldots, y_m=y_{(m-1)}$

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

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

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

( 12.10)

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

( 12.11)

Пример

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

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

$Описание: y''=\frac{1}{x}y' - \frac{x^2-1}{x^2}y.$

( 12.12)

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

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

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

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

$Описание: \left\{ \begin{array}{l} y_1=y'\\ y'_1=f(x,y,y_1). \end{array} \right.$

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

$Описание: \left\{ \begin{array}{l} y'=y_1\\ y_1^{'}=\frac{1}{x}y_1-\frac{x^2-1}{x^2}y \end{array} \right.$

( 12.13)

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

где

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

Дана система

$Описание: \left\{ \begin{array}{l} y'=f_1(x,y,y_1),\\ y_1^{'}=f_2(x,y,y_1). \end{array} \right.$

( 12.14)

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

$Описание: y(x_0) = y_0,\\ y_1(x_0)=(y_1)_0.$

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

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

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

и $Описание: (y_1)_i(x_i), i=\overline{1,n}$ , т.е. найти таблицу

i	x	y	y1
0	x0	y0	(y1)0
1	x1	y1	(y1)1
2	x2	y2	(y1)2
3	x3	y3	(y1)3
...	...	...	...
n	xn	yn	(y1)n

Здесь

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

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

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

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

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

При этом на каждом шаге, т.е. для каждого значения xi решением являются две узловые точки с координатами (xi, yi), (xi, (y1)i).

Описание:

Рис. 12.15.

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

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

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

Описание:

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

Описание: Схема алгоритма решения системы (12.6)

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

Здесь

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

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

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

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

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

$Описание: y_{(m)}=f(x,y,y',y'',\ldots,y^{(m-1)})$

( 12.15)

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

$Описание: \left\{ \begin{array}{l} y' = y_1\\ y'_1 = y_2\\ y'_2 = y_3\\ \ldots\\ y'_{(m-1)} = f(x,y_1,y_2,\ldots,y_{(m-1)} \end{array} \right.$

( 12.6)

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

$Описание: y, y_1=y', y_2=y'', \ldots, y_m=y^{(m-1)},$

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

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

Таблица 12.1.
i	x	y	y1=y',	y1=y''1,	:	y_m-1=y^(m-1)
0	x0	y0	(y1)0	(y2)0	:	(ym-1)0
1	x1	y1	(y1)1	(y2)1	:	(ym-1)1
2	x2	y2	(y1)2	(y2)2	:	(ym-1)2
3	x3	y3	(y1)3	(y2)3	:	(ym-1)3
:	:	:	:	:	:	:
n	xn	yn	(y1)n	(y2)n	:	(ym-1)n

где

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

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

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

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

На графике решением дифференциального уравнения m -го порядка (12.8) является совокупность $Описание: (n\times m)$ узловых точек. При этом каждому шагу интегрирования, т.е. каждому значению xi, $Описание: i=\overline{1,n}$ , соответствуют m узловых точек скоординатами

$Описание: (x_i,y_1), (x_i,(y_1)_i), (x_i,(y_2)_i), \ldots, (x_i,(y_{m-1})_i).$

При построении алгоритма задачи будем как и ранее расчет вести по шагам интегрирования, т.е. в цикле по $Описание: i=\overline{1,n}.$ При этом, как и ранее, на каждом i -ом шаге цикла будем рассчитывать решение дифференциального уравнения и тут же его печатать. Тогда нет необходимости формировать матрицу решений, а можно ограничиться формированием массива решений (12.15), который соответствует одной i-й строке матрицы:

$Описание: Y=[y(1), y(2), y(3), \ldots , y(m)],$

( 12.17)

где

$Описание: y(1)= y(x),\\ y(2)= y'(x),\\ y(3)= y"(x),\\ \ldots \ldots\\ y(m)= y^{m-1}(x).$

( 12.18)

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

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

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

$Описание: y''=\frac{1}{x}y' - \frac{x^2-1}{x^2}y$

( 12.19)

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

$Описание: y(1)=y(x),\\ y(2)=y'(x)$

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

$Описание: \left\{ \begin{array}{l} y'(1)=y(2),\\ y'(2)=\frac{1}{x}y(2)-\frac{x^2-1}{x^2}y(1). \end{array} \right.$

( 12.20)

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

$Описание: \left\{ \begin{array}{l} f(1)=y(2),\\ f(2)=\frac{1}{x}y(2)-\frac{x^2-1}{x^2}y(1). \end{array} \right.$

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

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)] - результат, массив значений правых частей уравнений системы (12.10) длиной m.

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

основной программы, в которой организуем циклический процесс по всем шагам интегрирования, $Описание: i=\overline{1,n}$ ;
подпрограммы RGK, в которой на каждом i -м шаге реализуется метод Рунге Кутта (4-го порядка) для системы дифференциальных уравнений m -го порядка. Здесь на каждом шаге интегрирования вычисляется массив решений Y длиной m. Для этого внутри подпрограммы RGK организуем циклы по j, где $Описание: j=\overline{1,m}$ ;
подпрограммы PRAV, обращение к которой осуществляется из подпрограммы RGK для вычисления массива F длиной m - значений правых частей уравнений системы.

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

Описание: Схема алгоритма основной программы

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

Здесь

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

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

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

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

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

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

Тогда

$Описание: T(1,j)=h\cdot f_ j(x, Y),\\ y1(j)=y(j)+T(1,j)/2, j=\overline{1,m},\\ T(2,j)=h\cdot f_ j(x+h/2, Y1),\\ y1(j)=y(j)+T(2,j)/2, j=\overline{1,m},\\ T(3,j)=h\cdot f_ j(x+h/2, Y1),\\ y1(j)=y(j)+T(3,j)/2, j=\overline{1,m},\\ T(4,j)=h\cdot f_ j(x+h, Y1),\\ y(j)=y(j)+\frac{1}{6}(T(1,j)+2T(2,j)+2T(3,j)+T(4,j)).$

Здесь

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

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

T - рабочая матрица порядка ( $Описание: 4\times m$ ).

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

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

Описание: Схема алгоритма подпрограммы RGK

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