Текст подпрограммы и версий
de11r_p.zip , de11e_p.zip 		Тексты тестовых примеров
tde11r_p.zip , tde11e_p.zip

Подпрограмма:  DE11R (модуль DE11R_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 ) ,
 методом Хойна.

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

И.С.Березин, Н.П.Жидков, Методы вычислений, т.2, Физматгиз, 1960.

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

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

Подпрограмма DE11R предназначена для численного решения дифференциальных уравнений и систем уравнений с правой частью, имеющей непрерывные частные производные вплоть до 4 порядка включительно.

Хотя заданная точность EPS не гарантируется в общем случае, большой опыт эксплуатации данной подпрограммы убедительно показывает, что вычисляемое ею численное решение достаточно близко приближает точное решение.

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

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

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

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

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

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

Unit TDE11R_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc, UtRes_p, FDE11R_p, DE11R_p;

function TDE11R: String;

implementation

function TDE11R: String;
var
M,_i,IERR :Integer;
XK,XN,H,EPS,P,HMIN :Real;
Y :Array [0..1] of Real;
YP :Array [0..1] of Real;
DELTY :Array [0..1] of Real;
RАВ :Array [0..3] of Real;
begin
Result := '';  { результат функции }
Y[0] := Sqrt(2.0)/2;
Y[1] := -Y[0];
ХК := 3.14159265359;
M := 2;
XN := 0.75*3.14159265359;
H := 0.01;
EPS := 3.E-5;
P := 1.E-6;
HMIN := 1.E-8;
DE11R(FDE11R,M,XN,Y,XK,HMIN,EPS,
     P,H,Y,YP,DELTY,RAB,IERR);
Result := Result + Format('%s',['       XK']);
Result := Result + Format('%s',['               YK(1)']);
Result := Result + Format('%s',['               YK(2)' + #$0D#$0A]);
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('TDE11R',Result);  { вывод результатов в файл TDE11R.res }
exit;
end;

end.

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

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

implementation

procedure fde11r(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.5928428474 * 10-7
          Y (2)  =  -0.9999983798

          IERR  =   0