Текст подпрограммы и версий ( Фортран )
de03r.zip , de03d.zip 		Тексты тестовых примеров ( Фортран )
tde03r.zip , tde03d.zip
Текст подпрограммы и версий ( Си )
de03r_c.zip , de03d_c.zip 		Тексты тестовых примеров ( Си )
tde03r_c.zip , tde03d_c.zip
Текст подпрограммы и версий ( Паскаль )
de03r_p.zip , de03e_p.zip 		Тексты тестовых примеров ( Паскаль )
tde03r_p.zip , tde03e_p.zip

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

Назначение

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

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

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

(1)                         Y ' (X)   =   A * Y(X) + U(X, Y)  ,

                                   Y    =   ( y1,..., yM ) , 
                                   A    =   ( ai j ) ,     i, j  =  1, ..., M ,
                           U(X, Y)   =   ( U1(X, Y), ... , UM(X, Y) ) ,

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

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

где A - постоянная числовая матрица.

Предполагается, что среди характеристических корней матрицы A имеются большие по модулю корни, а константа Липшица для функции U (X, Y), т.е. константа L из условия Липшица 

| Ui ( x, y1(1),..., yM(1) ) - Ui ( x, y1(2), ... , yM(2) ) |  ≤
                                                                                      M
                                                                             ≤  L  ∑   | yj(1) - yj(2)  |  , 
                                                                                     j=1

независящая от  i,  x,  y(1),  y(2) , невелика. Также предполагается, что нелинейный член U (x, y) является достаточно малым. Решение вычисляется в одной точке ХК, которая является концом интервала интегрирования. Для интегрирования системы применяется метод Лоусона.

Метод Лоусона является одношаговым и заключается в следующем. Допустим, что искомое решение системы (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 ) A ] Z(x)  ,

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

(2)   Z ' (x)  =  U1(x, z)  =  exp [ - ( x - xn ) A ]  U( x, exp [( x - xn ) A ] Z(x) )

Данное преобразование выполняется самой подпрограммой. Характеристические корни матрицы Якоби 

   ∂U1 / ∂Z   =   exp [ - ( x - xn ) A ]  (∂U / ∂y)  exp [ ( x - xn ) A ]

являются характеристическими корнями матрицы ∂U / ∂y и в силу малости константы Липшица функции U (x, y) невелики. Поэтому система (2) не жесткая и может быть решена традиционными методами численного интегрирования. Данная подпрограмма решает ее методом Рунге - Кутта  4 - ого порядка, при этом одновременно с решением (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  DE03R (FU, M, XN, YN, XK, A, HMIN, EPS, P, H, Y,
                                             R, IERR) 
   Параметры
FU - подпрограмма вычисления функции U (x, y) в правой части системы. Первый оператор подпрограммы должен иметь вид:
   SUBROUTINE  FU (X, Y, U, M).
Здесь Y и U - одномерные массивы длины M. В массив U помещается значение функции U (x, y), вычисленное при значении аргументов X и Y (тип параметров X, Y, U: вещественный);
M - количество уравнений в системе (тип: целый);
XN, YN - начальные значения аргумента и решения; в случае системы уравнений (т.е. M ≠ 1) YN представляет одномерный массив длины M (тип: вещественный);
XK - значение аргумента, при котором требуется вычислить решение задачи Коши (конец интервала интегрирования); XK может быть больше, меньше или равно XN (тип: вещественный);
A - вещественный двумерный массив размера M*M, содержащий элементы матрицы A;
HMIN - минимальное значение абсолютной величины шага, которое разрешается использовать при интегрировании данной системы уравнений (тип: вещественный);
EPS - допустимая мера погрешности, с которой требуется вычислить все компоненты решения (тип: вещественный);
P - граница перехода, используемая при оценке меры погрешности решения (тип: вещественный);
H - вещественная переменная, содержащая начальное значение шага интегрирования; может задаваться с учетом направления интегрирования, т.е. положительным, если XN < XK, отрицательным, если XN > XK, или без всякого учета в виде абсолютной величины;
Y - искомое решение задачи Коши, вычисленное подпрограммой при значении аргумента XK; для системы уравнений (когда M ≠ 1) задается одномерным массивом длины M. В случае совпадения значений параметров XN и XK значение Y полагается равным начальному значению YN (тип: вещественный);
R - вещественный одномерный рабочий массив длины 4*M*M+7*M+1;
IERR - целая переменная, значение которой в результате работы подпрограммы полагается равным 65, если какая - нибудь компонента решения не может быть вычислена с требуемой точностью EPS. В этом случае интегрирование системы прекращается; при желании интегрирование системы можно повторить обращением к подпрограмме с новыми значениями параметров HMIN и H. 
   Версии
DE03D - вычисление решения задачи Коши для квазилинейной системы обыкновенных дифференциальных уравнений первого порядка с большой константой Липшица в конце интервала интегрирования методом Лоусона с удвоенным числом значащих цифр. При этом параметры XN, YN, XK, A, HMIN, EPS, P, H, Y, R и параметры X, Y, R в подпрограмме FU должны иметь тип DOUBLE PRECISION. 
   Вызываемые подпрограммы
      DE02R -
      DE02D   		выполнение одного шага численного интегрирования квазилинейной системы обыкновенных дифференциальных уравнений первого порядка с большой константой Липшица методом Лоусона.
      UTDE20 -
      UTDE21   		подпрограмма выдачи диагностических сообщений; Подпрограмма DE02R, UTDE20 вызывается при работе подпрограммы DE03R, а подпрограммы DE02D, UTDE21 - при работе DE03D. 
   Замечания по использованию
 

Данная подпрограмма предназначена для интегрирования квазилинейных систем, имеющих малый нелинейный член U (x, y).

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

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

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

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

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

                                     
        y1'  =  - 500 y1  +  y2  +  sin ( y1 + y2 ) ( 1 + e - x/2 ) ,
        y1(0) = 50 
        y2'  =  - 1000 y1  +  y2  +  sin ( y1 - y2 ) ( 1 + e - x/2 ) ,
        y2(0) = 50 

        0 ≤ x ≤ 5

Приводятся подпрограмма вычисления функций 

        U1 (x, y)  =  sin ( y1 + y2 ) ( 1 + e - x/2 ) ,
        U2 (x, y)  =  sin ( y1 - y2 ) ( 1 + e - x/2 )

из правой части системы, фрагмент вызывающей программы и результаты счета. 

      SUBROUTINE  FU (X, Y, U, M) 
      DIMENSION Y(2), U(2) 
      T = 1. + EXP(-X/2.) 
      U(1) = T*SIN(Y(1) + Y(2)) 
      U(2) = T*SIN(Y(1) - Y(2)) 
      RETURN 
      END 

      DIMENSION A(2, 2), Y(2), YN(2), R(31) 
      EXTERNAL FU 
      M = 2 
      XN = 0. 
      YN(1) = 5.E1 
      YN(2) = 5.E1 
      HMIN = 1.E-10 
      EPS = 1.E-8 
      P = 1.E2 
      XK = 5. 
      H = 0.01 
      A(1, 1) = -5.E2 
      A(1, 2) = 1. 
      A(2, 1) = -1.E3 
      A(2, 2) = 1. 
      CALL  DE03R (FU, M, XN, YN, XK, A, HMIN, EPS, P, H, Y, R, IERR) 

 Результаты: 

             Y(1)                                   Y(2)                               H 
     -3.508894600484-04      -8.340807464992-02       1.250000000001-03 

      IERR = 0