Текст подпрограммы и версий ( Фортран ) de35r.zip , de35d.zip |
Тексты тестовых примеров ( Фортран ) tde35r.zip , tde35d.zip |
Текст подпрограммы и версий ( Си ) de35r_c.zip , de35d_c.zip |
Тексты тестовых примеров ( Си ) tde35r_c.zip , tde35d_c.zip |
Текст подпрограммы и версий ( Паскаль ) de35r_p.zip , de35e_p.zip |
Тексты тестовых примеров ( Паскаль ) tde35r_p.zip , tde35e_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' ) ,
отыскиваются необходимые при численном интегрировании этой системы многошаговым методом Штермера значение решения Y в первом (после начальной точки XN) узле интегрирования ХN + Н, его конечная разность первого порядка и конечные разности правой части системы до четвертого порядка включительно в том же узле XN + H. Четвертый порядок разностей соответствует шестому порядку точности используемого метода Штермера.
Указанная процедура вычисления необходимых для счета начальных значений называется разгоном.
Используемый в подпрограмме метод разгона является итерационным способом, опирающимся на пару формул, одна из которых представляет значение решения в узле XN + H через начальные условия (т.е. через его значение и значение его производной в начальной точке XN), а также через конечные разности правой части системы в точке XN, а другая является экстраполяционной формулой, применяемой для предсказания приближенного значения решения в используемом для интегрирования методе Штермера.
Березин И.С., Жидков Н.П. Методы вычислений, т.2. Физматгиз, М., 1960.
SUBROUTINE DE35R (F, M, IORDER, XN, YN, DYN, H, X, YX, DY, DF, RFN, RF, R)
Параметры
F - | имя подпрограммы вычисления значений правой части системы. Первый оператор подпрограммы должен иметь вид: |
SUBROUTINE F (X, Y, DY, M). Здесь: X, Y - значения независимой и зависимой переменных, соответственно; вычисленное значение правой части должно быть помещено в DY. B случае системы уравнений, т.е. когда M ≠ 1, параметры Y и DY представляют одномерные массивы длиной M (тип параметров X, Y и DY: вещественный); |
M - | количество уравнений в системе (тип: целый); |
IORDER - | порядок точности без единицы того метода Штермера, для которого выполняется разгон и который будет использоваться при интегрировании данной системы уравнений (тип: целый); |
XN, YN - DYN | начальные значения аргумента, решения и его производной; в случае системы уравнений (т.е. M ≠ 1) YN и DYN представляют одномерные массивы длиной M (тип: вещественный); |
H - | вещественное значение шага интегрирования; |
X - | вещественная переменная, значение которой на выходе из подпрограммы представляет первый (после начальной точки XN) узел интегрирования XN + H, в котоpом вычислены необходимые для интегрирования данной системы уравнений начальные значения; |
YX, DY - | одномерные вещественные массивы длиной M, в которых запоминаются значения решения и его первой разности, вычисленные в узле X, при этом погрешность решения и разности имеет (IОRDЕR + 1) - й порядок по H; |
DF - | двумерный вещественный массив размера M * IORDER, в котоpом запоминаются значения правой части системы и ее разностей до (IОRDЕR - 1) - го порядка включительно, вычисленные в узле X и умноженные на коэффициент H2 / 12, при этом элемент DF (I, 1) этого массива содержит значение правой части I - го уравнения системы, DF (I, J + 1) - ее J - ю разность, погрешность которой имеет (IОRDЕR + 1) - й порядок по H; |
RFN - RF, R | одномерные вещественные рабочие массивы длины M. |
Версии
DE35D - | построение начальных значений при интегрировании системы обыкновенных дифференциальных уpавнений второго порядка методом Штермера без контроля точности, при этом все вычисления с действительными числами выполняются с удвоенным числом значащих цифр; в этом случае параметры XN, YN, DYN, H, X, YX, DY, DF, N, RF, R и параметры X, Y и DY в подпрограмме F должны иметь тип DOUBLE PRECISION. |
Вызываемые подпрограммы
Подпрограммы DE35R и DE35D используют рабочие подпрограммы DE28RS и DE28DS, соответственно. |
Замечания по использованию
Значение параметра H, задаваемое при обращении к подпрограмме, должно быть таким, чтобы узел XN + (IORDER - 1) * H не выходил за конец интервала интегрирования. Значение параметра IORDER на входе в подпрограмму полагается равным 5. В дальнейшем предполагается расширить допустимое множество значений этого параметра. При работе подпрограммы и ее версии значения параметров M, IORDER, XN, YN, DYN, H сохраняются. |
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 YN(2), DYN(2), Y(2), DY(2), DELTY(2), DF(2, 5), * RFN(2), RF(2) EXTERNAL F M = 2 IORDER = 5 XN = 0.75 * 3.14159265359 YN(1) = SQRT(2.) / 2. YN(2) = - YN(1) DYN(1) = YN(2) DYN(2) = - YN(1) H = 0.01 CALL DE35R (F, M, IORDER, XN, YN, DYN, H, X, * Y, DY, DF, RFN, RF, DELTY) Результаты: X = 2.366194490 + 00 Y(1) = 7.000004762 - 01 Y(2) = - 7.141423761 - 01 DY(1) = - 7.106305006 - 03 DY(2) = - 7.035594917 - 03 DF(1, 1) = - 5.83333730151 - 06 DF(1, 2) = 5.92192083634 - 08 DF(1, 3) = 5.89251453187 - 10 DF(1, 4) = - 5.86490578325 - 12 DF(1, 5) = - 5.76760861293 - 14 DF(2, 1) = 5.95118646749 - 06 DF(2, 2) = 5.86299576250 - 08 DF(2, 3) = - 5.89250072347 - 10 DF(2, 4) = - 5.92374760355 - 12 DF(2, 5) = 6.01324545713 - 14