БЧА НИВЦ МГУ. DE45E_P. Вычисление решения задачи Коши для системы обыкновенных дифференциальных уравнений первого порядка в конце интервала интегрирования методом рядов Чебышёва

Текст подпрограммы и версий (Паскаль)
de45e_p.zip

Тексты тестовых примеров (Паскаль)
tde45e1_p.zip , tde45e2_p.zip , tde45e3_p.zip , tde45e4_p.zip

Подпрограмма: DE45E (модуль DE45E_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, которая является концом интервала интегрирования. Предполагается, что правая часть системы имеет непрерывные ограниченные частные производные по переменным X, Y. Тогда решение системы разлагается на промежутке интегрирования в равномерно сходящийся ряд по смещенным многочленам Чебышёва первого рода.

Интервал интегрирования разбивается на элементарные сегменты длины H: [x_s, x_s + H], x₀ = XN, s = 0,1,...,[(XK - XN) / H] - 1 (квадратные скобки означают целую часть числа). Если длина интервала интегрирования не является целым кратным H, то последний элементарный сегмент считается нестандартным; в этом случае s принимает значения s = 0,1,...,[(XK - XN)/H]. На каждом элементарном сегменте решение исходной задачи Коши приближенно представляется в виде (K + 1) - й частичной суммы смещенного ряда Чебышёва

                                    _K+1
           Y(x_s + αH)   ≈   ∑   a_i^*[Y] T_i^*(α) ,     0 ≤ α ≤ 1 ,
                                    ⁱ⁼⁰

здесь a_i^*[Y] - коэффициенты ряда, T_i^*(α) - смещенные многочлены Чебышёва первого рода на [0,1]. Значения H и K задаются пользователем при обращении к подпрограмме. Они выбираются таким образом, чтобы на каждом элементарном сегменте [x_s, x_s + H] ряд Чебышёва для решения задачи Коши (и его производной) являлся быстросходящимся рядом, а его сумма хорошо приближалась частичной суммой (K + 1) - го порядка (для производной - K - го порядка). Решение в конце каждого сегмента вычисляется по формуле

                                                _K+1
        Y(x_s + H)  =  Y(x_s+1)   ≈   ∑   a_i^*[Y] ,
                                                ⁱ⁼⁰

В качестве решения в конце интервала интегрирования принимается значение решения в конце последнего элементарного сегмента.

Коэффициенты a_i^*[Y] ряда Чебышёва для решения на сегменте [x_s, x_s + H] выражаются через коэффициенты a_i^*[Φ] ряда Чебышёва его производной Φ(α) = F(x_s + αH, Y(x_s + αH)), 0 ≤ α ≤ 1, на [x_s, x_s + H], которые, в свою очередь, вычисляются приближенно итерационным способом, исходя из некоторого начального приближения, с помощью квадратурной формулы Маркова на [x_s, x_s + H] с K + 1 узлом. При этом один из узлов квадратурной формулы совпадает с x_s, а остальные K узлов лежат внутри интервала (x_s, x_s + H). Количество итераций, которое предписывается выполнить в этом итерационном процессе, одинаково для всех сегментов и задается при обращении к подпрограмме. Если при выбранном H ряды Чебышёва для Y(X) = Y(x_s + αH), 0 ≤ α ≤ 1 , и его производной на элементарном сегменте [x_s, x_s + H] быстро сходятся, то для того, чтобы приближенное решение в конце одного такого сегмента имело максимальный порядок точности относительно H, необходимо выполнить не менее K итераций; при этом погрешность приближенного решения в конце элементарного сегмента является величиной порядка O(H^K + 2) при H --> 0. Если H подобрано достаточно малым, то хорошая точность приближенного решения может быть получена и при меньшем числе итераций. Вообще, число итераций зависит от K и H. С увеличением H число итераций может также возрастать.

Начальное приближение коэффициентов a_i^*[Φ] ряда Чебышёва для производной на сегменте [x_s, x_s + H] может вычисляться двумя способами. В первом способе начальное приближение определяется только с использованием значения решения Y(X) в узле x_s. При этом погрешность начального приближения для всех коэффициентов a₀^*[Φ], a₁^*[Φ], ..., a_K^*[Φ] является величиной O(H²) при H --> 0. Во втором способе начальное приближение определяется через коэффициенты ряда Чебышёва производной Φ(α) на предыдущем элементарном сегменте [x_s - 1, x_s]. В этом случае погрешности начального приближения для коэффициентов a₀^*[Φ], a₁^*[Φ], ..., a_K^*[Φ] имеют, соответственно, порядки O(H), O(H²), ..., O(H^K + 1). Второй способ определения начального приближения в некоторых случаях может привести к более быстрой сходимости итерационного процесса и, тем самым, к меньшему числу выполняемых итераций. Второй способ может быть применен только начиная со второго элементарного сегмента [x₀ + H, x₀ + 2H]. На начальном элементарном сегменте [x₀, x₀ + H] всегда применяется исключительно первый способ. Способ выбора начального приближения задается пользователем при обращении к подпрограмме.

В дальнейшем при описании параметров подпрограммы коэффициенты ряда Чебышёва будем называть коэффициентами Чебышёва.

О.Б.Арушанян, С.Ф.Залеткин. Приближенное решение задачи Коши для обыкновенных дифференциальных уравнений методом рядов Чебышёва. Вычислительные методы и программирование: Новые вычислительные технологии (Электронный научный журнал) (17), 121 - 131, 2016.

О.Б.Арушанян, Н.И.Волченскова, С.Ф.Залеткин. Вычисление коэффициентов разложения решения задачи Коши в ряд по многочленам Чебышёва. Вестник Московского университета. Серия 1: Математика. Механика. 5, 24 - 30, 2012.

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

procedure DE45E(F :Proc_F70E; var M :Integer; XN :Extended;
                var YN :Array of Extended; XK :Extended; var K :Integer;
                INIAPR :Integer; IMAX :Integer; var H :Extended;
                var Y :Array of Extended; var RAB :Array of Extended);

   Параметры

F -	имя подпрограммы вычисления значений правой части дифференциального уравнения. Первый оператор подпрограммы должен соответствовать процедурному типу:

 
       Procedure (X :Extended; var Y :Array of Extended; var Z :Array of Extended; M :Integer);

Здесь: X, Y - значения независимой и зависимой переменных, соответственно. Вычисленное значение правой части должно быть помещено в Z. Для скалярного случая при M = 1 параметры Y и Z представляют собой массивы длины 1, в случае системы уравнений, т.е. когда M ≠ 1 , параметры Y и Z представляют массивы длины M (тип параметров X, Y и Z: Extended);

M -	количество уравнений в системе (тип: целый);
XN, YN -	начальные значения аргумента и решения; в случае системы уравнений (т.е. когда M ≠ 1) YN представляет массив длины M (тип: Extended);
XK -	значение аргумента, при котором требуется вычислить решение задачи Коши (конец интервала интегрирования); XK может быть больше, меньше или равно XN (тип: Extended);
K -	порядок частичной суммы смещенного ряда Чебышёва, с помощью которой аппроксимируется производная решения задачи Коши на каждом элементарном сегменте разбиения интервала интегрирования; при этом само решение задачи Коши приближается на каждом элементарном сегменте частичной суммой (K + 1) - го порядка; K≥2 (см. "Математическое описание" и "Замечания по использованию"; тип: целый);
INIAPR -	целый указатель способа выбора начального приближения коэффициентов Чебышёва для производной решения на каждом элементарном сегменте:

INIAPR=1 -

для первого способа, когда начальное приближение определяется только с использованием значения решения в начале каждого элементарного сегмента;

INIAPR=2 -

для второго способа, когда начальное приближение коэффициентов Чебышёва на текущем элементарном сегменте (начиная со второго) определяется через коэффициенты Чебышёва, вычисленные на предыдущем элементарном сегменте, т.е. путем экстраполяции коэффициентов с предыдущего сегмента на следующий (см. "Математическое описание");

IMAX -	целая переменная, задающая количество итераций, которое предполагается выполнить в итерационном процессе вычисления коэффициентов Чебышёва для производной решения задачи Коши на каждом элементарном сегменте, исходя из некоторого начального приближения, способ определения которого задается параметром INIAPR; IMAX≥1. Для получения максимального порядка точности приближенного решения необходимо выполнить не менее K итераций (см. "Математическое описание" и "Замечания по использованию");
H -	переменная с двойной точностью, содержащая значение длины элементарных сегментов, на которые разбивается интервал интегрирования (диаметр разбиения промежутка интегрирования или аналог шага интегрирования для одношаговых методов). Предполагается, что на каждом элементарном сегменте ряды Чебышёва для решения и его производной являются быстросходящимися рядами. Может задаваться с учетом направления интегрирования, т.е. положительным, если XK > XN, отрицательным, если XK < XN , или без такого учета в виде абсолютной величины;
Y -	искомое решение задачи Коши, вычисленное подпрограммой при значении аргумента XK; для системы уравнений (когда M ≠ 1) задается массивом длины M. В случае совпадения значений параметров XN и XK значение Y полагается равным начальному значению YN (тип: Extended);
RAB -	одномерный рабочий массив длины 2 * K² + 7 * K + 5 * M * K + 8 * М + 1. (тип: Extended).

   Версии:  нет

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

DE44E -

выполнение одного шага приближенного интегрирования системы обыкновенных дифференциальных уравнений первого порядка методом рядов Чебышёва

Кроме того, используются рабочие подпрограммы DE70EK, DE70EH, DE70E0, DE70EI, DE70EF, DE70EQ, DE71EE, DE70EP, DE71ET, DE71EP, DE71EI, DE71EF, DE71ES, DE70EA, DE70EC.

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

Разбиение промежутка интегрирования на элементарные сегменты (шаги) выполняется для того, чтобы на каждом таком сегменте ряды Чебышёва для решения и его производной были быстросходящимися рядами. Другими словами, длина элементарных сегментов, задаваемая параметром H, подбирается таким образом, чтобы убывание коэффициентов этих рядов Чебышёва на элементарном сегменте происходило достаточно быстро, вследствие чего можно было бы считать частичные суммы этих рядов близкими к многочленам наилучшего равномерного приближения на элементарном сегменте для решения и его производной. Порядок этих частичных сумм задается параметром K.

Если начальное приближение для коэффициентов Чебышёва функции Φ(α) определяется первым способом (т.е. при INIAPR = 1), то для получения максимального порядка точности приближенного решения в конце элементарного сегмента необходимо выполнить в итерационном процессе не менее K итераций; тогда IMAX≥K. Если начальное приближение коэффициентов Чебышёва функции Φ(α) определяется вторым способом (т.е. при INIAPR = 2), то для получения максимального порядка точности приближенного решения необходимо выполнить в итерационном процессе не менее K + 1 итераций; в этом случае IMAX≥K + 1. Однако в некоторых случаях при втором способе определения начального приближения итерационный процесс может сойтись за значительно меньшее число итераций.

Если диаметр разбиения H подобран достаточно малым, то хорошая точность приближенного решения может быть получена и с существенно меньшим числом итераций при любом способе выбора начального приближения. Вообще, число итераций зависит от K и H. С увеличением H число итераций может также возрастать.

Если правая часть дифференциального уравнения не зависит от переменной Y, т.е. дифференциальное уравнение имеет вид Y' = F(X), то число итераций можно положить равным 1 при любых H и K, удовлетворяющих описанным выше условиям. В этом случае параметр IMAX = 1.

Как следует из вышеописанного, управлять точностью приближенного решения задачи Коши можно с помощью четырех параметров H, K, IMAX, INIAPR, подбирая для каждой конкретной задачи наиболее подходящий набор их значений.

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

Так как при интегрировании системы уравнений с помощью подпрограммы DE45E используется глобальная запись (структура данных) с именем _COM70D для хранения промежуточных значений, пользователь не должен портить элементы этой записи.

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

Вычисления во всех четырех примерах на Паскале проводились с 19-20 значащими цифрами.

 1)         y₁'  =   x/y₂ ,    y₁(0)  =  3  ,
             y₂'  =  -x/y₁ ,     y(0)  =  1/6,            0 ≤ x ≤ 3√2 .
     
        Точное решение системы:
                 y₁(x)  =  3*e^x*x  ,     y₂(x)  =  (1/6)*e^-x*x .

