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

Текст подпрограммы и версий ( Фортран ) de76d.zip	Тексты тестовых примеров ( Фортран ) tde76d1.zip , tde76d2.zip , tde76d3.zip , tde76d4.zip
Текст подпрограммы и версий ( Си ) de76d_c.zip	Тексты тестовых примеров ( Си ) tde76d1_c.zip , tde76d2_c.zip , tde76d3_c.zip , tde76d4_c.zip
Текст подпрограммы и версий ( Паскаль ) de76e_p.zip	Тексты тестовых примеров ( Паскаль ) tde76e1_p.zip , tde76e2_p.zip , tde76e3_p.zip , tde76e4_p.zip

Подпрограмма: DE76D

Назначение

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

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

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

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

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

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

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

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

Коэффициенты a_i^*[U_K+1] ряда Чебышёва для решения U_K+1(x_s+αH) на s - м сегменте [x_s, x_s + H] выражаются через коэффициенты a_i^*[Φ] ряда Чебышёва его производной Φ(α) = F(x_s + αH, U_K+1(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 подобрано достаточно малым, то хорошая точность приближенного решения может быть получена и при меньшем числе итераций.

Начальное приближение коэффициентов a_i^*[Φ] ряда Чебышёва для производной на сегменте [x_s, x_s + H] может вычисляться двумя способами. В первом способе начальное приближение определяется только с использованием значения решения U_K+1(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₁, x₂]. На начальном элементарном сегменте [x₀, x₁] всегда применяется исключительно первый способ. Способ выбора начального приближения задается пользователем при обращении к подпрограмме.

На каждом частичном сегменте [x_s, x_s + H] = [x_s, x_s+1] для оценки погрешности решения U_K+1(x_s + αH), полученного при некотором заданном значении K = K₁, вычисляется на этом же сегменте специальным образом второе приближенное решение при K = K₂ > K₁, которое также представляется указанными выше частичными суммами при K = K₂. Значение K₂ задается пользователем при обращении к подпрограмме. Второе решение U_K₂+1, как имеющее более высокий порядок точности O(H^K₂ + 2), используется для оценки погрешности первого решения U_K₁+1 на каждом элементарном сегменте. Заданные значения K₁ и K₂ используются для всех сегментов, на которые автоматически разбивается интервал интегрирования [XN, XK]. Второе приближенное решение на сегменте [x_s, x_s+1] вычисляется аналогично первому решению итерационным способом. Единственным отличием при этом является то, что начальным приближением в этом дополнителном итерационном процессе служит уже вычисленное на отрезке [x_s, x_s+1] первое решение U_K₁+1(x_s + αH). Количество итераций, которое потребуется выполнить в этом дополнительном итерационном процессе, одинаково для всех сегментов и задается при обращении к подпрограмме.

Длина элементарного сегмента H выбирается в подпрограмме из условия, чтобы на каждый элементарный сегмент [x_s, x_s + H] приходилась приблизительно одинаковая погрешность. Поэтому длина элементерного сегмента H - величина переменная, определяемая автоматически подпрограммой. Она изменяется от одного элементарного сегмента к другому элементарному сегменту. При обращении к подпрограмме задается только начальное значение длины первого элементерного сегмента [XN, XN + H] = [x₀, x₀ + H] которое может быть уменьшено подпрограммой, если на нем не будет достигнута заданная точность.

В качестве решения Y(x_s+1) задачи на сегменте [x_s, x_s+1] выбирается второе приближенное решение U_K₂+1(x_s+1), поскольку оно имеет более высокий порядок точности O(H^K₂ + 2). В качестве решения исходной задачи Коши в конце интервала интегрирования XK принимается значение решения U_K₂+1(XK) в конце последнего элементарного сегмента.

Для системы уравнений (M > 1) с проверкой на точность могут вычисляться либо все компоненты решения, либо некоторые из них (в частности, одна компонента). В последнем случае номера проверяемых на точность компонент задаются при обращении к подпрограмме.

Для оценки уклонения приближенного решения от точного (т.е. для оценки абсолютной погрешности приближенного решения) в подпрограмме используются два способа, или две формулы. Первый способ опирается на первую формулу для абсолютной погрешности приближенного решения. Она является асимптотической (т.е. справедливой при H --> 0) и представляет собою разность двух частичных сумм смещенного ряда Чебышёва на сегменте [x_s, x_s + H] с приближенными коэффициентами, а именно разность приближенных частичных сумм (K₂+1)-го и (K₁+1)-го порядков в точке x_s + H. Второй способ опирается на несколько завышенную оценку абсолютной погрешности, которая представляет сумму абсолютных величин (модулей) разностей приближенных коэффициентов этих частичных сумм и поэтому несколько завышена по сравнению с первой. Назовем эту формулу второй. Более того, при использовании этой, второй, формулы в том случае, когда точность решения оценивается по относительной погрешности или по мере погрешности, то дополнительно используется также завышенная оценка сверху для относительной погрешности и завышенная оценка сверху для меры погрешности. Второй способ оценки погрешности, хотя и является более точным, накладывает более сильное ограничение на размер элементарного сегмента и приводит к меньшей длине сегментов. Используемый способ оценки абсолютной погрешности задается при обращении к подпрограмме.

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

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

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

О.Б.Арушанян, С.Ф.Залеткин. Об оценке погрешности приближенного решения обыкновенных дифференциальных уравнений, определенного с помощью рядов Чебышёва// Вычислительные методы и программирование, 2020, т. 21, с. 241 - 250.

О.Б.Арушанян, С.Ф.Залеткин. О вычислении приближенного решения обыкновенных дифференциальных уравнений методом рядов Чебышёва и оценка его погрешности// Вестн. Моск. ун-та. Сер. 1, Математика. Механика. 2020, № 5 с. 22 - 26.

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

    SUBROUTINE  DE76D(F, M, XN, YN, XK, K, INIAPR, IMAX, IU, EPS, THRESH, MEXACT, NUMBES, MESTER, H, HMIN, K2, IMAX2, NATTEM, Y, RAB, IERR)

   Параметры

F -	имя подпрограммы вычисления значений правой части дифференциального уравнения. Первый оператор подпрограммы должен иметь вид: SUBROUTINE F (X, Y , DY, M). Здесь: X, Y - значения независимой и зависимой переменных, соответственно. Вычисленное значение правой части должно быть помещено в DY. B случае системы уравнений, т.е. когда M ≠ 1 , параметры Y и DY представляют массивы длины M (тип параметров X, Y и DY: с двойной точностью);
M -	количество уравнений в системе (тип: целый);
XN, YN -	начальные значения аргумента и решения; в случае системы уравнений (т.е. когда M ≠ 1) YN представляет массив длины M (тип: с двойной точностью);
XK -	значение аргумента, при котором требуется вычислить решение задачи Коши (конец интервала интегрирования); XK может быть больше, меньше или равно XN (тип: с двойной точностью);
K -	порядок частичной суммы смещенного ряда Чебышёва, с помощью которой аппроксимируется производная первого приближенного решения задачи Коши на каждом элементарном сегменте разбиения интервала интегрирования; при этом само приближенное решение задачи Коши приближается на каждом элементарном сегменте частичной суммой (K + 1) - го порядка; K≥2. Этот параметр относится только к первому приближенному решению, имеющему меньший порядок точности O(H^K + 2), и его производной (см. "Математическое описание" и "Замечания по использованию"; тип: целый);
INIAPR -	целый указатель способа выбора начального приближения коэффициентов Чебышёва для производной первого приближенного решения на каждом элементарном сегменте:

INIAPR=1 -

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

INIAPR=2 -

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

IMAX -

целая переменная, задающая количество итераций, которое предполагается выполнить в итерационном процессе вычисления коэффициентов Чебышёва для производной первого приближенного решения задачи Коши на каждом элементарном сегменте, исходя из некоторого начального приближения, способ определения которого задается параметром INIAPR; IMAX≥1. Для получения максимального порядка точности первого приближенного решения необходимо выполнить не менее K итераций. Этот параметр относится только к первому приближенному решению, имеющему меньший порядок точности O(H^K + 2), и его производной. Если правая часть дифференциального уравнения не зависит от переменной Y, т.е. дифференциальное уравнение имеет вид Y' = F(X), то число итераций при вычислении первого приближенного решения можно положить равным 1. В этом случае параметр IMAX=1 (см. "Математическое описание" и "Замечания по использованию");

IU -

целый указатель типа погрешности численного решения:

IU=1 -	компоненты решения проверяются на точность по абсолютной погрешности;
IU=2 -	компоненты решения проверяются на точность по относительной погрешности;
IU=3 -	компоненты решения проверяются на точность по мере погрешности.

	Контроль точности по мере погрешности заключается в следующем. Если некоторая компонента приближенного решения по абсолютной величине не меньше определенной наперед заданной положительной константы THRESH (называемой границей перехода), то контроль точности для этой компоненты ведется по относительной погрешности, иначе - по абсолютной. Таким образом, контроль точности по мере погрешности состоит в том, что на тех участках интервала интегрирования, где абсолютная величина компоненты решения меньше значения THRESH, контроль точности ведется для нее по абсолютной погрешности, а там, где абсолютная величина этой компоненты решения равна значению THRESH или превосходит значение THRESH, контроль точности для нее ведется по относительной погрешности;
EPS -	допустимая погрешность, с которой требуется вычислить проверяемые на точность компоненты решения; тип погрешности специфицируется с помощью параметра IU (тип: с двойной точностью);
THRESH -	граница перехода, используемая при оценке меры погрешности решения (тип: с двойной точностью);
MEXACT -	количество компонент численного решения, проверяемых на точность: 1 ≤ MEXACT ≤ M. Для скалярного дифференциального уравнения значение этого параметра в подпрограмме не используется, т.к. в этом случае численное решение обязательно проверяется на точность (тип: целый);
NUMBES -	массив, содержащий номера проверяемых на точность компонент решения. Для скалярного дифференциального уравнения, а также для системы уравнений в том случае, если все компоненты решения проверяются на точность, значение этого массива в подпрограмме не используется. Таким образом, массив NUMBES используется в подпрограмме только в том случае, когда число проверяемых на точность компонент решения меньше числа уравнений; NUMBES представляет массив длины MEXACT (тип: целый);
MESTER -	целый указатель способа оценки абсолютной погрешности компоненты приближенного решения:

MESTER=1 -	для оценки абсолютной погрешности используется первая формула (асимптотическая);
MESTER=2 -	для оценки абсолютной погрешности используется вторая формула (завышенная оценка);

	(См. "Математическое описание");
H -	переменная с двойной точностью, содержащая начальное значение длины первого элементарного сегмента (аналог шага интегрирования для разностных методов). Может задаваться с учетом направления интегрирования, т.е. положительным, если XK > XN, отрицательным, если XK < XN , или без такого учета в виде абсолютной величины;
HMIN -	минимальное значение длины элементарного сегмента, которое разрешается использовать при разбиении интервала интегрирования на элементарные (частичные) сегменты (тип: с двойной точностью);
K2 -	порядок частичной суммы смещенного ряда Чебышёва производной решения, с помощью которой вычисляется частичная сумма порядка (K2+1) для оценивающего решения более высокого порядка точности на каждом элементарном сегменте разбиения области интегрирования. Этот параметр относится только ко второму приближенному решению более высокого порядка точности O(H^K2 + 2) , K2 > K (т.е. к оценивающему решению) и его производной (см. "Математическое описание" и "Замечания по пользованию"; тип: целый);
IMAX2 -	целая переменная, задающая количество итераций, которое предполагается выполнить в дополнительном итерационном процессе для вычисления второго, оценивающего, решения более высокого порядка точности. Если правая часть дифференциального уравнения не зависит от переменной Y, т.е. дифференциальное уравнение имеет вид Y' = F(X), то число итераций при вычислении второго приближенного решения можно положить равным 1. В этом случае параметр IMAX2=1 (см. "Математическое описание" и "Замечания по использованию");
NATTEM -	целая переменная, значение которой ограничивает число последовательных сокращенипй длины элементарного сегмента, совершаемых в одной точке, если приближенное решение на этом элементарном сегменте не достигает требуемой точности;
Y -	искомое решение задачи Коши, вычисленное подпрограммой при значении аргумента XK; в качестве такого решения принимается значение второго, оценивающего, приближенного решения, поскольку оно имеет более высокий порядок точности по сравнению с первым приближенным решением. Для системы уравнений (когда M ≠ 1) задается массивом длины M. В случае совпадения значений параметров XN и XK значение Y полагается равным начальному значению YN (тип: с двойной точностью);
RAB -	одномерный рабочий массив длины 2 * K² + 5 * K + 5 * M * K + K2 * K2 +5 * K2 + 4 * K2 * M + 18 * М + 1.(тип: с удвоенной точностью);
IERR -	целая переменная, служащая для сообщения об ошибках, обнаруженных в процессе работы подпрограммы; при этом:

IERR=65 -
IEER=66

когда какая-нибудь компонента решения не может быть вычислена с требуемой точностью EPS; при этом IERR=65 указывает, что требуемая точность не может быть достигнута при минимальной длине частичного сегмента, равной HMIN, а IERR=66 показывает, что требуемая точность не может быть достигнута, т.к. исчерпано заданное число сокращений NATTEM длины элементарного сегмента в одной и той же точке интервала интегрирования.

Версии: нет

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

DE75D -

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

Кроме того, используются рабочие подпрограммы DE70DK, DE70DH, DE70D0, DE70DI, DE70DF, DE70DQ, DE71DE, DE70DP, DE71DT, DE71DP, DE71DI, DE71DF, DE71DS, DE70DA, DE70DC, DE750, DE75D1, DE75DE, DE75DK, DE75DN.

UTDE75 - подпрограмма выдачи диагностических сообщений.

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

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

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

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

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

Для получения максимального порядка точности второго приближенного (оценивающего) решения необходимо выполнить в дополнительном итерационном процессе не менее (K2 - K) итераций. Указанное в данном разделе число итераций IMAX и IMAX2 носит асимптотический характер (когда все наши оценки справедливы при H --> 0). На практике же количество итераций в каждом из двух итерационных процессов зависит от длины сегмента, от заданной точности, от используемых порядков частичных сумм и от конкретной системы дифференциальных уравнений. Поэтому число итераций в обоих итерационных процессах, т.е. значения параметров IMAX и IMAX2, могут быть как меньше, так и больше рекомендуемых здесь значений.

При работе подпрограммы значения параметров M, XN, YN, XK, K, INIAPR, IMAX, IU, EPS, THRESH, MEXACT, NUMBES, MESTER, HMIN, K2, IMAX2, NATTEM сохраняются. Если после работы подпрограммы нет необходимости иметь начальное значение решения, то параметры YN и Y при обращении к ней можно совместить.

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

При работе подпрограммы DE76D используется общий блок

COMMON/STAT76/NACCEP, NREJEC

для сбора полезной статистической информации. Целая переменная NACCEP после выхода из подпрограммы содержит число элементарных сегментов, на которые автоматически разбивался весь промежуток интегрирования [XN, XK], а целая переменная NREJEC содержит число отклоненных элементарных сегментов, а именно тех сегментов, на которых вычисленная оценка погрешности превышает заданное значение EPS, т.е. число таких сегментов, на которых не достигается заданная точность и длины которых подвергались сокращению. Когда приближенное решение задачи Коши не может быть вычислено в конце интервала интегрирования XK и работа подпрограммы DE76D прерывается в некоторой точке интервала интегрирования при значении параметра IERR=65 или IERR=66, то значение переменной NACCEP содержит число элементарных сегментов, считая от начала промежутка интегрирования XN до этой точки прерывания, на которых приближенное решение вычислено с требуемой точностью. Значение переменной NREJEC содержит в этом случае число отклоненных сегментов, считая от начала промежутка интегрирования и включая число отклоненных сегментов с началом в этой точке прерывания.

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

Данные примеры не только иллюстрируют правила использования подпрогрвммы DE76D, но и показывают способность реализованного в подпрограмме метода рядов Чебышёва вычислять приближенное решение задачи Коши со стабильно высокой точностью. Вычисления во всех четырех примерах на Фортране проводились с 15-16 значащими цифрами.

 1)         y'  =  qy,    y(0)  =  e^q,    q = 4,   0 ≤ x ≤ 7.

Точное решение задачи Коши y(x) = e^q(1+x). Решение вычисляется в точке XK = 7. Приводятся фрагмент вызывающей программы, подпрограмма вычисления значений правой части, а также результаты счета, включая точное значение решения YT, относительную погрешность приближенного решения (YT - Y) / Y, значения переменных NACCEP и NREJEC из общего блока и число обращений к правой части N.

DOUBLE PRECISION XN,YN,XK,EPS,THRESH,H,HMIN,Y,RAB,Q,YT
      PARAMETER (K=18,K2=25)
      DIMENSION RAB(2*K*K+5*K+5*K+K2*K2+5*K2+4*K2+18+1)
      EXTERNAL F
      COMMON/STAT76/NACCEP,NREJEC
      COMMON/BLOCK/N,Q
      M=1
      Q=4.0D0
      XN=0.
      YN=EXP(Q)
      XK=7.D0
      YT =EXP(Q*(XK+1.D0))
      K0=K
      INIAPR=1
      IMAX=28
      IU=2
      EPS=0.5D-11
      THRESH=1.D0
      MEXACT=1
      MESTER=1
      H=1.D0
      HMIN=1.D-3
      K20=K2
      IMAX2=3
      NATTEM=3
      N=0  
      CALL DE76D(F,M,XN,YN,XK,K0,INIAPR,IMAX,IU,EPS,THRESH,MEXACT,
     1           M,MESTER,H,HMIN,K20,IMAX2,NATTEM,Y,RAB,IERR)
      PRINT 50,IERR
      PRINT 49,N
      PRINT 52,NACCEP,NREJEC
      PRINT 51,Y,YT,(YT-Y)/Y
      STOP
   49 FORMAT(1X,'N=',I12)
   50 FORMAT(1X,'IERR=',I2)
   51 FORMAT(1X,'Y =',D25.16/1X,'YT=',D25.16/1X,'(YT-Y)/Y=',D25.16)
   52 FORMAT(1X, 'NACCEP=',I4,5X,'NREJEC=',I4)
      END
C
      SUBROUTINE F(X,Y,Z,M)
      DOUBLE PRECISION X,Y,Z,Q
      COMMON/BLOCK/N,Q
      Z=Q*Y
      N=N+1
      RETURN
      END

Результаты:

----------------------------------------------------------

IERR= 0
N=        3330
NACCEP=   5     NREJEC=   0
Y =    .7896296018267691D+14
YT=    .7896296018268069D+14
(YT-Y)/Y=    .4788637598251465D-13

-----------------------------------------------------------

2) Решается задача Коши из примера 1 со значением параметра MESTER=2. Приводятся результаты счета, включая те же данные, что и в примере 1.


