|
Текст подпрограммы и версий ( Фортран ) de16r.zip , de16d.zip |
Тексты тестовых примеров ( Фортран ) tde16r.zip , tde16d.zip |
|
Текст подпрограммы и версий ( Си ) de16r_c.zip , de16d_c.zip |
Тексты тестовых примеров ( Си ) tde16r_c.zip , tde16d_c.zip |
|
Текст подпрограммы и версий ( Паскаль ) de16r_p.zip , de16e_p.zip |
Тексты тестовых примеров ( Паскаль ) tde16r_p.zip , tde16e_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 ) ,
методом Рунге - Кутта - Фельберга пятого порядка точности. Решение вычисляется в одной точке XK, которая является концом интервала интегрирования. Каждая компонента решения вычисляется с контролем точности по относительной погрешности на тех участках интервала интегрирования, на которых модуль этой компоненты больше некоторого наперед заданного числа Р, и по абсолютной погрешности на остальных участках, т.е. там, где модуль компоненты меньше этого числа. B качестве абсолютной погрешности решения используется оценка главного члена асимптотического разложения погрешности метода на одном шаге, получаемая при вычитании двух выражений, представляющих приближенные значения решения пятого и четвертого порядка точности. При этом на каждом шаге интегрирования для определения решения и его погрешности используются всего шесть вычислений правой части системы.
Дж.Форсайт, М.Малькольм, К.Моулер, Машинные методы математических вычислений, "Мир", M., 1980.
SUBROUTINE DE16R (F, M, XN, YN, XK, HMIN, EPS,
P, H, Y, RAB, IERR)
Параметры
| F - |
имя подпрограммы вычисления значений правой
части дифференциального уравнения. Первый
оператоp подпрограммы должен иметь вид: SUBROUTINE F (X, Y , DY, M). Здесь: X, Y - значения независимой и зависимой переменных, соответственно. Вычисленное значение правой части должно быть помещено в DY. B случае системы уравнений, т.е. когда M ≠ 1 , параметры Y и DY представляют массивы длины M (тип параметров X, Y и DY: вещественный); |
| M - | количество уравнений в системе (тип: целый); |
| XN, YN - | начальные значения аргумента и решения; в случае системы уравнений (т.е. M ≠ 1) YN представляет одномерный массив длины M (тип: вещественный); |
| XK - | значение аргумента, при котоpом требуется вычислить решение задачи Коши (конец интервала интегрирования); XK может быть больше, меньше или pавно XN (тип: вещественный); |
| HMIN - | минимальное значение абсолютной величины шага, котоpое разрешается использовать при интегрировании данной системы уравнений (тип: вещественный); |
| EPS - | допустимая меpа погрешности, с которой требуется вычислить все компоненты решения (тип: вещественный); |
| P - | граница перехода, используемая при оценке погрешности решения (тип: вещественный); |
| H - | вещественная переменная, содержащая начальное значение шага интегрирования; может задаваться с учетом направления интегрирования, т.е. положительным, если XN < XK, отрицательным, если XN > XK, или без всякого учета в виде абсолютной величины; на выходе из подпрограммы содержит значение последнего шага интегрирования; |
| Y - | искомое решение задачи Коши, вычисленное подпрограммой при значении аргумента XK; для системы уравнений (когда M ≠ 1) задается одномерным массивом длины M; в случае совпадения значений параметров XN и XK значение Y полагается равным начальному значению YN (тип: вещественный); |
| RAB - | одномерный рабочий массив вещественного типа длиной 6*M + 1; |
| IERR - | целая переменная, значение которой в результате работы подпрограммы полагается равным 65, если какая - нибудь компонента решения не может быть вычислена с требуемой точностью EPS. B этом случае интегрирование системы прекращается; при желании интегрирование системы можно повторить обращением к подпрограмме с новыми значениями параметров HMIN и H. |
Версии
| DE16D - | вычисление решения задачи Коши для системы обыкновенных дифференциальных уравнений первого порядка в конце интервала интегрирования методом Рунге - Кутта - Фельберга с удвоенным числом значащих цифр. При этом параметры XN, YN, XK, HMIN, EPS, P, H, Y, RAB и параметры X, Y и DY в подпрограмме F должны иметь тип DOUBLE PRECISION. |
Вызываемые подпрограммы
|
DE15R - DE15D | выполнение одного шага численного интегрирования системы обыкновенных дифференциальных уравнений первого порядка методом Рунге - Кутта - Фельберга пятого порядка точности. |
|
UTDE16 - UTDE17 | подпрограммы выдачи диагностических сообщений. |
| Подпрограммы DE15R и UTDE16 вызываются при работе подпрограммы DE16R, а подпрограммы DE15D и UTDE17 - при работе подпрограммы DE16D. |
Замечания по использованию
|
B общем случае заданая точность не гарантируется. При работе подпрограммы значения параметров M, XN, YN, XK, HMIN, EPS, P сохраняются. Если после работы подпрограммы нет необходимости иметь начальное значение решения YN, то параметры YN и Y при обращении к ней можно совместить. При работе подпрограммы счета правой части F значения параметров X, Y и M не должны изменяться. |
y1' = -y2 - 0.1 x - 0.9
y2' = y1 - 0.1 x - 1.1 0 ≤ x ≤ π
y1 (0) = 1 , y2 (0) = - 2
Точное решение системы:
y1 = 0.1 x + sin x + 1.
y2 = - 0.1 x - cos x - 1.
Приводятся подпрограммы вычисления значений правой части и фрагмент вызывающей программы, а также результаты счета.
SUBROUTINE FFELB (X, Y, Z, M)
DIMENSION Y(2), Z(2)
R = 0.1*X
Z(1) = -Y(2) - R - 0.9
Z(2) = Y(1) - R - 1.1
RETURN
END
DIMENSION YN(2), Y(2), RAB(13)
EXTERNAL FFELB
M = 2
XN = 0.
YN(1) = 1.
YN(2) = -2.
XK = 3.14159265358979 D0
HMIN = 1.E-14
EPS = 1.E-8
P = 100.
H = 0.01
CALL DE16R (FFELB, M, XN, YN, XK, HMIN, EPS, P, H, Y, RAB,
* IERR)
Результаты:
IERR = 0
Y(1) = 1.314159265929 + 00
Y(2) = -3.141592563520 - 01