Решение вычисляется в точке XK = 3√2. Приводятся подпрограмма вычисления значений правой части и фрагмент вызывающей программы, а также результаты счета, включая абсолютную DELY и относительную RELERR погрешности приближенного решения,точные значения решения YT, значения параметров H, K, IMAX, INIAPR при двух способах определения начального приближения коэффициентов Чебышёва производной решения.

Unit Fde45e1_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc;
procedure fde45e1(X :Extended; var Y :Array of Extended; var Z :Array of Extended;
                M :Integer);
implementation
procedure fde45e1(X :Extended; var Y :Array of Extended; var Z :Array of Extended;
                M :Integer);
begin
Z[0] := X/Y[1];
Z[1] := -X/Y[0];
end;
end.

function tde45e1: String;
var
M,K,IMAX,INIAPR,_i :Integer;
XN,H,XK,H0 :Extended;
YN :Array [0..1] of Extended;
Y :Array [0..1] of Extended;
RAB :Array [0..1691] of Extended;
YT :Array [0..1] of Extended;
DELY :Array [0..1] of Extended;
RELERR :Array [0..1] of Extended;
begin
M := 2;
XN := 0.0;
YN[0] := 3.e0;
YN[1] := 1.e0/6.e0;
XK := 3.e0*Sqrt(2.e0);
YT[0] := 3.e0*Exp(XK*XK);
YT[1] := 1.e0/6.e0*Exp(-XK*XK);
H0 := 0.1e0;
K := 25;
IMAX := 31;
for INIAPR:=1 to 2 do
 begin
  H := H0;
  DE45E(Fde45e1,M,XN,YN,XK,K,INIAPR,IMAX,H,Y,RAB);
  DELY[0] := YT[0]-Y[0];
  DELY[1] := YT[1]-Y[1];
  RELERR[0] := DELY[0]/Y[0];
  RELERR[1] := DELY[1]/Y[1];
