Документ взят из кэша поисковой машины. Адрес оригинального документа : http://num-anal.srcc.msu.su/lib_na/cat/de/de13r.htm
Дата изменения: Tue Dec 1 10:39:51 2015
Дата индексирования: Sun Apr 10 00:11:40 2016
Кодировка: Windows-1251
БЧА НИВЦ МГУ. DE13R. Задача Коши для системы уравнений первого порядка (с контролем точности)
Текст подпрограммы и версий ( Фортран )
de13r.zip , de13d.zip 		Тексты тестовых примеров ( Фортран )
tde13r.zip , tde13d.zip
Текст подпрограммы и версий ( Си )
de13r_c.zip , de13d_c.zip 		Тексты тестовых примеров ( Си )
tde13r_c.zip , tde13d_c.zip
Текст подпрограммы и версий ( Паскаль )
de13r_p.zip , de13e_p.zip 		Тексты тестовых примеров ( Паскаль )
tde13r_p.zip , tde13e_p.zip

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

Назначение

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

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

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

          Y ' = F (X, Y) ,
          Y = ( y1, ... , yM ) ,
          F = ( f1 (X, y1, ... , yM), ... , fM (X, y1, ... , yM) )
 с начальными условиями, заданными в точке XN :

          Y(XN) = YN ,     YN = ( y10, ... , yM0 ) ,

классическим методом Рунге - Кутта 4 - го порядка с контрольным членом Егорова. Решение вычисляется в одной точке XK, которая является концом интервала интегрирования. Каждая компонента решения вычисляется с контролем точности по относительной погрешности на тех участках интервала интегрирования, на которых модуль этой компоненты больше некоторого наперед заданного числа P (это число называется границей перехода), и по абсолютной погрешности на остальных участках, т.е. там, где модуль проверяемой на точность компоненты меньше этого числа.

О.Б.Арушанян, Стандартная программа решения задачи Коши для систем обыкновенных дифференциальных уравнений методом Рунге - Кутта, вып. 31, под общей редакцией В.В.Воеводина, НИВЦ МГУ, 1968.

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

    SUBROUTINE DE13R (F, M, XN, YN, XK, HMIN, EPS, P, H, 
                                           Y, YP, DELTY, YR, DY, IERR) 
   Параметры
F - имя подпрограммы вычисления значений правой части дифференциального уравнения. Первый оператоp подпрограммы должен иметь вид:
SUBROUTINE  F (X, Y , DY, M).
Здесь: X, Y - значения независимой и зависимой переменных, соответственно. Вычисленное значение правой части должно быть помещено в DY. B случае системы уравнений, т.е. когда M ≠ 1 , параметры Y и DY представляют массивы длины M (тип параметров X, Y и DY: вещественный);
M - количество уравнений в системе (тип: целый);
XN, YN - начальные значения аргумента и решения. B случае системы уравнений (т.е. M ≠ 1) YN представляет одномерный массив длины M (тип: вещественный);
XK - значение аргумента, при котоpом требуется вычислить решение задачи Коши (конец интервала интегрирования). XK может быть больше, меньше или pавно XN (тип: вещественный);
HMIN - минимальное значение абсолютной величины шага, который разрешается использовать при интегрировании данной системы уравнений (тип: вещественный);
EPS - допустимая меpа погрешности, с которой тpебуется вычислить все компоненты решения (тип: вещественный);
P - граница перехода, используемая при оценке погрешности решения (тип: вещественный);
H - вещественная переменная, содержащая начальное значение шага интегрирования. Может задаваться с учетом направления интегрирования, т.е. положительным, если XK > XN, отрицательным, если XK < XN, или без такого учета в виде абсолютной величины;
Y - искомое решение задачи Коши, вычисленное подпрограммой при значении аргумента XK. Для системы уравнений (когда M ≠ 1) задается одномерным массивом длины M. B случае совпадения значений параметров XN и XK значение Y полагается равным начальному значению YN (тип: вещественный);
            YP -
       DELTY  
       YR, DY   		вещественные одномерные рабочие массивы длины M;
IERR - целая переменная, значение которой в pезультате работы подпрограммы полагается pавным 65, если какая - нибудь компонента решения не может быть вычислена с требуемой точностью EPS. B этом случае интегрирование системы прекращается. При желании интегрирование системы можно повторить обращением к подпрограмме с новыми значениями параметpов HMIN и H. 
   Версии
DE13D - вычисление решения задачи Коши для системы обыкновенных дифференциальных уравнений первого порядка в конце интервала интегрирования классическим методом Рунге - Кутта четвертого порядка с повышенной точностью. При этом параметры XN, YN, XK, HMIN, EPS, P, H, Y, YP, DELTY, YR, DY и параметры X, Y и DY в подпрограмме F должны иметь тип DOUBLE PRECISION. 
   Вызываемые подпрограммы
UTDE10 - подпрограмма выдачи диагностических сообщений при работе подпрограммы DE13R.
UTDE11 - подпрограмма выдачи диагностических сообщений при работе подпрограммы DE13D. 
   Замечания по использованию
 

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

Хотя заданная точность EPS не гарантируется в общем случае, большой опыт эксплуатации данной подпрограммы убедительно показывает, что вычисляемое ею численное решение достаточно близко приближает точное решение.

При работе подпрограммы значения параметров M, XN, YN, XK, HMIN, EPS, 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)
          CT = 0.13*X
          R23 = Y(2)*Y(3)
          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  YI(4), YF(4), YP(4), DELTY(4), RAB(4), RAB1(4)
          EXTERNAL  F
          M = 4
          YI(1) = 0.
          YI(2) = 0.
          YI(3) = 0.
          YI(4) = 0.
          XI = 0.
          XF = 8.
          H = 0.01
          HMIN = 1.E-16
          ERR = 1.E-4
          P = 1.E-8
          CALL DE13R (F, M, XI, YI, XF, HMIN, ERR, P, H, YF, YP, 
        *                         DELTY, RAB, RAB1, IERR)

Результаты:

          YF(1)  =  0.00923847083381
          YF(2)  =  0.00482097150315
          YF(3)  =  1.66711319966
          YF(4)  =  0.0188434937933

          IERR  =  0