|
Текст подпрограммы и версий ( Фортран ) de50r.zip , de50d.zip |
Тексты тестовых примеров ( Фортран ) tde50r.zip , tde50d.zip |
|
Текст подпрограммы и версий ( Си ) de50r_c.zip , de50d_c.zip |
Тексты тестовых примеров ( Си ) tde50r_c.zip , tde50d_c.zip |
|
Текст подпрограммы и версий ( Паскаль ) de50r_p.zip , de50e_p.zip |
Тексты тестовых примеров ( Паскаль ) tde50r_p.zip , tde50e_p.zip |
Вычисление решения двухточечной краевой задачи для линейного дифференциального уравнения второго порядка методом конечных разностей.
Решается двухточечная краевая задача для линейного дифференциального уравнения второго порядка
(1) y '' + f(x) y ' + g(x) y = r(x) , xN≤ x ≤ xK
c заданными на концах отрезка интегрирования [ xN, xК ] граничными условиями вида
(2) aN y '(xN) + bN y(xN) = cN , a2N + b2N ≠ 0 (3) aK y '(xK) + bK y(xK) = cK , a2K + b2K ≠ 0
- методом конечных разностей.
Решение вычисляется на заданной на отрезке [ xN, xK ] равномерной сетке с заданной абсолютной погрешностью EPS.
L.Fox, Proc. Roy. Soc. A 190, 31 - 59, 1947.
SUBROUTINE DE50R (NX, XN, XK, CN, CK, F, G, R, EPS,
IMAX, Y, IK, DELTY, RAB, IERR)
Параметры
| NX - | число узлов равномерной сетки на отрезке [ xN, xK ], (включая xN и xK), в которых требуется вычислить решение задачи. Предполагается, что узлы расположены в порядке возрастания: xN = x1 < x2< ... < xNX = xK; NX ≥ 3 (тип: целый); |
| XN, XK - | концы отрезка интегрирования, в которых заданы граничные условия (2) и (3), соответственно; XN < XK (тип: вещественный); |
| CN, CK - |
одномерные вещественные массивы длины 3, в
которых помещаются коэффициенты граничных условий
(2) и (3), соответственно, т.е.
CN (1) = aN,
CN (2) = bN,
CN (3) = cN,
CK (1) = aK,
CK (2) = bK,
CK (3) = cK. Коэффициенты граничных условий
должны удовлетворять условию a2N + b2N ≠ 0 и a2K + b2K ≠ 0; |
| F, G, R - |
одномерные вещественные масивы длины NX, в
которых перед обращением к подпрограмме должны
быть помещены значения коэффициентов уравнения
(1), вычисленные в узлах равномерной сетки, на
которой требуется получить решение краевой задачи, а именно:
F (1) = f ( x1), F (2) = f ( x2 ) ,..., F (NX) = f ( xNX );
G (1) = g ( x1), G (2) = g ( x2 ) ,..., G (NX) = g ( xNX );
R (1) = r ( x1), R (2) = r ( x2 ) ,..., R (NX) = r ( xNX );
|
| EPS - | абсолютная погрешность, с которой требуется вычислить решение краевой задачи (тип: вещественный); |
| IMAX - | максимальное допустимое число разностных итераций, котоpое разрешается использовать в процессе работы подпрограммы. Если эти итерации не нужны, то параметр IMAX при обращении к подпрограмме задается равным 0; в этом случае значение параметра EPS при работе подпрограммы не используется (тип: целый); |
| Y - | одномерный вещественный массив длины NX, в котоpом в результате работы подпрограммы получаются приближеннные значения решения краевой задачи, вычисленные на сетке (эти значения помещаются в Y в любом случае, даже и тогда, когда они не могут быть вычислены с требуемой точностью); |
| IK - | число разностных итераций, использованных в действительности подпрограммой для достижения требуемой точности EPS решения. B случае, когда данная точность не может быть достигнута за IMAX итераций, значение параметра IK в результате работы подпрограммы полагается равным - 1. Если IMAX = 0, то и значение IK в результате выполнения программы тоже полагается равным 0 (тип: целый); |
| DELTY - | одномерный вещественный массив длины NX, значение которого в результате работы подпрограммы полагается равным разности между начальным приближением решения и значением его последней итерации, вычисленной во всех NX узлах сетки. B случае IMAX = 0 все элементы массива DELTY полагаются равными нулю; |
| RAB - | одномерный вещественный рабочий массив длины 8 * NX; |
| IERR - | целая переменная, служащая для сообщения об ошибках, обнаруженных в процессе работы подпрограммы; при этом: |
| IERR=65 - | если в граничном условии (2) коэффициенты aN = bN = 0; |
| IERR=66 - | если в граничном условии (3) коэффициенты aK = bK = 0; |
| IERR = 1 - | если приближенное решение задачи, вычисленное за IMAX итераций, не достигает требуемой точности EPS; |
| В каждом из этих случаев решение задачи прекращается. При IERR = 1 интегрирование можно повторить обращением к подпрограмме с новыми значениями параметров NX, F, G, R или IMAX. |
Версии
| DE50D - | вычисление решения двухточечной краевой задачи для линейного дифференциального уравнения второго порядка методом конечных разностей с повышенной точностью. При этом параметры XN, XK, CN, CK, F, G, R, EPS, Y, DELTY, RAB должны иметь тип DOUBLE PRECISION. |
Вызываемые подпрограммы
| AS09R - | вычисление решения систем линейных алгебраических уравнений с ленточной матрицей. |
| AS09D - | вычисление решения систем линейных алгебраических уравнений с ленточной матрицей с двойной точностью. |
| I I10R - | подпрограмма вычисления центральных разностей таблично заданной функции. |
| I I10D - | подпрограмма вычисления центральных разностей таблично заданной функции с двойной точностью. |
| Подпрограммы AS09R и I I10R вызываются подпрограммой DE50R, а подпрограммы AS09D и I I10D - подпрограммой DE50D. |
| UTDE10 - | подпрограмма выдачи диагностических сообщений при работе подпрограммы DE50R. |
| UTDE11 - | подпрограмма выдачи диагностических сообщений при работе подпрограммы DE50D. |
Замечания по использованию
|
При задании значений коэффициентов уравнения в узлах равномерной сетки последняя должна быть достаточно густой для того, чтобы обеспечить сходимость метода конечных разностей. В общем случае заданная точность EPS приближенного pешения не гарантируется. Если для двух последовательных итераций Yi ( j ) и Yi ( j - 1 ) разность ei ( j ) = Yi ( j ) - Yi ( j - 1 ) удовлетворяет критерию | ei( j ) | ≤ EPS для i = 1, ..., NX , то yi( j ) принимается в качестве приближенного решения краевой задачи. Если же после выполнения IMAX итераций точность EPS не достигается, то, тем не менее, значение приближенного решения, полученное на последней итерации, помещается в массив Y. B этом случае параметр IERR = 1, а IK = - 1. Значения параметров NX, XN, XK, CN, CK, F, G, R, EPS, IMAX сохраняются. |
Решается краевая задача:
y '' + x2 y ' + x y = 10e- x*x ( 2x2 - 1 ) ( 2 - x ) , - 1≤ x ≤ 1
y ' - 2y = 0 при x = - 1
y ' + 2y = 0 при x = 1
(точное решение задачи y (x) = 10e- x*x )
DIMENSION CN(3), CK(3), F(200), G(200), R(200), YX(200),
* DELTY(200), RAB(1600)
XN = - 1.
XK = 1.
NX = 200
CN(1) = 1.
CN(2) = - 2.
CN(3) = 0.
CK(1) = 1.
CK(2) = 2.
CK(3) = 0.
EPS = 0.00001
IMAX = 4
N = NX - 1
H = (XK - XN) / N
X = XN - H
DO 1 I = 1, NX
X = X + H
F(I) = X**2
G(I) = X
R(I) = 10. * EXP(- X**2) * (2. * X**2 - 1.) * (2. - X)
1 CONTINUE
CALL DE50R (NX, XN, XK, CN, CK, F, G, R, EPS, IMAX, YX,
* IK, DELTY, RAB, IERR)
Результаты: (приводятся только для первых 10 узлов)
YX(1) = 3.678799304
YX(2) = 3.753113958
YX(3) = 3.828156408
YX(4) = 3.903910587
YX(5) = 3.980359671
YX(6) = 4.057486083
YX(7) = 4.135271487
YX(8) = 4.213696790
YX(9) = 4.292742144
YX(10) = 4.372386948
IK = 2
IERR = 0