|
Текст подпрограммы и версий 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.
procedure DE59R(var CK :Array of Real; var CG :Array of Real;
var CR :Array of Real; NX :Integer; XN :Real;
XK :Real; GN1 :Boolean; var GN :Array of Real;
GK1 :Boolean; var GK :Array of Real;
var Y :Array of Real; var ALP :Array of Extended;
var BET :Array of Extended; var IERR :Integer);
Параметры
| 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 | одномерные вещественные рабочие массивы расширенной (Extended) точности длины 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 - | вычисление решения двухточечной краевой задачи для обыкновенного линейного самосопряженного дифференциального уравнения второго порядка с непрерывными коэффициентами на неравномерной сетке. |
| DE57E - | то же, что и DE57R, при этом все арифметические операции с вещественными числами выполняются с удвоенным числом значащих цифр. |
| DE58E - | то же, что и DE58R, при этом все арифметические операции с вещественными числами выполняются с удвоенным числом значащих цифр. |
| DE59E - | то же, что и DE59R, при этом все арифметические операции с вещественными числами выполняются с удвоенным числом значащих цифр. |
| DE60E - | то же, что и DE60R, при этом все арифметические операции с вещественными числами выполняются с удвоенным числом значащих цифр. |
|
Первые операторы этих подпрограмм имеют следующий вид:
procedure DE60R(var CK :Array of Real; var CG :Array of Real;
var CR :Array of Real; NX :Integer;
var X :Array of Real; GN1 :Boolean;
var GN :Array of Real; GK1 :Boolean;
var GK :Array of Real; var Y :Array of Real;
var ALP :Array of Extended; var BET :Array of Extended;
var IERR :Integer);
Параметр X в подпрограммах DE60R и DE58R - имя одномерного вещественного массива длины NX, в последовательных элементах которого, начиная с первого, размещаются в возрастающем порядке узлы неравномерной сетки, причем начало и конец отрезка интегрирования, в которых заданы граничные условия (4) - (7), хранятся, соответственно, в первом и последнем элементах X (1), X (NX) этого массива. Параметры FK, FG, FR в подпрограммах DE57R и DE58R - имена вещественных подпрограмм - функций function FK (X :Real): Real; function FG (X :Real): Real; function FR (X :Real): Real; вычисляющих значения коэффициентов уравнения 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 имеет версию DE57E, DE58E, DE59E и DE60E, соответственно, с теми же самыми эффектом и списком параметpов, выполняющую все операции над действительными числами с удвоенным числом значащих цифр. При этом параметры FK, FG, FR (в подпрограммах DE57E,
DE58E), являются именами подпрограмм - функций |
|
CK, CG - CR | (в подпрограммах DE59E, DE60E), |
| XN, XK - | (в подпрограммах DE57E, DE59E), |
| X - | (в подпрограммах DE58E, DE60E), |
|
GN, GK - Y, ALP BET | (во всех подпрограммах) |
|
должны иметь тип Extended. Тип остальных параметров не изменяется. В версиях подпрограмм DE57R, DE58R, DE59R и DE60R, формальные параметры ALP и BET имеют тип Extended. |
Вызываемые подпрограммы
| UTDE14 - | подпрограмма выдачи диагностических сообщений при работе подпрограмм DE57R, DE58R, DE59R и DE60R. |
| UTDE15 - | подпрограмма выдачи диагностических сообщений при работе подпрограмм DE57E, DE58E, DE59E и DE60E. |
Замечания по использованию
|
Значение параметра NX, задаваемое при обращении к подпрограмме, должно быть достаточно большим, чтобы получить хорошее приближение численного решения к точному решению задачи. При интегрировании дифференциального уравнения (1) с непрерывными коэффициентами с помощью подпрограмм DE57R, DE58R, DE57E, DE58E значения коэффициента 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 сохраняются. |
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 .
Unit TDE59R_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc, UtRes_p, DE59R_p;
function TDE59R: String;
implementation
function TDE59R: String;
var
NX,I,_i,IERR :Integer;
XN,XK,GN,GK :Real;
GN1,GK1 :Boolean;
ALP :Array [0..19] of Extended;
ВЕТ :Array [0..19] of Extended;
Y :Array [0..20] of Real;
СК :Array [0..19] of Real;
CG :Array [0..39] of Real;
CR :Array [0..39] of Real;
label
_200,_5,_15;
begin
Result := ''; { результат функции }
NX := 21;
XN := 0.0;
ХК := 1.0;
GN1 := True;
GN := 1.0;
GK1 := True;
GK := 0.0;
for I:=1 to 40 do
begin
CG[I-1] := 0.0;
CR[I-1] := 0.0;
_200:
end;
for I:=1 to 10 do
begin
_5:
CK[I-1] := 1.0;
end;
for I:=11 to 20 do
begin
_15:
CK[I-1] := 10.0;
end;
DE59R(CK,CG,CR,NX,XN,XK,GN1,GN,GK1,GK,Y,ALP,BET,IERR);
Result := Result + Format('%s',[' IERR=']);
Result := Result + Format('%2d ',[IERR]) + #$0D#$0A;
Result := Result + Format('%s',[' XN=']);
Result := Result + Format('%20.16f ',[XN]);
Result := Result + Format('%s',[' XK=']);
Result := Result + Format('%20.16f ',[XK]);
Result := Result + Format('%s',[' NX=']);
Result := Result + Format('%5d ',[NX]) + #$0D#$0A;
Result := Result + #$0D#$0A;
for _i:=0 to 20 do
begin
Result := Result + Format('%20.16f ',[Y[_i]]);
if ( ((_i+1) mod 4)=0 )
then Result := Result + #$0D#$0A;
end;
Result := Result + #$0D#$0A;
UtRes('TDE59R',Result); { вывод результатов в файл TDE59R.res }
exit;
end;
end.
Результаты:
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 .
Решение вычисляется на неравномерной сетке.
Unit tde58r_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc, UtRes_p, FFKDE58R_p, FFGDE58R_p, FFRDE58R_p, DE58R_p;
function tde58r: String;
implementation
function tde58r: String;
var
NX,J,K,I,NXM1,IERR :Integer;
XN,XK,H,GK :Real;
GN1,GK1 :Boolean;
ALP :Array [0..281] of Extended;
BET :Array [0..281] of Extended;
Y :Array [0..282] of Real;
YT :Array [0..282] of Real;
Z :Array [0..282] of Real;
X :Array [0..282] of Real;
GN :Array [0..1] of Real;
label
_1,_100,_20;
begin
Result := ''; { результат функции }
NX := 283;
XN := 0.0;
XK := 1.57079632679;
H := 0.001;
X[0] := XN;
X[281] := XK-0.02;
X[282] := XK;
J := 0;
for K:=1 to 28 do
begin
for I:=1 to 10 do
begin
J := J+1;
_1:
X[J+0] := X[J-1]+(I)*H;
end;
end;
GN1 := False;
GN[0] := 1000.0;
GN[1] := -4.0;
GK1 := True;
GK := 2.0;
DE58R(FFKDE58R,FFGDE58R,FFRDE58R,NX,X,GN1,GN,GK1,GK,Y,ALP,BET,IERR);
Result := Result + Format('%s',[' IERR=']);
Result := Result + Format('%2d ',[IERR]) + #$0D#$0A;
Result := Result + Format('%s',[' XN=']);
Result := Result + Format('%20.16f ',[XN]);
Result := Result + Format('%s',[' XK=']);
Result := Result + Format('%20.16f ',[XK]);
Result := Result + Format('%s',[' NX=']);
Result := Result + Format('%5d ',[NX]);
Result := Result + Format('%s',[' Y']);
Result := Result + Format('%s',[' YT']);
Result := Result + Format('%s',[' Z' + #$0D#$0A]);
NXM1 := NX-1;
for I:=1 to NX do
begin
YT[I-1] := 2.0*Sin(X[I-1]);
Z[I-1] := Y[I-1]-YT[I-1];
_100:
end;
for I:=1 to NXM1 do
begin
Result := Result + Format('%20.16f %20.16f %20.16f ',
[Y[I-1],YT[I-1],Z[I-1]]) + #$0D#$0A;
_20:
end;
Result := Result + Format('%20.16f %20.16f %20.16f ',
[Y[NX-1],YT[NX-1],Z[NX-1]]) + #$0D#$0A;
UtRes('tde58r',Result); { вывод результатов в файл tde58r.res }
exit;
end;
end.
Unit ffkde58r_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc;
function ffkde58r(X :Real): Real;
implementation
function ffkde58r(X :Real): Real;
begin
{ Result - прототип имени функции FK на FORTRANe }
Result := 2.0+X;
exit;
end;
end.
Unit ffgde58r_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc;
function ffgde58r(X :Real): Real;
implementation
function ffgde58r(X :Real): Real;
begin
{ Result - прототип имени функции FG на FORTRANe }
Result := X;
exit;
end;
end.
Unit ffrde58r_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc;
function ffrde58r(X :Real): Real;
implementation
function ffrde58r(X :Real): Real;
begin
{ Result - прототип имени функции FR на FORTRANe }
Result := 4.0*Sin(X)*(1.0+X)-2.0*Cos(X);
exit;
end;
end.
Результаты: (приводятся только для первых 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.
Unit tde57r_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc, UtRes_p, FFKDE57R_p, FFGDE57R_p, FFRDE57R_p, DE57R_p;
function tde57r: String;
implementation
function tde57r: String;
var
NX,NXM1,I,IERR :Integer;
XN,XK,GN,GK,H,X :Real;
GN1,GK1 :Boolean;
Y :Array [0..314] of Real;
ALP :Array [0..313] of Extended;
BET :Array [0..313] of Extended;
YT :Array [0..314] of Real;
label
_100,_20;
begin
Result := ''; { результат функции }
NX := 315;
XN := 0.0;
XK := 1.57079632679;
GN1 := True;
GN := 0.0;
GK1 := True;
GK := 2.0;
DE57R(FFKDE57R,FFGDE57R,FFRDE57R,NX,XN,XK,GN1,GN,GK1,GK,Y,ALP,BET,IERR);
Result := Result + Format('%s',[' IERR=']);
Result := Result + Format('%2d ',[IERR]) + #$0D#$0A;
NXM1 := NX-1;
H := (XK-XN)/NXM1;
X := XN-H;
for I:=1 to NX do
begin
X := X+H;
YT[I-1] := 2.0*Sin(X);
_100:
end;
Result := Result + Format('%s',[' XN=']);
Result := Result + Format('%20.16f ',[XN]);
Result := Result + Format('%s',[' XK=']);
Result := Result + Format('%20.16f ',[XK]);
Result := Result + Format('%s',[' NX=']);
Result := Result + Format('%5d ',[NX]) + #$0D#$0A;
for I:=1 to NXM1 do
begin
Result := Result + Format(' %20.16f %20.16f ',
[Y[I-1],YT[I-1]]) + #$0D#$0A;
_20:
end;
Result := Result + Format(' %20.16f %20.16f ',
[Y[NX-1],YT[NX-1]]) + #$0D#$0A;
UtRes('tde57r',Result); { вывод результатов в файл tde57r.res }
exit;
end;
end.
Unit ffkde57r_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc;
function ffkde57r(X :Real): Real;
implementation
function ffkde57r(X :Real): Real;
begin
{ Result - прототип имени функции FK на FORTRANe }
Result := 2.0+X;
exit;
end;
end.
Unit ffgde57r_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc;
function ffgde57r(X :Real): Real;
implementation
function ffgde57r(X :Real): Real;
begin
{ Result - прототип имени функции FG на FORTRANe }
Result := X;
exit;
end;
end.
Unit ffrde57r_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc;
function ffrde57r(X :Real): Real;
implementation
function ffrde57r(X :Real): Real;
begin
{ Result - прототип имени функции FR на FORTRANe }
Result := 4.0*Sin(X)*(1.0+X)-2.0*Cos(X);
exit;
end;
end.
Результаты: (приводятся только для первых 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