|
Текст подпрограммы и версий ( Фортран ) de40r.zip , de40d.zip |
Тексты тестовых примеров ( Фортран ) tde40r.zip , tde40d.zip |
|
Текст подпрограммы и версий ( Си ) de40r_c.zip , de40d_c.zip |
Тексты тестовых примеров ( Си ) tde40r_c.zip , tde40d_c.zip |
|
Текст подпрограммы и версий ( Паскаль ) de40r_p.zip , de40e_p.zip |
Тексты тестовых примеров ( Паскаль ) tde40r_p.zip , tde40e_p.zip |
Выполнение одного шага численного интегрирования системы обыкновенных дифференциальных уравнений второго порядка методом Штермера без контроля точности.
Выполняется один шаг численного интегрирования системы М обыкновенных дифференциальных уравнений второго порядка
Y '' = F (X, Y) ,
Y = ( y1, ..., yM ) ,
F = ( f1 (X, y1, ..., yM), ..., fM (X, y1, ..., yM) )
многошаговым предсказывающе - исправляющим методом пятого порядка точности. Суть метода состоит в том, что сначала по явной формуле Адамса - Штермера строятся предсказанные значения приближенного решения, которые затем исправляются по неявной формуле Адамса - Штермера.
По значению решения YX в узле Xn вычиляется значение этого решения в узле Xn + H. При этом, если Xn является началом ХN интервала интегрирования, то для вычисления решения в XN + H требуется еще и знание производной решения в XN.
Bместе с решением подпрограмма определяет также всю информацию, необходимую для вычисления приближенного решения в любом следующем узле. Значение H шага интегрирования задается пользователем при обращении к подпрограмме.
Березин И.С., Жидков Н.П. Методы вычислений, т.2. Физматгиз, М., 1960.
Жоголев E.A. Программа интегрирования систем обыкновенных дифференциальных уравнений 2 - го порядка методом Штермера. Сб. "Вычислительные методы и программирование", вып. 1 , Изд - во МГУ, 1962.
SUBROUTINE DE40R (F, M, JSTART, YX, X, DY, H, DELTY, DF,
RFN, RF, YP, DYP, IERR)
Параметры
| F - |
имя подпрограммы вычисления значений
правой части системы дифференциальных
уpавнений. Первый оператор подпрограммы должен иметь вид: SUBROUTINE F (X, Y, DY, M). Здесь: X, Y - значения независимой и зависимой переменных, соответственно; вычисленное значение правой части должно быть помещено в DY; в случае системы уравнений, т.е. когда M ≠ 1, параметры Y и DY представляют одномерные массивы длиной M (тип параметров X, Y, DY: вещественный); |
| M - | количество уравнений в системе (тип: целый); |
| JSTART - | целый указатель режима использования подпрограммы, имеющий следующие значения: |
| 0 - | первое обращение к подпрограмме должно быть выполнено с нулевым значением JSTART; |
| +1 - | выполнить один шаг интегрирования системы дифференциальных уравнений для значений независимой и зависимой переменных и шага интегрирования, заданных параметрами X, YX и H, соответственно; |
| -1 - | повторить последний шаг интегрирования с новым значением параметра H; |
| на выходе из подпрограмы JSTART полагается равным + 1; |
| X, YX - | заданные вещественные значения аргумента и решения в нем; в результате работы подпрограммы в X получается новое значение аргумента, а в YX - соответствующее значение решения; в случае системы уравнений, т.е. когда M ≠ 1, YX задается одномерным массивом длины M; |
| DY - | одномерный вещественный рабочий массив длины M; при первом обращении к подпрограмме (т.е. со значением JSTART = 0) в этом массиве задается значение производной pешения в начальном узле интервала интегрирования; |
| H - | вещественное значение шага интегрирования; |
|
DELTY - RFN, RF YP, DYP | одномерные вещественные рабочие массивы длины M; |
| DF - | двумерный вещественный рабочий массив размеpа M * 5, в котоpом запоминаются значения правой части системы и ее разностей до четвертого порядка включительно, умноженные на H 2 / 12; |
| IERR - | целая переменная, значение которой в pезультате работы подпрограммы полагается равным 65, если при заданном значении шага интегрирования не может быть достигнута сходимость численного решения к точному значению решения; в этом случае интегрирование системы можно повторить обращением к подпрограмме с новым значением параметра H и со значением JSTART = + 1 . |
Версии
| DE40D - | выполнение одного шага численного интегрирования системы обыкновенных дифференциальных уpавнений второго порядка методом Штермера без контроля точности, при этом все вычисления над действительными числами ведутся с удвоенным числом значащих цифр. В этом случае параметры YX, X, DY, H, DELTY, DF, RFN, RF, YP, DYP и параметры X, Y и DY в подпрограмме F должны иметь тип DOUBLE PRECISION. |
Вызываемые подпрограммы
| DE35R - | построение начальных значений при интегрировании системы обыкновенных дифференциальных уравнений второго порядка методом Штермера без контроля точности; вызывается при работе подпрограммы DE40R. |
| DE35D - | построение начальных значений при интегрировании системы обыкновенных дифференциальных уравнений второго порядка методом Штермера без контроля точности, при этом все вычисления над действительными числами ведутся с удвоенным числом значащих цифр; вызывается при работе подпрограммы DE40D. |
| UTDE12 - | подпрограмма выдачи диагностических сообщений при работе подпрограммы DE40R. |
| UTDE13 - | подпрограмма выдачи диагностических сообщений при работе подпрограммы DE40D. |
| Kpоме того, DE40R и DE40D используют рабочие подпрограммы DE42RP, DE28RS и DE42DP, DE28DS, соответственно. |
Замечания по использованию
|
При работе подпрограммы и ее версии значения параметров M и H сохраняются. При многократном использовании подпрограммы и ее версии для вычисления решения на отрезке значения параметров M, YX, X, DY, DF, YP, DYP не должны изменяться в вызывающей программе между последовательными обращениями к ним. |
y1'' = - y1 ,
y2'' = - y2 ,
y1 (3/4 π) = √2 / 2 , y2 (3/4 π) = - √2 / 2 ,
y1' (3/4 π) = - √2 / 2 , y2' (3/4 π) = - √2 / 2 .
Приводятся подпрограмма вычисления значений правой части системы и фрагмент вызывающей программы, выполняющей несколько шагов из одной и той же точки, а также результаты счета.
SUBROUTINE F (X, Y, DY, M)
DIMENSION Y(2), DY(2)
DY(1) = - Y(1)
DY(2) = - Y(2)
RETURN
END
DIMENSION Y(2), DY(2), DELTY(2), DF(10), RFN(2), RF(2), YP(2),
* DYP (2)
EXTERNAL F
M = 2
JSTART = 0
Y(1) = SQRT(2.) / 2.
Y(2) = - Y(1)
X = 0.75 * 3.14159265359
DY(1) = Y(2)
DY(2) = - Y(1)
H = 0.01
IH = 0
6 CONTINUE
CALL DE40R (F, M, JSTART, Y, X, DY, H, DELTY, DF, RFN, RF,
* YP, DYP, IERR)
IH = IH + 1
Y1 = SIN(X)
Y2 = COS(X)
C ПEЧATЬ TOЧHOГO И ПPИБЛИЖEHHOГO ЗHAЧEHИЙ PEШEHИЯ
PRINT 1, X, Y, H, Y1, Y2
IF (IH .EQ. 1) GO TO 6
IF (IH .EQ. 5) GO TO 20
JSTART = - 1
IF (IH .EQ. 3) GO TO 7
IF (IH .EQ. 4) GO TO 8
H = - 0.005
GO TO 6
7 H = 0.02
GO TO 6
8 H = 0.01
GO TO 6
20 STOP
1 FORMAT (3E18.9)
Результаты:
после первого обращения к подпрограмме -
X Y(1) Y(2)
2.366194490 + 00 7.000004762 - 01 - 7.141423761 - 01
H Y1 Y2
1.000000000 - 02 7.000004762 - 01 - 7.141423761 - 01
после второго обращения к подпрограмме -
X Y(1) Y(2)
2.376194490 + 00 6.928241717 - 01 - 7.211065574 - 01
H Y1 Y2
1.000000000 - 02 6.928241717 - 01 - 7.211065574 - 01
после третьего обращения к подпрограмме -
X Y(1) Y(2)
2.361194490 + 00 7.035624232 - 01 - 7.106334615 - 01
H Y1 Y2
- 5.000000000 - 03 7.035624232 - 01 - 7.106334615 - 01
после четвертого обращения к подпрограмме -
X Y(1) Y(2)
2.386194490 + 00 6.855785854 - 01 - 7.279986286 - 01
H Y1 Y2
2.000000000 - 02 6.855785854 - 01 - 7.279986286 - 01
после пятого обращения к подпрограмме -
X Y(1) Y(2)
2.376194490 + 00 6.928241717 - 01 - 7.211065574 - 01
H Y1 Y2
1.000000000 - 02 6.928241717 - 01 - 7.211065574 - 01