// Операторы вывода на печать: H, K, INIAPR, IMAX, Y, YT, DELY, RELERR 
. . . . . . . . . . . . . . . . . . . . . . . . . . 
end;

Результаты:

 H0=  0,1000000000000000    K=25    INIAPR=1    IMAX= 31 
 Y= 
      1,96979907411991536E+008       2,53832995745210471E-009 
 YT=
      1,96979907411991533E+008       2,53832995745210474E-009 
 DELY  =
     -2,59024091064929962E-009       3,39273298117452757E-026 
 RELERR   =
     -1,31497722010382753E-017       1,33660045701073693E-017 
 ----------------------------------------------------------------
 H0=  0,1000000000000000    K=25    INIAPR=2    IMAX= 31 
 Y= 
      1,96979907411991535E+008       2,53832995745210471E-009 
 YT=
      1,96979907411991533E+008       2,53832995745210474E-009 
 DELY  =
     -2,15368345379829407E-009       2,88786200183308002E-026 
 RELERR   =
     -1,09335184592902514E-017       1,13770157947937726E-017 
 ----------------------------------------------------------------

Видно, что при втором способе выбора начального приближения коэффициентов Чебышёва производной (т.е. при INIAPR = 2) приближенное решение вычислено с лучшей точностью, чем при INIAPR = 1, т.е. с меньшей абсолютной и относительной погрешностью.

