|
Текст подпрограммы и версий ( Фортран ) de00r.zip , de00d.zip , de84r.zip , de84d.zip |
Тексты тестовых примеров ( Фортран ) tde00r.zip , tde00d.zip , tde84r.zip , tde84d.zip |
|
Текст подпрограммы и версий ( Си ) de00r_c.zip , de00d_c.zip , de84r_c.zip , de84d_c.zip |
Тексты тестовых примеров ( Си ) tde00r_c.zip , tde00d_c.zip , tde84r_c.zip , tde84d_c.zip |
|
Текст подпрограммы и версий ( Паскаль ) de00r_p.zip , de00e_p.zip , de84r_p.zip , de84e_p.zip |
Тексты тестовых примеров ( Паскаль ) tde00r_p.zip , tde00e_p.zip , tde84r_p.zip , tde84e_p.zip |
Выполнение одного шага численного интегрирования устойчивой системы обыкновенных дифференциальных уравнений первого порядка с большой константой Липшица методом Лоусона.
Выполняется один шаг численного интегрирования системы M обыкновенных дифференциальных уравнений первого порядка
(1) Y ' = F(X, Y) ,
Y = ( y1,..., yM ) , F = ( f1( X, y1,..., yM ), ... , fM( X, y1,..., yM ) )
Предполагается, что среди характеристических корней матрицы Якоби ∂F / ∂y функции F имеются большие по модулю корни. По заданному значению решения YX в узле xn вычисляется значение этого решения в узле xn + H . Предполагается также, что система (1) на [ xn, xN + H ] может быть достаточно хорошо аппроксимирована линейной однородной системой с постоянными коэффициентами, т.е. системой вида
Y ' = AY ,
где A - некоторая постоянная матрица. Вычисление производится по методу Лоусона. Метод Лоусона заключается в следующем. Исходная система уравнений с помощью замены искомой функции y (x) на [ xn, 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 DE00R (F, FJ, M, JSTART, HMIN, EPS, P, YX, X, H,
BUL, XP, YP, A, E1, E2, E3, E4, R, R0, R1, R2,
R3, YD, 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 - | количество уравнений в системе (тип: целый); |
| JSTART - | целый указатель режима использования подпрограммы, имеющий следующие значения: |
| 0 - | первое обращение к подпрограмме должно быть выполнено с нулевым значением JSTART; при этом значение параметра в качестве матрицы A0 принимается матрица ∂f / ∂y в точке [ x = xn, y = yn ]; вычисленное значение матрицы A0 будет использоваться в дальнейшем для значений JSTART = 1, 2, 3; |
| 1 - | выполнить один шаг интегрирования системы дифференциальных уравнений для значений независимой и зависимой переменных, заданных параметрами YX и X соответственно; при этом само значение величины шага интегрирования H, на который предполагается продолжить решение системы, должно совпадать с тем значением шага, который был фактически выполнен подпрограммой при предыдущем обращении к ней, а матрица A0 остается той же, что и на предыдущем шаге; |
| 2 - | то же, что и для JSTART = 1, но только с той разницей, что величина шага интегрирования, на который предполагается продолжить решение, в два раза больше того значения шага, который был фактически выполнен при предыдущем обращении к подпрограмме; матрица A0 остается той же, что и на предыдущем шаге; |
| 3 - | то же, что и для JSTART = 1, но только с той разницей, что шаг интегрирования может принимать произвольное значение, матрица A0 остается той же, что и на предыдущем шаге. |
| -1 - | повторить последний шаг интегрирования с новыми значениями параметров H и/или HMIN; при этом значении JSTART в качестве матрицы A0 принимается матрица ∂f / ∂y в точке [ x = XP, y = YP ] таким образом, до тех пор, пока JSTART принимает значения 1, 2 или 3, сохраняется одно и то же значение матрицы A0 , равное значению матрицы Якоби, вычисленному при последнем обращении к подпрограмме с параметром JSTART = 0 или JSTART = - 1; если требуется изменить матрицу A0 , то следует обратиться к подпрограмме со значением JSTART = 0 или - 1; при выходе из подпрограммы параметр JSTART полагается равным 2, если выполнив заданный в H шаг, подпрограмма рекомендовала для интегрирования вдвое большеe его значение, и 1 в противном случае, т.е. в том случае, когда рекомендованное значение шага равно только что выполненному шагу; рекомендованное значение шага на выходе из подпрограммы запоминается в переменной H; |
| HMIN - | минимальное значение абсолютной величины шага, которое разрешается использовать при интегрировании данной системы уравнений (тип: вещественный); |
| EPS - | допустимая мера погрешности, с которой требуется вычислить все компоненты решения (тип: вещественный); |
| P - | граница перехода, используемая при оценке меры погрешности решения (тип: вещественный); |
| YX, X - | заданные вещественные значения решения и соответствующее ему значение аргумента; в результате работы подпрограммы в X получается новое значение аргумента, в YX - соответствующее значение решения; в случае системы уравнений т.е. когда M ≠ 1, YX задается одномерным массивом длины M; |
| H - | вещественная переменная, содержащая значение шага интегрирования; если для этого значения шага точность приближенного решения достигается, то именно он и реализуется подпрограммой, иначе этот шаг уменьшается подпрограммой до тех пор, пока не будет достигнута заданная точность EPS; на выходе из подпрограммы H содержит рекомендуемое подпрограммой значение следующего шага, определяемое ею с целью достижения более экономного способа интегрирования; |
| BUL - | логическая переменная, значение которой при обращении к подпрограмме полагается равным .TRUE., если заданный в H шаг выводит в конец интервала интегрирования, и .FALSE. в противном случае; в результате работы подпрограммы BUL равно .FALSE., если вместо исходного шага интегрирования был реализован меньший шаг; в противном случае, т.е. когда был выполнен именно заданный при обращении в H шаг, значение параметра BUL не меняется; |
| XP, YP - | вещественная рабочая переменная и одномерный рабочий массив длины M соответственно; значения параметров XP, YP на выходе из подпрограммы равны тем значениям, которые имели параметры X, YX при входе в нее (т.е. предыдущий узел и решение в нем); |
|
A, E1, E2 - E3, E4 | вещественные двумерные рабочие массивы размера M*M; |
|
R, R0 - R1, R2 R3, YD | вещественные одномерные рабочие массивы длины M; |
| IERR - | целая переменная, значение которой в результате работы подпрограммы полагается равным 65, если какая - нибудь компонента решения не может быть вычислена с требуемой точностью EPS; в этом случае последний шаг интегрирования системы можно повторить обращением к подпрограмме с новыми значениями параметров H, HMIN и значением JSTART = - 1 . |
Версии
| DE00D - | выполнение одного шага численного интегрирования системы обыкновенных дифференциальных уравнений первого порядка с большой константой Липшица методом Лоусона с удвоенным числом значащих цифр. При этом параметры HMIN, EPS, P, YX, X, H, XP, YP, A, E1, E2, E3, E4, R, R0, R1, R2, R3, YD и параметры X, Y, DY в подпрограмме F и X, Y, DF в подпрограмме FJ должны иметь тип DOUBLE PRECISION. |
| DE84R - |
выполнение одного шага численного интегрирования
устойчивой системы обыкновенных дифференциальных
уравнений первого порядка с большой константой Липшица
методом Лоусона. В отличии от подпрограммы DE00R
первый оператор подпрограммы DE84R не содержит в
списке формальных параметров имени подпрограммы
вычисления матрицы Якоби и имеет следующий вид: SUBROUTINE DE84R (F, M, JSTART, HMIN, EPS, P, YX, X, H, BUL, XP, YP, A, E1, E2, E3, E4, R, R0, R1, R2, R3, YD, IERR). При использовании этой подпрограммы пользователь не составляет подпрограмму вычисления матрицы Якоби. Все параметры подпрограммы DE84R имеют тот же смысл, что и одноименные параметры подпрограммы DE00R. |
| DE84D - | выполнение одного шага численного интегрирования устойчивой системы обыкновенных дифференциальных уравнений первого порядка с большой константой Липшица методом Лоусона с удвоенным числом значащих цифр. Первый оператор подпрограммы имеет тот же вид, что и в подпрограмме DE84R, при этом параметры HMIN, EPS, P, YX, X, H, XP, YP, A, E1, E2, E3, E4, R, R0, R1, R2, R3, YD и параметры X, Y, DY в подпрограмме F должны иметь тип DOUBLE PRECISION. |
Вызываемые подпрограммы
| AME2R - | подпрограмма вычисления матричной экспоненты; вызывается при работе подпрограммы DE00R. |
| AME2D - | подпрограмма вычисления матричной экспоненты; вызывается при работе подпрограммы DE00D. |
| UTDE20 - | подпрограмма выдачи диагностических сообщений при работе подпрограммы DE00R. |
| UTDE21 - | подпрограмма выдачи диагностических сообщений при работе подпрограммы DE00D. |
| Кроме того, при работе подпрограмм DE00R, DE84R и DE00D, DE84D вызываются рабочие подпрограммы DE00RS, DE04RM, DE00RD и DE00DS, DE04DM, DE00DD соответственно. | |
| При работе подпрограмм DE84R и DE84D дополнительно вызываются также рабочие подпрограммы DE84RJ и DE84DJ соответственно. |
Замечания по использованию
|
Данная подпрограмма предназначена для интегрирования таких
нелинейных систем, которые могут быть достаточно хорошо
аппроксимированы линейной однородной системой с постоянными
коэффициентами, т.е. системой вида Y ' = AY , где A - некоторая постоянная матрица. В общем случае заданная точность не гарантируется. При работе подпрограммы значения параметров M, HMIN, EPS, P сохраняются. Между последовательными обращениями к подпрограмме со значениями параметра JSTART = 1, 2, 3 пользователь не должен изменять содержимое массивов A, E1, E2. При работе подпрограмм F и FJ значения параметров M, X, Y не должны изменяться. При обращении к подпрограмме со значением JSTART = - 1 в качестве исходных значений аргумента и решения принимаются значения параметров XP и YP соответственно, т.е. те значения, которые эти параметры получили после самого последнего обращения к подпрограмме с неотрицательным значением JSTART. После работы подпрограммы в массиве R1 содержится значение оценки абсолютной погрешности на шаге, вычисленной по правилу Рунге. Данная подпрограмма и ее версия предназначены также для интегрирования жестких дифференциальных уравнений (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 YX(2), YP(2), A(2, 2), E1(2, 2), E2(2, 2), E3(2, 2),
* E4(2, 2), R(2), R0(2), R1(2), R2(2), R3(2), YD(2)
LOGICAL BUL
EXTERNAL F, FJ
M = 2
X = 1.
YX(1) = 49.
YX(2) = 50.
HMIN = 1.E-10
EPS = 1.E-8
P = 1.E2
JSTART = 0
H = 0.01
IH = 0
10 IH = IH + 1
CALL DE00R (F, FJ, M, JSTART, HMIN, EPS, P, YX, X, H, BUL, XP,
* YP, A, E1, E2, E3, E4, R, R0, R1, R2, R3, YD, IERR)
GO TO (11, 12, 13, 14, 15, 16), IH
11 H = 1.
JSTART = 3
GO TO 10
12 EPS = 1.E-2
GO TO 10
13 JSTART = -1
EPS = 1.E-8
GO TO 10
14 EPS = 1.E-2
H = 1.E-8
JSTART = 3
GO TO 10
15 H = -1.E-3
EPS = 1.E-8
JSTART = 0
GO TO 10
16 CONTINUE
Результаты:
После первого обращения к подпрограмме
X YX(1) YX(2)
1.000039062499+00 4.895766071531+00 4.808757932455+00
H R1(1) R1(2)
3.906250000002-05 9.597251282538-10 -7.233271996176-09
После второго обращения к подпрограмме
X YX(1) YX(2)
1.000069580077+00 4.733408100787+01 4.661901490984+01
H R1(1) R1(2)
3.051757812500-05 1.319373647371-10 6.946114202335-02
После третьего обращения к подпрограмме
X YX(1) YX(2)
1.000100097655+00 4.662140584405+01 4.517257692956+01
H R1(1) R1(2)
6.103515625000-05 -6.208817164108-11 1.830362084864-09
После четвертого обращения к подпрограмме
X YX(1) YX(2)
1.000100097655+00 4.662140584411+01 4.517257692973+01
H R1(1) R1(2)
3.051757812500-05 -5.820766091346-11 1.851003617047-09
После пятого обращения к подпрограмме
X YX(1) YX(2)
1.000100107654+00 4.662117408356+01 4.517210655950+01
H R1(1) R1(2)
2.000000000004-08 0.000000000000 + 00 0.000000000000 + 00
После шестого обращения к подпрограмме
X YX(1) YX(2)
1.000068857653+00 4.735108317091+01 4.665352269082+01
H R1(1) R1(2)
-3.125000000004-05 3.880510727560-11 -2.250696221985-09