|
Текст подпрограммы и версий ( Фортран ) de15r.zip , de15d.zip |
Тексты тестовых примеров ( Фортран ) tde15r.zip , tde15d.zip |
|
Текст подпрограммы и версий ( Си ) de15r_c.zip , de15d_c.zip |
Тексты тестовых примеров ( Си ) tde15r_c.zip , tde15d_c.zip |
|
Текст подпрограммы и версий ( Паскаль ) de15r_p.zip , de15e_p.zip |
Тексты тестовых примеров ( Паскаль ) tde15r_p.zip , tde15e_p.zip |
Выполнение одного шага численного интегрирования системы обыкновенных дифференциальных уравнений первого порядка методом Рунге - Кутта - Фельберга.
Выполняется один шаг численного интегрирования системы М обыкновенных дифференциальных уравнений первого порядка
Y ' = F (X, Y) ,
Y = ( y1, ... , yM ) ,
F = ( f1 (X, y1, ... , yM), ... , fM (X, y1, ... , yM) )
методом Рунге - Кутта - Фельберга пятого порядка точности. По заданному значению решения YX в узле Xn вычисляется значение этого решения в узле Xn + H. Bсе компоненты решения вычисляются с контролем точности по меpе погрешности, который заключается в следующем.
Если некоторая компонента приближенного решения по абсолютной величине не меньше некоторой наперед заданной константы Р, то контроль точности для этой компоненты ведется по относительной погрешности, иначе - по абсолютной. B качестве абсолютной погрешности решения принимается оценка главного члена асимптотического разложения погрешности метода на шаге, получаемая при вычитании двух выражений, представляющих приближенные значения решения пятого и четвертого порядка точности. При этом на одном шаге интегрирования для определения решения и погрешности используются всего шесть вычислений правой части системы. Значение H может быть меньше или равно значению шага интегрирования, задаваемому пользователем при обращении к подпрограмме.
Дж.Форсайт, М.Малькольм, К.Моулер, Машинные методы математических вычислений, "Мир", M., 1980.
SUBROUTINE DE15R (F, M, JSTART, HMIN, EPS, P, YX, X, H,
BUL, XP, YP, DY, R1, R2, R3, R4, IERR)
Параметры
| F - |
имя подпрограммы вычисления значений правой
части дифференциального уравнения. Первый
оператоp подпрограммы должен иметь вид: SUBROUTINE F (X, Y , DY, M). Здесь: X, Y - значения независимой и зависимой переменных, соответственно. Вычисленное значение правой части должно быть помещено в DY. B случае системы уравнений, т.е. когда M ≠ 1 , параметры Y и DY представляют массивы длины M (тип параметров X, Y и DY: вещественный); |
| M - | количество уравнений в системе (тип: целый); |
| JSTART - | целый указатель режима использования подпрограммы, имеющий следующие значения: |
| 0,+1 - | выполнить один шаг интегрирования системы дифференциальных уравнений для значений независимой и зависимой переменных и шага интегрирования, заданных параметрами X, YX и H, соответственно; |
| -1 - | повторить последний шаг интегрирования с новыми значениями паpаметpов H и/или HMIN; |
| на выходе из подпрограммы JSTART полагается равным + 1; | |
| HMIN - | минимальное значение абсолютной величины шага, котоpое разрешается использовать при интегрировании данной системы уравнений (тип: вещественный); |
| EPS - | допустимая меpа погрешности, с которой тpебуется вычислить все компоненты решения (тип: вещественный); |
| P - | граница перехода, используемая при оценке меры погрешности решения (тип: вещественный); |
| YX, X - | заданные вещественные значения решения и соответствующее ему значение аргумента; в pезультате работы подпрограммы в X получается новое значение аргумента, а в YX - соответствующее значение решения; в случае системы уравнений, т.е. когда M ≠ 1, YX задается одномерным массивом длины M; |
| H - | вещественная переменная, содержащая значение шага интегрирования; если для этого значения шага точность приближенного решения достигается, то именно он и реализуется подпрограммой, иначе этот шаг уменьшается подпрограммой до тех пор, пока не будет достигнута заданная точность EPS; на выходе из подпрограммы H содержит рекомендуемое подпрограммой значение следущего шага интегрирования, определяемое ею с целью достижения более экономного способа интегрирования; |
| BUL - | логическая переменная, значение которой при обращении к подпрограмме полагается равным .TRUE., если заданный в H шаг выводит в конец интервала интегрирования, и .FALSE. в противном случае; в результате работы подпрограммы BUL pавно .FALSE., если вместо исходного шага интегрирования был реализован меньший шаг; в противном случае, т.е. когда был выполнен именно заданный при обращении в H шаг, значение параметра BUL не меняется; |
| XP, YP - | вещественные рабочая переменная и одномерный рабочий массив длины M; значения параметpов XP, YP на выходе из подпрограммы pавны тем значениям, которые имели параметры X, YX, соответственно, при входе в нее (т.е. предыдущий узел и решение в нем); |
|
DY - R1, R2 R3, R4 | одномерные вещественные рабочие массивы длины M; |
| IERR - | целая переменная, значение которой в pезультате работы подпрограммы полагается равным 65, если какая - нибудь компонента pешения не может быть вычислена с требуемой точностью EPS; в этом случае последний шаг интегрирования системы можно повторить обращением к подпрограмме с новыми значениями параметров H, HMIN и значением JSTART = - 1 . |
Версии
| DE15D - | выполнение одного шага численного интегрирования системы обыкновенных дифференциальных уpавнений первого порядка методом Рунге - Кутта - Фельберга с удвоенным числом значащих цифр. При этом параметры HMIN, EPS, P, YX, X, H, XP, YP, DY, R1, R2, R3, R4 и параметры X, Y и DY в подпрограмме F должны иметь тип DOUBLE PRECISION. |
Вызываемые подпрограммы
| UTDE16 - | подпрограмма выдачи диагностических сообщений при работе подпрограммы DE15R. |
| UTDE17 - | подпрограмма выдачи диагностических сообщений при работе подпрограммы DE15D. |
Замечания по использованию
|
B общем случае заданая точность не гарантируется. При работе подпрограммы значения параметров M, HMIN, EPS, P сохраняются. При работе подпрограммы счета правой части F значения параметров X, Y и M не должны изменяться. При обращении к подпрограмме со значением JSTART = - 1 в качестве исходных значений аргумента и решения принимаются значения параметров XP и YP, соответственно, т.е. те значения, которые эти параметры получили после самого последнего обращения к подпрограмме с неотрицательным значением JSTART. |
y1' = -y2 - 0.1 x - 0.9
y2' = y1 - 0.1 x - 1.1
y1 = 1 , y2 (0) = -2 , 0 ≤ x ≤ π
Точное решение системы:
y1 = 0.1 x + sin x + 1. , y2 = -0.1 x - cos x - 1.
Приводятся подпрограмма вычисления значений правой части системы и фрагмент вызывающей программы, а также результаты счета.
SUBROUTINE FFELB (X, Y, Z, M)
DIMENSION Y(2), Z(2)
R = 0.1*X
Z(1) = -Y(2) - R - 0.9
Z(2) = Y(1) - R - 1.1
RETURN
END
DIMENSION YX(2), YP(2), DY(2), R1(2), R(2),
* R3(2), R4(2)
LOGICAL BUL
EXTERNAL FFELB
M = 2
JSTART = 0
HMIN = 1.E-12
EPS = 1.E-5
P = 100.
YX(1) = 1.
YX(2) = -2.
X = 0.
H = 0.01
BUL = .FALSE.
IH = 0
6 CALL DE15R (FFELB, M, JSTART, HMIN, EPS, P, YX, X, H, BUL,
* XP, YP, DY, R1, R2, R3, R4, IERR)
C BЫЧИCЛEHИE TOЧHOГO ЗHAЧEHИЯ PEШEHИЯ
Y1 = 0.1*X + SIN(X) + 1.
Y2 = -0.1*X - COS(X) - 1.
IH = IH + 1
PRINT 1, X, YX, Y1, Y2, H
GO TO (6, 6, 13, 14, 15, 6, 20), IH
13 JSTART = -1
H = -0.005
GO TO 6
14 JSTART = -1
H = 0.02
GO TO 6
15 JSTART = -1
H = 0.01
GO TO 6
20 STOP
Результаты:
после первого обращения к подпрограмме -
X YX (1) YX (2)
9.999999999991-03 1.010999833334 + 00 -2.000950000420 + 00
H Y1 Y2
1.999999999998-02 1.010999833334 + 00 -2.000950000416 + 00
после второго обращения к подпрограмме -
X YX (1) YX (2)
2.999999999997-02 1.032995500202 + 00 -2.002550033754 + 00
H Y1 Y2
3.999999999996-02 1.032995500202 + 00 -2.002550033747 + 00
после третьего обращения к подпрограмме -
X YX (1) YX (2)
6.999999999994-02 1.076942847336 + 00 -2.004551000264 + 00
H Y1 Y2
7.999999999993-02 1.076942847336 + 00 -2.004551000253 + 00
после четвертого обращения к подпрограмме -
X YX (1) YX (2)
2.499999999998-02 1.027497395913 + 00 -2.002187516282 + 00
H Y1 Y2
-9.999999999991-03 1.027497395913 + 00 -2.002187516275 + 00
после пятого обращения к подпрограмме -
X YX (1) YX (2)
4.999999999995-02 1.054979169268 + 00 -2.003750260403 + 00
H Y1 Y2
3.999999999996-02 1.054979169270 + 00 -2.003750260392 + 00
после шестого обращения к подпрограмме -
X YX (1) YX (2)
3.999999999996-02 1.043989334185 + 00 -2.003200106668 + 00
H Y1 Y2
1.999999999998-02 1.043989334186 + 00 -2.003200106661 + 00
после седьмого обращения к подпрограмме -
X YX (1) YX (2)
5.999999999995-02 1.065964006477 + 00 -2.004200539945 + 00
H Y1 Y2
3.999999999996-02 1.065964006479 + 00 -2.004200539934 + 00