|
Текст подпрограммы и версий ( Фортран ) de38r.zip , de38d.zip |
Тексты тестовых примеров ( Фортран ) tde38r.zip , tde38d.zip |
|
Текст подпрограммы и версий ( Си ) de38r_c.zip , de38d_c.zip |
Тексты тестовых примеров ( Си ) tde38r_c.zip , tde38d_c.zip |
|
Текст подпрограммы и версий ( Паскаль ) de38r_p.zip , de38e_p.zip |
Тексты тестовых примеров ( Паскаль ) tde38r_p.zip , tde38e_p.zip |
Построение начальных значений при интегрировании системы обыкновенных дифференциальных уравнений второго порядка с правой частью, зависящей от производной, методом Штермера с контролем точности.
Для системы M обыкновенных дифференциальных уравнений второго порядка
Y '' = F( X, Y, Y ' ) ,
Y = ( y1, ..., yM ) ,
F = ( f1( X, y1, ..., yM, y1', ..., yM' ), ..., fM( X, y1, ..., yM, y1', ..., yM' ) ) ,
с начальными условиями, заданными в точке XN:
Y(XN) = YN, YN = ( y10, ..., yM0 ) ,
Y ' (XN) = DYN, DYN = ( y10', ..., yM0' ) ,
отыскиваются необходимые при численном интегрировании этой системы методом Штермера значения решения Y и производной Y ' в первом (после начальной точки XN) узле интегрирования ХN + Н, первая разностная производная назад и конечные разности правой части системы до четвертого порядка включительно в том же узле XN + H.
Четвертый порядок разностей соответствует пятому порядку точности используемого метода Штермера.
Указанная процедура вычисления необходимых для счета начальных значений называется разгоном.
Используемый в подпрограмме метод разгона является итерационным способом, опирающимся на три формулы.
Первая формула представляет значение решения в узле ХN + Н через начальные условия (т.е. через его значение и значение его производной в начальной точке XN), а также через конечные расности правой части системы в точке XN.
Вторая формула является предсказывающей формулой Адамса, по которой вычисляется производная решения.
Третья формула является экстраполяционной формулой, применяемой для предсказания приближенного значения решения в используемом для интегрирования методе Штермера.
При этом приближенные значения вычисляются через значения его разностной производной назад.
B качестве критерия точности при разгоне принята оценка первого отброшенного члена экстраполяционной формулы Штермера, при этом все компоненты решения проверяются на точность по мере погрешности, который заключается в следующем: если некоторая компонента приближенного решения по абсолютной величине не меньше некоторой наперед заданной константы P, то контроль точности для этой компоненты ведется по относительной погрешности, иначе - по абсолютной.
Березин И.С., Жидков Н.П. Методы вычислений, т.2. Физматгиз, М., 1960.
Бахвалов H.C., Численные методы, т.1, "Hаука", 1975.
SUBROUTINE DE38R (F, M, IORDER, XN, YN, DYN, HMIN, EPS, P,
H, X, YX, DYX, Z, DF, RFN, RF, R, 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 - | количество уравнений в системе (тип: целый); |
| IORDER - | порядок точности без единицы того метода Штермера, для которого выполняется разгон и который будет использоваться при интегрировании данной системы уравнений (тип: целый); |
|
XN, YN - DYN | начальные значения аргумента, решения и его производной; в случае системы уравнений (т.е. M ≠ 1) YN и DYN представляют одномерные массивы длиной M (тип: вещественный); |
| HMIN - | минимальное значение абсолютной величины шага, котоpое разрешается использовать при интегрировании данной системы уравнений (тип: вещественный); |
| EPS - | допустимая меpа погрешности, с которой тpебуется вычислить все компоненты решения (тип: вещественный); |
| P - | граница перехода, используемая при оценке меры погрешности решения (тип: вещественный); |
| H - | вещественная переменная, содержащая значение шага интегрирования. Если для этого значения шага точность при разгоне достигается, то именно он и реализуется на разгоне, иначе этот шаг уменьшается подпрограммой до тех пор, пока не будет достигнута заданная точность; |
| X - | вещественная переменная, значение которой на выходе из подпрограммы представляет первый (после начальной точки XN) узел интегрирования XN + H, в котоpом вычислены необходимые для интегрирования данной системы уравнений начальные значения; |
|
YX - DYX, Z | одномерные вещественные массивы длины M, в которых запоминаются значения решения, его производной и первой разностной производной назад, вычисленные в узле X; при этом погрешность решения имеет (IОRDЕR + 2) - й порядок по H, а погрешность DYX и Z - (IОRDЕR + 1) - й порядок по H; |
| DF - | двумерный вещественный массив размера M * IORDER, в котоpом запоминаются значения правой части системы и ее разностей до (IОRDЕR - 1) - го порядка включительно, вычисленные в узле X и умноженные на коэффициент H / 12; при этом элемент DF (I, 1) этого массива содержит значение правой части I - ого уравнения системы, а DF (I, J + 1) - ее J - ю разность, погрешность которой имеет (IОRDЕR + 1) - й порядок по H; |
|
RFN - RF, R | одномерные вещественные рабочие массивы длины M. |
| IERR - | целая переменная, значение которой в pезультате работы подпрограммы полагается равным 65, если на разгоне не может быть достигнута требуемая точность EPS; в этом случае разгон можно начать сначала обращением к подпрограмме с новыми значениями паpаметpов H и HMIN. |
Версии
| DE38D - | построение начальных значений при интегрировании системы обыкновенных дифференциальных уpавнений второго порядка, с правой частью, зависящей от производной, методом Штермера с контpолем точности, при этом все промежуточные вычисления выполняются с удвоенным числом значащих цифр. B этом случае параметры XN, YN, DYN, HMIN, EPS, P, H, X, YX, DYX, Z, DF, RFN, RF, R и параметры X, Y, DY и D2Y в подпрограмме F должны иметь тип DOUBLE PRECISION. |
Вызываемые подпрограммы
| UTDE16 - | подпрограмма выдачи диагностических сообщений при работе подпрограммы DE38R. |
| UTDE17 - | подпрограмма выдачи диагностических сообщений при работе подпрограммы DE38D. |
Замечания по использованию
|
B общем случае требуемая точность при разгоне не гарантируется. Значение параметра H, задаваемое при обращении к подпрограмме, должно быть таким, чтобы узел XN + (IORDER - 1) * H не выходил за конец интервала интегрирования. Значение параметра IORDER на входе в подпрограмму полагается равным 5. B дальнейшем предполагается расширить допустимое множество значений этого параметра. При работе подпрограммы и ее версии значения параметров M, IORDER, XN, YN, DYN, HMIN, EPS, P сохраняются. |
y1'' = 0.5 ( y2 - y1 + y1' + y2' ) - 0.5 , y1(0) = 1 , y1' (0) = 1.5 ,
y2'' = y1 - 0.5x , y2(0) = - 1 , y2' (0) = - 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 YN(2), DYN(2), YX(2), DYX(2), Z(2), DF(2, 5), RFN(2),
* RF(2), R(2)
EXTERNAL F2
YN(1) = 1.
YN(2) = - 1.
DYN(1) = 1.5
DYN(2) = - 0.5
XN = 0.
P = 100.
H = 0.01
HMIN = 1.E - 11
EPS = 0.0001
CALL DE38R (F2, 2, 5, XN, YN, DYN, HMIN, EPS, P, H, X,
* YX, DYX, Z, DF, RFN, RF, R, IERR)
Результаты:
IERR = 0
X = 1.000000000001 - 02
YX(1) = 1.014949833751 YX(2) = - 1.004949833752
DYX(1) = 1.489950167079 DYX(2) = - 4.899501670834 - 01
Z(1) = 1.494983375078 Z(2) = - 4.949833750843 - 01
DF(1, 1) = - 8.416248614562 - 04 DF(1, 2) = - 8.291528123650 - 06
DF(1, 3) = 8.333255330228 - 08 DF(1, 4) = 8.377316618179 - 10
DF(1, 5) = - 8.498979298110 - 12
DF(2, 1) = 8.416248614518 - 04 DF(2, 2) = 8.291528119209 - 06
DF(2, 3) = - 8.333252932146 - 08 DF(2, 4) = - 8.377654125979 - 10
DF(2, 5) = 8.512301974406 - 12
Точные значения решения и производной такие:
y1 = 1.014949833745 , y2 = - 1.004949833748
y1' = 1.489950167079 , y2' = - 4.899501670798 - 01