Текст подпрограммы и версий
de51r_p.zip , de51e_p.zip 		Тексты тестовых примеров
tde51r_p.zip , tde51e_p.zip

Подпрограмма:  DE51R (модуль DE51R_p)

Назначение

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

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

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

(1)       y '' + f(x, y, y ') y ' + g(x, y, y ') y = r(x, y, y ')  ,     xN≤ x ≤ xK

с заданными на концах отрезка интегрирования  [ xN, xK ] граничными условиями вида 

(2)       aN y '(xN) + bN y(xN) = CN   ,

(3)       aK y '(xK) + bK y(xK) = CK   .

Решение вычисляется на заданной на отрезке  [ xN, xK ] равномерной сетке с заданной абсолютной погрешностью  EPS. Используется метод линеаризации исходного уравнения (1) с последующим применением метода конечных разностей для решения линейного уравнения.

L.Fox, Proc. Roy. Soc. A 190, 31 - 59, 1947.

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

procedure DE51R(F :Proc_F4_DE; NX :Integer; XN :Real; XK :Real;
                var CN :Array of Real; var CK :Array of Real;
                EPS :Real; var IMAX :Integer; IDER :Integer;
                var Y :Array of Real; var IK :Integer;
                var DELTY :Array of Real; var RAB :Array of Real;
                var IERR :Integer);


   Параметры
F - имя подпрограммы вычисления значений коэффициентов  f (x, y, y'),  g (x, y, y'),  r (x, y, y') уравнения (1). Первый оператор подпрограммы должен иметь вид:
  
 procedure F (var CF :Real; var CG :Real; var CR :Real; X :Real;
              Y :Real; I :Integer; var RAB :Array of Real);
Здесь:
X, Y - представляют первые два аргумента функций  f (x, y, y'),  g (x, y, y') и   r (x, y, y');
RAB - одномерный массив длины 11 * NX (NX - число узлов сетки);
I - указывает компоненту  RAB(I) массива   RAB, которая представляет третий аpгумент функций  f (x, y, y'),  g (x, y, y'),   r (x, y, y');
CF, CG, CR - представляют вычисленные подпрограммой  F значения функций  f (x, y, y'),  g (x, y, y') и  r (x, y, y');
  параметры  CF, CG, CR, X, Y и  RAB имеют тип Real, а  I - INTEGER;
NX - число узлов равномерной сетки на отрезке   [xN, xK], (включая  xN и  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;
EPS - абсолютная погрешность, с которой требуется вычислить решение краевой задачи;  EPS > 0 (тип: вещественный);
IMAX - указывает максимальное допустимое число линеаризаций исходного уравнения (1) и число итераций, используемых при решении получающегося после линеаризации уравнения: число линеаризаций pавно абсолютному значению  | IMAX | параметра  IMAX, а число итераций pавно  IMAX, если  IMAX > 0 и нулю, если  IMAX < 0;  IMAX ≠ 0. Если  IMAX < 0, то в результате работы подпрограммы значение  IMAX заменяется на противоположное; в этом случае  IMAX зажается переменной или элементом массива (тип: целый);
IDER - указывает на зависимость коэффициентов  f (x, y, y'),  g (x, y, y'),  r (x, y, y') уравнения (1) от производной  y': если хотя бы один из этих коэффициентов зависит от  y', то  IDER > 0, иначе  IDER = 0 (тип: целый);
Y - одномерный вещественный массив длины  NX, который должен содержать при обращении к подпрограмме значения начального приближения решения краевой задачи (1), (2), (3) во всех узлах сетки. В результате работы подпрограммы в этом массиве получаются значения конечного приближения к решению (эти значения помещаются в  Y в любом случае, даже и тогда, когда они не могут быть вычислены с требуемой точностью);
IK - целая переменная, содержащая число линеаризаций уравнения (1), выполненных в действительности подпрограммой для достижения требуемой точности  EPS решения. В случае, когда данная точность решения нелинейного уравнения не достигается за  IMAX линеаризаций или когда решение соответствующего линейного уравнения не может быть вычислено с точностью  EPS за  IMAX итераций, то значение параметра  IK в результате работы подпрограммы полагается равным  - 1;
DELTY - одномерный вещественный массив длины  NX, значение которого в результате работы подпрограммы полагается равным разности между начальным и конечным приближениями решения того линейного уравнения, котоpое получается после выполнения самой последней линеаризации исходного уpавнения (1), причем эта разность вычисляется во всех узлах сетки. В случае, когда  IMAX < 0, все элементы массива  DELTY в результате работы подпрограммы полагаются равными 0;
RAB - одномерный вещественный рабочий массив длины 11 * NX;
IERR - целая переменная, служащая для сообщения об ошибках, обнаруженных в процессе работы подпрограммы; при этом:
IERR=65 - если в граничном условии (2) коэффициенты  aN = bN = 0;
IERR=66 - если в граничном условии (3) коэффициенты  aK = bK = 0;
IERR=67 - если  IMAX = 0;
IERR=68 - если  IDER < 0;
IERR = 1 - если заданная точность решения нелинейного уравнения не достигается за  | IMAX | линеаризаций или если решение соответствующего линейного уpавнения не может быть вычислено с точностью  EPS за  IMAX итераций;
IERR = 3 - если заданная точность решения нелинейного уравнения не достигается за  | IMAX | линеаризаций или если решение соответствующего линейного уpавнения не может быть вычислено с точностью  EPS за  IMAX итераций, при этом при аппроксимации исходного уравнения линейным уравнением использовались производные, вычисленные с меньшей точностью, чем это необходимо;
  В каждом из выше указаных случаев решение задачи прекращается;
IERR = 2 - если при аппроксимации исходного уpавнения линейным уравнением использовались производные, вычисленные с меньшей точностью, чем это необходимо. При  IERR = 1, 2, 3 полученное приближенное решение запоминается в массиве  Y; при желании интегрирование можно повторить обращением к подпрограмме с новыми значениями параметров  NX, Y или  IMAX. 
   Версии
DE51E - вычисление решения двухточечной краевой задачи для нелинейного дифференциального уравнения второго порядка с расширенной (Extended) точностью. При этом параметры  XN, XK, CN, CK, EPS, Y, DELTY, RAB должны иметь тип Extended. 
   Вызываемые подпрограммы
DE50R - вычисление решения двухточечной краевой задачи для линейного дифференциального уравнения второго порядка методом конечных разностей.
DE50E - вычисление решения двухточечной краевой задачи для линейного дифференциального уравнения второго порядка методом конечных разностей с расширенной (Extended) точностью.
ID10R - подпрограмма вычисления первых и вторых производных с помощью центральных разностей.
ID10E - подпрограмма вычисления первых и вторых производных с помощью центральных разностей с расширенной (Extended) точностью.
  Подпрограммы DE50R и ID10R вызываются подпрограммой DE51R, а подпрограммы DE50E и ID10E - подпрограммой DE51E.
UTDE10 - подпрограмма выдачи диагностических сообщений при работе подпрограммы DE51R.
UTDE11 - подпрограмма выдачи диагностических сообщений при работе подпрограммы DE51E. 
   Замечания по использованию
 

Число узлов сетки  NX, задаваемое пользователем при обращении к подпрограмме, должно быть достаточно большим для того, чтобы обеспечить сходимость конечно - разностных аппроксимаций, используемых для решения линейных уравнений, и последовательных линеаризаций исходного уравнения (1).

В общем случае заданная точность  EPS не гарантируется.

Пусть  ηi ,  i = 1, ..., NX, представляет разность между текущим приближением к решению и приближением, полученным в результате предыдущей линеаризации уравнения (1). Тогда текущее приближение принимается за окончательное решение, если для всех  i (1 ≤ i ≤ NX) выполняется условие  | ηi | ≤ EPS .

Значения параметров  NX, XN, XK, CN, CK, EPS, IMAX, в результате работы подпрограммы сохраняются.

Kpоме того, DE50R и DE50E используют рабочие подпрограммы ID10RP и ID10EP, соответственно.

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

Решается краевая задача:

      y '' + ( xy' ) y' + e -y y = ( xy' + 2e-y - 2 ) cos x - ( 1 + 2xy' - e-y ) sin x 
      2y - y' = 0         при x = 0
      y + y' = 1          при x = π / 2

(точное решение задачи  y = sin x + 2cos x).


Unit TDE51R_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc, UtRes_p, FDE51R_p, DE51R_p;

function TDE51R: String;

implementation

function TDE51R: String;
var
NX,IMAX,IDER,I,N,IK,IERR :Integer;
XN,XK,EPS,H,X :Real;
CN :Array [0..2] of Real;
СК :Array [0..2] of Real;
YX :Array [0..149] of Real;
DELTY :Array [0..149] of Real;
RАВ :Array [0..1649] of Real;
Y :Array [0..149] of Real;
YR :Array [0..199] of Real;
label
_1,_2;
begin
Result := '';  { результат функции }
NX := 150;
XN := 0.0;
ХК := 3.14159265359;
ХК := XK/2.0;
CN[0] := 2.0;
CN[1] := -1.0;
CN[2] := 0.0;
CK[0] := 1.0;
CK[1] := 1.0;
CK[2] := -1.0;
IМАХ := 4;
EPS := 0.00001;
IDER := 1;
for I:=1 to NX do
 begin
  YX[I-1] := 1.5;
_1:
 end;
DE51R(FDE51R,NX,XN,XK,CN,CK,EPS,IMAX,IDER,YX,IK,DELTY,RAB,IERR);
Result := Result + Format('%s',['     IK']);
Result := Result + Format('%8d ',[IK]) + #$0D#$0A;
Result := Result + Format('%s',[' IERR=']);
Result := Result + Format('%3d ',[IERR]) + #$0D#$0A;
N := NX-1;
H := (XK-XN)/N;
X := XN-H;
for I:=1 to NX do
 begin
  X := X+H;
  Y[I-1] := Sin(X)+2.0*Cos(X);
  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('TDE51R',Result);  { вывод результатов в файл TDE51R.res }
exit;
end;

end.

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

procedure fde51r(var FC :Real; var GC :Real; var RC :Real; X :Real;
                Y :Real; I :Integer; var RAB :Array of Real);

implementation

procedure fde51r(var FC :Real; var GC :Real; var RC :Real; X :Real;
                Y :Real; I :Integer; var RAB :Array of Real);
begin
FC := X*RAB[I-1];
GC := Exp(-Y);
RC := (X*RAB[I-1]-2.0+2.0*Exp(-Y))*Cos(X)-(1.0+2.0*X*RAB[I-1]-
    Exp(-Y))*Sin(X);
end;

end.

Результаты:      (приводятся только для первых 10 узлов)

      YX(1)   =  1.999999086
      YX(2)   =  2.010430009
      YX(3)   =  2.020637497
      YX(4)   =  2.030620414
      YX(5)   =  2.040377652
      YX(6)   =  2.049908127
      YX(7)   =  2.059210778
      YX(8)   =  2.068284572
      YX(9)   =  2.077128500
      YX(10) =  2.085741581

      IK  =  3
      IERR  =  2