2) Решается задача Коши из примера 1 с другими значениями параметров H, IMAX: H = 0,4, IMAX = 39. Так как значение H здесь увеличено по сравнению с примером 1, то и число итераций также возрастает. Приводятся результаты счета, включающие те же данные, что и в примере 1.

 H0=  0,4000000000000000    K=25    INIAPR=1    IMAX= 39 
 Y= 
      1,96979907411991521E+008       2,53832995745210494E-009 
 YT=
      1,96979907411991533E+008       2,53832995745210474E-009 
 DELY  =
      1,20780896395444870E-008      -1,91850972149750071E-025 
 RELERR   =
      6,13163535216953334E-017      -7,55815734619166647E-017 
 ----------------------------------------------------------------
 H0=  0,4000000000000000    K=25    INIAPR=2    IMAX= 39 
 Y= 
      1,96979907411991530E+008       2,53832995745210479E-009 
 YT=
      1,96979907411991533E+008       2,53832995745210474E-009 
 DELY  =
      3,23052518069744110E-009      -4,64481300994131751E-026 
 RELERR   =
      1,64002776889353776E-017      -1,82986967328850883E-017 
 ----------------------------------------------------------------

Здесь так же, как и в примере 1, при втором способе выбора начального приближения (т.е. при INIAPR = 2) приближенное решение вычислено с лучшей точностью, чем при INIAPR = 1, т.е. с меньшей абсолютной и относительной погрешностью.

