|
Текст подпрограммы и версий 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.
procedure DE50R(NX :Integer; XN :Real; XK :Real;
var CN :Array of Real; var CK :Array of Real;
var F :Array of Real; var G :Array of Real;
var R :Array of Real; EPS :Real; var IMAX :Integer;
var Y :Array of Real; var IK :Integer;
var DELTY :Array of Real; var RAB :Array of Real;
var IERR :Integer);
Параметры
| 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. |
Версии
| DE50E - | вычисление решения двухточечной краевой задачи для линейного дифференциального уравнения второго порядка методом конечных разностей с расширенной (Extended) точностью. При этом параметры XN, XK, CN, CK, F, G, R, EPS, Y, DELTY, RAB должны иметь тип Extended. |
Вызываемые подпрограммы
| AS09R - | вычисление решения систем линейных алгебраических уравнений с ленточной матрицей. |
| AS09E - | вычисление решения систем линейных алгебраических уравнений с ленточной матрицей с расширенной (Extended) точностью. |
| I I10R - | подпрограмма вычисления центральных разностей таблично заданной функции. |
| I I10E - | подпрограмма вычисления центральных разностей таблично заданной функции с расширенной (Extended) точностью. |
| Подпрограммы AS09R и I I10R вызываются подпрограммой DE50R, а подпрограммы AS09E и I I10E - подпрограммой DE50E. |
| UTDE10 - | подпрограмма выдачи диагностических сообщений при работе подпрограммы DE50R. |
| UTDE11 - | подпрограмма выдачи диагностических сообщений при работе подпрограммы DE50E. |
Замечания по использованию
|
При задании значений коэффициентов уравнения в узлах равномерной сетки последняя должна быть достаточно густой для того, чтобы обеспечить сходимость метода конечных разностей. В общем случае заданная точность 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 )
Unit TDE50R_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc, UtRes_p, DE50R_p;
function TDE50R: String;
implementation
function TDE50R: String;
var
NX,IMAX,N,I,IK,IERR :Integer;
XN,XK,EPS,H,X :Real;
CN :Array [0..2] of Real;
СК :Array [0..2] of Real;
F :Array [0..199] of Real;
G :Array [0..199] of Real;
R :Array [0..199] of Real;
YX :Array [0..199] of Real;
DELTY :Array [0..199] of Real;
RАВ :Array [0..1599] of Real;
Y :Array [0..199] of Real;
YR :Array [0..199] of Real;
label
_1,_2;
begin
Result := ''; { результат функции }
XN := -1.0;
ХК := 1.0;
NX := 200;
CN[0] := 1.0;
CN[1] := -2.0;
CN[2] := 0.0;
CK[0] := 1.0;
CK[1] := 2.0;
CK[2] := 0.0;
IМАХ := 4;
EPS := 0.00001;
N := NX-1;
H := (XK-XN)/N;
X := XN-H;
for I:=1 to NX do
begin
X := X+H;
F[I-1] := IntPower(X,2);
G[I-1] := X;
R[I-1] := 10.0*Exp(-IntPower(X,2))*(2.0*IntPower(X,2)-1.0)*(2.0-X);
_1:
end;
DE50R(NX,XN,XK,CN,CK,F,G,R,EPS,IMAX,YX,IK,DELTY,RAB,IERR);
Result := Result + Format('%s',[' IK']);
Result := Result + Format('%8d ',[IK]) + #$0D#$0A;
X := XN-H;
for I:=1 to NX do
begin
X := X+H;
Y[I-1] := 10.0*Exp(-IntPower(X,2));
YR[I-1] := YX[I-1]-Y[I-1];
Result := Result + Format(' %20.16f %20.16f %20.16f %20.16f ',
[X,YX[I-1],Y[I-1],YR[I-1]]) + #$0D#$0A;
_2:
end;
UtRes('TDE50R',Result); { вывод результатов в файл TDE50R.res }
exit;
end;
end.
Результаты: (приводятся только для первых 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