Текст подпрограммы и версий ( Фортран )
de01r.zip , de01d.zip , de85r.zip , de85d.zip 		Тексты тестовых примеров ( Фортран )
tde01r.zip , tde01d.zip , tde85r.zip , tde85d.zip
Текст подпрограммы и версий ( Си )
de01r_c.zip , de01d_c.zip , de85r_c.zip , de85d_c.zip 		Тексты тестовых примеров ( Си )
tde01r_c.zip , tde01d_c.zip , tde85r_c.zip , tde85d_c.zip
Текст подпрограммы и версий ( Паскаль )
de01r_p.zip , de01e_p.zip , de85r_p.zip , de85e_p.zip 		Тексты тестовых примеров ( Паскаль )
tde01r_p.zip , tde01e_p.zip , tde85r_p.zip , tde85e_p.zip

Подпрограмма:  DE01R (версия DE85R)

Назначение

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

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

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

(1)                             Y '  =  F(X, Y) ,

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

с начальными условиями, заданными в точке XN: 

                         Y(XN)  =  YN ,   YN  =  ( y1,..., yM )  .

Предполагается, что среди характеристических корней матрицы Якоби ∂F / ∂y функции F имеются большие по модулю корни. Решение вычисляется в одной точке XK, которая является концом интервала интегрирования. Предполагается также, что система (1) может быть достаточно хорошо аппроксимирована линейной однородной системой с кусочнопостоянными коэффициентами, т.е. на каждом из [ xn, xn + H ], на которые разбивается весь интервал интегрирования, она может быть аппроксимирована системой вида 

                                     Y '  =   AY  ,

где A - постоянная матрица, вообще говоря, своя для каждого отрезка. Для интегрирования системы применяется метод Лоусона.

Метод Лоусона является одношаговым методом и заключается в следующем. Допустим, что искомое решение системы (1) уже вычислено в некоторой точке  x = xn  интервала интегрирования, т.е. известно  yn ≈ y(xn). Для отыскания решения  y(xn + 1) = y(xn + H) в следующем узле  xn + 1 = xn + H выполняются такие действия. Исходная система уравнений с помощью замены искомой функции  y (x) на  xn ≤ x ≤ xn + H  по формуле 

                             y(x)  =  exp [ ( x - xn ) A0 ] Z(x)  ,

где A0 - некоторая постоянная матрица, преобразуется в систему уравнений относительно новой неизвестной функции Z(X): 

(2) Z ' (x)  =  F1(x, z)  =  exp [ - ( x - xn ) A0 ] { F(x, exp [ ( x - xn ) A0 ] Z(x)) -
                                      - A0 exp [ ( x - xn ) A0 ] Z(x) }

Матрица Якоби ∂ F1 / ∂Z системы (2) и матрица Якоби ∂F / ∂y системы (1) связаны между собой соотношением 

         ∂ F1 / ∂Z  =  exp [ - (x - xn) A0] { ∂F / ∂y - A0 } exp [ ( x - xn ) A0 ]

Указанное преобразование выполняется самой подпрограммой. В качестве матрицы A0 подпрограмма выбирает матрицу ∂F / ∂y, вычисленную либо в точке ( xn, yn ), либо в некоторой предыдущей точке ( xn - k, yn - k ). Если шаг H достаточно мал, то данное преобразование позволяет уменьшить характеристические корни матрицы ∂ F1 / ∂Z по сравнению с характеристическими корнями матрицы ∂F / ∂y. А это приводит к уменьшению константы Липшица системы (2) по сравнению с константой Липшица системы (1).

Для решения системы (2) применяются классический метод Рунге - Кутта четвертого порядка точности, причем одновременно с решением (2) производится обратное преобразование от функции Z (x) к функции Y (x).

Все компоненты решения вычисляются с контролем точности по мере погрешности, который заключается в следующем. Если некоторая компонента приближенного решения по абсолютной величине не меньше некоторой наперед заданной константы P, то контроль точности для этой компоненты ведется по относительной погрешности, иначе - по абсолютной. Абсолютная погрешность приближенного решения оценивается по правилу Рунге.

J.Douglas Lowson, Generalized Runge - Kutta processes for stable systems with large lipshitz constants, SIAM Journal on Numerical Analisys, 1967, Vol.4, No 3.

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

    SUBROUTINE  DE01R (F, FJ, M, XN, YN, XK, HMIN, EPS, P, N, H, Y,
                                             R, IERR) 
   Параметры
F - подпрограмма вычисления правой части системы в точке X. Первый оператор подпрограммы должен иметь вид:
   SUBROUTINE  F (X, Y, DY, M).
