Текст подпрограммы и версий
de54r_p.zip , de54e_p.zip , de56r_p.zip , de56e_p.zip 		Тексты тестовых примеров
tde54r_p.zip , tde54e_p.zip , tde56r_p.zip , tde56e_p.zip

Подпрограмма:  DE54R (модуль DE54R_p) (версия DE56R)

Назначение

Вычисление решения двухточечной краевой задачи для линейного дифференциального уравнения второго порядка на равномерной сетке методом прогонки A.A.Абрамова с контролем точности.

Математическое описание

Решается двухточечная краевая задача для линейного дифференциального уравнения второго порядка 

(1)       y '' + f(x) y ' + g(x) y = r(x)

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

- методом прогонки А.А.Абрамова. Решение вычисляется на заданной на отрезке интегрирования равномерной сетке с заданной мерой погрешности  EPS. Контроль точности по меpе погрешности заключается в следующем: если численное решение по абсолютной величине не меньше некоторой заданной константы  P, то контроль точности ведется по относительной погрешности, иначе - по абсолютной.

А.А.Абрамов. Вариант метода прогонки. Ж. вычисл. матем. и матем. физ., 1961, 1, N 2, 349 - 351.

Использование 

procedure DE54R(FF :Func_F_DE; FG :Func_F_DE; FR :Func_F_DE; XN :Real;
                XK :Real; var CN :Array of Real;
                var CK :Array of Real; NX :Integer; X0 :Real;
                HX :Real; BULDY :Boolean; HMIN :Real;
                EPS :Real; P :Real; var H :Real;
                var Y :Array of Real; var DY :Array of Real;
                var NX1 :Integer; var IERR :Integer);


   Параметры
         FF, FG -
         FR   		имена вещественных подпрограмм - функций вычисления значений коэффициентов уравнения (1)  f (x), g (x) и  r (x), соответственно. Эти подпрограммы - функции должны иметь по одному паpаметpу, представляющему вещественное значение независимой переменной  x;
XN, XK - концы отрезка интегрирования, в которых заданы граничные условия (2) и (3), соответственно;  XK > XN или  XK < XN (тип: вещественный);
CN, CK - одномерные вещественные массивы длины 3, в которых помещаются коэффициенты граничных условий (2) и (3), т.е.
CN (1) = aN , CN (2) = bN , CN (3) = cN ,
CK (1) = aK , CK (2) = bK , CK (3) = cK ;
коэффициенты граничных условий должны удовлетворять неpавенствам:
a2N + b2N ≠ 0  и  a2K + b2K ≠ 0 ;
NX - число узлов равномерной сетки на отрезке интегрирования, в которых требуется вычислить решение задачи;  NX ≥ 1 (тип: целый);
XO - начальный узел равномерной сетки;  XO должен находиться между  XN и  XK или совпадать с каким-нибудь концом (тип: вещественный);
HX - шаг равномерной сетки;  HX > 0, если  XK > XN и  HX < 0, если  XK < XN (тип: вещественный);
BULDY - логическая переменная или константа, принимающая при входе в подпрограмму значение TRUE, если вместе с решением задачи требуется вычислить и его производную, и FALSE - в противном случае, т.е. когда необходимо иметь только решение задачи;
HMIN - минимальное значение абсолютной величины шага интегрирования, котоpое разрешается использовать при численном решении воспомагательных задач Коши, к которым сводится данная краевая задача (тип: вещественный);
EPS - допустимая меpа погрешности, с которой тpебуется вычислить решение краевой задачи (тип: вещественный);
P - граница перехода, используемая при оценке меры погрешности решения (тип: вещественный);
H - вещественная переменная, содержащая начальное значение шага интегрирования, используемое при численном решении воспомагательных задач Коши;
NX1 - целая переменная, указывающая на выходе из подпрограммы число узлов сетки, в которых вычислено решение данной краевой задачи;  NX1 ≤ NX;
Y, DY - одномерные вещественные массивы длины не более  NX1, в которых запоминаются значения pешения и его производной, вычисленные в  NX1 узлах сетки, при этом элементы  Y (I) и  DY (I) содержат значения решения и производной в узле  XO + (I - 1) * H, причем производная вычисляется только тогда, когда  BULDY=TRUE (таким образом, предполагается, что узлы сетки занумерованы в направлении от XN до  XK);
IERR - целая переменная, служащая для сообщения об ошибках, обнаруженных в процессе работы подпрограммы; при этом:
IERR=1 - когда последний узел сетки  XO + (N - 1) * H выходит за конец  XK интервала интегрирования, а предпоследний узел  XO + (N - 2) * H является внутренним узлом, т.е. лежит внутри интервала интегрирования;
IERR=3 - когда несколько узлов сетки (больше одного) выходят за конец  XK интервала интегрирования, а самый последний узел, принадлежащий отрезку интегрирования, является внутренним узлом;
IERR=2 - когда один из узлов сетки совпадает с концом  XK интервала интегрирования, а все последующие узлы не попадают в интервал интегрирования;
  во всех указанных выше случаях решение и его производная вычисляются только в тех узлах сетки, которые принадлежат отрезку интегрирования, а в случае  IERR = 1 и  IERR = 3 они вычисляются также в конце  XK отрезка и запоминаются в  Y (NX1) и  DY (NX1), соответственно;
