|
Текст подпрограммы и версий ( Фортран ) de04r.zip , de04d.zip |
Тексты тестовых примеров ( Фортран ) tde04r.zip , tde04d.zip |
|
Текст подпрограммы и версий ( Си ) de04r_c.zip , de04d_c.zip |
Тексты тестовых примеров ( Си ) tde04r_c.zip , tde04d_c.zip |
|
Текст подпрограммы и версий ( Паскаль ) de04r_p.zip , de04e_p.zip |
Тексты тестовых примеров ( Паскаль ) tde04r_p.zip , tde04e_p.zip |
Выполнение одного шага численного интегрирования устойчивой системы обыкновенных дифференциальных уравнений первого порядка с большой константой Липшица методом Лоусона.
Выполняется один шаг численного интегрирования линейной устойчивой системы M обыкновенных дифференциальных уравнений первого порядка с переменными коэффициентами
(1) Y ' (X) = A(X) * Y(X) + φ(X) ,
Y = ( y1,..., yM ) ,
A(X) = ( ai j(X) ) , i, j = 1, ..., M ,
φ(X) = ( φ1(X), ... , φM(X) )
Предполагается, что среди характеристических корней матрицы A(X) имеются большие по модулю корни, а функция φ(x) является достаточно малой. По заданному значению решения YX в узле xn вычисляется значение этого решения в узле xn + H . Вычисление производится по методу Лоусона.
Метод Лоусона заключается в следующем. Исходная система уравнений с помощью замены искомой функции Y (x) на [ xn, xn + H ] по формуле
Y(x) = exp [ ( x - xn ) A0 ] Z(x),
где A0 - некоторая постоянная матрица, преобразуется в систему уравнений относительно новой неизвестной функции Z (X):
(2) Z ' (x) = A1(x) Z(x) + φ1(x) = exp [ - ( x - xn ) A0 ] { A(x) - A0 }
exp [ ( x - xn ) A0 ] Z(x) + exp [ ( - ( x - xn ) A0 ] φ(x)
Данное преобразование выполняется самой подпрограммой. В качестве матрицы A0 подпрограмма выбирает матрицу A0 = A (xn + H /2). Если шаг H достаточно мал, то преобразование позволяет уменьшить характеристические корни матрицы A1 (x) по сравнению с характеристическими корнями исходной матрицы A (x). Это приводит к уменьшению константы Липшица системы (2) по сравнению с константой Липшица системы (1). Для решения системы (2) применяются формулы классического метода Рунге - Кутта четвертого порядка точности, причем одновременно с решением (2) производится обратное преобразование от функции Z (x) к функции Y (x).
Значение H может быть меньше или равно значению шага интегрирования, задаваемому пользователем при обращении к подпрограмме. Все компоненты решения вычисляются с контролем точности по мере погрешности, который заключается в следующем. Если некоторая компонента приближенного решения по абсолютной величине не меньше некоторой наперед заданной константы P, то контроль точности для этой компоненты ведется по относительной погрешности, иначе - по абсолютной. Абсолютная погрешность приближенного решения оценивается по правилу Рунге.
J.Douglas Louson. Generalized Runge - Kutta processes for stable systems with large Lipshitz constants // SIAM Journal on Numerical Analisys. - 1967. - Vol 4, No 3.
SUBROUTINE DE04R (FA, FI, M, JSTART, HMIN, EPS, P, YX, X, H,
BUL, XP, YP, A, E1, E2, E3, W1,
R1, R2, R3, R4, R5, R6, YD, IERR)
Параметры
| FA - |
подпрограмма вычисления матрицы системы A(X) в
точке X. Первый оператор подпрограммы должен
иметь вид: SUBROUTINE FA (A, X, M). Здесь A - двумерный массив размера M*M, в котором помещается матрица системы, вычисленная при значении аргумента X (тип параметров A, X: вещественный); |
| FI - |
подпрограмма вычисления неоднородности правой части системы
φ (X) в любой точке X. Первый
оператор подпрограммы должен иметь вид: SUBROUTINE FI (G, X, M). Здесь G - одномерный массив длины M, в который помещается неоднородность правой части системы, вычисленная при значении аргумента X (тип параметров G, X: вещественный); |
| M - | количество уравнений в системе (тип: целый); |
| JSTART - | целый указатель режима использования подпрограммы, имеющий следующие значения: |
| 0,+1 - | выполнить один шаг интегрирования системы дифференциальных уравнений для значений независимой и зависимой переменных и шага интегрирования, заданных параметрами X, YX, и H, соответственно; |
| -1 - | повторить последний шаг интегрирования с новыми значениями параметров H и/или HMIN; на выходе из подпрограммы JSTART полагается равным + 1; |
| 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, W1 | вещественные двумерные рабочие массивы размера M*M; |
|
R1, R2 - R3, R4, R5, R6, YD | вещественные одномерные рабочие массивы длины M; |
| IERR - | целая переменная, значение которой в результате работы подпрограммы полагается равным 65, если какая - нибудь компонента решения не может быть вычислена с требуемой точностью EPS; в этом случае последний шаг интегрирования системы можно повторить обращением к подпрограмме с новыми значениями параметров H, HMIN и значением JSTART = - 1 . |
Версии
| DE04D - | выполнение одного шага численного интегрирования линейной устойчивой системы обыкновенных дифференциальных уравнений первого порядка с большой константой Липшица методом Лоусона с удвоенным числом значащих цифр. При этом параметры HMIN, EPS, P, YX, X, H, XP, YP, A, E1, E2, E3, W1, R1, R2, R3, R4, R5, R6, YD и параметры A, G, X в подпрограммах FA и FI должны иметь тип DOUBLE PRECISION. |
Вызываемые подпрограммы
| AME2R - | подпрограмма вычисления матричной экспоненты; вызывается при работе подпрограммы DE04R. |
| AME2D - | подпрограмма вычисления матричной экспоненты; вызывается при работе подпрограммы DE04D. |
| UTDE20 - | подпрограмма выдачи диагностических сообщений при работе подпрограммы DE04R. |
| UTDE21 - | подпрограмма выдачи диагностических сообщений при работе подпрограммы DE04D. |
|
Кроме того при работе подпрограмм DE04R и DE04D вызываются подпрограммы DE04RD, DE04RM, DE04RP, DE04RQ и DE04DD, DE04DM, DE04DP, DE04DQ соответственно. |
Замечания по использованию
|
Данная подпрограмма предназначена для интегрирования линейных систем, имеющих малую неоднородность φ (x). В общем случае заданная точность не гарантируется. При работе подпрограммы значения параметров M, HMIN, EPS, P сохраняются. При работе подпрограмм FA и FI значения параметров X и M не должны изменяться. При обращении к подпрограмме со значением JSTART = - 1 в качестве исходных значений аргумента и решения принимаются значения параметров XP, YP, соответственно, т.е. те значения, которые эти параметры получили после самого последнего обращения к подпрограмме с неотрицательным значением JSTART. После работы подпрограммы в массиве R1 содержится значение оценки погрешности на шаге, вычисленной по правилу Рунге. Данная подпрограмма и ее версия предназначены также для интегрирования жестких дифференциальных уравнений (1). |
y1' = - ( 2 + x ) y1 /( 1 + x ) + 20 x y2 ,
y1(0) = 0 ,
y2' = -20 x y1 + ( 2 + x ) y2 /( 1 + x ) ,
y2(0) = 1 .
Приводятся подпрограммы вычисления матрицы системы и неоднородной части, фрагмент вызывающей программы, выполняющей несколько шагов из одной точки, а также результаты счета.
SUBROUTINE FA (A, X, M)
DIMENSION A(2, 2)
A(1, 1) = -(2. + X)/(1. + X)
A(1, 2) = 20.*X
A(2, 1) = -A(1, 2)
A(2, 2) = A(1, 1)
RETURN
END
SUBROUTINE FI (R1, X, M)
DIMENSION R1(2)
R1(1) = 0.
R1(2) = 0.
RETURN
END
DIMENSION A(2, 2), E1(2, 2), E2(2, 2), E3(2, 2), W1(2, 2), R1(2),
* R2(2), R3(2), R4(2), R5(2), R6(2), YD(2), YP(2), YX(2)
LOGICAL BUL
EXTERNAL FA, FI
M = 2
X = 0.
YX(1) = 0.
YX(2) = 1.
HMIN = 1.E-10
EPS = 1.E-8
P = 100.
JSTART = 1
H = 0.01
IH = 0
10 IH = IH + 1
CALL DE04R (FA, FI, M, JSTART, HMIN, EPS, P, YX, X, H, BUL,
* XP, YP, A, E1, E2, E3, W1, R1, R2, R3,
* R4, R5, R6, YD, IERR)
GO TO (11, 12, 13, 14, 15, 16, 17), IH
11 EPS = 1.E-10
GO TO 10
12 EPS = 1.0
GO TO 10
13 JSTART = -1
H = 0.01
EPS = 1.E-5
GO TO 10
14 JSTART = -1
H = -0.0001
EPS = 1.E-2
GO TO 10
15 H = 1.
EPS = 1.E-10
GO TO 10
16 EPS = 1.E-4
GO TO 10
17 CONTINUE
Результаты:
После первого обращения к подпрограмме
X YX(1) YX(2)
1.000000000001-02 9.802471967628-04 9.802468700200-01
H R1(1) R1(2)
2.000000000001-02 0.000000000000+00 6.063298011814-14
после второго обращения к подпрограмме
X YX(1) YX(2)
3.000000000000-02 8.479506692396-03 9.421419716109-01
H R1(1) R1(2)
4.000000000002-02 1.231607408649-14 1.333925562599-12
после третьего обращения к подпрограмме
X YX(1) YX(2)
6.999999999994-02 4.268132414620-02 8.703501915807-01
H R1(1) R1(2)
8.000000000004-02 1.970571853839-12 3.565219230945-11
после четвертого обращения к подпрограмме
X YX(1) YX(2)
3.999999999996-02 1.478074532271-02 9.237177506848-01
H R1(1) R1(2)
2.000000000001-02 0.000000000000+00 6.063298011814-14
после пятого обращения к подпрограмме
X YX(1) YX(2)
2.990000000000-02 8.424732657545-03 9.423281849631-01
H R1(1) R1(2)
-2.000000000004-04 0.000000000000+00 6.063298011814-14
после шестого обращения к подпрограмме
X YX(1) YX(2)
6.115000000000-02 3.314040588361-02 8.858545422663-01
H R1(1) R1(2)
3.125000000000-02 4.395891058563-13 1.097456940137-11
после седьмого обращения к подпрограмме
X YX(1) YX(2)
9.240000000000-02 7.117143016819-02 8.315812963056-01
H R1(1) R1(2)
6.250000000000-02 7.80649190203-13 8.852415097246-12