Донецкий национальный технический университет
Опубликован: 15.03.2007 | Доступ: свободный | Студентов: 5956 / 2116 | Оценка: 4.11 / 3.78 | Длительность: 12:32:00
Специальности: Математик
Лекция 12:

Компьютерное моделирование и решение нелинейных уравнений

< Лекция 11 || Лекция 12: 12345678910

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

Как уже было сказано, любое дифференциальное уравнение 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-х программных модулей:

  1. основной программы, в которой организуем циклический процесс по всем шагам интегрирования, i=\overline{1,n} ;
  2. подпрограммы RGK, в которой на каждом i -м шаге реализуется метод Рунге Кутта (4-го порядка) для системы дифференциальных уравнений m -го порядка. Здесь на каждом шаге интегрирования вычисляется массив решений Y длиной m. Для этого внутри подпрограммы RGK организуем циклы по j, где j=\overline{1,m} ;
  3. подпрограммы 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
< Лекция 11 || Лекция 12: 12345678910
Равиль Султанов
Равиль Султанов

В уравнениях движения кривошипно-шатунного механизма вместо обозначения радиуса кривошипа "r" ошибочно записан символ "γ" (гамма).

P.S. Может быть это слишком очевидно, но не упомянуто, что угол поворота кривошипа φ считается малым.

Александр Никитин
Александр Никитин

Добрый день.

В расчете параметра Т4 xi суммируется с величиной h/2 ?