|
Текст подпрограммы и версий ( Фортран ) de59r.zip de59d.zip de60r.zip de60d.zip de57r.zip de57d.zip de58r.zip de58d.zip |
Тексты тестовых примеров ( Фортран ) tde59r.zip tde59d.zip tde60r.zip tde60d.zip tde57r.zip tde57d.zip tde58r.zip tde58d.zip |
|
Текст подпрограммы и версий ( Си ) de59r_c.zip de59d_c.zip de60r_c.zip de60d_c.zip de57r_c.zip de57d_c.zip de58r_c.zip de58d_c.zip |
Тексты тестовых примеров ( Си ) tde59r_c.zip tde59d_c.zip tde60r_c.zip tde60d_c.zip tde57r_c.zip tde57d_c.zip tde58r_c.zip tde58d_c.zip |
|
Текст подпрограммы и версий ( Паскаль ) de59r_p.zip de59e_p.zip de60r_p.zip de60e_p.zip de57r_p.zip de57e_p.zip de58r_p.zip de58e_p.zip |
Тексты тестовых примеров ( Паскаль ) tde59r_p.zip tde59e_p.zip tde60r_p.zip tde60e_p.zip tde57r_p.zip tde57e_p.zip tde58r_p.zip tde58e_p.zip |
Вычисление решения двухточечной краевой задачи для обыкноновенного линейного самосопряженного дифференциального уравнения второго порядка с разрывными коэффициентами на равномерной сетке с помощью однородной консервативной разностной схемы второго порядка точности.
Решается двухточечная краевая задача для линейного самосопряженного уравнения второго порядка
(1) d (k(x) dy/dx) / dx - g(x) y = - r(x) , xN < x < xK
с разрывными коэффициентами
(2) k(x) , g(x) , r(x) ∈ Q( 2 ) [xN, xK]
( Q ( 2 ) [xN, xK] - класс функций, кусочно - непрерывных на [xN, xK] вместе со вторыми производными ), удовлетворяющими неpавенствам
(3) k(x) ≥ c > 0 , g(x) ≥ 0 .
Hа концах отрезка [xN, xK] задаются произвольные комбинации следующих краевых условий
при x = xN (4) 1) y( xN ) = GN , (5) 2) k dy/dx - ( σN y - μN ) = 0 , σN ≥ 0 , при x = xK (6) 1) y( xK ) = GK , (7) 2) - k dy/dx - ( σK y - μK ) = 0 , σK ≥ 0 ;
при σN > 0 условие (5) (а при σK > 0 условие (7)) является условием третьего рода, при σN = 0 (или при σK = 0) - условием второго рода.
Решение вычисляется на равномерной сетке узлов, заданной на [xN, xK], с помощью однородной консервативной разностной схемы, имеющей на указанном классе коэффициентов (2) второй порядок точности, при этом предполагается, что точки разрыва коэффициентов k (x), g (x) и r (x) уравнения (1) совпадают с узлами этой сетки.
Самарский A.A., Теория разностных схем, "Hаука", M., 1977.
SUBROUTINE DE59R (CK, CG, CR, NX, XN, XK, GN1, GN, GK1,
GK, Y, ALP, BET, IERR)
Параметры
| NX - | число узлов равномерной сетки, включая концы отрезка интегрирования, в которых требуется вычислить решение краевой задачи; NX > 1 (тип: целый); |
| XN, XK - | концы отрезка интегрирования, в которых задаются граничные условия (4) или (5) и (6) или (7), соответственно; XK > XN (тип: вещественный); |
| CK - |
одномерный вещественный массив длины NX - 1,
содержащий на входе в подпрограмму значения
коэффициента k (x) уравнения, а именно:
I - й элемент CK (I) массива содержит значение
функции k (x) в точке x = XN + (I - 1/2) * (XK - XN)/(NX - 1); функция k (x) должна принимать положительные значения; |
| CG, CR - |
одномерные вещественные массивы длины 2*NX - 2,
содержащие на входе в подпрограмму значения коэффициентов
g (x) и r (x), соответственно, а именно: CG (1) содержит предел справа функции g (x) в начале XN отрезка интегрирования, (2*I) - й и (2*I + 1) - й элементы CG (2*I), CG (2*I + 1) массива CG содержат, соответственно, пределы слева и справа функции g (x) в узле сетки x = XN + I*(XK - XN)/(NX - 1) (I = 1, 2, ..., NX - 2), (2*NХ - 2) - й элемент CG (2*NX - 2) массива CG содержит предел слева функции g (x) в конце XK отрезка интегрирования, элементы массива CR содержат аналогичные значения, вычисленные для функции r (x); |
| GN1 - | логическая переменная (или логическая константа), имеющая на входе в подпрограмму значение .TRUE., если в начале XN отрезка интегрирования задано граничное условие первого рода (4), и .FALSE. в противном случае, т.е. когда в XN задано условие 2 - го или 3 - его рода (5); |
| GN - |
если GN1 = .TRUE., то GN представляет значение
решения в начале XN отрезка интегрирования; если GN1 = .FALSE., то GN является именем одномерного массива длины 2, первый элемент GN (1) которого содержит значение σN, а второй элемент GN (2) - значение μN, входящие в граничное условие (5) (тип: вещественный); |
| GK1 - | логическая переменная (или логическая константа), имеющая на входе в подпрограмму значение .TRUE., если в конце XK отрезка интегрирования задано граничное условие первого рода (6), и .FALSE. в противном случае, т.е. когда в XK задано условие 2 - го (σK = 0) или 3 - го (σK > 0) рода (7); |
| GK - |
если GK1 = .TRUE., то GK представляет значение
решения в конце отрезка интегрирования; если GK1 = .FALSE., то GK является именем одномерного массива длины 2, первый элемент GK (1) которого содержит значение σK, а второй элемент GK (2) - значение μK, входящие в граничное условие (7) (тип: вещественный); |
| Y - |
одномерный вещественный массив длины NX,
содержащий на выходе из подпрограммы решение
краевой задачи, вычисленное на равномерной
сетке; при этой I - й элемент Y (I) этого
массива содержит приближенное значение решения в узле XN + (I - 1)*(XK - XN)/(NX - 1); |
|
ALP - BET | одномерные вещественные рабочие массивы двойной точности длины NX - 1; |
| IERR - | целая переменная, служащая для сообщения об ошибках, обнаруженных в процессе работы подпрограммы; при этом: |
| IERR=65 - | если в граничном условии (5) коэффициент σN < 0; |
| IERR=66 - | если в граничном условии (7) коэффициент σK < 0; |
| IERR=67 - | если какое - нибудь значение коэффициента k (x), заданного в массиве CK, меньше или pавно 0; |
| IERR=68 - | если какое - нибудь значение коэффициента g (x), заданного в массиве CG, меньше 0; |
| IERR=69 - | если значение NX ≤ 1; |
| IERR=70 - | если XK ≤ XN; |
| во всех указанных случаях интегрирование уравнения прекращается. |
Версии
| DE60R - | вычисление решения двухточечной краевой задачи для обыкновенного линейного самосопряженного дифференциального уравнения второго порядка с разрывными коэффициентами на неравномерной сетке. |
| DE57R - | вычисление решения двухточечной краевой задачи для обыкновенного линейного самосопряженного дифференциального уравнения второго порядка с непрерывными коэффициентами на равномерной сетке. |
| DE58R - | вычисление решения двухточечной краевой задачи для обыкновенного линейного самосопряженного дифференциального уравнения второго порядка с непрерывными коэффициентами на неравномерной сетке. |
| DE57D - | то же, что и DE57R, при этом все арифметические операции с вещественными числами выполняются с удвоенным числом значащих цифр. |
| DE58D - | то же, что и DE58R, при этом все арифметические операции с вещественными числами выполняются с удвоенным числом значащих цифр. |
| DE59D - | то же, что и DE59R, при этом все арифметические операции с вещественными числами выполняются с удвоенным числом значащих цифр. |
| DE60D - | то же, что и DE60R, при этом все арифметические операции с вещественными числами выполняются с удвоенным числом значащих цифр. |
|
Первые операторы этих подпрограмм имеют следующий вид: SUBROUTINE DE60R (CK, CG, CR, NX, X, GN1, GN, GK1, GK, Y, ALP, BET, IERR) SUBROUTINE DE57R (FK, FG, FR, NX, XN, XK, GN1, GN, GK1, GK, Y, ALP, BET, IERR) SUBROUTINE DE58R (FK, FG, FR, NX, X, GN1, GN, GK1, GK, Y, ALP, BET, IERR) Параметр X в подпрограммах DE60R и DE58R - имя одномерного вещественного массива длины NX, в последовательных элементах которого, начиная с первого, размещаются в возрастающем порядке узлы неравномерной сетки, причем начало и конец отрезка интегрирования, в которых заданы граничные условия (4) - (7), хранятся, соответственно, в первом и последнем элементах X (1), X (NX) этого массива. Параметры FK, FG, FR в подпрограммах DE57R и DE58R - имена вещественных подпрограмм - функций REAL FUNCTION FK (X), REAL FUNCTION FG (X), REAL FUNCTION FR (X), вычисляющих значения коэффициентов уравнения k (x), g (x) и r (x), соответственно, при вещественном значении аргумента x = X. Остальные параметры указанных версий имеют такой же смысл, что и одноименные параметры подпрограммы DE59R. При этом значение параметра IERR = 70 в подпрограммах DE60R и DE58R указывает на то, что узлы сетки расположены в массиве X не в возрастающем порядке, значения IERR = 67, 68 в п/п DE58R и DE57R - что значения коэффициентов k (x) и g (x) уравнения, вычисленные с помощью подпрограмм - функций FK и FG, не удовлетворяют требуемое неpавенство (3). B подпрограмме DE60R в массиве CK размещаются значения коэффициента k (x) уравнения, вычисленные в "потоковых" точках, т.е. в середине между узлами сетки, а в массивах CG и CR - пределы слева и справа функций g (x) и r (x) в узлах неравномерной сетки. Kpоме того, каждая из подпрограмм DE57R, DE58R, DE59R и DE60R имеет версию DE57D, DE58D, DE59D и DE60D, соответственно, с теми же самыми эффектом и списком параметpов, выполняющую все операции над действительными числами с удвоенным числом значащих цифр. При этом параметры FK, FG, FR (в подпрограммах DE57D,
DE58D), являются именами подпрограмм - функций |
|
CK, CG - CR | (в подпрограммах DE59D, DE60D), |
| XN, XK - | (в подпрограммах DE57D, DE59D), |
| X - | (в подпрограммах DE58D, DE60D), |
|
GN, GK - Y, ALP BET | (во всех подпрограммах) |
|
должны иметь тип DOUBLE PRECISION. Тип остальных параметров не изменяется. В версиях подпрограмм DE57R, DE58R, DE59R и DE60R, формальные параметры ALP и BET имеют тип DOUBLE PRECISION. |
Вызываемые подпрограммы
| UTDE14 - | подпрограмма выдачи диагностических сообщений при работе подпрограмм DE57R, DE58R, DE59R и DE60R. |
| UTDE15 - | подпрограмма выдачи диагностических сообщений при работе подпрограмм DE57D, DE58D, DE59D и DE60D. |
Замечания по использованию
|
Значение параметра NX, задаваемое при обращении к подпрограмме, должно быть достаточно большим, чтобы получить хорошее приближение численного решения к точному решению задачи. При интегрировании дифференциального уравнения (1) с непрерывными коэффициентами с помощью подпрограмм DE57R, DE58R, DE57D, DE58D значения коэффициента k (x) вычисляются в "потоковых" точках, т.е. в середине между узлами, а значения коэффициентов g (x) и r (x) - в узлах сетки последовательно слева направо (начало и конец XN, XK отрезка интегрирования также считаются узлами сетки). Поэтому эти продпрограммы можно использовать также и в случае, когда какой - нибудь из коэффициентов k (x), g (x) и r (x) (или, быть может, все) задан в виде массива значений в указанных точках. Эти значения должны быть расположены в последовательных элементах массива, начиная с первого. При этом: если в XN задано граничное условие 2 - го или 3 - го рода, то значения коэффициентов g (x) и r (x) должны быть заданы, начиная с первого узла XN сетки, а если 1 - го рода, то начиная со второго узла XN + (XK - XN)/(NX - 1); если в конце промежутка XK задано условие 1 - го рода, то значения коэффициентов g (x) и r (x) используются только до (NХ - 1) - го узла включительно. B этом случае каждое обращение к соответствующей подпрограмме - функции вычисления коэффициента сводится к выборке значения очередного элемента массива (см. примеp использования N 3). Значения параметров CK, CG, CR, NX, XK, X, GN1, GN, GK1, GK сохраняются. |
1). d (k(x) dy/dx) / dx = 0 , y(0) = 1 , y(1) = 0 ,
k(x) = 1 , x ≤ 0.55 ,
k(x) = 10 , x > 0.55 ;
Точное решение:
y = 1 - x/0.595 , x ≤ 0.55 , (1 - x) / 5.95 , x > 0.55 .
DOUBLE PRECISION ALP(20), BET(20)
REAL Y(21), CK(20), CG(40), CR(40)
LOGICAL GN1, GK1
NX = 21
XN = 0.
XK = 1.
GN1 = .TRUE.
GN = 1.
GK1 = .TRUE.
GK = 0.
DO 10 I = 1, 40
CG(I) = 0.
10 CR(I) = 0.
DO 20 I = 1, 10
20 CK(I) = 1.
DO 30 I = 11, 20
30 CK(I) = 10.
CALL DE59R (CK, CG, CR, NX, XN, XK, GN1, GN, GK1, GK, Y,
* ALP, BET, IERR)
Результаты:
IERR = 0
Y(1) = 1.00000000000 + 00 Y(11) = 9.09090909100 - 02
Y(2) = 9.09090909089 - 01 Y(12) = 8.18181818191 - 02
Y(3) = 8.18181818179 - 01 Y(13) = 7.27272727281 - 02
Y(4) = 7.27272727271 - 01 Y(14) = 6.36363636372 - 02
Y(5) = 6.36363636362 - 01 Y(15) = 5.45454545463 - 02
Y(6) = 5.45454545453 - 01 Y(16) = 4.54545454552 - 02
Y(7) = 4.54545454545 - 01 Y(17) = 3.63636363643 - 02
Y(8) = 3.63636363636 - 01 Y(18) = 2.72727272732 - 02
Y(9) = 2.72727272727 - 01 Y(19) = 1.81818181822 - 02
Y(10) = 1.81818181818 - 01 Y(20) = 9.09090909109 - 03
Абсолютная погрешность во всех узлах сетки pавна 0.00000000000
2). d ((2 + x) dy/dx) / dx - xy = - 4sin x - 4x sin x + 2cos x , 0 < x < π / 2 ;
2 dy/dx - (1000y + 4) = 0 при x = 0 ,
y = 2 при x = π / 2 .
Решение вычисляется на неравномерной сетке.
REAL FUNCTION FK (X)
FK = 2. + X
RETURN
END
REAL FUNCTION FG (X)
FG = X
RETURN
END
REAL FUNCTION FR (X)
FR = 4. * SIN(X) * (1. + X) - 2. * COS(X)
RETURN
END
DOUBLE PRECISION ALP(282), BET(282)
REAL Y(283), X(283), GN(2)
LOGICAL GN1, GK1
EXTERNAL FK, FG, FR
NX = 283
XK = 1.57079632679
H = 0.001
X(1) = 0.
X(282) = XK - 0.02
X(283) = XK
J = 0
DO 1 K = 1, 28
DO 1 I = 1, 10
J = J + 1
1 X(J + 1) = X(J) + FLOAT(I) * H
GN1 = .FALSE.
GN(1) = 1000.
GN(2) = - 4
GK1 = .TRUE.
GK = 2.
CALL DE58R (FK, FG, FR, NX, X, GN1, GN, GK1, GK, Y, ALP,
* BET, IERR)
Результаты: (приводятся только для первых 10 узлов)
IERR = 0
Y YT - точное решение
- 4.11615097562 - 08 0.00000000000 + 00
1.99993826288 - 03 1.99999966667 - 03
5.99989114401 - 03 5.99999100000 - 03
1.19997755539 - 02 1.19999280001 - 02
1.99994536197 - 02 1.99996666683 - 02
2.99986015503 - 02 2.99988750126 - 02
4.19965896832 - 02 4.19969130680 - 02
5.59923325737 - 02 5.59926829533 - 02
7.19841091670 - 02 7.19844490075 - 02
3). d ((2 + x) dy/dx) / dx - xy = - 4sin x - 4x sin x + 2cos x , 0 < x < π / 2 ,
y = 0 при x = 0 ,
y = 2 при x = π / 2 ;
Подпрограммы - функции для вычисления значений коэффициентов k (x) = 2 + x, g (x) = x и r (x) = - 4sin x - 4x sin x + 2cos x ; те же, что и в примере N 2.
DIMENSION ALP(314), BET(314), Y(315)
DOUBLE PRECISION ALP(314), BET(314)
LOGICAL GN1, GK1
EXTERNAL FK, FG, FR
NX = 315
XN = 0.
XK = 1.57079632679
GN1 = .TRUE.
GN = 0.
GK1 = .TRUE.
GK = 2.
CALL DE57R (FK, FG, FR, NX, XN, GN1, GN, GK1, GK, Y, ALP,
* BET, IERR)
Результаты: (приводятся только для первых 10 узлов)
IERR = 0
Y YT - точное решение
0.00000000000 + 00 0.00000000000 + 00
1.00050273414 - 02 1.00050304151 - 02
2.00098044150 - 02 2.00098104511 - 02
3.00140808468 - 02 3.00140897354 - 02
4.00176062755 - 02 4.00176179076 - 02
5.00201303585 - 02 5.00201446263 - 02
6.00214027785 - 02 6.00214195751 - 02
7.00211732491 - 02 7.00211924689 - 02
8.00191915218 - 02 8.00192130604 - 02