Здесь X и Y - значения независимой и зависимой переменных, соответственно. Вычисленное значение правой части должно быть помещено в DY. В случае системы уравнений, т.е. M ≠ 1, параметры Y и DY представляют одномерные массивы длиной M (тип параметров X, Y и DY: вещественный);
FJ - подпрограмма вычисления матрицы Якоби правой части системы. Первый оператор подпрограммы должен иметь вид:
   SUBROUTINE  FJ (X, Y, DF, M).
Здесь X и Y - значения независимой и зависимой переменных, соответственно, причем Y представляет одномерный массив длины M. DF - двумерный массив размера M*M. Вычисленное значение матрицы Якоби ∂F / ∂y должно быть помещено в массив DF, при этом частная производная от правой части I - ого уравнения F по J - ой переменной Y (J) запоминается в DF (I, J), т.е. DF (I, J) = ∂Fi / ∂yj (тип параметров X, Y и DF: вещественный);
M - количество уравнений в системе (тип: целый);
XN, YN - начальные значения аргумента и решения; в случае системы уравнений (т.е. M ≠ 1) YN представляет одномерный массив длины M;
XK - значение аргумента, при котором требуется вычислить решение задачи Коши (конец интервала интегрирования); XK может быть больше, меньше или равно XN (тип: вещественный);
HMIN - минимальное значение абсолютной величины шага, которое разрешается использовать при интегрировании данной системы уравнений (тип: вещественный);
EPS - допустимая мера погрешности, с которой требуется вычислить все компоненты решения (тип: вещественный);
P - граница перехода, используемая при оценке меры погрешности решения (тип: вещественный);
N - целый параметр, управляющий режимом вычисления матрицы A0 и имеющий следующие значения:
Если N ≤ 0 , (т.е. имеет произвольное отрицательное или нулевое значение), то матрица A0 вычисляется только один раз и полагается равной значению A0 = ∂F / ∂y в точке [x = x0, y = y0]  на протяжении всего интервала интегрировария;
если N > 0, то матрица A0 изменяется периодически, при этом она остается постоянной на протяжении каждого ряда из N последовательных шагов и вычисляется вновь на первом шаге каждого ряда, т.е. матрица A0 = ∂F / ∂y вычисляется только в следующих точках:
[x = x0, y = y0]   при x0 ≤ x ≤ xN;
[x = xN, y = yN]   при xN < x ≤ x2N;
[x = x2N, y = y2N]  при x2N < x≤ x3N;
Чем сильнее меняется матрица Якоби ∂f / ∂y системы при изменении  x, тем меньше должно быть значение параметра N. В частности, если N = 1, то матрица Якоби системы ∂f / ∂y будет вычисляться в каждой точке  xn,  yn .
H - вещественная переменная, содержащая начальное значение шага интегрирования; может задаваться с учетом направления интегрирования, т.е. положительным, если XN < XK, отрицательным, если XN > XK, или без всякого учета в виде абсолютной величины;
Y - искомое решение задачи Коши, вычисленное подпрограммой при значении аргумента XK; для системы уравнений (когда M ≠ 1) задается одномерным массивом длины M. В случае совпадения значений параметров XN и XK значение Y полагается равным начальному значению YN (тип: вещественный);
R - одномерный рабочий массив вещественного типа длины 5*M*M+7*M+1;
IERR - целая переменная, значение которой в результате работы подпрограммы полагается равным 65, если какая - нибудь компонента решения не может быть вычислена с требуемой точностью EPS. В этом случае интегрирование системы прекращается; при желании интегрирование системы можно повторить обращением к подпрограмме с новыми значениями параметров HMIN, H и N. 
   Версии