3) Решается система обыкновенных дифференциальных уравнений из примера 1 на отрезке -3√2 ≤ x ≤ 0. Начальные условия задаются в точке t = -3√2:

     y₁(t)  =  3*e^t*t ,    y₂(t)  =  (1/6)*e^-t*t  .

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

Unit Fde45e3_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc;
procedure fde45e3(X :Extended; var Y :Array of Extended; var Z :Array of Extended;
                M :Integer);
implementation
procedure fde45e3(X :Extended; var Y :Array of Extended; var Z :Array of Extended;
                M :Integer);
begin
Z[0] := X/Y[1];
Z[1] := -X/Y[0];
end;
end.

function tde45e3: String;
var
M,K,IMAX,INIAPR,_i :Integer;
XN,H,XK,H0 :Extended;
YN :Array [0..1] of Extended;
Y :Array [0..1] of Extended;
RAB :Array [0..1691] of Extended;
YT :Array [0..1] of Extended;
DELY :Array [0..1] of Extended;
RELERR :Array [0..1] of Extended;
begin
M := 2;
XN := -3.e0*Sqrt(2.e0);
YN[0] := 3.e0*Exp(XN*XN);
YN[1] := 1.e0/6.e0*Exp(-XN*XN);
YT[0] := 3.e0;
YT[1] := 1.e0/6.e0;
XK := 0.0;
H0 := 0.1e0;
K := 25;
IMAX := 30;
for INIAPR:=1 to 2 do
 begin
  H := H0;
  DE45E(Fde45e3,M,XN,YN,XK,K,INIAPR,IMAX,H,Y,RAB);
  DELY[0] := YT[0]-Y[0];
  DELY[1] := YT[1]-Y[1];
  RELERR[0] := DELY[0]/Y[0];
  RELERR[1] := DELY[1]/Y[1];
// Операторы вывода на печать: H, K, INIAPR, IMAX, Y, YT, DELY, RELERR 
. . . . . . . . . . . . . . . . . . . . . . . . . . 
end;

