|
Текст подпрограммы и версий ( Фортран ) de95r.zip de95d.zip de91r.zip de91d.zip de93r.zip de93d.zip de97r.zip de97d.zip |
Тексты тестовых примеров ( Фортран ) tde95r.zip tde95d.zip tde91r.zip tde91d.zip tde93r.zip tde93d.zip tde97r.zip tde97d.zip |
|
Текст подпрограммы и версий ( Си ) de95r_c.zip de95d_c.zip de91r_c.zip de91d_c.zip de93r_c.zip de93d_c.zip de97r_c.zip de97d_c.zip |
Тексты тестовых примеров ( Си ) tde95r_c.zip tde95d_c.zip tde91r_c.zip tde91d_c.zip tde93r_c.zip tde93d_c.zip tde97r_c.zip tde97d_c.zip |
|
Текст подпрограммы и версий ( Паскаль ) de95r_p.zip de95e_p.zip de91r_p.zip de91e_p.zip de93r_p.zip de93e_p.zip de97r_p.zip de97e_p.zip |
Тексты тестовых примеров ( Паскаль ) tde95r_p.zip tde95e_p.zip tde91r_p.zip tde91e_p.zip tde93r_p.zip tde93e_p.zip tde97r_p.zip tde97e_p.zip |
Вычисление решения задачи Коши для жесткой системы обыкновенных дифференциальных уравнений первого порядка в конце интервала интегрирования методом типа Розенброка четвертого порядка.
Решается задача Коши для системы M обыкновенных дифференциальных уравнений первого порядка
Y ' = F(X,Y) ,
Y = ( y1,...,yM ) , F = ( f1( X, y1,..., yM ),..., fM( X, y1,..., yM ) )
методом типа Розенброка четвертого порядка.
Решение вычисляется в одной точке XK, которая является концом интервала интегрирования. Для всех компонент решения осуществляется контроль точности по мере погрешности, который заключается в следующем. Если некоторая компонента приближенного решения по абсолютной величине не меньше некоторой наперед заданной константы P (называемой границей перехода), то контроль точности для этой компоненты ведется по относительной погрешности, иначе - по абсолютной.
Артемьев С.С., Демидов Г.В. A - устойчивый метод типа Розенброка четвертого порядка точности решения задачи Коши для жестких систем обыкновенных дифференциальных уравнений// Некоторые проблемы вычислительной и прикладной математики. Новосибирск: "Наука", 1975.
Современные численные методы решения обыкновенных дифференциальных уравнений/ Под ред. Дж.Холла и Дж.Уатта. М.: "Мир", 1979.
SUBROUTINE DE95R (F, FJ, FX, M, XN, YN, XK, HMIN, EPS, P, H, Y,
R, IERR)
Параметры
| F - |
имя подпрограммы вычисления значений правой
части дифференциальных уравнений. Первый оператор
подпрограммы должен иметь вид: SUBROUTINE F (X, Y, DY, M). Здесь X, Y - значения независимой и зависимой переменных, соответственно. Вычисленное значение правой части должно быть помещено в DY; в случае системы уравнений, т.е. когда M ≠ 1, параметры Y и DY представляют одномерные массивы длиной M (тип параметров X, Y и DY: вещественный); |
| FJ - |
имя подпрограммы вычисления значений элементов матрицы Якоби
∂f / ∂y
правой части системы. Первый оператор подпрограммы должен иметь вид: SUBROUTINE FJ (X, Y, Z, M). Здесь X, Y - значения независимой и зависимой переменных, соответственно. В случае системы уравнений, т.е. когда M ≠ 1, параметр Y представляет собой одномерный массив длины M, а параметр Z - двумерный массив размера M*M. Значения элементов матрицы Якоби ∂f / ∂y правой части системы должны быть помещены в массив Z, при этом частная производная от правой части I - го уравнения по J - ой переменной Y (J) запоминается в элементе Z (I,J) (тип параметров X, Y и Z: вещественный); |
| FX - |
имя подпрограммы вычисления значений частных
производных по X правой части системы
∂f / ∂x. Первый
оператор подпрограммы должен иметь вид: SUBROUTINE FX (X, Y, Z1, M). Здесь: X, Y - значения независимой и зависимой переменных, соответственно. В случае системы уравнений, т.е. когда M ≠ 1, параметры Y и Z1 представляют собой одномерные массивы длиной M. Значения частных производных по X правой части системы ∂f / ∂x должны быть помещены в массив Z1; при этом частная производная от правой части I - го уравнения запоминается в элементе Z1 (I) (тип параметров X, Y и Z1: вещественный); |
| M - | количество уравнений в системе (тип: целый); |
| XN, YN - | начальные значения аргумента и решения; в случае системы уравнений (т.е. M ≠ 1) YN представляет одномерный массив длины M (тип: вещественный); |
| XK - | значение аргумента, при котором требуется вычислить решение задачи Коши (конец интервала интегрирования); XK может быть больше, меньше или равно XN (тип: вещественный); |
| HMIN - | минимальное значение абсолютной величины шага, которое разрешается использовать при интегрировании данной системы уравнений (тип: вещественный); |
| EPS - | допустимая мера погрешности, с которой требуется вычислить все компоненты решения (тип: вещественный); |
| P - | граница перехода, используемая при оценке меры погрешности решения (тип: вещественный); |
| H - | вещественная переменная, содержащая начальное значение шага интегрирования; может задаваться с учетом направления интегрирования, т.е. положительным, если XN < XK, отрицательным, если XN > XK, или без всякого учета в виде абсолютной величины; |
| Y - | искомое решение задачи Коши, вычисленное подпрограммой при значении аргумента XK; для системы уравнений (когда M ≠ 1) задается одномерным массивом длины M. В случае совпадения значений параметров XN и XK значение Y полагается равным начальному значению YN (тип: вещественный); |
| R - | одномерный рабочий массив вещественного типа длины 3*M*M + 11*M + 1; |
| IERR - | целая переменная, значение которой в результате работы подпрограммы полагается равным 65, если какая - нибудь компонента решения не может быть вычислена с требуемой точностью EPS; в этом случае интегрирование системы можно повторить обращением к подпрограмме с новыми значениями параметров H и HMIN. |
Версии
| DE97R - |
вычисление решения задачи Коши для жесткой системы обыкновенных дифференциальных уравнений первого порядка в конце интервала интегрирования методом типа Розенброка четвертого порядка. Отличие подпрограммы DE97R от подпрограммы DE95R состоит в том, что при обращении к подпрограмме DE97R не требуется задавать подпрограмму вычисления матрицы Якоби системы и подпрограмму вычисления частных производных правой части системы по x. Все частные производные от правой части вычисляются в подпрограмме DE97R с помощью разностных отношений. Первый оператор подпрограммы DE97R имеет вид: SUBROUTINE DE97R (F, M, XN, YN, XK, HMIN, EPS,
P, H, Y, R, IERR)
Список формальных параметров подпрограммы DE97R
отличается от списка параметров подпрограммы DE95R
отсутствием параметров FJ и FX.
Параметры подпрограммы DE97R имеют тот же смысл, что
и одноименные параметры подпрограммы DE95R.
|
| DE91R - |
вычисление решения задачи Коши для автономной
жесткой системы обыкновенных дифференциальных уравнений
первого порядка в конце интервала интегрирования
методом типа Розенброка четвертого порядка.
Данная подпрограмма предназначена для интегрирования
систем уравнений, в правые части которых не входит
независимая переменная x, т.е. систем вида
y ' = f (y).
Первый оператор подпрограммы DE91R имеет вид:
SUBROUTINE DE91R (F, FJ, M, XN, YN, XK, HMIN, EPS,
P, H, Y, R, IERR)
Список формальных параметров подпрограммы DE91R отличается от списка параметров DE95R отсутствием параметра FX. Параметры подпрограммы DE91R имеют тот же смысл, что и одноименные параметры подпрограммы DE95R, кроме параметра R. В подпрограмме DE91R параметр R представляет одномерный вещественный рабочий массив длины 3M2 + 8M + 1 . |
| DE93R - |
вычисление решения задачи Коши для автономной
жесткой системы обыкновенных дифференциальных уравнений
первого порядка в конце интервала интегрирования
методом Розенброка четвертого порядка. Данная
подпрограмма предназначена для интегрирования систем
уравнений, в правые части которых не входит
независимая переменная x, т.е. систем вида
y ' = f (y).
Отличие подпрограммы DE93R от подпрограммы DE95R состоит
в том, что при обращении к подпрограмме DE93R не
требуется задавать подпрограмму вычисления матрицы Якоби системы
∂f / ∂y
и подпрограмму вычисления частных производных
∂f / ∂x.
Все частные производные от правой
части вычисляются в подпрограмме DE93R с помощью
разностных отношений. Первый оператор подпрограммы DE93R имеет вид:
SUBROUTINE DE93R (F, M, XN, YN, XK, HMIN, EPS,
P, H, Y, R, IERR)
Список формальных параметров подпрограммы DE93R отличается от списка параметров подпрограммы DE95R отсутствием параметров FJ и FX. Параметры подпрограммы DE93R имеют тот же смысл, что и одноименные параметры подпрограммы DE95R, кроме параметра R. В подпрограмме DE93R параметр R представляет одномерный вещественный рабочий массив длины 3M2 + 8M + 1 . |
| DE95D - | вычисление решения задачи Коши для жесткой системы обыкновенных дифференциальных уравнений первого порядка в конце интервала интегрирования методом типа Розенброка четвертого порядка с повышенной точностью. Первый оператор подпрограммы имеет тот же вид, что и в подпрограмме DE95R; при этом параметры XN, YN, XK, HMIN, EPS, P, H, Y, R и параметры X, Y, DY, Z, Z1 в подпрограммах F, FJ и FX должны иметь тип DOUBLE PRECISION. |
| DE97D - | вычисление решения задачи Коши для жесткой системы обыкновенных дифференциальных уравнений первого порядка в конце интервала интегрирования методом типа Розенброка четвертого порядка с повышенной точностью. Первый оператор подпрограммы имеет тот же вид, что и в подпрограмме DE97R; при этом параметры XN, YN, XK, HMIN, EPS, P, H, Y, R и параметры X, Y, DY в подпрограмме F должны иметь тип DOUBLE PRECISION. |
| DE91D - | вычисление решения задачи Коши для автономной жесткой системы обыкновенных дифференциальных уравнений первого порядка в конце интервала интегрирования методом типа Розенброка четвертого порядка с повышенной точностью. Первый оператор подпрограммы имеет тот же вид, что и в подпрограмме DE91R; при этом параметры XN, YN, XK, HMIN, EPS, P, H, Y, R и параметры X, Y, DY, Z в подпрограммах F и FJ должны иметь тип DOUBLE PRECISION. |
| DE93D - | вычисление решения задачи Коши для автономной жесткой системы обыкновенных дифференциальных уравнений первого порядка в конце интервала интергирования методом типа Розенброка четвертого порядка с повышенной точностью. Первый оператор подпрограммы имеет тот же вид, что и в подпрограмме DE93R; при этом параметры XN, YN, XK, HMIN, EPS, P, H, Y, R и параметры X, Y, DY в подпрограмме F должны иметь тип DOUBLE PRECISION. |
Вызываемые подпрограммы
|
DE94R - DE94D | выполнение одного шага численного интегрирования жесткой системы обыкновенных дифференциальных уравнений первого порядка методом типа Розенброка с обыкновенной и повышенной точностью. Вызываются при работе подпрограмм DE95R и DE95D соответственно. |
|
DE96R - DE96D | выполнение одного шага численного интегрирования жесткой системы обыкновенных дифференциальных уравнений первого порядка методом типа Розенброка с обыкновенной и повышенной точностью. Вызываются при работе подпрограмм DE97R и DE97D соответственно. |
|
DE90R - DE90D | выполнение одного шага численного интегрирования автономной жесткой системы обыкновенных дифференциальных уравнений первого порядка методом типа Розенброка с обыкновенной и повышенной точностью. Вызываются при работе подпрограмм DE91R и DE91D соответственно. |
|
DE92R - DE92D | выполнение одного шага численного интегрирования автономной жесткой системы обыкновенных дифференциальных уравнений первого порядка методом типа Розенброка с обыкновенной и повышенной точностью. Вызываются при работе подпрограмм DE93R и DE93D соответственно. |
|
UTDE20 - UTDE21 | подпрограммы выдачи диагностических сообщений. Подпрограмма UTDE20 вызывается при работе подпрограмм DE91R, DE93R, DE95R, DE97R; подпрограмма UTDE21 вызывается при работе подпрограмм DE91D, DE93D, DE95D, DE97D. |
Замечания по использованию
|
В общем случае заданная точность не гарантируется. При работе подпрограммы и ее версий значения параметров M, XN, YN, XK, HMIN, EPS, P сохраняются. При работе подпрограмм F, FJ и FX значения параметров M, X, Y не должны изменяться. Если после работы подпрограммы нет необходимости иметь начальное значение решения YN, то параметры YN и Y при обращении к ней можно совместить. При этом следует иметь в виду, что в случае аварийного выхода из подпрограммы, т.е. со значением IERR = 65, значение параметра YN будет испорчено. В подпрограмме DE95R и ее версиях используется одношаговый метод типа Розенброка четвертого порядка точности. На каждом шаге h интегрирования, выполняемом из текущего узла интегрирования xn, вычисляются четыре значения правой части системы дифференциальных уравнений. Среди этих значений есть одно, которое должно вычисляться при значении независимой переменной x, равном xn - h. Это значение аргумента находится слева от точки xn при h > 0 и справа от точки xn при h < 0. В частности, при выполнении первого шага h из начальной точки xN (напомним, что абсолютная величина первого шаг h берется равной абсолютной величине значения, заданного параметром H при обращении к подпрограмме, а знак первого шага определяется соотношением значений параметров XN и XK) значение xN - h не будет принадлежать интервалу интегрирования, ограниченному значениями, заданными параметрами XN и XK при обращении к подпрограмме DE95R (или ее версиями). Это следует учитывать при составлении подпрограммы F вычисления правой части системы. Если правая часть системы не определена для x, не принадлежащих интервалу интегрирования, то попытка вычислить правую часть для указанного значения аргумента x может привести к аварийному прерыванию. Если при решении системы диффиренциальных уравнений не задаются подпрограммы FJ вычисления матрицы Якоби и FX вычисления частных производных по x правой части системы (т.е. производится обращение к подпрограммам DE93R, DE97R, DE93D, DE97D), то все частные производные по y и по x в этих подпрограммах аппроксимируются центральными разностными отношениями. Замена точных значений производных разностными аппроксимациями может привести к росту погрешности приближенного решения системы дифференциальных уравнений по сравнению с тем, что было бы, если бы все частные производные вычислялись точно с помощью подпрограмм FJ и FX. Даже если будет мала погрешность аппроксимации частных производных (например, если правая часть системы является линейной функцией своих аргументов, то погрешность аппроксимации частных производных равна нулю), может оказаться значительной вычислительная погрешность, особенно, когда вычисления выполняются с одинарной точностью. При этом вычислительная погрешность может даже превосходить верхний предел погрешности приближенного решения, заданный при обращении к подпрограмме параметром EPS. В этом случае целесообразно использовать не подпрограммы DE93R, DE97R, а их версии, выполняющие вычисления с удвоенным числом значащих цифр, т.е. подпрограммы DE93D, DE97D. |
Пример 1.
Решается задача Коши
y1' = -100 y1
y2' = -100 y1 - 2 y2 + 20 e - 100 x + 2 e - x cos x
y3' = -100 y1 + 9998 y2 - 9990 y3 - 10 y4 + 20 e - 100 x + 2 e - x cos x
y4' = -100 y1 + 9988 y2 + 20 y3 - 10010 y4 + 20 e - 100 x + 2 e - x cos x
0 ≤ x ≤ 10 , y1(0) = 10 , y2(0) = 11 , y3(0) = 111 , y4(0) = 111.
Точное решение задачи имеет вид:
y1 = 10 e -100 x
y2 = 10 e - 100 x + e - x cos x + e - x sin x
y3 = y2 + 100 e - 10000 x cos 10x
y4 = y3 + 100 e - 10000 x sin 10x
Матрица Якоби правой части системы имеет вид:
| -100 0 0 0
| -100 -2 0 0
| -100 9998 -9990 -10
| -100 9988 20 -10010
Частные производные по x правой части системы имеют вид:
| 0
| -2000 e - 100 x - 2 e - x (cos x + sin x)
| -2000 e - 100 x - 2 e - x (cos x + sin x)
| -2000 e - 100 x - 2 e - x (cos x + sin x)
Ниже приводятся фрагмент вызывающей программы для DE95R, подпрограммы F, FJ, FX и результаты счета, полученные после нескольких обращений к подпрограмме DE95R.
DIMENSION YN(4), Y(4), R(93)
EXTERNAL F, FJ, FX
M = 4
XN = 0.
YN(1) = 10.
YN(2) = 11.
YN(3) = 111.
YN(4) = 111.
HMIN = 1.E-10
EPS = 1.E-2
P = 1.E3
XK = 10.
C BЫЧИCЛEHИE TOЧHOГO PEШEHИЯ CИCTEMЫ:
Y1 = 10.*EXP(-100.*XK)
Y2 = Y1 + EXP(-XK)*(COS(XK) + SIN(XK))
Y3 = Y2 + 100.*EXP(-1.E4*XK)*COS(10.*XK)
Y4 = Y3 + 100.*EXP(-1.E4*XK)*SIN(10.*XK)
IH = 0
10 IH = IH + 1
H = 0.01
CALL DE95R (F, FJ, FX, M, XN, YN, XK, HMIN, EPS, P, H, Y, R,
* IERR)
EPS = EPS*1.E-2
IF (IH .LT. 4) GO TO 10
SUBROUTINE F (X, Y, DY, M)
DIMENSION Y(4), DY(4)
DY(1) = -100.*Y(1)
T1 = 20.*EXP(-100.*X) + 2.*EXP(-X)*COS(X)
DY(2) = DY(1) - 2.*Y(2) + T1
DY(3) = DY(1) + 9998.*Y(2) - 9990.*Y(3) - 10.*Y(4) + T1
DY(4) = DY(1) + 9988.*Y(2) + 20.*Y(3) - 10010.*Y(4) + T1
RETURN
END
SUBROUTINE FJ (X, Y, DF, M)
DIMENSION Y(4), DF(4, 4)
DF(1, 1) = -100.
DF(1, 2) = 0.
DF(1, 3) = 0.
DF(1, 4) = 0.
DF(2, 1) = -100.
DF(2, 2) = -2.
DF(2, 3) = 0.
DF(2, 4) = 0.
DF(3, 1) = -100.
DF(3, 2) = 9998.
DF(3, 3) = -9990.
DF(3, 4) = -10.
DF(4, 1) = -100.
DF(4, 2) = 9988.
DF(4, 3) = 20.
DF(4, 4) = -10010.
RETURN
END
SUBROUTINE FX (X, Y, DX, M)
DIMENSION Y(4), DX(4)
DX(1) = 0.
DX(2) = -2.E3*EXP(-1.E2*X) - 2.*EXP(-X)*(COS(X) + SIN(X))
DX(3) = DX(2)
DX(4) = DX(2)
RETURN
END
Результаты:
Y1 Y2
0.000000000000+00 -6.279230870976-05
Y3 Y4
-6.279230870976-05 -6.279230870976-05
после первого обращения к подпрограмме -
EPS = 1.0-02
Y(1) Y(2)
7.146873880164-13 -6.764660892966-05
Y(3) Y(4)
-6.764660892122-05 -6.764660892966-05
H = 2.560000000001+00
после второго обращения к подпрограмме -
EPS = 1.0-04
Y(1) Y(2)
6.621227792960-20 -6.330159900469-05
Y(3) Y(4)
-6.330159900392-05 -6.330159900414-05
H = 2.560000000001+00
после третьего обращения к подпрограмме -
EPS = 1.0-06
Y(1) Y(2)
0.000000000000+00 -6.286382905407-05
Y(3) Y(4)
-6.286382905418-05 -6.286382905407-05
H = 1.280000000001+00
после четвертого обращения к подпрограмме -
EPS = 1.0-08
Y(1) Y(2)
0.000000000000+00 -6.279451115976-05
Y(3) Y(4)
-6.279451115976-05 -6.279451116020-05
H = 3.200000000002-01
Пример 2.
Решается система из примера 1. Приводятся фрагмент вызывающей программы для DE97R, подпрограмма F и результаты счета, полученные после нескольких обращений к подпрограмме DE97R.
DIMENSION YN(4), Y(4), R(93)
EXTERNAL F
M = 4
XN = 0.
YN(1) = 10.
YN(2) = 11.
YN(3) = 111.
YN(4) = 111.
HMIN = 1.E-10
EPS = 1.E-2
P = 1.E3
XK = 10.
C BЫЧИCЛEHИE TOЧHOГO PEШEHИЯ CИCTEMЫ:
Y1 = 10.*EXP(-100.*XK)
Y2 = Y1 + EXP(-XK)*(COS(XK) + SIN(XK))
Y3 = Y2 + 100.*EXP(-1.E4*XK)*COS(10.*XK)
Y4 = Y3 + 100.*EXP(-1.E4*XK)*SIN(10.*XK)
IH = 0.
10 IH = IH + 1
H = 0.01
CALL DE97R (F, M, XN, YN, XK, HMIN, EPS, P, H, Y, R, IERR)
EPS = EPS*1.E-2
IF (IH .LT. 3) GOTO 10
SUBROUTINE F (X, Y, DY, M)
DIMENSION Y(4), DY(4)
DY(1) = -100.*Y(1)
T1 = 20.*EXP(-100.*X) + 2.*EXP(-X)*COS(X)
DY(2) = DY(1) - 2.*Y(2) + T1
DY(3) = DY(1) + 9998.*Y(2) - 9990.*Y(3) - 10.*Y(4) + T1
DY(4) = DY(1) + 9988.*Y(2) + 20.*Y(3) - 10010.*Y(4) + T1
RETURN
END
Результаты:
Y1 Y2
0.000000000000+00 -6.279230870976-05
Y3 Y4
-6.279230870976-05 -6.279230870976-05
после первого обращения к подпрограмме -
EPS = 1.0-02
Y(1) Y(2)
7.146874056701-13 -6.764660658287-05
Y(3) Y(4)
-6.764662603498-05 -6.764660998348-05
H = 2.560000000001+00
после второго обращения к подпрограмме -
EPS = 1.0-04
Y(1) Y(2)
6.614556599143-02 -6.330157681200-05
Y(3) Y(4)
-6.330157682632-05 -6.330157680712-05
H = 2.560000000001+00
после третьего обращения к подпрограмме -
EPS = 1.0-06
Y(1) Y(2)
2.733180996095-20 -6.286373127617-05
Y(3) Y(4)
-6.286373127695-05 -6.286373127762-05
H = 1.280000000001+00
Пример 3.
Решается задача Коши
y1' = - 10000 y1 + 100 y2 - 10 y3 + y4
y2' = - 1000 y2 + 10 y3 -10 y4
y3' = - y3 + 10 y4
y4' = - 0.1 y4
0 ≤ x ≤ 20 , y1(0) = 1 , y2(0) = 1 , y3(0) = 1 , y4(0) = 1
Точное решение задачи имеет вид:
y1 = (1 - d - f - g) t3 + d t2 + f t1 + g y4
y2 = (1 - a - b) t2 + a t1 + b y4
y3 = (- 9.1 / 0.9) t1 + (10 / 0.9) y4
y4 = e - 0.1 x
где t1 = e - x , t2 = e - 1000 x , t3 = e - 10000 x ,
a = - 91 / (0.9*999) , b = 91 / (0.9*999.9) , d = ( 1 - a - b ) / 90 ,
f = [ 100 a + (91 / 0.9) ] / 9999 , g = [ 100 b - (100 / 0.9) + 1 ] / 9999.9 .
Матрица Якоби правой части системы имеет вид:
-10000 100 -10 1
0 -1000 10 -10
0 0 -1 10
0 0 0 -0.1
Ниже приводятся фрагмент вызывающей программы для DE91R, подпрограммы F и FJ и результаты счета, полученные после нескольких обращений к подпрограмме DE91R.
DIMENSION YN(4), Y(4), R(81)
EXTERNAL F, FJ
M = 4
XN = 0.
YN(1) = 1.
YN(2) = 1.
YN(3) = 1.
YN(4) = 1.
HMIN = 1.E-10
EPS = 1.E-2
P = 100.
XK = 20.
C BЫЧИCЛEHИE TOЧHOГO PEШEHИЯ CИCTEMЫ
Y4 = EXP(-0.1*XK)
T1 = EXP(-XK)
Y3 = -(9.1/0.9)*T1 + (10./0.9)*Y4
T2 = EXP(-1.E3*XK)
A = -91./(0.9*999.)
B = 91./(0.9*999.9)
C = 1. - A - B
Y2 = C*T2 + A*T1 + B*Y4
T3 = EXP(-1.E4*XK)
DD = C/90.
FF = (100.*A + 91./0.9)/9999.
GG = (100.*B - 100./0.9 + 1.)/9999.9
QQ = 1. - DD - FF- GG
Y1 = QQ*T3 + DD*T2 + FF*T1 + GG*Y4
IH = 0
10 IH = IH + 1
H = 0.01
CALL DE91R (F, FJ, M, XN, YN, XK, HMIN, EPS, P, H, Y, R, IERR)
EPS = EPS*1.E-2
IF (IH .LT. 4) GOTO 10
SUBROUTINE F (X, Y, DY, M)
DIMENSION Y(4), DY(4)
DY(1) = -1.E4*Y(1) + 1.E2*Y(2) - 10.*Y(3) + Y(4)
DY(2) = -1.E3*Y(2) + 10.*Y(3) - 10.*Y(4)
DY(3) = -Y(3) + 10.*Y(4)
DY(4) = -0.1*Y(4)
RETURN
END
SUBROUTINE FJ (X, Y, DF, M)
DIMENSION Y(4), DF(4, 4)
DF(1, 1) = -1.E4
DF(1, 2) = 1.E2
DF(1, 3) = -10.
DF(1, 4) = 1.
DF(2, 1) = 0.
DF(2, 2) = -1.E3
DF(2, 3) = 10.
DF(2, 4) = -10.
DF(3, 1) = 0.
DF(3, 2) = 0.
DF(3, 3) = -1.
DF(3, 4) = 10.
DF(4, 1) = 0.
DF(4, 2) = 0.
DF(4, 3) = 0.
DF(4, 4) = -0.1
RETURN
END
Результаты:
Y1 Y2
-1.353352661869-03 1.368526917891-02
Y3 Y4
1.503725348452+00 1.353352832364-01
после первого обращения к подпрограмме -
EPS = 1.0-02
Y(1) Y(2)
-1.352902557240-03 1.368071768447-02
Y(3) Y(4)
1.503225232165+00 1.352902708961-01
H = 1.024000000001+01
после второго обращения к подпрограмме -
EPS = 1.0-04
Y(1) Y(2)
-1.353311838306-03 1.368485637501-02
Y(3) Y(4)
1.503679988920+00 1.353311999785-01
H = 2.560000000001+00
после третьего обращения к подпрограмме -
EPS = 1.0-06
Y(1) Y(2)
-1.353350127546-03 1.368524355264-02
Y(3) Y(4)
1.503722532527+00 1.353350296827-01
H = 1.280000000001+00
после четвертого обращения к подпрограмме -
EPS = 1.0-08
Y(1) Y(2)
-1.353352646026-03 1.368526901869-02
Y(3) Y(4)
1.503725330851+00 1.353352816498-01
H = 3.200000000002-01
Пример 4.
Решается задача Коши из примера 3. Приводятся фрагмент вызывающей программы для DE93R, подпрограмма F и результаты счета, полученные после нескольких обращений к подпрограмме DE93R.
DIMENSION YN(4), Y(4), R(81)
EXTERNAL F
M = 4
XN = 0.
YN(1) = 1.
YN(2) = 1.
YN(3) = 1.
YN(4) = 1.
HMIN = 1.E-10
EPS = 1.E-8
P = 100
XK = 20.
C BЫЧИCЛEHИE TOЧHOГO PEШEHИЯ CИCTEMЫ:
Y4 = EXP(-0.1*XK)
T1 = EXP(-XK)
Y3 = -(9.1/0.9)*T1 + (10./0.9)*Y4
T2 = EXP(-1.E3*XK)
A = -91./(0.9*999.)
B = 91./(0.9*999.9)
C = 1. - A - B
Y2 = C*T2 + A*T1 + B*Y4
T3 = EXP(-1.E4*XK)
DD = C/90.
FF = (100.*A + 91./0.9)/9999.
GG = (100.*B - 100./0.9 + 1.)/9999.9
QQ = 1. - DD - FF - GG
Y1 = QQ*T3 + DD*T2 + FF*T1 + GG*Y4
IH = 0
10 IH = IH + 1
H = 0.01
CALL DE93R (F, M, XN, YN, XK, HMIN, EPS, P, H, Y, R, IERR)
EPS = EPS*1.E-2
IF (IH .LT. 4) GO TO 10
SUBROUTINE F (X, Y, DY, M)
DIMENSION Y(4), DY(4)
DY(1) = -1.E4*Y(1) + 1.E2*Y(2) - 10.*Y(3) + Y(4)
DY(2) = -1.E3*Y(2) + 10.*Y(3) - 10.*Y(4)
DY(3) = -Y(3) + 10.*Y(4)
DY(4) = -0.1*Y(4)
RETURN
END
Результаты:
Y1 Y2
-1.353352661869-03 1.368526917891-02
Y3 Y4
1.503725348452+00 1.353352832364-01
после первого обращения к подпрограмме -
EPS = 1.0-02
Y(1) Y(2)
-1.352902572352-03 1.368071785960-02
Y(3) Y(4)
1.503225234967+00 1.352902709414-01
H = 1.024000000001+01
после второго обращения к подпрограмме -
EPS = 1.0-04
Y(1) Y(2)
-1.353311843749-03 1.368485643758-02
Y(3) Y(4)
1.503679988458+00 1.353311999721-01
H = 2.560000000001+00
после третьего обращения к подпрограмме -
EPS = 1.0-06
Y(1) Y(2)
-1.353350112996-03 1.368524337161-02
Y(3) Y(4)
1.503722531763+00 1.353350295963-01
H = 1.280000000001+00
после четвертого обращения к подпрограмме -
EPS = 1.0-08
Y(1) Y(2)
-1.353351314448-03 1.368525627954-02
Y(3) Y(4)
1.503725330012+00 1.353352828808-01
H = 3.200000000002-01
Пример 5.
Решается задача Коши
y1' = - a2 y1 y2
y2' = 2 a1 y5 x92 - a2 y1 y2 - a5 x6 y2 + a6 x7
y3' = a3 x6 x8 - a7 y3
y4' = a6 x7 + a7 y3 + a4 x6
y5' = - a1 y5 x92
где
x6 = 1 - a8 - 2 y1 + y2 - y3 - y4 + 2y5
x7 = a8 + y1 - y2 - 2 y5
x8 = a9 - y3
x9 = - 0.8745 - a8 + y3 + y4 + 2 y5
a1 = 80, a2 = 29, a3 = 1, a4 = 0.288*10-3 , a5 = 166 ,
a6 = 0.959 * 10-4 , a7 = 0.232 * 10-3 , a8 = 0.0477 , a9 = 0.602 ,
0 ≤ x ≤ 3000 .
y1(0) = 1 , y2(0) = 0.0477 , y3(0) = 0 , y4(0) = 0 , y5(0) = 0.5
Этот пример взят из статьи С.С.Артемьева и Г.В.Демидова, которая приведена в списке литературы, данном в разделе "Математическое описание".
Ниже приводятся фрагмент вызывающей программы для DE91R, подпрограммы F и FJ и результаты счета, полученные после двух обращений к подпрограмме DE91R со значениями параметра EPS, равными 10-6 и 10-7.
DIMENSION YN(5), Y(5), R(116)
EXTERNAL F, FJ
M = 5
XN = 0.
YN(1) = 1.
YN(2) = 0.0477
YN(3) = 0.
YN(4) = 0.
YN(5) = 0.5
HMIN = 1.E-10
EPS = 1.E-6
P = 100.
XK = 3.E3
IH = 0
10 IH = IH + 1
H = 0.01
CALL DE91R (F, FJ, M, XN, YN, XK, HMIN, EPS, P, H, Y, R, IERR)
EPS = EPS*1.E-1
IF (IH .LT. 2) GO TO 10
SUBROUTINE F (X, Y, DY, M)
DIMENSION Y(5), DY(5)
A1 = 80.
A2 = 29.
A3 = 1.
A4 = 0.288E-3
A5 = 166.
A6 = 0.959E-4
A7 = 0.232E-3
A8 = 0.0477
A9 = 0.602
X6 = 1. - A8 - 2.*Y(1) + Y(2) - Y(3) - Y(4) + 2.*Y(5)
X7 = A8 + Y(1) - Y(2) - 2*Y(5)
X8 = A9 - Y(3)
X9 = -0.8745 - A8 + Y(3) + Y(4) + 2.*Y(5)
T1 = A1*Y(5)*X9*X9
T2 = A6*X7
T3 = A7*Y(3)
DY(1) = -A2*Y(1)*Y(2)
DY(2) = 2.*T1 + DY(1) - A5*X6*Y(2) + T2
DY(3) = A3*X6*X8 - T3
DY(4) = T2 + T3 + A4*X6
DY(5) = -T1
RETURN
END
SUBROUTINE FJ (X, Y, DF, M)
DIMENSION Y(5), DF(5, 5)
A1 = 80.
A2 = 29.
A3 = 1.
A4 = 0.288E-3
A5 = 166.
A6 = 0.959E-4
A7 = 0.232E-3
A8 = 0.0477
A9 = 0.602
X6 = 1. - A8 - 2.*Y(1) + Y(2) - Y(3) - Y(4) + 2.*Y(5)
X8 = A9 - Y(3)
X9 = -0.8745 - A8 + Y(3) + Y(4) + 2.*Y(5)
T1 = 2.*A1*Y(5)*X9
T2 = A5*Y(2)
T3 = A1*X9*X9
DF(1, 1) = -A2*Y(2)
DF(1, 2) = -A2*Y(1)
DF(1, 3) = 0.
DF(1, 4) = 0.
DF(1, 5) = 0.
DF(2, 1) = 2.*T2 + DF(1, 1) + A6
DF(2, 2) = DF(1, 2) - A5*X6 - T2 - A6
DF(2, 3) = 2.*T1 + T2
DF(2, 4) = DF(2, 3)
DF(2, 5) = 2.*T3 + 4.*T1 - 2.*T2 - 2.*A6
DF(3, 2) = A3*X8
DF(3, 1) = -2.*DF(3, 2)
DF(3, 3) = -DF(3, 2) - A3*X6 - A7
DF(3, 4) = -DF(3, 2)
DF(3, 5) = 2.*DF(3, 2)
DF(4, 1) = A6 - 2.*A4
DF(4, 2) = -A6 + A4
DF(4, 3) = A7 - A4
DF(4, 4) = -A4
DF(4, 5) = 2.*(A4 - A6)
DF(5, 1) = 0.
DF(5, 2) = 0.
DF(5, 3) = -T1
DF(5, 4) = DF(5, 3)
DF(5, 5) = -T3 + 2.*DF(5, 3)
RETURN
END
Результаты:
после первого обращения к подпрограмме -
EPS = 1.0-06
Y(1) Y(2) Y(3)
7.743161591600-02 3.837857646993-05 5.035858808897-01
Y(4) Y(5)
3.578708481186-01 3.240421442632-02
H = 1.638400000001+02
после второго обращения к подпрограмме -
EPS = 1.0-07
Y(1) Y(2) Y(3)
7.743302086658-02 3.837918184429-05 5.035845024640-01
Y(4) Y(5)
3.578709420071-01 3.240497334838-02
H = 8.192000000004+01
Пример 6.
Решается задача Коши из примера 5. Приводятся фрагмент вызвающей программы для DE93R, подпрограмма F и результаты счета, полученные после двух обращений к подпрограмме DE93R со значениями параметра EPS, равными 10-6 и 10-7.
DIMENSION YN(5), Y(5), R(116)
EXTERNAL F
M = 5
XN = 0.
YN(1) = 1.
YN(2) = 0.0477
YN(3) = 0.
YN(4) = 0.
YN(5) = 0.5
HMIN = 1.E-10
EPS = 1.E-6
P = 100.
XK = 3.E3
IH = 0
10 IH = IH + 1
H = 0.01
CALL DE93R (F, M, XN, YN, XK, HMIN, EPS, P, H, Y, R, IERR)
EPS = EPS*1.E-1
IF (IH .LT. 2) GO TO 10
SUBROUTINE F (X, Y, DY, M)
DIMENSION Y(5), DY(5)
A1 = 80.
A2 = 29.
A3 = 1.
A4 = 0.288E-3
A5 = 166.
A6 = 0.959E-4
A7 = 0.232E-3
A8 = 0.0477
A9 = 0.602
X6 = 1. - A8 - 2.*Y(1) + Y(2) - Y(3) - Y(4) + 2.*Y(5)
X7 = A8 + Y(1) - Y(2) - 2*Y(5)
X8 = A9 - Y(3)
X9 = -0.8745 - A8 + Y(3) + Y(4) + 2.*Y(5)
T1 = A1*Y(5)*X9*X9
T2 = A6*X7
T3 = A7*Y(3)
DY(1) = -A2*Y(1)*Y(2)
DY(2) = 2.*T1 + DY(1) - A5*X6*Y(2) + T2
DY(3) = A3*X6*X8 - T3
DY(4) = T2 + T3 + A4*X6
DY(5) = -T1
RETURN
END
Результаты:
после первого обращения к подпрограмме -
EPS = 1.0-06
Y(1) Y(2) Y(3)
7.743161056910-02 3.837859005229-05 5.035858775518-01
Y(4) Y(5)
3.578708481446-01 3.240421646177-02
H = 1.638400000001+02
после второго обращения к подпрограмме -
EPS = 1.0-07
Y(1) Y(2) Y(3)
7.743301273331-02 3.837911001064-05 5.035844456870-01
Y(4) Y(5)
3.578709422377-01 3.240497096766-02
H = 8.192000000004+01