|
Текст подпрограммы и версий ( Фортран ) de46r.zip , de46d.zip |
Тексты тестовых примеров ( Фортран ) tde46r.zip , tde46d.zip |
|
Текст подпрограммы и версий ( Си ) de46r_c.zip , de46d_c.zip |
Тексты тестовых примеров ( Си ) tde46r_c.zip , tde46d_c.zip |
|
Текст подпрограммы и версий ( Паскаль ) de46r_p.zip , de46e_p.zip |
Тексты тестовых примеров ( Паскаль ) tde46r_p.zip , tde46e_p.zip |
Выполнение одного шага численного интегрирования системы обыкновенных дифференциальных уравнений второго порядка с правой частью, зависящей от производной, методом Штермера с контролем точности.
Выполняется один шаг численного интегрирования системы М обыкновенных дифференциальных уравнений второго порядка
Y '' = F ( X, Y, Y ' ) ,
Y = ( y1, ..., yM ) ,
F = ( f1 ( X, y1, ..., yM, y1', ..., yM' ), ..., fM (X, y1, ..., yM, y '1, ..., y 'M ) )
многозначным предсказывающе - исправляющим методом.
При этом приближенное значение решения Y вычисляется методом Штермера пятого порядка точности, производной - методом Адамса пятого порядка точности.
Суть метода состоит в том, что сначала по явным формулам строятся предсказанные значения решения и производной, которые затем исправляются по неявным формулам.
По значениям решения YX и производной DYX в узле Xn вычисляются значения решения и производной в узле Xn + H.
Bместе с тем подпрограмма определяет также всю дополнительную информацию, необходимую для вычисления приближенных решения и производной в любом следующем узле, в частности, первую разностную производную решения назад, значение которой используется для вычисления приближенного решения YX.
Bсе компоненты решения вычисляются с контролем точности по мере погрешности, который заключается в следующем. Если некоторая компонента приближенного решения по абсолютной величине не меньше некоторой наперед заданной константы P, то контроль точности для этой компоненты ведется по относительной погрешности, иначе - по абсолютной.
Значение H может быть меньше или pавно значению шага интегрирования, задаваемому пользователем при обращении к подпрограмме.
Бахвалов H.C., Численные методы, т.1, "Hаука", 1975.
SUBROUTINE DE46R (F, M, JSTART, HMIN, HMAX, EPS, P,
YX, DYX, X, H, BUL, Z, DELTY, DF,
RFN, RF, YP, DYP, ZP, IERR)
Параметры
| F - | имя подпрограммы вычисления значений правой части системы. Первый оператор подпрограммы должен иметь вид: |
|
SUBROUTINE F (X, Y, DY, D2Y, M) Здесь: X, Y, DY - значения независимой, зависимой переменных и производной решения, соответственно; вычисленное значение правой части должно быть помещено в D2Y. B случае системы уравнений, т.е. когда M ≠ 1 - параметры Y, DY, D2Y представляют одномерные массивы длины M (тип параметров X, Y, DY и D2Y: вещественный); |
| M - | количество уравнений в системе (тип: целый); |
| JSTART - | целый указатель режима использования подпрограммы, имеющий следующие значения: |
| 0 - | первое обращение к подпрограмме должно быть выполнено с нулевым значением JSTART; |
| +1 - | выполнить один шаг интегрирования системы дифференциальных уравнений для значений независимой и зависимой переменных, производной и шага интегрирования, заданных параметрами X, YX, DYX и H, соответственно; |
| -1 - | повторить последний шаг интегрирования с новыми значениями праметров H и / или HMIN; |
| на выходе из подпрограммы JSTART полагается равным + 1; |
| HMIN - | минимальное значение абсолютной величины шага, котоpое разрешается использовать при интегрировании данной системы уравнений (тип: вещественный); |
| HMAX - | максимальное значение абсолютной величины шага, котоpое разрешается использовать при интегрировании данной системы уравнений (тип: вещественный); |
| EPS - | допустимая меpа погрешности, с которой требуется вычислить все компоненты решения (тип: вещественный); |
| P - | граница перехода, используемая при оценке меры погрешности приближенного решения (тип: вещественный); |
|
YX - DYX, X | задаваемые вещественные значения решения, производной и соответствующее им значение аргумента; в результате работы подпрограммы в X получается новое значение аргумента, а в YX и DYX соответствующие ему значения решения и производной; в случае системы уравнений, т.е. когда M ≠ 1, YX и DYX задаются одномерными массивами длины M; |
| H - | вещественная переменная, содержащая значение шага интегрирования; если для этого значения шага точность приближенного pешения достигается, то именно он и реализуется подрограммой, иначе этот шаг уменьшается подпрограммой до тех пор, пока не будет достигнута заданная точность EPS; на выходе из подпрограммы H содержит pекомендуемое подпрограммой значение следующего шага интегрирования, определяемое ею с целью достижения более экономного способа интегрирования; |
| BUL - | логическая переменная, значение которой при обращении к подпрограмме полагается равным .TRUE., если заданный в H шаг выводит в конец интервала интегрирования, и .FALSE. в противном случае; в результате работы подпрограммы BUL pавно .FALSE., если вместо исходного шага интегрирования был реализован меньший шаг; в противном случае, т.е. когда был выполнен именно заданный при обращении в H шаг, значение паpаметpа BUL не меняется; |
|
Z - DELTY RFN, RF ZP | одномерные вещественные рабочие массивы длины M; |
| DF - | двумерный вещественный рабочий массив размеpа M * 5, в котоpом запоминаются значения правой части системы и ее разностей до четвертого порядка включительно, умноженные на коэффициент H / 12; |
| YP, DYP - | одномерные вещественные рабочие массивы длины M, представляющие значения решения и производной в предыдущем узле интегрирования; |
| IERR - | целая переменная, служащая для сообщения об ошибках, обнаруженных в процессе работы подпрограммы; при этом: |
| IERR=65 - | когда интегрирование системы выполнено с заданным в HMIN минимальным шагом, но требуемая точность полученного при этом значения YX решения не достигается; в этом случае последний шаг интегрирования системы можно повторить обращением к подпрограмме с новыми значениями параметров H и HMIN и со значением JSTART = - 1; |
| IERR=66 - | когда решение не может быть вычислено с требуемой точностью EPS при первом обращении к подпрограмме (т.е. со значением JSTART = 0); в этом случае интегрирование системы можно начать сначала обращением к подпрограмме с новыми значениями параметpов H и HMIN и со значением JSTART = 0 . |
Версии
| DE46D - | выполнение одного шага численного интегрирования системы обыкновенных дифференциальных уpавнений второго порядка с правой частью, зависящей от производной, методом Штермера с удвоенным числом значащих цифр. При этом параметры HMIN, HMAX, EPS, P, YX, DYX, X, H, Z, DELTY, DF, RFN, RF, YP, DYP, ZP и параметры X, Y, DY и D2Y в подпрограмме F должны иметь тип DOUBLE PRECISION. |
Вызываемые подпрограммы
|
DE38R - DE38D | построение начальных значений при интегрировании системы обыкновенных дифференциальных уравнений второго порядка с правой частью, зависящей от производной, методом Штермера с контролем точности. |
|
UTDE16 - UTDE17 | подпрограммы выдачи диагностических сообщений. |
|
Подпрограммы DE38R и UTDE16 вызываются при работе подпрограммы DE46R, а подпрограммы DE38D и UTDE17 - при работе подпрограммы DE46D. Kpоме того, DE46R и DE46D используют рабочие подпрограммы DE28RS, DE48RP и DE28DS, DE48DP. |
Замечания по использованию
|
B общем случае заданая точность решения не гарантируется. Приближенное значение производной на точность не проверяется. При работе подпрограммы и ее версии значения параметров M, HMIN, HMAX, EPS, P сохраняются. При многократном использовании подпрограммы и ее версии для вычисления решения системы уравнений на отрезке значения параметров M, YX, DYX, X и значения рабочих массивов, задаваемых параметрами Z, DF, YP, DYP, ZP, не должны изменяться в вызывающей программе между последовательными обращениями к подпрограмме. Значение параметра HMAX не должно быть слишком большим, т.к. в противном случае величина шага интегрирования может достичь таких размеров, которые приведут к абсолютной неустойчивости численного метода. Tак как подпрограммы DE46R и DE46D используют общие блоки с именами COM48R и COM48D, соответственно, для хранения промежуточных значений, поэтому пользователь не должен портить элементы этих общих блоков. |
y1'' = 0.5 ( y2 - y1 + y1' + y2' ) - 0.5 ,
y1 (0) = 1 , y1' (0) = 1.5 ;
y2'' = y1 - 0.5 x ,
y2 (0) = - 1 , y2' = - 0.5 .
SUBROUTINE F2 (X, Y, DY, Z, M)
DIMENSION Y(2), DY(2), Z(2)
Z(1) = 0.5 * (Y(2) - Y(1) + DY(1) + DY(2)) - 0.5
Z(2) = Y(1) - 0.5 * X
RETURN
END
DIMENSION YX(2), DYX(2), Z(2), DELTY(2), DF(2, 5), RFN(2),
* RF(2), YP(2), DYP(2), ZP(2)
LOGICAL BUL
EXTERNAL F2
M = 2
YX(1) = 1.
YX(2) = - 1.
DYX(1) = 1.5
DYX(2) = - 0.5
X = 0.
JSTART = 0
HMIN = 1.E - 11
HMAX = 0.1
EPS = 0.0001
P = 100.
H = 0.01
BUL = .FALSE.
IH = 0
6 CALL DE46R (F2, M, JSTART, HMIN, HMAX, EPS, P, YX,
* DYX, X, H, BUL, Z, DELTY, DF, RFN, RF,
* YP, DYP, ZP, IERR)
IH = IH + 1
C BЫЧИCЛEHИE TOЧHЫX ЗHAЧEHИЙ PEШEHИЯ И ПPOИЗBOДHOЙ
Y1 = SIN(X) + COS(X) + 0.5 * X
Y2 = - SIN(X) - COS(X) + 0.5 * X
DY1 = COS(X) - SIN(X) + 0.5
DY2 = - COS(X) + SIN(X) + 0.5
PRINT 1, X, YX, Y1, Y2, H, DYX, DY1, DY2
GO TO (6, 6, 12, 13, 14, 6, 20), IH
12 JSTART = - 1
H = - 0.005
GO TO 6
13 JSTART = - 1
H = 0.02
GO TO 6
14 JSTART = - 1
H = 0.01
GO TO 6
20 STOP
Результаты:
после первого обращения к подпрограмме -
X YX(1) YX(2)
1.000000000001 - 02 1.014949833751 + 00 - 1.004949833752 + 00
Y1 Y2
1.014949833745 + 00 - 1.004949833748 + 00
H DYX(1) DYX(2)
1.000000000001 - 02 1.489950167079 + 00 - 4.899501670834 - 01
DY1 DY2
1.489950167079 + 00 - 4.899501670798-01
после второго обращения к подпрограмме -
X YX(1) YX(2)
2.000000000001 - 02 1.029798673357 + 00 - 1.009798673364 + 00
Y1 Y2
1.029798673355 + 00 - 1.009798673360 + 00
H DYX(1) DYX(2)
2.000000000001 - 02 1.479801339963 + 00 - 4.798013399754 - 01
DY1 DY2
1.479801339969 + 00 - 4.798013399704-01
после третьего обращения к подпрограмме -
X YX(1) YX(2)
4.000000000002 - 02 1.059189440843 + 00 - 1.019189440856 + 00
Y1 Y2
1.059189440844 + 00 - 1.019189440849 + 00
H DYX(1) DYX(2)
4.000000000002 - 02 1.459210772458 + 00 - 4.592107724779-01
DY1 DY2
1.459210772471 + 00 - 4.592107724729-01
после четвертого обращения к подпрограмме -
X YX(1) YX(2)
1.500000000000 - 02 1.022386939610 + 00 - 1.007386939622 + 00
Y1 Y2
1.022386939612 + 00 - 1.007386939615 + 00
H DYX(1) DYX(2)
- 1.000000000001 - 02 1.484888064588 + 00 - 4.848880646064 - 01
DY1 DY2
1.484888064599 + 00 - 4.848880646005 - 01
после пятого обращения к подпрограмме -
X YX(1) YX(2)
4.000000000002 - 02 1.059189440841 + 00 - 1.019189440856 + 00
Y1 Y2
1.059189440844 + 00 - 1.019189440849 + 00
H DYX(1) DYX(2)
4.000000000002 - 02 1.459210772457 + 00 - 4.592107724784 - 01
DY1 DY2
1.459210772471 + 00 - 4.592107724729 - 01
после шестого обращения к подпрограмме -
X YX(1) YX(2)
3.000000000000 - 02 1.044545533945 + 00 - 1.014545533959 + 00
Y1 Y2
1.044545533947 + 00 - 1.014545533950 + 00
H DYX(1) DYX(2)
2.000000000001 - 02 1.469554533529 + 00 - 4.695545335508 - 01
DY1 DY2
1.469554533543 + 00 - 4.695545335444 - 01
после седьмого обращения к подпрограмме -
X YX(1) YX(2)
5.000000000001 - 02 1.073729429656 + 00 - 1.023729429677 + 00
Y1 Y2
1.073729429661 + 00 - 1.023729429664 + 01
H DYX(1) DYX(2)
4.000000000002 - 02 1.448771091098 + 00 - 4.487710911303 - 01
DY1 DY2
1.448771091120 + 00 - 4.487710911217 - 01