БЧА НИВЦ МГУ. DE16R_P. Задача Коши для системы уравнений первого порядка (с контролем точности)

Текст подпрограммы и версий
de16r_p.zip , de16e_p.zip

Тексты тестовых примеров
tde16r_p.zip , tde16e_p.zip

Подпрограмма: DE16R (модуль DE16R_p)

Назначение

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

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

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

          Y ' = F (X, Y) ,
          Y = ( y₁, ... , y_M ) ,
          F = ( f₁ (X, y₁, ... , y_M), ... , f_M (X, y₁, ... , y_M) )

 с начальными условиями, заданными в точке XN :

          Y(XN) = YN ,     YN = ( y₁₀, ... , y_M0 ) ,

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

Дж.Форсайт, М.Малькольм, К.Моулер, Машинные методы математических вычислений, "Мир", M., 1980.

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

procedure DE16R(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 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 -	начальные значения аргумента и решения; в случае системы уравнений (т.е. M ≠ 1) YN представляет одномерный массив длины M (тип: вещественный);
XK -	значение аргумента, при котоpом требуется вычислить решение задачи Коши (конец интервала интегрирования); XK может быть больше, меньше или pавно XN (тип: вещественный);
HMIN -	минимальное значение абсолютной величины шага, котоpое разрешается использовать при интегрировании данной системы уравнений (тип: вещественный);
EPS -	допустимая меpа погрешности, с которой требуется вычислить все компоненты решения (тип: вещественный);
P -	граница перехода, используемая при оценке погрешности решения (тип: вещественный);
H -	вещественная переменная, содержащая начальное значение шага интегрирования; может задаваться с учетом направления интегрирования, т.е. положительным, если XN < XK, отрицательным, если XN > XK, или без всякого учета в виде абсолютной величины; на выходе из подпрограммы содержит значение последнего шага интегрирования;
Y -	искомое решение задачи Коши, вычисленное подпрограммой при значении аргумента XK; для системы уравнений (когда M ≠ 1) задается одномерным массивом длины M; в случае совпадения значений параметров XN и XK значение Y полагается равным начальному значению YN (тип: вещественный);
RAB -	одномерный рабочий массив вещественного типа длиной 6*M + 1;
IERR -	целая переменная, значение которой в результате работы подпрограммы полагается равным 65, если какая - нибудь компонента решения не может быть вычислена с требуемой точностью EPS. B этом случае интегрирование системы прекращается; при желании интегрирование системы можно повторить обращением к подпрограмме с новыми значениями параметров HMIN и H.

   Версии

DE16E -

вычисление решения задачи Коши для системы обыкновенных дифференциальных уравнений первого порядка в конце интервала интегрирования методом Рунге - Кутта - Фельберга с удвоенным числом значащих цифр. При этом параметры XN, YN, XK, HMIN, EPS, P, H, Y, RAB и параметры X, Y и DY в подпрограмме F должны иметь тип Extended.

   Вызываемые подпрограммы

DE15R - DE15E	выполнение одного шага численного интегрирования системы обыкновенных дифференциальных уравнений первого порядка методом Рунге - Кутта - Фельберга пятого порядка точности.
UTDE16 - UTDE17	подпрограммы выдачи диагностических сообщений.

Подпрограммы DE15R и UTDE16 вызываются при работе подпрограммы DE16R, а подпрограммы DE15E и UTDE17 - при работе подпрограммы DE16E.

   Замечания по использованию

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

При работе подпрограммы значения параметров M, XN, YN, XK, HMIN, EPS, P сохраняются.

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

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

          y₁'  =  -y₂ - 0.1 x - 0.9
          y₂'  =   y₁ - 0.1 x - 1.1     0 ≤ x ≤ π
    
          y₁ (0)  =   1 ,     y₂ (0)  =  - 2
   
 Точное решение системы:

         y₁  =    0.1 x  +  sin x  +  1.
         y₂  =  - 0.1 x - cos x -

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

Unit TDE16R_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc, UtRes_p, FDE16R_p, DE16R_p;

function TDE16R: String;

implementation

function TDE16R: String;
var
M,_i,IERR :Integer;
XN,XK,HMIN,EPS,P,H,Y1,Y2 :Real;
YN :Array [0..1] of Real;
Y :Array [0..1] of Real;
RАВ :Array [0..12] of Real;
begin
Result := '';  { результат функции }
M := 2;
XN := 0.0;
YN[0] := 1.0;
YN[1] := -2.0;
ХК := 3.141592653589793238;
HMIN := 1.E-14;
EPS := 1.E-5;
P := 100.0;
H := 0.01;
DE16R(FDE16R,M,XN,YN,XK,HMIN,EPS,P,H,Y,RAB,IERR);
Result := Result + Format('%s',[' IERR=']);
Result := Result + Format('%4d ',[IERR]) + #$0D#$0A;
Y1 := 0.1*XK+Sin(XK)+1.0;
Y2 := -0.1*XK-Cos(XK)-1.0;
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;
Result := Result + Format('%20.16f  %20.16f  %20.16f ',
 [H,Y1,Y2]) + #$0D#$0A;
UtRes('TDE16R',Result);  { вывод результатов в файл TDE16R.res }
exit;
end;

end.

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

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

implementation

procedure fde16r(X :Real; var Y :Array of Real; var Z :Array of Real;
                M :Integer);
var
R :Real;
begin
R := 0.1*X;
Z[0] := -Y[1]-R-0.9;
Z[1] := Y[0]-R-1.1;
end;

end.

Результаты:

          IERR  =   0

          Y(1)  =   1.314159265929 + 00
          Y(2)  =  -3.141592563520 - 01