|
Текст подпрограммы и версий ( Фортран ) de08r.zip , de08d.zip |
Тексты тестовых примеров ( Фортран ) tde08r.zip , tde08d.zip |
|
Текст подпрограммы и версий ( Си ) de08r_c.zip , de08d_c.zip |
Тексты тестовых примеров ( Си ) tde08r_c.zip , tde08d_c.zip |
|
Текст подпрограммы и версий ( Паскаль ) de08r_p.zip , de08e_p.zip |
Тексты тестовых примеров ( Паскаль ) tde08r_p.zip , tde08e_p.zip |
Выполнение одного шага численного интегрирования системы обыкновенных дифференциальных уравнений первого порядка методом Инглэнда.
Выполняется один шаг численного интегрирования системы M обыкновенных дифференциальных уравнений первого порядка
Y ' = F(X, Y) ,
Y = ( y1, ... , yM ) , F = ( f1(X, y1,..., yM), ... , fM(X, y1,..., yM) )
методом Инглэнда пятого порядка точности. По заданному значению решения YX в узле Xn вычисляется значение этого решения в узле Xn + H. Все компоненты решения вычисляются с контролем точности по мере погрешности следующим образом.
Если некоторая компонента приближенного решения по абсолютной величине не меньше некоторой наперед заданной константы P, то контроль точности для этой компоненты ведется по относительной погрешности, иначе - по абсолютной. В качестве абсолютной погрешности решения принимается оценка главного члена асимптотического разложения погрешности метода на шаге, получаемая при вычитании двух выражений, представляющих приближенные значения решения пятого и четвертого порядка точности. При этом на одном шаге интегрирования для определения решения и погрешности используются всего шесть вычислений правой части системы.
Значение H может быть меньше или равно значению шага интегрирования, задаваемому пользователем при обращении к подпрограмме.
| 1. | Бахвалов Н.С., Жидков Н.П., Кобельков Г.М. Численные методы. - М.: Наука, 1987. |
| 2. | England R. Error estimates for Runge - Kutta type solutions to systems of ordinary differential equations. The Computer Journal. 1969 - v.12, No.2. |
SUBROUTINE DE08R (F, M, JSTART, HMIN, EPS, P, YX, X, H, BUL,
XP, YP, R1, R2, R3, IERR)
Параметры
| F - |
имя подпрограммы вычисления значений правой части
системы дифференциальных уравнений. Первый оператор
подпрограммы должен иметь вид SUBROUTINE F (X, Y, DY, M). Здесь X, Y - значения независимой и зависимой переменных, соответственно; вычисленное значение правой части должно быть помещено в DY. В случае системы уравнений, т.е. когда M ≠ 1, параметры Y и DY представляют одномерные массивы длины M (тип параметров X, Y и DY: вещественный); |
| 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, соответственно, при входе в нее (т.е. предыдущий узел и решение в нем); |
| R1 - | двумерный вещественный рабочий массив размера М*5; |
| R2, R3 - | одномерные вещественные рабочие массивы длины M; |
| IERR - | целая переменная, значение которой в результате работы подпрограммы полагается равным 65, если какая - нибудь компонента решения не может быть вычислена с требуемой точностью EPS; в этом случае последний шаг интегрирования системы можно повторить обращением к подпрограмме с новыми значениями параметров H, HMIN и значением JSTART = - 1 . |
Версии
| DE08D - | выполнение одного шага численного интегрирования системы обыкновенных дифференциальных уравнений первого порядка методом Инглэнда с удвоенным числом значащих цифр. При этом параметры HMIN, EPS, P, YX, X, H, XP, YP, R1, R2, R3 и параметры X, Y и DY в подпрограмме F должны иметь тип DOUBLE PRECISION. |
Вызываемые подпрограммы
| UTDE20 - | подпрограмма выдачи диагностических сообщений при работе подпрограммы DE08R. |
| UTDE21 - | подпрограмма выдачи диагностических сообщений при работе подпрограммы DE08D. |
Замечания по использованию
|
В общем случае заданная точность не гарантируется. При работе подпрограммы значения параметров M, HMIN, EPS, P сохраняются. При работе подпрограммы счета правой части F значения параметров X, Y и M не должны изменяться. При обращении к подпрограмме со значением JSTART = - 1 в качестве исходных значений аргумента и решения принимаются значения параметров XP и YP, соответственно, т.е. те значения, которые эти параметры получили после самого последнего обращения к подпрограмме с неотрицательным значением JSTART. На выходе из подпрограммы в массиве R3 содержится оценка локальной погрешности метода на выполненном шаге. |
y1' = - y1 - 5y2 , y1(0) = 1
y2' = y1 + y2 , y2(0) = 1
Точное решение системы:
y1 = cos2x - 3 sin2x ,
y2 = cos2x + sin2x
Приводятся подпрограмма вычисления значений правой части системы и фрагмент вызывающей программы, а также результаты счета.
SUBROUTINE FENGL (X, Y, Z, M)
DIMENSION Y(2), Z(2)
Z(1) = -Y(1) - 5.*Y(2)
Z(2) = Y(1) + Y(2)
RETURN
END
DIMENSION YX(2), YP(2), R1(2, 5), R2(2), R3(2)
LOGICAL BUL
EXTERNAL FENGL
M = 2
P = 100.
HMIN = 1.E-12
H = 1.E-2
EPS = 1.E-10
JSTART = 1
X = 0.
YX(1) = 1.
YX(2) = 1.
BUL = .FALSE.
IH = 0
2 CALL DE08R (FENGL, M, JSTART, HMIN, EPS, P, YX, X, H, BUL,
* XP, YP, R1, R2, R3, IERR)
C BЫЧИCЛEHИE TOЧHOГO ЗHAЧEHИЯ PEШEHИЯ
Y(1) = COS(2.*X) - 3.*SIN(2.*X)
Y(2) = COS(2.*X) + SIN(2.*X)
IH = IH + 1
GO TO (2, 2, 13, 14, 15, 2, 10), IH
13 JSTART = -1
H = -0.005
GO TO 2
14 JSTART = -1
H = 0.02
GO TO 2
15 EPS = 1.E-12
H = 0.01
GO TO 2
10 STOP
END
Результаты:
после первого обращения к подпрограмме -
X YX(1) YX(2)
1.000000000001-02 9.398040065844-01 1.019798673355 + 00
H Y1 Y2
1.000000000001-02 9.398040065835-01 1.019798673356 + 00
контрольный член: -7.983658179000-11 2.673061771929-11
после второго обращения к подпрограмме -
X YX(1) YX(2)
2.000000000001-02 8.792321040974-01 1.039189440839 + 00
H Y1 Y2
1.000000000001-02 8.792321040992-01 1.039189440846 + 00
контрольный член: -8.037659426918-11 2.623323780426-11
после третьего обращения к подпрограмме -
X YX(1) YX(2)
3.00000000000-02 8.183085204928-01 1.058164546403 + 00
H Y1 Y2
1.00000000001-02 8.183085204937-01 1.058164546412 + 00
контрольный член: -8.088818503892-11 2.565069276093-11
после четвертого обращения к подпрограмме -
X YX(1) YX(2)
1.500000000000-02 9.095635331360-01 1.029545533938 + 00
H Y1 Y2
-1.000000000001-02 9.095635331387-01 1.029545533947 + 00
контрольный член: 2.458477865729-12 -8.100187187665-13
после пятого обращения к подпрограмме -
X YX(1) YX(2)
3.000000000000-02 8.183085204928-01 1.058164546403 + 00
H Y1 Y2
1.000000000001-02 8.183085204937-01 1.058164546412 + 00
контрольный член: -8.088818503892-11 2.565059276093-11
после шестого обращения к подпрограмме -
X YX(1) YX(2)
3.249999999997-02 8.030255271615-01 1.062842482492+00
H Y1 Y2
2.500000000001-03 8.030255271660-01 1.062842482504+00
контрольный член: -1.136868377216-13 1.953992523340-14
после седьмого обращения к подпрограмме -
X YX(1) YX(2)
3.499999999997-02 7.877224582326-01 1.067493847573+00
H Y1 Y2
2.500000000001-03 7.877224582389-01 1.067493847588+00
контрольный член: -9.237055564881-14 1.065814103640-14