DE01D - вычисление решения задачи Коши для устойчивой системы обыкновенных дифференциальных уравнений первого порядка с большой константой Липшица в конце интервала интегрирования методом Лоусона с удвоенным числом значащих цифр. При этом параметры XN, YN, XK, HMIN, EPS, P, H, Y, R и параметры X, Y, DY в подпрограмме F и X, Y, DF в подпрограмме FJ должны иметь тип DOUBLE PRECISION.
DE85R - вычисление решения задачи Коши для устойчивой системы обыкновенных дифференциальных уравнений первого порядка с большой константой Липшица в конце интервала интегрирования методом Лоусона. В отличии от подпрограммы DE01R первый оператор подпрограммы DE85R не содержит в списке формальных параметров имени подпрограммы вычисления матрицы Якоби и имеет следующий вид:
SUBROUTINE  DE85R (F, M, XN, YN, XK, HMIN, EPS, P, N, H, Y, R, IERR)
При использовании этой подпрограммы пользователь не составляет подпрограмму вычисления матрицы Якоби. Все параметры подпрограммы DE85R имеют тот же смысл, что и одноименные параметры подпрограммы DE01R.
DE85D - вычисление решения задачи Коши для устойчивой системы обыкновенных дифференциальных уравнений первого порядка с большой константой Липшица в конце интервала интервала интегрирования методом Лоусона с удвоенным числом значащих цифр. Первый оператор подпрограммы имеет тот же вид, что и в подпрограмме DE85R, при этом параметры XN, YN, XK, HMIN, EPS, P, H, Y, R и параметры X, Y, DY в подпрограмме F должны иметь тип DOUBLE PRECISION. 
   Вызываемые подпрограммы
      DE00R -
      DE00D   		выполнение одного шага интегрирования системы обыкновенных дифференциальных первого порядка с большой константой Липшица методом Лоусона.
      UTDE20 -
      UTDE21   		подпрограммы выдачи диагностических сообщений.
  Подпрограммы DE00R, UTDE20 вызываются при работе подпрограммы DE01R, а подпрограммы DE00D, UTDE21 - при работе DE01D. 
   Замечания по использованию
 

Данная подпрограмма предназначена для интегрирования таких нелинейных систем, которые могут быть достаточно хорошо аппроксимированы линейной однородной системой с кусочнопостоянными коэффициентами, т.е. на каждом из отрезков [x, x + H] , на которые разбивается весь интервал интегрирования, они могут быть аппроксимированы системой вида

                     Y '  =  AY

где A - постоянная матрица, вообще говоря, своя для каждого отрезка.

В общем случае заданная точность не гарантируется.

При работе подпрограммы значения параметров M, XN, YN, XK, HMIN, EPS, P, N сохраняются.

При работе подпрограмм F и FJ значения параметров X, Y, M не должны изменяться.

Если после работы подпрограммы нет необходимости иметь начальное значение решения YN, то параметры YN и Y при обращении к ней можно совместить. При этом следует иметь в виду, что в случае аварийного выхода из подпрограммы, т.е. со значением IERR = 65, значение параметра YN будет испорчено.

Подпрограммы DE01R и DE01D предназначены также для решения задачи Коши для жестких дифференциальных уравнений (1).

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

                                     
     y1'  =  (- 500 + 2 x) y1  +  (1 - x /10 ) y2  +  (1 + 1/(x + 1) ) sin y1 ,
     y1(1) = 49 
     y2'  =  (- 1000 - 10 x) y1  +  (1 + x /10 ) y2  +  (2  + 1/(x + 2) ) cos y2 ,
     y2(1) = 50 

              1 ≤ X ≤ 6

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

      SUBROUTINE  F (X, Y, DY, M) 
      DIMENSION Y(2), DY(2) 
      DY(1) = (2.*X - 5.E2)*Y(1) + (1. - X*1.E-1)*Y(2) + (1. + 1./(1. + X))* 
     *              SIN(Y(1)) 
      DY(2) = (-1.E3 - 1.E1*X)*Y(1) + (1. + X*1.E-1)*Y(2) +
     *              (2. + 1./(2. + X))*COS(Y(2)) 
      RETURN 
      END 

      SUBROUTINE  FJ (X, Y, DF, M) 
      DIMENSION DF(2, 2), Y(2) 
      DF(1, 1) = 2.*X - 5.E2 + (1. + 1./(1. + X))*COS(Y(1)) 
      DF(1, 2) = 1. - X*1.E-1 
      DF(2, 1) = -1.E3 - 1.E1*X 
      DF(2, 2) = 1. + X*1.E-1 - (2. + 1./(2. + X))SIN(Y(2)) 
      RETURN 
      END 

      DIMENSION YN(2), Y(2), R(33) 
      EXTERNAL F, FJ 
      M = 2 
      XN = 1. 
      YN(1) = 49. 
      YN(2) = 50. 
      HMIN = 1.E-10 
      EPS = 1.E-8 
      P = 100 
      H = 0.01 
      XK = 0 
      N = 2 
      CALL  DE01R (F, FJ, M, XN, YN, XK, HMIN, EPS, P, N, H, Y, R,
     *                          IERR) 

Результаты: 

             Y(1)                                    Y(2)                              H 
     -4.175740488392-02       -5.087148371158+01      2.500000000001-03 
      IERR = 0

Результаты решения этой же задачи, полученные с помощью подпрограммы DE85R 

     -4.175740472033-02       -5.087148363417+01      5.000000000003-03 
      IERR = 0

Результаты решения этой же задачи, полученные с помощью подпрограммы DE25R 

     -4.175740429849-02       -5.087148203468+01      1.756831014061-03 
      IERR = 0