Документ взят из кэша поисковой машины. Адрес
оригинального документа
: http://num-anal.srcc.msu.su/lib_na/cat/de/de03r.htm
Дата изменения: Mon Nov 30 15:27:40 2015 Дата индексирования: Sun Apr 10 00:10:53 2016 Кодировка: Windows-1251 |
Текст подпрограммы и версий ( Фортран ) de03r.zip , de03d.zip |
Тексты тестовых примеров ( Фортран ) tde03r.zip , tde03d.zip |
Текст подпрограммы и версий ( Си ) de03r_c.zip , de03d_c.zip |
Тексты тестовых примеров ( Си ) tde03r_c.zip , tde03d_c.zip |
Текст подпрограммы и версий ( Паскаль ) de03r_p.zip , de03e_p.zip |
Тексты тестовых примеров ( Паскаль ) tde03r_p.zip , tde03e_p.zip |
Вычисление решения задачи Коши для квазилинейной устойчивой системы обыкновенных дифференциальных уравнений первого порядка с большой константой Липшица методом Лоусона.
Решается задача Коши для квазилинейной системы M обыкновенных дифференциальных уравнений первого порядка:
(1) Y ' (X) = A * Y(X) + U(X, Y) , Y = ( y1,..., yM ) , A = ( ai j ) , i, j = 1, ..., M , U(X, Y) = ( U1(X, Y), ... , UM(X, Y) ) ,
с начальными условиями, заданными в точке XN:
Y(XN) = YN , YN = ( y1, ... , yM ) .
где A - постоянная числовая матрица.
Предполагается, что среди характеристических корней матрицы A имеются большие по модулю корни, а константа Липшица для функции U (X, Y), т.е. константа L из условия Липшица
| Ui ( x, y1(1),..., yM(1) ) - Ui ( x, y1(2), ... , yM(2) ) | ≤ M ≤ L ∑ | yj(1) - yj(2) | , j=1
независящая от i, x, y(1), y(2) , невелика. Также предполагается, что нелинейный член U (x, y) является достаточно малым. Решение вычисляется в одной точке ХК, которая является концом интервала интегрирования. Для интегрирования системы применяется метод Лоусона.
Метод Лоусона является одношаговым и заключается в следующем. Допустим, что искомое решение системы (1) уже вычислено в некоторой точке x = xn интервала интегрирования, т.е. известно yn ≈ y(xn). Для отыскания решения Y(xn + 1) = Y(xn + H) в следующем узле xn + 1 = xn + H выполняются такие действия. Исходная система уравнений с помощью замены искомой функции Y (x) на xn ≤ x ≤ xn + H по формуле
y(x) = exp [ ( x - xn ) A ] Z(x) ,
преобразуется в систему уравнений относительно новой неизвестной функции Z (X):
(2) Z ' (x) = U1(x, z) = exp [ - ( x - xn ) A ] U( x, exp [( x - xn ) A ] Z(x) )
Данное преобразование выполняется самой подпрограммой. Характеристические корни матрицы Якоби
∂U1 / ∂Z = exp [ - ( x - xn ) A ] (∂U / ∂y) exp [ ( x - xn ) A ]
являются характеристическими корнями матрицы ∂U / ∂y и в силу малости константы Липшица функции U (x, y) невелики. Поэтому система (2) не жесткая и может быть решена традиционными методами численного интегрирования. Данная подпрограмма решает ее методом Рунге - Кутта 4 - ого порядка, при этом одновременно с решением (2) производится обратное преобразование от функции Z (x) к функции y (x) .
Все компоненты решения вычисляются с контролем точности по мере погрешности, который заключается в следующем. Если некоторая компонента приближенного решения по абсолютной величине меньше некоторой наперед заданной константы P, то контроль точности для этой компоненты ведется по относительной погрешности, иначе - по абсолютной. Абсолютная погрешность приближенного решения оценивается по правилу Рунге.
J.Douglas Lowson, Generalized Runge - Kutta processes for stable systems with large lipshitz constants // SIAM Journal on Numerical Analisys - 1967 - Vol 4, No 3.
SUBROUTINE DE03R (FU, M, XN, YN, XK, A, HMIN, EPS, P, H, Y, R, IERR)
Параметры
FU - |
подпрограмма вычисления функции U (x, y) в правой
части системы. Первый оператор подпрограммы
должен иметь вид: SUBROUTINE FU (X, Y, U, M). Здесь Y и U - одномерные массивы длины M. В массив U помещается значение функции U (x, y), вычисленное при значении аргументов X и Y (тип параметров X, Y, U: вещественный); |
M - | количество уравнений в системе (тип: целый); |
XN, YN - | начальные значения аргумента и решения; в случае системы уравнений (т.е. M ≠ 1) YN представляет одномерный массив длины M (тип: вещественный); |
XK - | значение аргумента, при котором требуется вычислить решение задачи Коши (конец интервала интегрирования); XK может быть больше, меньше или равно XN (тип: вещественный); |
A - | вещественный двумерный массив размера M*M, содержащий элементы матрицы A; |
HMIN - | минимальное значение абсолютной величины шага, которое разрешается использовать при интегрировании данной системы уравнений (тип: вещественный); |
EPS - | допустимая мера погрешности, с которой требуется вычислить все компоненты решения (тип: вещественный); |
P - | граница перехода, используемая при оценке меры погрешности решения (тип: вещественный); |
H - | вещественная переменная, содержащая начальное значение шага интегрирования; может задаваться с учетом направления интегрирования, т.е. положительным, если XN < XK, отрицательным, если XN > XK, или без всякого учета в виде абсолютной величины; |
Y - | искомое решение задачи Коши, вычисленное подпрограммой при значении аргумента XK; для системы уравнений (когда M ≠ 1) задается одномерным массивом длины M. В случае совпадения значений параметров XN и XK значение Y полагается равным начальному значению YN (тип: вещественный); |
R - | вещественный одномерный рабочий массив длины 4*M*M+7*M+1; |
IERR - | целая переменная, значение которой в результате работы подпрограммы полагается равным 65, если какая - нибудь компонента решения не может быть вычислена с требуемой точностью EPS. В этом случае интегрирование системы прекращается; при желании интегрирование системы можно повторить обращением к подпрограмме с новыми значениями параметров HMIN и H. |
Версии
DE03D - | вычисление решения задачи Коши для квазилинейной системы обыкновенных дифференциальных уравнений первого порядка с большой константой Липшица в конце интервала интегрирования методом Лоусона с удвоенным числом значащих цифр. При этом параметры XN, YN, XK, A, HMIN, EPS, P, H, Y, R и параметры X, Y, R в подпрограмме FU должны иметь тип DOUBLE PRECISION. |
Вызываемые подпрограммы
DE02R - DE02D | выполнение одного шага численного интегрирования квазилинейной системы обыкновенных дифференциальных уравнений первого порядка с большой константой Липшица методом Лоусона. |
UTDE20 - UTDE21 | подпрограмма выдачи диагностических сообщений; Подпрограмма DE02R, UTDE20 вызывается при работе подпрограммы DE03R, а подпрограммы DE02D, UTDE21 - при работе DE03D. |
Замечания по использованию
Данная подпрограмма предназначена для интегрирования квазилинейных систем, имеющих малый нелинейный член U (x, y). В общем случае заданная точность не гарантируется. При работе подпрограммы значения параметров M, A, XN, YN, XK, HMIN, EPS, P сохраняются. При работе подпрограммы FU значения параметров X, Y и M не должны изменяться. Если после работы подпрограммы нет необходимости иметь начальное значение решения YN, то параметры YN и Y при обращении к ней можно совместить. При этом следует иметь в виду, что в случае аварийного выхода из подпрограммы, т.е. со значением IERR = 65, значение параметра YN будет испорчено. Подпрограммы DE03R и DE03D предназначены также для решения задачи Коши для жестких дифференциальных уравнений (1). |
y1' = - 500 y1 + y2 + sin ( y1 + y2 ) ( 1 + e - x/2 ) , y1(0) = 50 y2' = - 1000 y1 + y2 + sin ( y1 - y2 ) ( 1 + e - x/2 ) , y2(0) = 50 0 ≤ x ≤ 5
Приводятся подпрограмма вычисления функций
U1 (x, y) = sin ( y1 + y2 ) ( 1 + e - x/2 ) , U2 (x, y) = sin ( y1 - y2 ) ( 1 + e - x/2 )
из правой части системы, фрагмент вызывающей программы и результаты счета.
SUBROUTINE FU (X, Y, U, M) DIMENSION Y(2), U(2) T = 1. + EXP(-X/2.) U(1) = T*SIN(Y(1) + Y(2)) U(2) = T*SIN(Y(1) - Y(2)) RETURN END DIMENSION A(2, 2), Y(2), YN(2), R(31) EXTERNAL FU M = 2 XN = 0. YN(1) = 5.E1 YN(2) = 5.E1 HMIN = 1.E-10 EPS = 1.E-8 P = 1.E2 XK = 5. H = 0.01 A(1, 1) = -5.E2 A(1, 2) = 1. A(2, 1) = -1.E3 A(2, 2) = 1. CALL DE03R (FU, M, XN, YN, XK, A, HMIN, EPS, P, H, Y, R, IERR) Результаты: Y(1) Y(2) H -3.508894600484-04 -8.340807464992-02 1.250000000001-03 IERR = 0