IERR=65 - когда в граничном условии (2) на конце  XN коэффициенты  aN = bN = 0;
IERR=66 - когда в граничном условии (3) на конце  XK коэффициенты  aK = bK = 0;
IERR=67 - когда начальный узел  XO сетки не принадлежит отрезку интегрирования;
IERR=68 и IERR=70 - когда решение краевой задачи не может быть вычислено с требуемой точностью  EPS при заданном минимальном значении  HMIN шага интегрирования, при этом  IERR = 68 указывает, что точность не достигается при прямой прогонке Абрамова, а IERR=70 - при обратной прогонке;
  при  IERR = 65, 66, 67, 68, 70 интегрирование уравнения прекращается; в случае  IERR = 68 и  IERR = 70 интегрирование можно повторить обращением к подпрограмме с новыми значениями  H и  HMIN;
IERR=69 - когда данная краевая задача не может быть решена методом прогонки А.А.Абрамова. 
   Версии
DE54E - вычисление решения двухточечной краевой задачи для линейного дифференциального уравнения второго порядка на равномерной сетке методом прогонки А.А.Абрамова с расширенной (Extended) точностью. При этом параметры  FF, FG, FR, XN, XK, CN, CK, XO, HX, HMIN, EPS, P, H, Y, DY и формальный параметp в подпрограммах - функциях FF, FG и  FR должны иметь тип Extended.
DE56R - вычисление решения двухточечной краевой задачи для линейного дифференциального уравнения второго порядка на неравномерной сетке методом прогонки А.А.Абрамова с контролем точности. Первый оператор подпрограммы имеет вид:
  
procedure DE56R(FF :Func_F_DE; FG :Func_F_DE; FR :Func_F_DE; XN :Real;
                XK :Real; var CN :Array of Real;
                var CK :Array of Real; NX :Integer;
                var X :Array of Real; BULDY :Boolean; HMIN :Real;
                EPS :Real; P :Real; var H :Real;
                var Y :Array of Real; var DY :Array of Real;
                var NX1 :Integer; var IERR :Integer);

Здесь:  X - одномерный вещественный массив длины  NX, представляющий узлы неравномерной сетки, в которых требуется вычислить решение задачи, при этом узлы сетки располагаются в  X в направлении от  XN до  XK.

Остальные параметры подпрограммы DE56R имеют тот же смысл, что и одноименные параметры подпрограммы DE54R. При этом выход из подпрограммы DE56R со значением  IERR = 67 указывает на то, что не принадлежит интервалу интегрирования.
DE56E - вычисление решения двухточечной краевой задачи для линейного дифференциального уравнения второго порядка на неравномерной сетке методом прогонки А.А.Абрамова с расширенной (Extended) точностью. Первый оператор подпрограммы имеет тот же вид, что и в подпрограмме DE56R; при этом параметры  FF, FG, FR, XN, XK, CN, CK, X, HMIN, EPS, P, H, Y, DY и формальный параметр в подпрограммах - функциях  FF, FG, FR должны иметь тип Extended. 
   Вызываемые подпрограммы
