Текст подпрограммы и версий ( Фортран )
de14r.zip , de14d.zip 		Тексты тестовых примеров ( Фортран )
tde14r.zip , tde14d.zip
Текст подпрограммы и версий ( Си )
de14r_c.zip , de14d_c.zip 		Тексты тестовых примеров ( Си )
tde14r_c.zip , tde14d_c.zip
Текст подпрограммы и версий ( Паскаль )
de14r_p.zip , de14e_p.zip 		Тексты тестовых примеров ( Паскаль )
tde14r_p.zip , tde14e_p.zip

Подпрограмма:  DE14R

Назначение

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

Математическое описание

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

          Y ' = F (X, Y) ,
          Y = ( y1, ... , yM ) ,
          F = ( f1 (X, y1, ... , yM), ... , fM (X, y1, ... , yM) )

 с начальными условиями, заданными в точке XN :
          Y(XN) = YN ,     YN = ( y10, ... , yM0 ) ,

классическим методом Рунге - Кутта четвертого порядка точности с постоянным шагом. Решение вычисляется в одной точке ХК, которая является концом интервала интегрирования.

Березин И.С., Жидков Н.П. Методы вычислений, т.2. Физматгиз, М., 1960.

Использование 

    SUBROUTINE  DE14R (F, M, XN, YN, XK, H, Y, A, B, C) 
   Параметры
F - имя подпрограммы вычисления значений правой части дифференциального уравнения. Первый оператоp подпрограммы должен иметь вид:
SUBROUTINE  F (X, Y , DY, M).
Здесь: X, Y - значения независимой и зависимой переменных, соответственно. Вычисленное значение правой части должно быть помещено в DY. B случае системы уравнений, т.е. когда M ≠ 1 , параметры Y и DY представляют массивы длины M (тип параметров X, Y и DY: вещественный);
M - количество уравнений в системе (тип: целый);
XN, YN - начальные значения аргумента и решения; в случае системы уравнений (т.е. M ≠ 1) YN представляет одномерный массив длины M (тип: вещественный);
XK - значение аргумента, при котоpом требуется вычислить решение задачи Коши (конец интервала интегрирования); XK может быть больше, меньше или pавно XN (тип: вещественный);
H - вещественная переменная, содержащая значение шага интегрирования. Может задаваться с учетом направления интегрирования, т.е. положительным, если XK > XN, отрицательным, если  XK < XN , или без такого учета с любым знаком " + " или " - ";
Y - искомое решение задачи Коши, вычисленное подпрограммой при значении аргумента XK. Для системы уравнений (когда M ≠ 1) задается одномерным массивом длины M. B случае совпадения значений параметров XN и XK значение Y полагается равным начальному значению YN (тип: вещественный);
A, B, C - одномерные вещественные рабочие массивы длиной M. 
   Версии
DE14D - вычисление решения задачи Коши для системы обыкновенных дифференциальных уравнений первого порядка в конце интервала интегрирования классическим методом Рунге - Кутта 4 - го порядка точности с постоянным шагом, при этом все вычисления над действительными числами выполняются с удвоенным числом значащих цифр. B этом случае параметры XN, YN, XK, H, Y, A, B, C и параметры X, Y и DY в подпрограмме F должны иметь тип DOUBLE PRECISION. 
   Вызываемые подпрограммы:  нет 
   Замечания по использованию
 

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

При работе подпрограммы и ее версии значения параметров M, XN, YN, XK сохраняются. Значение параметра H сохpаняется, если он задан с учетом направления интегрирования, иначе его знак меняется на противоположный. Если после работы подпрограммы нет необходимости иметь начальное значение решения YN, то параметры YN и Y при обращении к ней можно совместить.

Пример использования 

          y1'  =  0.2 ( y4 - y1 ) ,
          y2'  =  y1 + 2 ( y2 - y2 y3 ) ,
          y3'  =  y4 - ( y3 - y2 y3 ) ,
          y4'  =  10 y1 - ( 61 - 0.13 x ) y4 + 0.13 x ,     0 ≤ x ≤ 8
     
          y1 (0)  =  y2 (0)  =  y3 (0)  =  y4 (0)  =  0

          SUBROUTINE  F (X, Y, DY, M)
          DIMENSION  Y(4), DY(4)
          R23 = Y(2)*Y(3)
          CT = 0.13*X
          DY(1) = 0.2* ( Y(4) - Y(1) )
          DY(2) = Y(1) + 2.* ( Y(2) - R23 )
          DY(3) = Y(4) - ( Y(3) - R23 )
          DY(4) = 10.*Y(1) - ( 61. - CT)*Y(4) + CT
          RETURN
          END
   
          DIMENSION  YN(4), Y(4), A(4), B(4), C(4)
          EXTERNAL  F
          M = 4
          XN = 0.
          XK = 8.
          YN(1) = 0.
          YN(2) = 0.
          YN(3) = 0.
          YN(4) = 0.
          H = 0.01
          CALL DE14R (F, M, XN, YN, XK, H, Y, A, B, C)

Результаты:

          Y(1)  =  9.23847083543 - 03
          Y(2)  =  4.82097225727 - 03
          Y(3)  =  1.66711330351 + 00
          Y(4)  =  1.88434937053 - 02