IERR= 0
N=        4662
NACCEP=   7     NREJEC=   0
Y =    .7896296018267920D+14
YT=    .7896296018268069D+14
(YT-Y)/Y=    .1879837073693702D-13

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

3) Решается задача Коши из примера 1 со значением EPS = 0.5D^-13. Приводятся результаты счета, включая те же данные, что и в примере 1.

IERR= 0
N=        3996
NACCEP=   6     NREJEC=   0
Y =    .7896296018268075D+14
YT=    .7896296018268069D+14
(YT-Y)/Y=   -.7915103468183855D-15

Видно, что уменьшение параметра EPS приводит к увеличению числа сегментов и повышению точности приближенного решения.

4) Решается задача Коши из примера 1 со значением параметра INIAPR=2 и со значениями параметров IMAX=25, IMAX2=6. Приводятся фрагмент вызывающей программы и результаты счета, включая те же данные, что и в примере 1.

      DOUBLE PRECISION XN,YN,XK,EPS,THRESH,H,HMIN,Y,RAB,Q,YT
      PARAMETER (K=18,K2=25)
      DIMENSION RAB(2*K*K+5*K+5*K+K2*K2+5*K2+4*K2+18+1)
      EXTERNAL F
      COMMON/STAT76/NACCEP,NREJEC
      COMMON/BLOCK/N,Q
      M=1
      Q=4.0D0
      XN=0.
      YN=EXP(Q)
      XK=7.D0
      YT =EXP(Q*(XK+1.D0))
      K0=K
      INIAPR=2
      IMAX=25
      IU=2
      EPS=0.5D-11
      THRESH=1.D0
      MEXACT=1
      MESTER=1
      H=1.D0
      HMIN=1.D-3
      K20=K2
      IMAX2=6
      NATTEM=3
      N=0  
      CALL DE76D(F,M,XN,YN,XK,K0,INIAPR,IMAX,IU,EPS,THRESH,MEXACT,
     1           M,MESTER,H,HMIN,K20,IMAX2,NATTEM,Y,RAB,IERR)
      PRINT 50,IERR
      PRINT 49,N
      PRINT 52,NACCEP,NREJEC
      PRINT 51,Y,YT,(YT-Y)/Y
      STOP
   49 FORMAT(1X,'N=',I12)
   50 FORMAT(1X,'IERR=',I2)
   51 FORMAT(1X,'Y =',D25.16/1X,'YT=',D25.16/1X,'(YT-Y)/Y=',D25.16)
   52 FORMAT(1X, 'NACCEP=',I4,5X,'NREJEC=',I4)
      END
C
      SUBROUTINE F(X,Y,Z,M)
      DOUBLE PRECISION X,Y,Z,Q
      COMMON/BLOCK/N,Q
      Z=Q*Y
      N=N+1
      RETURN
      END


Результаты:

----------------------------------------------------------------------

IERR= 0
N=        3363
NACCEP=   5     NREJEC=   0
Y =    .7896296018268813D+14
YT=    .7896296018268069D+14
(YT-Y)/Y=   -.9418973127137907D-13

-----------------------------------------------------------------------

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