|
Текст подпрограммы и версий ( Фортран ) de43r.zip , de43d.zip , de49r.zip , de49d.zip |
Тексты тестовых примеров ( Фортран ) tde43r.zip , tde43d.zip , tde49r.zip , tde49d.zip |
|
Текст подпрограммы и версий ( Си ) de43r_c.zip , de43d_c.zip , de49r_c.zip , de49d_c.zip |
Тексты тестовых примеров ( Си ) tde43r_c.zip , tde43d_c.zip , tde49r_c.zip , tde49d_c.zip |
|
Текст подпрограммы и версий ( Паскаль ) de43r_p.zip , de43e_p.zip , de49r_p.zip , de49e_p.zip |
Тексты тестовых примеров ( Паскаль ) tde43r_p.zip , tde43e_p.zip , tde49r_p.zip , tde49e_p.zip |
Вычисление решения задачи Коши для системы обыкновенных дифференциальных уравнений второго порядка в конце интервала интегрирования методом Штермера с контролем точности.
Решается задача Коши для системы M обыкновенных дифференциальных уравнений второго порядка
Y '' = F (X, Y) ,
Y = ( y1, ..., yM ) ,
F = ( f1 ( X, y1, ..., yM ), ..., fM ( X, y1, ..., yM ) )
с начальными условиями, заданными в точке XN:
Y (XN) = YN , YN = ( y10, ..., yM0 ) ,
Y ' (XN) = DYN , DYN = ( y10', ..., yM0' ) , -
многошаговым предсказывающе - исправляющим методом пятого порядка точности.
Суть метода состоит в том, что сначала по явной формуле Адамса - Штермера строятся предсказанные значения приближенного решения, которые затем исправляются по неявной формуле Адамса - Штермера.
Решение вычисляется в одной точке XK, которая является концом интервала интегрирования.
Bсе компоненты решения проверяются на точность по меpе погрешности, т.е. если некоторая компонента приближенного решения по абсолютной величине не меньше некоторой заданной константы P (называемой границей перехода), то контроль точности для этой компоненты ведется по относительной погрешности, иначе - по абсолютной.
Березин И.С., Жидков Н.П. Методы вычислений, т.2, Физматгиз, М., 1960.
Бахвалов H.C., Численные методы, т.1, "Hаука", 1975.
Жоголев E.A. Программа интегрирования систем обыкновенных дифференциальных уравнений 2 - го порядка методом Штермера. Сб."Вычислительные методы и программирование", вып.1, Изд - во МГУ, 1962.
SUBROUTINE DE43R (F, M, XN, YN, DYN, XK, HMIN, HMAX, EPS, P,
H, Y, DY, DELTY, DF, RFN, RF, YP, DYP, IERR)
Параметры
| F - | имя подпрограммы вычисления значений правой части дифференциальных уравнений. Первый оператор подпрограммы должен иметь вид: |
|
SUBROUTINE F (X, Y, DY, M) Здесь: X, Y - значения независимой и зависимой переменных, соответственно; вычисленное значение правой части должно быть помещено в DY; в случае системы уравнений, т.е. когда M ≠ 1, параметры Y и DY представляют одномерные массивы длиной M (тип параметров X, Y и DY: вещественный); |
| M - | количество уравнений в системе (тип: целый); |
|
XN, YN - DYN | начальные значения аргумента, решения и его производной; в случае системы уравнений (т.е. M ≠ 1) YN и DYN представляют одномерные массивы длины M (тип вещественный); |
| XK - | значение аргумента, при котоpом требуется вычислить решение задачи Коши (конец интервала интегрирования); XK может быть больше, меньше или pавно XN (тип: вещественный); |
| HMIN - | минимальное значение абсолютной величины шага, котоpое разрешается использовать при интегрировании данной системы уравнений (тип: вещественный); |
| HMAX - | максимальное значение абсолютной величины шага, котоpое разрешается использовать при интегрировании данной системы уравнений (тип: вещественный); |
| EPS - | допустимая меpа погрешности, с которой тpебуется вычислить все компоненты решения (тип: вещественный); |
| P - | граница перехода, используемая при оценке меры погрешности решения (тип: вещественный); |
| H - | вещественная переменная, содержащая начальное значение шага интегрирования; может задаваться с учетом направления интегрирования, т.е. положительным, если XK > XN, отрицательным, если XK < XN, или без такого учета в виде абсолютной величины; |
| Y - | искомое решение задачи Коши, вычисленное подпрограммой для значения аргумента XK; для системы уравнений задается одномерным массивом длины M; в случае совпадения значений параметров XN и XK значение Y полагается равным начальному значению YN (тип: вещественный); |
|
DY - DELTY RFN, RF YP, DYP | одномерные вещественные рабочие массивы длины M; |
| DF - | двумерный вещественный рабочий массив размеpа M * 5; |
| IERR - | целая переменная, служащая для сообщения об ошибках, обнаруженных в процессе работы подпрограммы; при этом: |
| IERR=65 | и IERR=66 - когда какая - нибудь компонента решения не может быть вычислена с требуемой точностью EPS при заданных начальном шаге H и его минимальном значении HMIN, при этом IERR = 66 указывает, что тpебуемая точность не может быть достигнута при разгоне, а IERR=65 - после разгона; при желании интегрирование системы можно повторить обращением к подпрограмме с новыми значениями параметров H и HMIN. |
Версии
| DE43D - | вычисление решения задачи Коши для системы обыкновенных дифференциальных уравнений второго порядка в конце интервала интегрирования методом Штермера с повышенной точностью. При этом параметры XN, YN, DYN, XK, HMIN, HMAX, EPS, P, H, Y, DY, DELTY, DF, RFN, RF, YP, DYP и параметры X, Y и DY в подпрограмме F должны иметь тип DOUBLE PRECISION. |
| DE49R - | назначение этой подпрограммы такое же, как и подпрограммы DE43R, но схема организации вычислений в ней несколько иная, чем в DE43R, а именно: во - первых, исправленное и предсказанное значения решения вычисляются через значения некоторой промежуточной переменной, введение которой позволяет уменьшить вычислительную погрешность в приближенном значении решения; во - вторых, для вычисления приближенного решения на одном шаге используется на одно вычисление правой части уравнения больше, чем в подпрограмме DE43R. Таким образом, обеспечивается возможность вычислить решение уравнения более точно, чем в DE43R. |
| DE49D - | то же, что и DE49R, но все вычисления ведутся с удвоенным числом значащих цифр; при этом параметры XN, YN, DYN, XK, HMIN, HMAX, EPS, P, H, Y, DY, DELTY, DF, RFN, RF, YP, DYP и параметры X, Y и DY в подпрограмме F должны иметь тип DOUBLE PRECISION. |
Вызываемые подпрограммы
| DE42R - | выполнение одного шага численного интегрирования системы обыкновенных дифференциальных уравнений второго порядка методом Штермера с контролем точности. |
| DE42D - | выполнение одного шага численного интегрирования системы обыкновенных дифференциальных уравнений второго порядка методом Штермера с повышенной точностью. |
|
DE48R - DE48D | выполнение одного шага численного интегрирования системы обыкновенных дифференциальных уравнений второго порядка методом Штермера, преобразованным к виду, где влияние вычислительной погрешности меньше. |
|
UTDE12 - UTDE13 UTDE16 UTDE17 | подпрограммы выдачи диагностических сообщений. |
|
Подпрограммы DE42R и UTDE12 вызываются при работе подпрограммы DE43R, а подпрограммы DE42D и UTDE13 - при работе подпрограммы DE43D. Подпрограммы DE48R и UTDE16 вызываются при работе подпрограммы DE49R, а подпрограммы DE48D и UTDE17 - при работе подпрограммы DE49D. |
Замечания по использованию
|
B общем случае заданная точность не гарантируется. При работе подпрограммы и ее версии значения параметров M, XN, YN, DYN, XK, HMIN, HMAX, EPS, P сохраняются. Если после работы подпрограммы нет необходимости иметь начальные значения решения YN и / или его производной DYN, то параметры YN и Y, а в случае производной параметры DYN и DY, при обращении к ней можно совместить. При этом следует иметь в виду, что в случае аварийного выхода из подпрограммы, т.е. со значением IERR = 65 или 66, значения параметров YN и DYN будут испорчены. B том случае, когда важнее сократить время вычисления, целесообразно использовать подпрограмму DE43R; если требуется поточнее вычислить решение, то следует использовать подпрограмму DE49R. Значение параметра HMAX не должно быть слишком большим, т.к. в противном случае величина шага интегрирования может достичь таких размеров, которые приведут к абсолютной неустойчивости численного метода. Tак как при интегрировании уравнений с помощью подпрограммы DE49R и DE49D используются общие блоки с именами COM48R и COM48D, соответственно, поэтому пользователю не рекомендуется использовать для своих целей общие блоки с такими же именами. |
y1'' = - 6y1 - 7y2 ,
y2'' = - 3y1 - 2y2 + 2x ;
y1 (0) = 0 , y2 (0) = 0 ,
y1' (0) = - 1 / 9 , y2' (0) = - 3 .
Приводятся подпрограмма вычисления значений правой части и фрагмент вызывающей программы, выполняющей дважды интегрирование данной системы, сначала на отрезке [0, 5], затем на отрезке [-5, 0], справа налево, а также результаты счета.
SUBROUTINE F (X, Y, Z, M)
DIMENSION Y(2), Z(2)
Z(1) = - 6. * Y(1) - 7. * Y(2)
Z(2) = - 3. * Y(1) - 2. * Y(2) + 2. * X
RETURN
END
DIMENSION YN(2), DYN(2), Y(2), DY(2), DELTY(2), DF(10),
* RFN(2), RF(2), YP(2), DYP(2)
EXTERNAL F
M = 2
XN = 0.
YN(1) = 0.
YN(2) = 0.
DYN(1) = - 1. / 9.
DYN(2) = - 3.
XK = 5.
HMIN = 1.E - 10
HMAX = 0.1
EPS = 0.0001
P = 100.
H = 0.01
C Y1 И Y2 ЯBЛЯЮTCЯ TOЧHЫM ЗHAЧEHИEM PEШEHИЯ
Y3 = (EXP(5.) - EXP(- 5.)) / 3.
Y1 = Y3 - 7. * SIN(15.) / 9. + 70. / 9.
Y2 = - Y3 - SIN(15.) / 3. - 20. / 3.
DO 2000 K = 1, 2
CALL DE43R (F, M, XN, DYN, XK, HMIN, HMAX, EPS, P, H, Y, DY,
* DELTY, DF, RFN, RF, YP, DYP, IERR)
C ПEЧATЬ ПPИБЛИЖEHHOГO И TOЧHOГO ЗHAЧEHИЙ PEШEHИЯ
PRINT 1, XK, Y(1), Y(2), H, Y1, Y2
XK = - 5.
Y1 = - Y1
Y2 = - Y2
H = 0.01
2000 CONTINUE
STOP
1 FORMAT (3E18.9)
Результаты:
после первого обращения к подпрограмме -
XK Y(1) Y(2)
5.000000000 + 00 5.674103844 + 01 - 5.655217994 + 01
H Y1 Y2
1.000000000 - 01 5.674080540 + 01 - 5.635223633 + 01
после второго обращения к подпрограмме -
XK Y(1) Y(2)
- 5.000000000 + 00 - 5.674103838 + 01 5.635217988 + 01
H Y1 Y2
- 1.000000000 - 01 - 5.674080540 + 01 5.635223633 + 01