UTDE12 - подпрограмма выдачи диагностических сообщений при работе подпрограмм DE54R и DE56R.
UTDE13 - подпрограмма выдачи диагностических сообщений при работе подпрограмм DE54E и DE56E.
  Kpоме того, подпрограммы DE54R, DE56R и подпрограммы DE54E, DE56E используют рабочие подпрограммы DE54RU и DE54EU, соответственно. 
   Замечания по использованию
 

B общем случае заданная точность приближенного решения не гарантируется.

Значения параметров  XN, XK, CN, CK, NX, XO, HX, BULDY, HMIN, EPS, P, X сохраняются.

Пример использования 

    y '' + x2 y ' + xy = 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 TDE54R_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc, UtRes_p, FDE54R_p, FGDE54R_p, FRDE54R_p, DE54R_p;

function TDE54R: String;

implementation

function TDE54R: String;
var
NX,I,NX1,IERR :Integer;
XN,XK,X0,HX,P,HMIN,EPS,H :Real;
BULDY :Boolean;
Y :Array [0..10] of Real;
DY :Array [0..10] of Real;
const
CN :Array [0..2] of Real = ( +1.0,-2.0,0.0 );
СК :Array [0..2] of Real = ( 1.0,2.0,0.0 );
begin
Result := '';  { результат функции }
XN := -1.0;
ХК := 1.0;
NX := 11;
X0 := XN;
НХ := 0.01;
P := 11.0;
НХ := 2.0/199.0;
X0 := 1.0-10.0*HX;
BULDY := True;
HMIN := 1.E-10;
EPS := 0.00001;
H := 0.01;
DE54R(FDE54R,FGDE54R,FRDE54R,XN,XK,CN,CK,NX,X0,HX,BULDY,HMIN,EPS,P,
     H,Y,DY,NX1,IERR);
Result := Result + Format('%s',[' IERR=']);
Result := Result + Format('%3d ',[IERR]) + #$0D#$0A;
Result := Result + Format('%s',[' NX1=']);
Result := Result + Format('%5d ',[NX1]);
Result := Result + Format('%s',['       Y']);
Result := Result + Format('%s',['                 DY' + #$0D#$0A]);
Result := Result + #$0D#$0A;
for I:=1 to 11 do
 begin
  Result := Result + Format('%20.16f %20.16f ',
 [Y[I-1],DY[I-1]]) + #$0D#$0A;
 end;
Result := Result + #$0D#$0A;
UtRes('TDE54R',Result);  { вывод результатов в файл TDE54R.res }
exit;
end;

end.

Unit fde54r_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc;

function fde54r(X :Real): Real;

implementation

function fde54r(X :Real): Real;
begin
{ Result - прототип имени функции F на FORTRANe }
Result := X*X;
exit;
end;

end.

Unit fgde54r_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc;

function fgde54r(X :Real): Real;

implementation

function fgde54r(X :Real): Real;
begin
{ Result - прототип имени функции G на FORTRANe }
Result := X;
exit;
end;

end.

Unit frde54r_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc;

function frde54r(X :Real): Real;

implementation

function frde54r(X :Real): Real;
var
X2 :Real;
begin
{ Result - прототип имени функции R на FORTRANe }
X2 := X*X;
Result := 10.0*Exp(-X2)*(2.0*X2-1.0)*(2.0-X);
exit;
end;

end.

Результаты:

      IERR  =  0
      NX1   =  11

                    Y                                     DY
      4.45260617314 + 00       - 8.01021606015 + 00
      4.37238323319 + 00       - 7.95378250048 + 00
      4.29273839036 + 00       - 7.89518713214 + 00
      4.21369299696 + 00       - 7.83450451556 + 00
      4.13526765439 + 00       - 7.77180951660 + 00
      4.05748221040 + 00       - 7.70717722605 + 00
      3.98035575729 + 00       - 7.64068288024 + 00
      3.90390663078 + 00       - 7.57240178309 + 00
      3.82815240978 + 00       - 7.50240922954 + 00
      3.75310991678 + 00       - 7.43078043029 + 00
      3.67879521918 + 00       - 7.35759043838 + 00