В цикле 19 - 22 полученные приращения прибавляются к начальной точке и результат записывается в выходной массив.

В блоке 23 вычисляются производные в конечной точке шага.

4.3 Блок - схема алгоритма одного шага по методу Рунге - Кутта

Рисунок 2

4.4 Подпрограмма одного шага по методу Рунге-Кутта.

SUBROUTINE SH(TN,H,XN,XK,F,PRAV,N,R)

DIMENSION XN(N),XK(N),F(N),P(4,2),R(N,3)

DATA P/0.,.5,.5,1.,6.,3.,3.,6./

DO 1 K=1,N

R(K,2)=0.

1 R(K,3)=0.

DO 4 J=1,4 ! Начало внешнего цикла определения 4-х приращений

T=TN+P(J,1)*H ! Задание независимой переменной.

DO 2 K=1,N ! Цикл задания массива значений зависимых переменных.

2 R(K,1)=XN(K)+P(J,1)*R(K,2)

CALL PRAV(T,R,F,N) ! Вычисление производных в заданной точке.

DO 3 K=1,N ! Цикл вычисления и накопления частичных приращений.

R(K,2)=H*F(K)

3 R(K,3)=R(K,3)+R(K,2)/P(J,2)

4 CONTINUE

DO 5 K=1,N

5 XK(K)=XN(K)+R(K,3) ! Вычисление переменных в конечной точке.

CALL PRAV(TN+H,XK,F,N) ! Вычисление производных в конечной точке.

RETURN

END

4.5 Описание алгоритма метода Рунге - Кутта с автоматическим выбором шага

Блок -схема алгоритма приведена на Рисунке 3.

В алгоритме используются следующие идентификаторы

Таблица 5

Второй столбец массива Х должен содержать весовые коэффициенты погрешности, на которые умножаются найденные погрешности каждой интегральной переменной, чтобы после сложения этих произведений получить общий критерий погрешности системы и сравнить его с заданной погрешностью. Во втором столбце они задаются с целью удобства ввода. Первый столбец массива Х заполняется начальными условиями, а затем , подряд, вводятся весовые коэффициенты. Алгоритм начинается с переписывания весовых коэффициентов в четвертый столбец массива F(блоки 2,3). Номера столбцов обозначены нижним индексом, а номера строк - верхним. После задания начальных значений параметров (4) вызывается подпрограмма вычисления производных (5) в начальной точке интервала интегрирования и начальная точка с производными в ней передается подпрограмме вывода (40). Затем начинается основной цикл выполнения шагов интегрирования. Задается шаг, равный максимальному (6), и выполняются шаги из точки Т в точку Т1 и из точки Т1 в точку Т2. Результаты записываются, соответственно, во второй и третий столбцы массивов X и F. Затем, для проверки точности выполняется удвоенный шаг из точки Т в точку Т2. Результаты этого шага записываются в четвертый столбец массива Х и в первый столбец массива F . В цикле (13, 14) накапливается критерий погрешности ЕН, как сумма взятых с весами погрешностей по каждому из уравнений. Погрешность каждой переменной вычисляется как 1/15 модуля разности между значениями этой переменной, вычисленными с разными шагами. Далее выполняется анализ критерия (15) и в зависимости от его значения шаг увеличивается, уменьшается или остается прежним. Если текущая погрешность ЕН не больше заданной Е , то результаты шага выводятся (25). При этом, если выполнялось два малых шага (КLР=1), то выводятся и результаты предыдущего шага (23). Так случается в начале интервала интегрирования и тогда, когда предыдущий шаг оказался неудачным и из-за большой погрешности величина шага уменьшена вдвое. После вывода двух шагов признак KLP сбрасывается в ноль (24). Выполненный шаг может быть последним на интервале (КР=1), тогда осуществляется выход из подпрограммы(26, 27). В блоке 28 выполняется проверка, может ли быть выполнен удвоенный шаг без выхода за пределы интервала? Если нет, то в (29) «взводится» признак конца интервала и устанавливается величина удвоенного шага равной оставшейся части интервала. В блоке 33 и цикле 34-35 последняя вычисленная точка делается начальной для выполнения двух малых шагов Н и контрольного удвоенного НВ. Соответственно, в 36 устанавливается признак двух шагов (KLP=1) и осуществляется возврат на блок 6 . Если «дошагивание» не нужно, то в 30 проверяется, является ли точность расчетов завышенной и в 31 можно ли удвоить малый шаг? При завышенной точности шаг можно удвоить, если он не превзойдет максимального НМ и удвоенный шаг не выведет за пределы интервала интегрирования. Если увеличение шага допустимо, то блок 32 это выполняет и далее все производится как при дошагивании, но без взвода признака конца. Если увеличение шага недопустимо, то в цикле 37, 38 выполняется подготовка к продолжению расчетов с прежним шагом. Из трех последних точек средняя делается начальной для выполнения контрольного шага удвоенной величины НВ, а последняя , - начальной для очередного малого шага Н.