|
Текст подпрограммы и версий ( Фортран ) de02r.zip , de02d.zip |
Тексты тестовых примеров ( Фортран ) tde02r.zip , tde02d.zip |
|
Текст подпрограммы и версий ( Си ) de02r_c.zip , de02d_c.zip |
Тексты тестовых примеров ( Си ) tde02r_c.zip , tde02d_c.zip |
|
Текст подпрограммы и версий ( Паскаль ) de02r_p.zip , de02e_p.zip |
Тексты тестовых примеров ( Паскаль ) tde02r_p.zip , tde02e_p.zip |
Выполнение одного шага численного интегрирования квазиленийной устойчивой системы обыкновенных дифференциальных уравнений первого порядка с большой константой Липшица методом Лоусона.
Выполняется один шаг численного интегрирования квазилинейной системы 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) ) ,
где 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) является достаточно малым. По заданному значению решения YX в узле xn вычисляется значение этого решения в узле xn + H . Вычисление производится по методу Лоусона.
Метод Лоусона заключается в следующем. Исходная система уравнений с помощью замены искомой функции y (x) на [ xn , 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) .
Значение H может быть меньше или равно значению шага интегрирования, задаваемому пользователем при обращении к подпрограмме.
Все компоненты решения вычисляются с контролем точности по мере погрешности, который заключается в следующем. Если некоторая компонента приближенного решения по абсолютной величине меньше некоторой наперед заданной константы 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 DE02R (FU, M, JSTART, A, HMIN, EPS, P, YX, X, H,
BUL, XP, YP, E1, E2, E3, E4, R, R0, R1, R2, R3,
YD, IERR)
Параметры
| FU - |
подпрограмма вычисления функции U (x, y) в правой
части системы. Первый оператор подпрограммы
должен иметь вид: SUBROUTINE FU (X, Y, U, M). Здесь Y и U - одномерные массивы длины M. В массив U помещается значение функции U (x, y), вычисленное при значении аргументов X и Y (тип параметров X, Y, U: вещественный); |
| M - | количество уравнений в системе (тип: целый); |
| JSTART - | целый указатель режима использования подпрограммы, имеющий следующие значения: |
| 0 - | первое обращение к подпрограмме должно быть выполнено с нулевым значением JSTART; |
| +1 - | выполнить один шаг интегрирования системы дифференциальных уравнений для значений независимой и зависимой переменных и шага интегрирования, заданных параметрами X, YX и H соответственно; при этом само значение величины шага интегрирования H, на который предполагается продолжить решение системы, должно совпадать с тем значением шага, который был фактически выполнен подпрограммой при предыдущем обращении к ней. Если же требуется выполнить шаг интегрирования, отличный от реально сделанного подпрограммой во время предыдущего обращения к ней, то для этого необходимо обратиться к подпрограмме со значением JSTART = 0; |
| +2 - | то же, что и для JSTART = + 1, но только с той разницей, что величина шага интегрирования, на который предполагается продолжить решение, в два раза больше того значения шага, который был фактически выполнен при предыдущем обращении к подпрограмме; |
| -1 - | повторить последний шаг интегрирования с новыми значениями параметров H и/или HMIN; |
| при выходе из подпрограммы параметр JSTART полагается равным 2, если выполнив заданный в H шаг, подпрограмма рекомендовала использовать для интегрирования вдвое большее его значение, и 1 в противном случае, т.е. в том случае, когда рекомендованное значение шага равно только что выполненному шагу; рекомендованное значение шага на выходе из подпрограммы запоминается в переменной H; | |
| A - | вещественный двумерный массив размера M*M, содержащий элементы матрицы A; |
| 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 при входе в нее (т.е. предыдущий узел и решение в нем); |
|
E1, E2 - E3, E4 | вещественные двумерные рабочие массивы размера M*M; |
|
R, R0, R1 - R2, R3, YD | вещественные одномерные рабочие массивы длины M; |
| IERR - | целая переменная, значение которой в результате работы подпрограммы полагается равным 65, если какая - нибудь компонента решения не может быть вычислена с требуемой точностью EPS; в этом случае последний шаг интегрирования системы можно повторить обращением к подпрограмме с новыми значениями параметров H, HMIN и значением JSTART = - 1 . |
Версии
| DE02D - | выполнение одного шага численного интегрирования квазилинейной системы обыкновенных дифференциальных уравнений первого порядка с большой константой Липшица методом Лоусона с удвоенным числом значащих цифр. При этом параметры HMIN, EPS, P, A, YX, X, H, XP, YP, E1, E2, E3, E4, R, R0, R1, R2, R3, YD и параметры X, Y, DY в подпрограмме FU должны иметь тип DOUBLE PRECISION. |
Вызываемые подпрограммы
| AME2R - | подпрограмма вычисления матричной экспоненты; вызывается при работе подпрограммы DE02R; |
| AME2D - | подпрограмма вычисления матричной экспоненты; вызывается при работе подпрограммы DE02D. |
| UTDE20 - | подпрограмма выдачи диагностических сообщений при работе подпрограммы DE02R. |
| UTDE21 - | подпрограмма выдачи диагностических сообщений при работе подпрограммы DE02D. |
| Кроме того, при работе подпрограмм DE02R и DE02D вызываются рабочие подпрограммы DE02RS, DE04RQ и DE02DS, DE04DQ соответственно. |
Замечания по использованию
|
Данная подпрограмма предназначена для интегрирования квазилинейных систем, имеющих малый нелинейный член U (x, y). В общем случае заданная точность не гарантируется. При работе подпрограммы значения параметров M, HMIN, EPS, P, A сохраняются. Между последовательными обращениями к подпрограмме со значениями параметра JSTART = 1, 2 пользователь не должен изменять содержимое массивов A, E1, E2. При работе подпрограммы FU значения параметров X, Y и M не должны изменяться. При обращении к подпрограмме со значением JSTART = - 1 в качестве исходных значений аргумента и решения принимаются значения параметров XP и YP соответственно, т.е. те значения, которые эти параметры получили после самого последнего обращения к подпрограмме с неотрицательным значением JSTART. После работы подпрограммы в массиве R1 содержится значение оценки абсолютной погрешности на шаге, вычисленной по правилу Рунге. Подпрограммы DE02R и DE02D предназначены также для интегрирования жестких дифференциальных уравнений (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
x ≥ 0
Приводятся подпрограмма вычисления функций
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), E1(2, 2), E2(2, 2), E3(2, 2), E4(2, 2), R(2),
* R0(2), R1(2), R2(2), R3(2), YD(2), YP(2), YX(2)
LOGICAL BUL
EXTERNAL FU
M = 2
X = 0.
YX(1) = 5.E1
YX(2) = 5.E1
HMIN = 1.E-10
EPS = 1.E-8
P = 1.E2
JSTART = 0
H = 0.01
A(1, 1) = -5.E2
A(1, 2) = 1.
A(2, 1) = -1.E3
A(2, 2) = 1.
H = 0
10 IH = IH + 1
CALL DE02R (FU, M, JSTART, A, HMIN, EPS, P, YX, X, H, BUL,
* XP, YP, E1, E2, E3, E4, R, R0, R1, R2, R3, YD, IERR)
GO TO (11, 12, 13, 14, 15, 16), IH
11 H = 1.
GO TO 10
12 EPS = 0.1
GO TO 10
13 JSTART = -1
EPS = 1.E-10
GO TO 10
14 EPS = 1.E-4
H = -1.E-5
JSTART = -1
GO TO 10
15 H = -1.E-4
EPS = 1.E-8
GO TO 10
16 CONTINUE
Результаты:
После первого обращения к подпрограмме
X YX(1) YX(2)
1.953125000001 - 05 4.951502369251 + 01 4.902915763238 + 01
H R1(1) R1(2)
1.953125000001 - 05 3.523503740624 - 09 -1.435788969199 - 10
После второго обращения к подпрограмме
X YX(1) YX(2)
5.004882812498 - 05 4.876668035128 + 01 4.753104564862 + 01
H R1(1) R1(2)
3.051757812500 - 05 5.044663945828 - 10 1.331015179553 - 09
После третьего обращения к подпрограмме
X YX(1) YX(2)
8.056640624998 - 05 4.802966483223 + 01 4.605556652445 + 01
H R1(1) R1(2)
6.103515625000 - 05 -2.398155629634 - 08 -9.972912569839 - 10
После четвертого обращения к подпрограмме
X YX(1) YX(2)
5.767822265623 - 05 4.858137986262 + 01 4.716006872657 + 01
H R1(1) R1(2)
7.629394531250 - 06 -4.268561800317 - 11 -3.880510727560 - 12
После пятого обращения к подпрограмме
X YX(1) YX(2)
4.004882812492 - 05 4.901063239173 + 01 4.801943553990 + 01
H R1(1) R1(2)
-2.000000000005 - 05 8.537123600634 - 11 7.761021455121 - 12
После шестого обращения к подпрограмме
X YX(1) YX(2)
-5.995117187529 - 05 4.950228966336 + 01 4.900345166435 + 01
H R1(1) R1(2)
-1.000000000002 - 04 -7.190586378176 - 09 -8.769954244295 - 10