Текст подпрограммы и версий
de10r_p.zip , de10e_p.zip 		Тексты тестовых примеров
tde10r_p.zip , tde10e_p.zip

Подпрограмма:  DE10R (модуль DE10R_p)

Назначение

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

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

Решается задача Коши для системы M обыкновенных дифференциальных уравнений 

          Y ' = F (X, Y) ,
          Y = ( y1, ... , yM ) ,
          F = ( f1 (X, y1, ... , yM), ... , fM (X, y1, ... , yM) )
 с начальными условиями, заданными в точке XN :
          Y(XN) = YN ,     YN = ( y10, ... , yM0 ) ,
 методом Mеpсона.

Решение вычисляется в одной точке XK, которая является концом интервала интегрирования. Каждая компонента решения вычисляется с контролем точности по относительной погрешности на тех участках интервала интегрирования, на которых модуль этой компоненты больше некоторого наперед заданного числа Р (это число называется границей перехода), и по абсолютной погрешности на остальных участках, т.е. там, где модуль проверяемой на точность компоненты меньше этого числа.

Дж.Н.Ланс, Численные методы для быстродействующих вычислительных машин, Изд-во иностранной литературы, M., 1962.

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

procedure DE10R(F :Proc_F_DE; M :Integer; XN :Real;
                var YN :Array of Real; XK :Real; HMIN :Real;
                EPS :Real; P :Real; var H :Real; var Y :Array of Real;
                var RA :Array of Real; var IERR :Integer);


   Параметры
F - имя подпрограммы вычисления значений правой части дифференциального уравнения. Первый оператоp подпрограммы должен иметь вид:
procedure F (X :Real; var Y :Array of Real; var DY :Array of Real; M :Integer);
Здесь: X, Y - значения независимой и зависимой переменных, соответственно. Вычисленное значение правой части должно быть помещено в DY. B случае системы уравнений, т.е. когда M ≠ 1, параметры Y и DY представляют массивы длины M (тип параметров X, Y и DY: вещественный);
M - количество уравнений в системе (тип: целый);
XN, YN - начальные значения аргумента и решения. В случае системы уравнений (т.е. M ≠ 1) YN представляет одномерный массив длины M (тип: вещественный);
XK - значение аргумента, при котоpом требуется вычислить решение задачи Коши (конец интервала интегрирования). XK может быть больше, меньше или pавно XN (тип: вещественный);
HMIN - минимальное значение абсолютной величины шага, который разрешается использовать при интегрировании данной системы уравнений (тип: вещественный);
EPS - допустимая меpа погрешности, с которой требуется вычислить все компоненты решения (тип: вещественный);
P - граница перехода, используемая при оценке погрешности решения (тип: вещественный);
H - вещественная переменная, содержащая начальное значение шага интегрирования. Может задаваться с учетом направления интегрирования, т.е. положительным, если XK > XN, отрицательным, если XK < XN, или без всякого учета в виде абсолютной величины;
Y - искомое решение задачи Коши, вычисленное подпрограммой при значении аргумента XK. Для системы уравнений (когда M ≠ 1) задается одномерным массивом длины M. В случае совпадения значений параметров XN и XK значение Y полагается равным начадьному значению YN (тип: вещественный);
RA - одномерный рабочий массив вещественного типа длины 4*M.
IERR - целая переменная, значение которой в результате работы подпрограммы полагается равным 65, если какая-нибудь компонента решения не может быть вычислена с требуемой точностью EPS. В этом случае интегрирование системы прекращается. При желании интегрирование системы можно повторить обращением к подпрограмме с новыми значениями параметров HMIN и H. 
   Версии
DE10E - вычисление решения задачи Коши для системы обыкновенных дифференциальных уравнений первого порядка в конце интервала интегрирования методом Mеpсона с расширенной (Extended) точностью. При этом параметры XN, YN, XK, HMIN, EPS, P, H, Y, RA и параметры X, Y и DY в подпрограмме F должны иметь тип Extended. 
   Вызываемые подпрограммы
UTDE10 - подпрограмма выдачи диагностических сообщений при работе подпрограммы DE10R.
UTDE11 - подпрограмма выдачи диагностических сообщений при работе подпрограммы DE10E. 
   Замечания по использованию
 

Подпрограммы DE10R и DE10E предназначены для численного решения дифференциальных уравнений и систем уравнений с правой частью, имеющей непрерывные частные производные вплоть до 5 порядка включительно. Они являются эффективными для нежестких уравнений и систем уравнений с несложными правыми частями (т.е. не являющимися трудоемкими для вычислений). Хотя заданная точность EPS не гарантируется в общем случае, большой опыт эксплуатации данной подпрограммы убедительно показывает, что вычисляемое ею численное решение достаточно близко приближает точное решение.

При работе подпрограммы значения параметров M, XN, YN, XK, HMIN, EPS и P сохраняются. Если после работы подпрограммы нет необходимости иметь начальное значение решения YN, то параметры YN и Y при обращении к ней можно совместить. При работе подпрограммы счета правой части F значения параметров X, Y и DY не должны изменяться.

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

Использование подпрограммы иллюстрируется на примере: 

          y '1  =  y2 ,
          y '2  =  -y1 ,      3π/4 ≤ x ≤ π

          y1(3π/4)  =  √2 / 2 ,      y2(3π/4)  =  - √2 / 2

Приводятся подпрограмма вычисления значений правой части и фрагмент вызывающей (основной) программы, а также результаты счета. 

Unit TDE10R_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc, UtRes_p, FDE10R_p, DE10R_p;

function TDE10R: String;

implementation

function TDE10R: String;
var
M,_i,IERR :Integer;
XN,H,XK,HMIN,EPS,P :Real;
RA :Array [0..7] of Real;
Y :Array [0..1] of Real;
begin
Result := '';  { результат функции }
M := 2;
XN := 0.75*3.14159265359;
H := 0.01;
Y[0] := Sqrt(2.0)/2.0;
Y[1] := -Y[0];
ХК := 3.14159265359;
HMIN := 1.E-4;
EPS := 1.E-5;
P := 1.E-7;
DE10R(FDE10R,M,XN,Y,XK,HMIN,EPS,P,H,Y,RA,IERR);
Result := Result + Format('%s',['      XF']);
Result := Result + Format('%s',['                YF(1)']);
Result := Result + Format('%s',['           YF(2)']);
Result := Result + Format('%20.16f ',[XK]) + #$0D#$0A;
Result := Result + #$0D#$0A;
for _i:=0 to 1 do
 begin
  Result := Result + Format('%20.16f ',[Y[_i]]);
  if ( ((_i+1) mod 2)=0 )
   then Result := Result + #$0D#$0A;
 end;
Result := Result + #$0D#$0A;
UtRes('TDE10R',Result);  { вывод результатов в файл TDE10R.res }
exit;
end;

end.

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

procedure fde10r(X :Real; var Y :Array of Real; var Z :Array of Real;
                M :Integer);

implementation

procedure fde10r(X :Real; var Y :Array of Real; var Z :Array of Real;
                M :Integer);
begin
Z[0] := Y[1];
Z[1] := -Y[0];
end;

end.

Результаты:
       
          Y (1)  =   0.4495909707 * 10-8
          Y (2)  =  -0.9999999996

          IERR  =  0