Результаты:

 H0=  0,1000000000000000    K=25    INIAPR=1    IMAX= 30 
 Y= 
      3,00000000000000000E+000       1,66666666666666667E-001 
 YT=
      3,00000000000000000E+000       1,66666666666666667E-001 
 DELY  =
     -1,51788304147970621E-018      -1,62630325872825665E-019 
 RELERR   =
     -5,05961013826568736E-019      -9,75781955236953990E-019 
 ----------------------------------------------------------------
 H0=  0,1000000000000000    K=25    INIAPR=2    IMAX= 30 
 Y= 
      3,00000000000000000E+000       1,66666666666666667E-001 
 YT=
      3,00000000000000000E+000       1,66666666666666667E-001 
 DELY  =
     -4,11996825544491685E-018       2,71050543121376108E-020 
 RELERR   =
     -1,37332275181497228E-018       1,62630325872825665E-019 
 ----------------------------------------------------------------

4) Решается система обыкновенных дифференциальных уравнений из примера 1, но начальное условие задается в точке t = 3√2:

     y₁(t)  =  3*e^t*t ,    y₂(t)  =  (1/6)*e^-t*t  ,

а решение задачи Коши вычисляется в точке XK = 0. Интегрирование системы выполняется справа налево с отрицательным шагом H = -0.1. Приводятся фрагмент вызывающей программы и результаты счета, включающие такие же данные, как и в примере 1.

Unit Fde45e4_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc;
procedure fde45e4(X :Extended; var Y :Array of Extended; var Z :Array of Extended;
                M :Integer);
implementation
procedure fde45e4(X :Extended; var Y :Array of Extended; var Z :Array of Extended;
                M :Integer);
begin
Z[0] := X/Y[1];
Z[1] := -X/Y[0];
end;
end.

function tde45e4: String;
var
M,K,IMAX,INIAPR,_i :Integer;
XN,H,XK,H0 :Extended;
YN :Array [0..1] of Extended;
Y :Array [0..1] of Extended;
RAB :Array [0..1691] of Extended;
YT :Array [0..1] of Extended;
DELY :Array [0..1] of Extended;
RELERR :Array [0..1] of Extended;
begin
Result := '';  { результат функции }
M := 2;
XN := 3.e0*Sqrt(2.e0);
YN[0] := 3.e0*Exp(XN*XN);
YN[1] := 1.e0/6.e0*Exp(-XN*XN);
YT[0] := 3.e0;
YT[1] := 1.e0/6.e0;
XK := 0.0;
H0 := -0.1e0;
K := 25;
IMAX := 30;
for INIAPR:=1 to 2 do
 begin
  H := H0;
  DE45E(Fde45e4,M,XN,YN,XK,K,INIAPR,IMAX,H,Y,RAB);
  DELY[0] := YT[0]-Y[0];
  DELY[1] := YT[1]-Y[1];
  RELERR[0] := DELY[0]/Y[0];
  RELERR[1] := DELY[1]/Y[1];
// Операторы вывода на печать: H, K, INIAPR, IMAX, Y, YT, DELY, RELERR 
. . . . . . . . . . . . . . . . . . . . . . . . . . 
end;

Результаты:
 H0= -0,1000000000000000    K=25    INIAPR=1    IMAX= 30 
 Y= 
      3,00000000000000000E+000       1,66666666666666667E-001 
 YT=
      3,00000000000000000E+000       1,66666666666666667E-001 
 DELY  =
     -1,51788304147970621E-018      -1,62630325872825665E-019 
 RELERR   =
     -5,05961013826568736E-019      -9,75781955236953990E-019 
 ----------------------------------------------------------------
 H0= -0,1000000000000000    K=25    INIAPR=2    IMAX= 30 
 Y= 
      3,00000000000000000E+000       1,66666666666666667E-001 
 YT=
      3,00000000000000000E+000       1,66666666666666667E-001 
 DELY  =
     -4,11996825544491685E-018       2,71050543121376108E-020 
 RELERR   =
     -1,37332275181497228E-018       1,62630325872825665E-019 
 ----------------------------------------------------------------

Примеры 3 и 4 относятся к тому случаю, когда решение задачи Коши вычисляется при первом способе выбора начального приближения для коэффициентов Чебышёва производной (т.е. при INIAPR = 1) с лучшей точностью, чем при INIAPR = 2, т.е. с меньшей абсолютной и относительной погрешностью.