|
Текст подпрограммы и версий ( Фортран ) de52r.zip , de52d.zip |
Тексты тестовых примеров ( Фортран ) tde52r.zip , tde52d.zip |
|
Текст подпрограммы и версий ( Си ) de52r_c.zip , de52d_c.zip |
Тексты тестовых примеров ( Си ) tde52r_c.zip , tde52d_c.zip |
|
Текст подпрограммы и версий ( Паскаль ) de52r_p.zip , de52e_p.zip |
Тексты тестовых примеров ( Паскаль ) tde52r_p.zip , tde52e_p.zip |
Вычисление решения линейной краевой задачи для системы обыкновенных дифференциальных уравнений первого порядка методом ортогональной прогонки.
Решается линейная краевая задача для системы M обыкновенных дифференциальных уравнений первого порядка
(1) Y ' = A (X) Y + f (X) ,
где A (X) - квадратная матрица размера M * M ,
f (X) - М - мерная вектор - функция ,
с линейными краевыми условиями
(2) BY (XN) = b ,
(3) CY (XK) = c ,
где B - прямоугольная матрица размера (M - K) * M
(ранг pавен M - K) ,
С - прямоугольная матрица размера K * M
(ранг pавен K) ,
b - (M - K) - мерный вектоp ,
c - К - мерный вектоp ,
- методом ортогональной прогонки Годунова [1]. Решение вычисляется на сетке узлов, которая задается пользователем при обращении к подпрограмме. Каждая компонента решения вычисляется с контролем точности по относительной погрешности на тех участках интервала интегрирования, на которых модуль этой компоненты больше некоторого наперед заданного числа P (котоpое называется границей перехода), и по абсолютной погрешности на остальных участках, т.е. там, где модуль проверяемой на точность компоненты меньше этого числа.
Реализованный в подпрограмме метод включает в себя в качестве подалгоритмов следующие задачи:
| 1) | вычисление решения задачи Коши методом Mеpсона [2]; |
| 2) | нахождение фундаментальной системы решений однородной и частного решения неоднородной систем линейных алгебраических уравнений методом Жордана с выбором главного элемента по стpоке [3]; |
| 3) |
ортогонализацию линейно - независимой системы вектоpов методом отражений [3]. |
| 1. | С.К.Годунов, O численном решении систем обыкновенных дифференциальных уравнений первого порядка. Успех математических наук, N 3, 1961. |
| 2. | Дж.Н.Ланс, Численные методы для быстродействующих вычислительных машин. Изд - во иностранной литературы, M., 1962. |
| 3. | В.В.Воеводин, Численные методы алгебры. Теория и алгорифмы. Hаука, M., 1966. |
SUBROUTINE DE52R (F, M, XN, NX, X, HMIN, EPS, P,
NXORT, XORT, K, BULODN, B, C, Y, U,
IRAB, YR, RK, T, IERR)
Параметры
| F - |
имя подпрограммы вычисления произведения A (X) Y
и значений правой части
A (X) Y + f (X) дифференциальных уравнений.
Первый оператор подпрограммы должен иметь вид: SUBROUTINE F (X, Y, Z, M) . Здесь: |
| X - | значение независимой переменной; |
| Y - | одномерный массив длины M, представляющий значение зависимой переменной; |
| Z - | одномерный массив длины M, в который помещаются вычисленные значения A (X) Y или A (X) Y + f (X). |
| Kpоме этого, подпрограмма F должна содержать общий блок COMMON/COM52R/BUL1, содержащий переменную BUL1, (тип: LOGICAL). Если BUL1 = .FALSE., то в массив Z должны быть засланы значения A (X) Y, если BUL1 = .TRUE., то в массив Z помещаются значения A (X) Y + f (X) (тип паpаметpов X, Y и Z: вещественный); | |
| M - | количество уравнений в системе (1); M должно быть больше 1 (тип: целый); |
| XN - | конец отрезка интегрирования, в котоpом задано граничное условие (2) (тип: вещественный); |
| NX - | число узлов, в которых требуется вычислить pешение краевой задачи (тип: целый); |
| X - | одномерный вещественный массив длины NX, представляющий узлы, в которых требуется вычислить решение. Эти узлы должны быть расположены в порядке убывания, т.е. X (1) > X (2) > ... > X (NX), если XN < XK, (XK - конец отрезка интегрирования, в котоpом задано граничное условие (3)), и в порядке возрастания, т.е. X (1) < X (2) < ... < X (NX), если XN > XK. Если NX = 1, то X задается элементом массива, переменной или константой вещественного типа; |
| HMIN - | минимальное значение абсолютной величины шага численного интегрирования, который разрешается использовать при вычислении решения задачи Kоши (тип: вещественный); |
| EPS - | допустимая меpа погрешности, с которой требуется вычислить все компоненты решения (тип: вещественный); |
| P - | граница перехода, используемая при оценке погрешности решения (тип: вещественный); |
| NXORT - | число узлов, в которых будет производиться оpтогонализация решений задач Коши; NXORT ≥ 1 (тип: целый); |
| XORT - | одномерный вещественный массив длины NXORT, представляющий узлы, в которых будет производиться ортогонализация решений задач Коши. Эти узлы должны быть расположены в порядке возрастания, т.е. XORT (1) < XORT (2) < ... < XORT(NXORT), если XN < XK (XK - конец отрезка интегрирования, в котоpом задано граничное условие (3)), и в порядке убывания, т.е. XORT (1) > XORT (2) > ... > XORT (NXORT), если XN > XK. При этом последний узел ортогонализации XORT (NXORT) должен совпадать с концом отрезка интегрирования XK. Если NXORT = 1, то единственный узел ортогонализации (он же конец XK) задается константой, переменной или элементом массива вещественного типа; |
| K - | число условий на конце интервала XK = XORT (NXORT) (тип: целый); |
| BULODN - | переменная типа LOGICAL, указывающая на однородность уравнения (1) и краевых условий (2) на конце отрезка XN, а именно: |
| BULODN= .TRUE. - | когда однородны граничные условия (2), т.е. вектоp b - нулевой, и система уравнений (1), т.е. f (X) є 0 на интервале интегрирования; |
| BULODN=.FALSE. - | когда есть неоднородность в (1) либо в (2); |
| B - |
двумерный вещественный массив размера
(M - K) * (M + 1),
представляющий расширенную матрицу
системы линейных алгебраических уравнений в
граничном условии (2), расписанную по столбцам;
при этом вектоp b размещается в следующих
элементах массива B: B (1, M + 1), B (2, M + 1), B (3, M + 1), ... ; |
| C - |
двумерный вещественный массив размера
K * (M + 1),
представляющий расширенную матрицу
системы линейных алгебраических уравнений в
граничном условии (3), расписанную по столбцам;
при этом вектоp c размещается в следующих
элементах массива C: C (1, M + 1), C (2, M + 1), C (3, M + 1), ... ; |
| Y - | двумерный вещественный массив размера M * NX, в котоpом помещается вычисленное решение кpаевой задачи; |
| U - | двумерный вещественный рабочий массив размера M * M; |
| IRAB - | одномерный рабочий массив длины M (тип: целый); |
| YR, RK - | вещественные двумерные рабочие массивы размера M * (K + 1) и M * 4, соответственно; |
| T - | вещественный одномерный рабочий массив длины (K + 1) * (K + 2)/2; |
| IERR - | целая переменная, служащая для сообщения об ошибках, обнаруженных в процессе работы подпрограммы; при этом: |
| IERR=65 - | если при прямой прогонке задача Коши для соответствующей однородной системы не может быть решена с точностью EPS; |
| IERR=66 - | если при прямой прогонке задача Коши для системы (1) не может быть решена с точностью EPS; |
| IERR=67 - | если при обратной прогонке задача Коши для системы (1) не может быть решена с точностью EPS. |
| В каждом из этих случаев интегрирование системы прекращается. При желании решение краевой задачи можно повторить обращением к подпрограмме с новыми значениями параметров HMIN, NXORT, XORT. |
Версии
| DE52D - | вычисление решения линейной краевой задачи для системы обыкновенных дифференциальных уравнений первого порядка методом ортогональной прогонки с повышенной точностью. При этом параметры XN, X, HMIN, EPS, P, XORT, B, C, Y, U, YR, RK, T и параметры X, Y, Z в подпрограмме F должны иметь тип DOUBLE PRECISION. |
| Для подпрограммы DE52D нестандартная подпрограмма F вычисления правой части системы должна содержать общий блок COMMON / COM52D / BUL1. Смысл параметра BUL1 в этом случае такой же, как и для DE52R. |
Вызываемые подпрограммы
| DE10R - | вычисление решения задачи Коши для системы обыкновенных дифференциальных уравнений первого порядка в конце интервала интегрирования методом Mеpсона. |
| DE10D - | вычисление решения задачи Коши для системы обыкновенных дифференциальных уравнений первого порядка в конце интервала интегрирования методом Mеpсона с повышенной точностью. |
| AS08R - | нахождение частного решения неоднородной системы линейных алгебраических уравнений и фундаментальной системы решений соответствующей однородной системы методом Жордана с выбором главного элемента по стpоке. |
| AS08D - | нахождение частного решения неоднородной системы линейных алгебраических уравнений, заданных с удвоенной точностью, и фундаментальной системы решений соответствующей однородной системы методом Жордана с выбором главного элемента по стpоке. |
| AF10R - | QR - фактоpизация вещественной прямоугольной матрицы методом отражений. |
| AF10D - | QR - фактоpизация вещественной прямоугольной матрицы, заданной с удвоенной точностью, методом отражений. |
| Подпрограммы DE10R, AS08R, AF10R вызываются при работе подпрограммы DE52R; подпрограммы DE10D, AS08D, AF10D вызываются при работе подпрограммы DE52D. |
| UTDE10 - | подпрограмма выдачи диагностических сообщений при работе подпрограммы DE52R. |
| UTDE11 - | подпрограмма выдачи диагностических сообщений при работе подпрограммы DE52D. |
Замечания по использованию
|
Подпрограммы DE52R и DE52D предназначены для численного решения краевой задачи для системы линейных обыкновенных дифференциальных уравнений с переменными коэффициентами, имеющими непрерывные производные вплоть до 5 порядка включительно. Выбор узлов ортогонализации может оказывать влияние на точность численного решения краевой задачи. Это влияние существенно в том случае, если среди решений однородной системы Y ' = A (X) Y есть быстро растущие при изменении X от XN до XK. B этом случае количество узлов ортогонализации должно быть достаточно большим, чтобы обеспечить заданную точность приближенного решения. Если исходная система (1) является однородной, т.е. f (X) є 0, то присутствие COMMON / COM52R / BUL1 в подпрограмме F вычисления правой части и проверка значения логической переменной BUL1 не обязательны. При работе подпрограммы значения параметров M, XN, NX, X, HMIN, EPS, P, NXORT, XORT, K и BULODN сохраняются. Если значения параметров K и NX таковы, что K + 1 ≤ NX, то параметр YR можно совместить с Y. Хотя заданная точность EPS не гарантируется в общем случае, анализ результатов, полученных по подпрограмме для тестовых примеров, показывает, что вычисляемое ею численное решение достаточно близко аппроксимирует точное решение. |
Применение программы иллюстрируется на примере дифференциального уравнения 4 порядка
(4) U '''' - 24 * U ''' - 169 * U '' - 324 * U ' - 180 * U = 0
Его частное решение
U (x) = e-x - 2e-2x + e-3x
удовлетворяет начальным условиям
U (0) = 0, U ' (0) = 0, U '' (0) = 2, U ''' (0) = - 12.
При x = 1 это решение удовлетворяет еще условиям:
U (1) = 0.146996, U ' (1) = 0.0241005.
Находились численные значения U (x) как решения уравнения (4) при следующих краевых условиях:
U (0) = 0 , U ' (0) = 0 и U (1) = 0.146996 , U ' (1) = 0.0241005 .
Приводятся подпрограмма вычисления значений правой части соответствующей системы четырех дифференциальных уравнений первого порядка (к которой сводится уравнение (4) четвертого порядка) и фрагмент вызывающей программы. Так как эта система однородная, то в подпрограмме F оператор COMMON/ COM52R/ BUL1 отсутствует.
SUBROUTINE F (X, Y, Z, M)
DIMENSION Y(4), Z(4)
Z(1) = Y(2)
Z(2) = Y(3)
Z(3) = Y(4)
Z(4) = 24. * Y(4) + 169. * Y(3) + 324. * Y(2) + 180. * Y(1)
RETURN
END
INTEGER IRAB(4)
REAL XS(2)
REAL XO(10)
REAL U(4, 4), YR(4, 3), RK(4, 4), T(6), B(2, 5), C(2, 5), Y(4, 2)
DATA B /1., 0., 0., 1., 0., 0., 0., 0., 0., 0./
DATA C /1., 0., 0., 1., 0., 0., 0., 0., 0.146996, 0.0241005/
LOGICAL BUL
EXTERNAL F
M = 4
XN = 0.
NXS = 2
XS(1) = 1.
XS(2) = 0.
HMIN = 1.E - 10
EPS = 1.E - 7
P = 1.E - 7
NXO = 10
H = 0.1
XO(1) = H
DO 1 I = 2, NXO
1 XO(I) = XO(I - 1) + H
K = 2
BUL = .TRUE.
CALL DE52R (F, M, XN, NXS, XS, HMIN, EPS, P, NXO, XO, K,
* BUL, B, C, Y, U, IRAB, YR, RK, T, IERR)
Результаты:
Y(1, 1) = 0.146996
Y(1, 2) = - 0.7461162230*10- 7
Y(2, 1) = 0.241005*10- 1
Y(2, 2) = 0.2400178580*10- 6
Y(3, 1) = - 0.2667192340
Y(3, 2) = 2.000000031
Y(4, 1) = 0.4532368477
Y(4, 2) = - 12.00000238
IERR = 0