|
Текст подпрограммы и версий ( Фортран ) de06d.zip |
Тексты тестовых примеров ( Фортран ) tde06d1.zip , tde06d2.zip , tde06d3.zip , tde06d4.zip |
|
Текст подпрограммы и версий ( Паскаль ) de06e_p.zip |
Тексты тестовых примеров ( Паскаль ) tde06e1_p.zip , tde06e2_p.zip , tde06e3_p.zip , tde06e4_p.zip |
Подпрограмма DE06D является комплексом программ, предназначенным дл выполнения одного шага приближенного интегрирования системы обыкновенных дифференциальных уравнений первого порядка методом рядов Чебышёва с контролем точности или без контроля точности.
Решается задача Коши для системы M обыкновенных дифференциальных уравнений первого порядка
(1) Y ' = F (x, Y) ,
Y = (y1, ... , yM) , F = (f1 (x, y1, ... , yM), ... , fM (x, y1, ... , yM)) , XN ≤ x ≤ XK
с начальными условиями, заданными в точке XN :
(2) Y(XN) = YN , YN = (y10, ... , yM0) ,
при условии, что правая часть системы (1) имеет непрерывные ограниченные частные производные по переменным x, Y. Предполагается, что на отрезке [XN, XK] задача (1), (2) имеет единственное решение. Тогда решение задачи Коши и его производная Y'(XN + αΔ) = F(XN + αΔ, Y(XN + αΔ)) = Φ(α) разлагаются на промежутке интегрирования [XN, XK] в равномерно сходящиеся ряды по смещенным многочленам Чебышёва первого рода
∞ 1
(3) Y(XN + αΔ) = ∑'ai*[Y]Ti*(α), 0 ≤ α ≤ 1, ai*[Y] = 2/π ∫ Y(XN + αΔ)Ti*(α)/√(α(1 - α))dα,
i=0 0
∞ 1
(4) Φ(α) = ∑'ai*[Φ]Ti*(α), 0 ≤ α ≤ 1, ai*[Φ] = 2/π ∫ Φ(α)Ti*(α)/√(α(1 - α))dα.
i=0 0
Здесь: штрих у знака суммы означает, что слагаемое с индексом 0 берется с дополнительным множителем 1/2; Ti* (α) - смещенный многочлен Чебышёва первого рода на [0, 1]: Ti* (α) = Ti (2α-1) ; Ti (t) - многочлен Чебышёва первого рода i-го порядка на [-1, 1]; Δ = XK - XN. Если ряд Чебышёва (3) (и ряд (4)) является быстросходящимся, то его сумма на [XN, XK] (и сумма ряда (4)) хорошо приближается частичной суммой некоторого порядка. Эта частичная сумма принимается в качестве приближенного аналитического решения задачи (1), (2) на промежутке [XN, XK]. В противном случае, т.е. при медленной сходимости ряда (3) на интервале [XN, XK], получение аналитического решения в виде одной частичной суммы на всем отрезке интегрирования [XN, XK] может быть затруднительным. Поэтому целесообразно использовать разбиение промежутка интегрирования [XN, XK] на такие элементарные сегменты, вообще говоря, переменной длины H: [xs, xs+H], x0 = XN, s = 0, 1, ... , на каждом из которых ряды Чебышёва для решения Y(x) и его производной Y'(x) будут сходиться значительно быстрее. На каждом подобном сегменте решение исходной задачи Коши приближенно представляется в виде (K + 1) - й частичной суммы смещенного ряда Чебышёва
K+1 1
(5) Y(xs + αH) ≈ UK+1(xs+αH) = ∑'ai*[UK+1]Ti*(α), 0 ≤ α ≤ 1, ai*[UK+1] = 2/π ∫ UK+1(xs + αH)Ti*(α)/√(α(1 - α))dα,
i=0 0
а его производная - в виде частичной суммы K-го порядка
K 1
(6) Y'(xs + αH) ≈ U'K+1(xs+αH) = ∑'ai*[U'K+1]Ti*(α), 0 ≤ α ≤ 1, ai*[U'K+1] = 2/π ∫ U'K+1(xs + αH)Ti*(α)/√(α(1 - α))dα.
i=0 0
В этом случае аналитическое решение задачи (1), (2) будет состоять из совокупности частичных сумм UK+1(xs+αH) рядов Чебышёва, построенных на этих элементарных сегментах. Порядок частичных сумм задается пользователем при обращении к подпрограмме.
При обращении к комплексу программ (подпрограмме) DE06D задаются начало элементарного сегмента и его длина. По заданному значению решения в начале X = xs элементарного сегмента [X, X + H] подпрограмма DE06D вычисляет значение решения в конце элементарного сегмента, т.е. в узле X + H = xs + H. Одновременно вычисляются коэффициенты Чебышёва ai*[UK+1] (i = 0, 1, ... , K + 1) на элементарном сегменте [X, X + H] = [xs, xs + H] для решения задачи Коши Y(X + αH) ≈ UK+1(X + αH), 0 ≤ α ≤ 1, и коэффициенты Чебышёва ai*[U'K+1] (i = 0, 1, ... , K) его производной U'K+1(X + αH). Причем длина сегмента H здесь может быть меньше или равна заданной при обращении длине. При многократном обращении к подпрограмме задаваемое значение H может изменяться от сегмента к сегменту, в общем случае H = Hs.
Приводимое здесь описание метода решения преследует две цели. Во первых, дать пользователю общее представление об используемом подходе к интегрированию обыкновенных дифференциальных уравнений. Во-вторых, на основе данного знания помочь пользователю яснее усвоить смысл и назначение формальных параметров. С более подробным изложением различных аспектов математического метода решения можно ознакомиться по приводимому в конце данного раздела списку литературы.
При разбиении промежутка интегрирования на элементарные сегменты решение задачи на [XN, XK] сводится к определению нескольких наборов коэффициентов ai*[UK+1], i = 0, 1, ... , K + 1. Коэффициенты ai*[UK+1] ряда Чебышёва для решения на сегменте [X, X + H] = [xs, xs + H] выражаются через коэффициенты ai*[Φ] ряда Чебышёва его производной Φ(α) = F(xs + αH, UK+1(xs + αH)), 0 ≤ α ≤ 1, на [X, X + H] = [xs, xs + H], которые, в свою очередь, вычисляются приближенно итерационным способом, исходя из некоторого начального приближения. Вычисления выполняются с помощью квадратурной формулы Маркова на [X, X + H] с (K + 1) узлом. При этом один из узлов квадратурной формулы совпадает с X, а остальные K узлов лежат внутри интервала (X, X + H). Количество итераций, которое предписывается выполнить в этом итерационном процессе, задается при обращении к подпрограмме и может меняться от сегмента к сегменту. Если при выбранном H ряды Чебышёва для решения UK+1(X + αH), 0 ≤ α ≤ 1 , и его производной U'K+1(X + αH) на элементарном сегменте [X, X + H] = [xs, xs + H] быстро сходятся, то для того, чтобы приближенное решение в конце такого сегмента имело максимальный порядок точности относительно H, необходимо выполнить не менее K итераций; при этом погрешность приближенного решения в конце элементарного сегмента является величиной порядка O(HK + 2) при H --> 0. Если H подобрано достаточно малым, то хорошая точность приближенного решения может быть получена и при меньшем числе итераций.
Начальное приближение коэффициентов ai*[Φ] ряда Чебышёва для производной на сегменте [X, X + H] = [xs, xs + H] может вычисляться двумя способами. В первом способе начальное приближение определяется только с использованием значения решения в точке X = xs. При этом погрешность начального приближения для всех коэффициентов a0*[Φ], a1*[Φ], ..., aK*[Φ] является величиной O(H2) при H --> 0. Во втором способе начальное приближение определяется через коэффициенты ряда Чебышёва производной решения Φ(α) на предыдущем элементарном сегменте [X - H', X] = [xs - 1, xs], где H' - длина предыдущего элементарного сегмента. В этом случае погрешности начального приближения для коэффициентов a0*[Φ], a1*[Φ], ... , aK*[Φ] имеют, соответственно, порядки O(H), O(H2), ..., O(HK + 1). Заметим, что длина H текущего элементарного сегмента [X, X + H] может быть больше или меньше длины H' предыдущего элементарного сегмента [X - H', X] или равна ей. Второй способ определения начального приближения в некоторых случаях может привести к более быстрой сходимости итерационного процесса и, тем самым, к меньшему числу выполняемых итераций. Второй способ может быть применен только начиная со второго элементарного сегмента [x0 + H0, x1 + H], H0 -длина начального (первого) элементарного сегмента, x1 = x0 + H0. На начальном элементарном сегменте [x0, x0 + H0] всегда применяется исключительно первый способ. Способ выбора начального приближения задается пользователем при обращении к подпрограмме и может меняться от сегмента к сегменту..
Если шаг интегрирования нужно выполнять с контролем точности, то на заданном при обращении к комплексу программ (подпрограмме) DE06D частичном сегменте [X, X + H] = [xs, xs + H] для оценки погрешности решения UK+1(X + αH), полученного при некотором заданном значении K = K1, вычисляется на этом же сегменте специальным образом второе приближенное решение при K = K2 > K1, которое вместе со своей производной также представляется указанными выше частичными суммами (5), (6) при K = K2. Второе решение UK2+1, как имеющее более высокий порядок точности O(HK2+2), используется для оценки погрешности первого решения UK1+1 на элементарном сегменте [X, X + H]. Значения K1 и K2 задаются пользователем при обращении к подпрограмме и могут изменяться от сегмента к сегменту. Второе приближенное решение на сегменте [X, X + H] вычисляется аналогично первому решению итерационным способом. Единственным отличием при этом является то, что начальным приближением в этом дополнительном итерационном процессе служит уже вычисленное на отрезке [X, X + H] первое приближенное решение UK+1(X + αH). Количество итераций, которое потребуется выполнить в этом дополнительном итерационном процессе, задается при обращении к подпрограмме и может меняться от сегмента к сегменту.
Для системы уравнений (M > 1) с проверкой на точность могут вычисляться либо все компоненты решения, либо некоторые из них (в частности, одна компонента). В последнем случае номера проверяемых на точность компонент задаются при обращении к подпрограмме.
Если шаг интегрирования нужно выполнять с контролем точности, то для оценки уклонения приближенного решения от точного (т.е. для оценки абсолютной погрешности приближенного решения) в подпрограмме используются два способа, или две формулы. Первый способ опирается на первую формулу для абсолютной погрешности приближенного решения. Она является асимптотической (т.е. справедливой при H --> 0) и представляет собою разность двух частичных сумм смещенного ряда Чебышёва на сегменте [X, X + H] с приближенными коэффициентами, а именно разность приближенных частичных сумм (K2 + 1) - го и (K1 + 1) - го порядков в точке X + H. Второй способ опирается на несколько завышенную оценку абсолютной погрешности, которая представляет сумму абсолютных величин (модулей) разностей приближенных коэффициентов этих частичных сумм и поэтому несколько завышена по сравнению с первой. Назовем эту формулу второй. Более того, при использовании этой, второй, формулы в том случае, когда точность решения оценивается по относительной погрешности или по мере погрешности, то дополнительно используется также завышенная оценка сверху для относительной погрешности и завышенная оценка сверху для меры погрешности. Второй способ оценки погрешности, хотя и является более точным, накладывает более сильное ограничение на размер элементарного сегмента и приводит к меньшей длине сегментов. Используемый способ оценки абсолютной погрешности задается при обращении к подпрограмме и может меняться от сегмента к сегменту.
Когда интегрирование на шаге требуется вести с контролем точности, то используется следующий алгоритм достижения заданной точности. Если для заданной при обращении к подпрограмме DE06D длины H сегмента [X, X + H] точность приближенного решения UK1+1 не достигается, то данный сегмент исключается из рассмотрения и программа сокращает его длину H до тех пор, пока на новом, сокращенном, сегменте не будет достигнута заданная точность вновь вычисленного приближенного решения UK1+1. При выполнении условия достижения на сегменте заданной точности приближенного решения UK1+1 этот сегмент (заданной, первоначальной, длины или сокращенной длины) принимается в качестве элементарного (частичного) сегмента, на который продлено решение дифференциального уравнения. В качестве такого решения на сегменте [X, X + H] принимается частичная сумма порядка K1 + 1 с коэффициентами, равными коэффициентам Чебышёва ai*[UK2+1], i = 0, 1, ... , K1+1, второго приближенного решения UK2+1, поскольку эти коэффициенты имеют более высокий порядок точности по сравнению с коэффициентами Чебышёва ai*[UK1+1] первого приближенного решения. В качестве значения решения в конце X + H сегмента [X, X + H] принимается значение второго приближенного решения в точке X + H, т.е. UK2+1(X + H), как имеющее более высокий порядок точности O(HK2+2) по сравнению с первым приближенным решением UK1+1(X + H). В качестве коэффициентов Чебышёва производной решения на сегменте [X, X + H] принимаются по той же причине коэффициенты Чебышёва производной второго приближенного решения ai*[U'K2+1], i = 0, 1, ... , K1.
Таким образом, один шаг приближенного интегрирования по методу рядов Чебышёва, выполняемый комплексом программ (подпрограммой) DE06D, заключается в следующем. По заданному значению решения YX в начале X сегмента [X, X + H] вычисляется (при необходимости, то с контролем точности) значение этого решения в конце X + H сегмента, а также определяются значения коэффициентов Чебышёва решения и значения коэффициентов Чебышёва его производной на сегменте [X, X + H]; другими словами, находится ортогональное разложение решения и ортогональное разложение его производной на сегменте [X, X + H].
В дальнейшем при описании параметров подпрограммы коэффициенты ряда Чебышёва будем называть коэффициентами Чебышёва.
О.Б.Арушанян, С.Ф.Залеткин. Приближенное решение задачи Коши для обыкновенных дифференциальных уравнений методом рядов Чебышёва. Вычислительные методы и программирование: Новые вычислительные технологии (Электронный научный журнал) (17), 121 - 131, 2016.
О.Б.Арушанян, Н.И.Волченскова, С.Ф.Залеткин. Вычисление коэффициентов разложения решения задачи Коши в ряд по многочленам Чебышёва. Вестник Московского университета. Серия 1: Математика. Механика. 5, 24 - 30, 2012.
О.Б.Арушанян, С.Ф.Залеткин. Об оценке погрешности приближенного решения обыкновенных дифференциальных уравнений, определенного с помощью рядов Чебышёва // Вычислительные методы и программирование. 21, 241 - 250, 2020.
О.Б.Арушанян, С.Ф.Залеткин. О вычислении приближенного решения обыкновенных дифференциальных уравнений методом рядов Чебышёва и оценка его погрешности // Вестн. Моск. ун-та. Сер. 1, Математика. Механика. 2020, № 5, 22 - 26.
О.Б.Арушанян, С.Ф.Залеткин. Использование рядов Чебышёва для приближенного аналитического решения обыкновенных дифференциальных уравнений. Вестник Московского университета. Серия 1: Математика. Механика. 5, 52 - 56, 2016.
SUBROUTINE DE06D (F, M, K, INIAPR, IMAX, JSTART, YX, X, AU, AJK, H, XP, MEXACT, IU, EPS, THRESH, NUMBES,
MESTER, K2, IMAX2, NATTEM, HMIN, BUL, RAB, IERR)
Параметры
| F - |
имя подпрограммы вычисления значений правой
части дифференциального уравнения. Первый
оператор подпрограммы должен иметь вид: SUBROUTINE F (X, Y , DY, M). Здесь: X, Y - значения независимой и зависимой переменных, соответственно. Вычисленное значение правой части должно быть помещено в DY. B случае системы уравнений, т.е. когда M ≠ 1 , параметры Y и DY представляют массивы длины M (тип параметров X, Y и DY: с двойной точностью); |
| M - | количество уравнений в системе (тип: целый); |
| K - | порядок частичной суммы смещенного ряда Чебышёва, с помощью которой аппроксимируется производная первого приближенного решения задачи Коши на элементарном сегменте [X, X + H] разбиения интервала интегрирования; при этом само решение задачи Коши приближается на элементарном сегменте [X, X + H] частичной суммой (K + 1) - го порядка; K ≥ 2. Этот параметр относится только к первому приближенному решению, имеющему меньший порядок точности O(HK+2), и его производной (см. "Математическое описание", "Замечания по использованию" и "Примеры"; тип: целый); |
| INIAPR - |
целый указатель способа выбора начального
приближения коэффициентов Чебышёва для производной первого приближенного
решения на элементарном сегменте [X, X + H]: |
| INIAPR=1 - |
для первого способа, когда начальное
приближение определяется только с
использованием значения решения в
начале X элементарного сегмента; |
| INIAPR=2 - | для второго способа, когда начальное приближение коэффициентов Чебышёва на текущем элементарном сегменте [X, X + H] (начиная со второго) определяется через коэффициенты Чебышёва, вычисленные на предыдущем элементарном сегменте, т.е. путем экстраполяции коэффициентов с предыдущего сегмента на следующий (см. "Математическое описание"). |
| Значение этого параметра может меняться от сегмента к сегменту. | |
| IMAX - | целая переменная, задающая количество итераций, которое предполагается выполнить в итерационном процессе вычисления коэффициентов Чебышёва для производной первого приближенного решения задачи Коши на элементарном сегменте [X, X + H], исходя из некоторого начального приближения, способ определения которого задается параметром INIAPR; IMAX ≥ 1. Для получения максимального порядка точности первого приближенного решения необходимо выполнить не менее K итераций. Значение IMAX может изменяться от сегмента к сегменту. Этот параметр относится только к первому приближенному решению, имеющему меньший порядок точности O(HK+2), и его производной. Если правая часть дифференциального уравнения не зависит от переменной Y, т.е. дифференциальное уравнение имеет вид Y' = F(X), то число итераций при вычислении первого приближенного решения можно положить равным 1. В этом случае параметр IMAX = 1 (см. "Математическое описание", "Замечания по использованию" и "Примеры"); |
| JSTART - |
целый указатель режима использования подпрограммы, имеющий следующие значения: |
| 0 - |
первое обращение к подпрограмме должно быть исполнено с нулевым значением JSTART:
выполнить первый (начальный) шаг интегрирования для значений независимой и зависимой переменных и
длины элементарного сегмента, заданных параметрами X, YX, H,
соответственно. При данном значении параметра JSTART будет применен исключительно первый способ определения
начального приближения для коэффициентов Чебышёва производной первого приближенного решения
независимо от значения параметра INIAPR.
Нулевое значение параметра JSTART также может означать, что необходимо выполнить очередной шаг
интегрирования с измененным значением параметра K (см. выше) или с измененным значением
параметра K2 (см. ниже); |
| 1 - | выполнить следующий (очередной) шаг интегрирования системы дифференциальных уравнений для значений независимой и зависимой переменных и длины элементарного сегмента, заданных параметрами X, YX, H, соответственно. При данном значении параметра JSTART способ определения начального приближения для коэффициентов Чебышёва производной первого приближенного решения определяется параметром INIAPR. При обращении к подпрограмме DE06D со значением JSTART = 1 значения параметра K (см. выше) и параметра K2 (см. ниже) должны совпадать с их значениями при предыдущем обращении к подпрограмме. |
| На выходе из подпрограммы DE06D параметр JSTART всегда принимает значение, равное 1 (см. "Примеры"); | |
| X, YX - | начальное значение аргумента и решения (начало элементарного сегмента X и решение в нем YX); в результате работы подпрограммы в X получается новое значение аргумента, отстоящее от начального значения на величину выбранного подпрограммой элементарного сегмента, т.е. такого сегмента, на который продлено решение дифференциального уравнения (конец элементарного сегмента), а в YX - соответствующее значение решения; если шаг интегрирования выполняется с контролем точности, то в качестве такого решения принимается второе приближенное, оценивающее, решение, поскольку оно имеет более высокий порядок точности по сравнению с первым приближенным решением. В случае системы уравнений, т.е. когда M ≠ 1, YX задается одномерным массивом длины M (тип параметров X, YX: с двойной точностью); |
| AU - | двумерный массив размера M * (K + 2). На выходе из подпрограммы содержит коэффициенты Чебышёва ai*[Y] для решения Y (X + αH), 0 ≤ α ≤ 1, на элементарном сегменте [X, X + H]. При этом переменная с индексом AU (N, I + 1) представляет I-й коэффициент Чебышёва N-й компоненты решения yN(x)(I = 0, 1, ... , K + 1), (N = 1, ... , M) (тип: с удвоенной точностью); |
| AJK - | двумерный массив размера M * (K + 1). На выходе из подпрограммы содержит коэффициенты Чебышёва ai*[Y'] производной решения Y'(X + αH), 0 ≤ α ≤ 1, на элементарном сегменте [X, X + H]. При этом переменная с индексом AJK(N, I + 1) представляет I-й коэффициент Чебышёва N-й компоненты производной решения y'N(x)(I = 0, 1, ... , K). Если при JSTART = 1 (т.е. при повторных обращениях к подпрограмме) параметр INIAPR = 2, то содержащиеся в массиве AJK при входе в подпрограмму коэффициенты Чебышёва производной решения, вычисленные на предыдущем элементарном сегменте [X - H', X] во время предыдущего обращения к подпрограмме DE06D, используются в подпрограмме при вычислении коэффициентов Чебышёва производной решения на текущем элементарном сегменте [X, X + H]. Заметим, что длина H текущего элементарного сегмента может быть больше или меньше длины H' предыдущего сегмента [X - H', X] или равна ей. |
| H - | переменная, содержащая задаваемое при обращении к подпрограмме значение длины текущего элементарного сегмента (тип: с двойной точностью); |
| XP - | переменная, содержащая на выходе из подпрограммы начальное значение аргумента, т.е. то значение, которое имел параметр X на входе в подпрограмму. Таким образом, на выходе из подпрограммы DE06D границы элементарного сегмента, к которому относятся вычисленные подпрограммой коэффициенты Чебышёва решения и коэффициенты Чебышёва его производной и содержащиеся в массивах AU и AJK соответственно, показываются параметрами XP и X, а именно: параметр XP содержит начало сегмента, а параметр X содержит конец этого сегмента (тип: с двойной точностью); |
| MEXACT - | количество компонент численного решения, проверяемых на точность: 1 ≤ MEXACT ≤ M. Если решение вычисляется без контроля точности, то значение MEXACT полагается равным нулю. При этом значения нижеследующих параметров IU, EPS, THRESH, NUMBES, MESTER, K2, IMAX2, NATTEM, HMIN, BUL, IERR в комплексе (подпрограмме) DE06D не используются. В этом случае конкретные значения соответствующих им фактических параметров при обращении к DE06D можно не задавать, однако их типы, количество и порядок записи должны соответствовать указанным формальным параметрам (тип: целый); |
| IU - | целый указатель типа погрешности численного решения: |
| IU=1 - | компоненты решения проверяются на точность по абсолютной погрешности; |
| IU=2 - | компоненты решения проверяются на точность по относительной погрешности; |
| IU=3 - | компоненты решения проверяются на точность по мере погрешности. |
| Контроль точности по мере погрешности заключается в следующем. Если некоторая компонента приближенного решения по абсолютной величине не меньше определенной наперед заданной положительной константы THRESH (называемой границей перехода), то контроль точности для этой компоненты ведется по относительной погрешности, иначе - по абсолютной. Таким образом, контроль точности по мере погрешности состоит в том, что на тех участках интервала интегрирования, где абсолютная величина компоненты решения меньше значения THRESH, контроль точности ведется для нее по абсолютной погрешности, а там, где абсолютная величина этой компоненты решения равна значению THRESH или првосходит значение THRESH, контроль точности для нее ведется по относительной погрешности; | |
| EPS - | допустимая погрешность, с которой требуется вычислить проверяемые на точность компоненты решения; тип погрешности специфицируется с помощью параметра IU (тип: с двойной точностью); |
| THRESH - | граница перехода, используемая при оценке меры погрешности решения (тип: с двойной точностью); |
| NUMBES - | массив, содержащий номера проверяемых на точность компонент решения. Для скалярного дифференциального уравнения, а также для системы уравнений в том случае, если все компоненты решения проверяются на точность, значение этого массива в подпрограмме не используется. Таким образом, массив NUMBES используется в подпрограмме только в том случае, когда число проверяемых на точность компонент решения меньше числа уравнений; NUMBES представляет массив длины MEXACT (тип: целый); |
| MESTER - | целый указатель способа оценки абсолютной погрешности компоненты приближенного решения: |
| MESTER=1 - | для оценки абсолютной погрешности используется первая формула (асимптотическая); |
| MESTER=2 - | для оценки абсолютной погрешности используется вторая формула (завышенная оценка). |
| (См. "Математическое описание"); | |
| K2 - | порядок частичной суммы смещенного ряда Чебышёва производной решения, с помощью которой вычисляется частичная сумма порядка (K2 + 1) для оценивающего решения более высокого порядка точности на элементарном сегменте [X, X + H] разбиения области интегрирования. Этот параметр относится только ко второму приближенному решению более высокого порядка точности O(HK2+2), K2 > K (т.е. к оценивающему решению) и его производной. Если обращение к подпрограмме DE06D осуществляется со значением параметра JSTART = 1 (см. выше), то значение параметра K2 должно совпадать с его значением при предыдущем обращении к подпрограмме. Если же при очередном обращении к подпрограмме DE06D значение параметра K2 необходимо изменить, то это обращение должно выполняться при нулевом значении параметра JSTART (JSTART = 0) (см. "Математическое описание", "Замечания по использованию" и "Примеры"; тип: целый); |
| IMAX2 - | целая переменная, задающая количество итераций, которое предполагается выполнить в дополнительном итерационном процессе для вычисления второго, оценивающего, приближенного решения более высокого порядка точности. Для получения максимального порядка точности оценивающего решения необходимо выполнить не менее (K2 - K) итераций. Значение IMAX2 может изменяться от сегмента к сегменту. Если правая часть дифференциального уравнения не зависит от переменной Y, т.е. дифференциальное уравнение имеет вид Y' = F(X), то число итераций при вычислении второго приближенного решения можно положить равным 1. В этом случае параметр IMAX2 = 1 (см. "Математическое описание" и "Замечания по использованию"); |
| NATTEM - | целая переменная, значение которой ограничивает число последовательных сокращений длины элементарного сегмента [X, X + H], если приближенное решение на этом элементарном сегменте не достигает требуемой точности; |
| HMIN - | минимальное значение длины элементарного сегмента, которое разрешается использовать при разбиении интервала интегрирования на элементарные (частичные) сегменты (тип: с двойной точностью); |
| BUL - | логическая переменная, значение которой при обращении к подпрограмме полагается равным .TRUE., если заданная в H длина сегмента выводит в конец интервала интегрирования, и .FALSE. в противном случае; в результате работы подпрограммы BUL равно .FALSE., если вместо исходного сегмента был выбран меньший сегмент, в противном случае, т.е. когда был выбран именно заданный при обращении в H сегмент, значение параметра BUL не меняется; |
| RAB - | одномерный рабочий массив. Длина массива равна 2 * K * K + 5 * K + 3 * K* M + K2 * K2 + 5 * K2 + 4 * K2 * M + 15 * M, если шаг интегрирования выполняется с контролем точности, и 2 * K * K + 7 * K + 3 * K * M + 5 * M, если шаг интегрирования выполняется без контроля точности (тип: с двойной точностью); |
| IEER - | целая переменная, служащая для сообщения об ошибках, обнаруженных в процессе работы подпрограммы; при этом: |
|
IERR=65 - IERR=66 | когда какая-нибудь компонента решения не может быть вычислена с требуемой точностью EPS; при этом IERR=65 указывает, что требуемая точность не может быть достигнута на сегменте [X, X + H] при минимальной длине частичного сегмента H, равной HMIN, а IERR=66 показывает, что требуемая точность не может быть достигнута на сегменте [X, X + H], т.к. исчерпано заданное число сокращений NATTEM длины H элементарного сегмента. |
Версии: нет
Вызываемые подпрограммы: DE75D, DE44D
|
DE06D использует рабочие подпрограммы DE70DH, DE70DI, DE70DF, DE70DQ, DE71DE, DE70DP, DE71DT, DE71DP, DE71DI, DE71DF, DE71DS, DE70DA, DE70DC, DE75D0, DE75D1, DE75DE, DE75DK, DE75DN. UTDE75 - подпрограмма выдачи диагностических сообщений. |
Замечания по использованию
|
В ообщем случае заданная точность не гарантируется. Разбиение промежутка интегрирования на элементарные сегменты [X, X + H] выполняется для того, чтобы на каждом таком сегменте ряды Чебышёва для решения и его производной были быстросходящимися рядами; другими словами, чтобы убывание коэффициентов этих рядов Чебышёва на элементарном сегменте [X, X + H] происходило достаточно быстро, вследствие чего можно было бы считать частичные суммы этих рядов близкими к многочленам наилучшего равномерного приближения на элементарном сегменте [X, X + H] для решения и его производной. Порядок этих частичных сумм задается параметром K и параметром K2. Если при вычислении первого приближенного решения начальное приближение для коэффициентов Чебышёва правой части системы (т.е. функции Φ(α)) определяется первым способом (т.е. при INIAPR = 1), то для получения максимального порядка точности первого приближенного решения в конце элементарного сегмента [X, X + H] необходимо выполнить в итерационном процессе не менее K итераций; тогда IMAX≥K. Если начальное приближение коэффициентов Чебышёва функции Φ(α) определяется вторым способом (т.е. при INIAPR = 2), то для получения максимального порядка точности первого приближенного решения необходимо выполнить в итерационном процессе не менее K + 1 итераций; в этом случае IMAX ≥ K + 1. Однако в некоторых случаях при втором способе определения начального приближения итерационный процесс может сойтись за значительно меньшее число итераций. Если длина H элементарного сегмента [X, X + H] подобрана достаточно малой, то хорошая точность первого приближенного решения может быть получена и с существенно меньшим числом итераций при любом способе выбора начального приближения. Для получения максимального порядка точности второго приближенного (оценивающего) решения необходимо выполнить в дополнительном итерационном процессе не менее (K2 - K) итераций. Указанное в данном разделе число итераций IMAX и IMAX2 носит асимптотический характер (когда все наши оценки справедливы при H --> 0). На практике же количество итераций в каждом из двух итерационных процессов зависит от длины сегмента, от требуемой точности, от заданных порядков частичных сумм и от конкретной системы дифференциальных уравнений. Поэтому число итераций в обоих итерационных процессах, т.е. значения параметров IMAX и IMAX2, может быть как меньше, так и больше рекомендуемых здесь значений. Если правая часть дифференциального уравнения не зависит от переменной Y, т.е. дифференциальное уравнение имеет вид Y'= F(X), то число итераций в обоих итерационных процессах можно положить равным 1. В этом случае параметры IMAX = 1 и IMAX2 = 1. При работе подпрограммы значения параметров M, K, INIAPR, IMAX, IU, EPS, THRESH, MEXACT, NUMBES, MESTER, HMIN, K2, IMAX2, NATTEM сохраняются. При многократном использовании подпрограммы DE06D со значением параметра JSTART = 1 для вычисления коэффициентов Чебышёва решения задачи Коши (1), (2) и его производной на последовательности элементарных сегментов, образующей промежуток интегрирования [XN, XK] системы дифференциальных уравнений, значения параметров M, K, YX, X, AU, AJK, K2, RAB не должны изменяться в вызывающей программе между последовательными обращениями к подпрограмме. Подпрограмма DE06D использует общие блоки с именами COM70D, COM75D, COM06D. Поэтому, если вызывающая (главная) программа содержит хотя бы одно обращение к подпрограмме DE06D со значением параметра JSTART = 1 (и, следовательно, число обращений к DE06D более одного), то данные общие блоки должны быть описаны в вызывающей (главной) программе в операторе COMMON, например, таким способом: COMMON / COM70D / WC1, WC2, WC3, LASN, HD2, HD4, HOLD COMMON / COM75D / LASN0, LASN2, WC30, WC12, WC22, WC32 COMMON / COM06D / I(22) Глобальные (общие) переменные, составляющие содержание общих блоков, описываются как переменные с двойной точностью: DOUBLE PRECISION WC1, WC2, WC3, HD2, HD4, HOLD, WC30, WC12, WC22, WC32 а массив I и переменные LASN, LASN0, LASN2 - как целый массив и как целые переменные (явно или неявно): INTEGER LASN, LASN0, LASN2, I(22) Пользователю не рекомендуется применять для своих целей общие блоки с указанными именами (пользователь не должен изменять содержимое указанных переменных из данных блоков). Если вызывающая программа содержит одно обращение к подпрограмме DE06D со значением параметра JSTART = 0 и при этом других обращений к DE06D нет или если вызывающая программа содержит несколько обращений к DE06D и при этом все обращения осуществляются с нулевым значением параметра JSTART, то указанные общие блоки описывать в вызывающей программе вовсе не обязательно (см. "Примеры использования"). При работе подпрограммы DE06D используется общий блок COMMON/STAT76/NACCEP, NREJEC для сбора полезной статистической информации, если шаг интегрирования выполняется с контролем точности. Такой же блок должен содержаться в вызывающей (главной) программе. Если подпрограмма DE06D заканчивает свою работу со значением IERR = 0, то на выходе из подпрограммы переменная NACCEP содержит значение, увеличенное на 1 по сравнению со значением этой переменной на входе в подпрограмму. Если же подпрограмма DE06D заканчивает свою работу со значением IERR = 65 или IERR = 66, то значение переменной NACCEP остается без изменения. Переменная NREJEC учитывает, сколько раз сокращалась длина сегмента [X, X + H] при вычислении решения с требуемой точностью. Если заданная точность была достигнута при первоначальной длине сегмента (задаваемой параметром H при входе в подпрограмму DE06D), то значение переменной NREJEC на выходе из подпрограммы равно значению этой переменной на входе в подпрограмму; если же первоначально заданная длина H подвергалась сокращению, то значение переменной NREJEC (которое эта переменная имела на входе в подпрограмму) увеличивается на столько, сколько раз сокращалась длина сегмента. Таким образом, переменная NACCEP подсчитывает число выбранных (или принятых) подпрограммой DE06D сегментов, а переменная NREJEC подсчитывает число отклоненных сегментов при решении дифференциальных уравнений на всем промежутке интегрирования (см. примеры 3 и 4). Если шаг интегрирования выполняется без контроля точности, то значения переменных NACCEP, NREJEC из общего блока STAT76 при работе DE06D не используются и описывать общий блок STAT76 в вызывающей программе не обязательно (см. примеры 1 и 2). |
Данные примеры не только иллюстрируют правила использования комплекса (подпрограммы) DE06D; они также показывают, как можно применять подпрограмму DE06D для решения обыкновенных дифференциальных уравнений на заданном промежутке интегрирования. Вычисления во всех примерах проводились на Фортране с 15-16 значащими цифрами. В первых двух примерах решаются задачи Коши для системы обыкновенных дифференциальных уравнений.
(7) y'1 = -2q (y2-1) + (1 -e1-y1+cos q(2x-1) ) / (x+1),
y'2 = 2q (y1-1) + ( 1 -e1-y2+sin q(2x-1) ) / (x+1), q = 0,5
без контроля точности.
Частное решение системы (7) имеет вид
y1(x) = 1 + cos q(2x - 1), y2(x) = 1 + sin q(2x - 1).
Функции y1(x), y2(x) разлагаются на отрезке [0, 1] в смещенные ряды Чебышёва, коэффициенты которых выражаются через цилиндрические функции:
∞
(8) y1(x) = 1 + J0 (q) + 2 ∑ (-1)i J2i (q) T2i* (x),
i=1
∞
(9) y2(x) = 1 + 2 ∑ (-1)i J2i+1 (q) T2i+1* (x).
i=0
Задачи Коши отличаются начальными условиями.
1) Начальное условие задачи Коши для системы (7)
XN = 0, y1(0) = 1 + cos q, y2(0) = 1 - sin q.
Рассматривается один элементарный сегмент [0, 1]. Выполняется одно обращение к подпрограмме DE06D из начальной точки X = 0 с шагом H = 1 без контроля точности. Вычисляются решение YX в конце данного сегмента, т.е. при x = 1, коэффициенты Чебышёва AU решения на данном сегменте и коэффициенты Чебышёва AJK его производной. Приводятся вызывающая программа, подпрограмма вычисления правой части системы (7), результаты счета, включая точное значение решения YT в точке x = 1, абсолютную погрешность DELY приближенного решения YX, вычисленного в точке x = 1. Кроме вышеперечисленного приводятся: значения параметров XP, X на выходе из подпрограммы, которые представляют элементарный сегмент [0, 1]; значения параметров H, K, INIAPR, IMAX, при которых были вычислены приближенное решение YX и коэффициенты Чебышёва. Даются также точные значения коэффициентов Чебышёва на данном сегменте для компонент решения y1(x), y2(x) в (8), (9) и абсолютные погрешности их приближенных значений, вычисленных подпрограммой DE06D.
DOUBLE PRECISION RAB,YT,DELY(2),
1 H,XP,XK,Q,JK(32),EXACT(32),
2 ERROR(42),V0,V1,W1,W2,
3 YX,X,AU,AJK,R
PARAMETER (K=11)
DIMENSION RAB(13),YT(2),
1 YX(2),AU(2,K+2),AJK(2,K+1),
2 R(2*K*K+7*K+3*K*2+5*2)
C__________________________________________________
C
C
LOGICAL BUL
EXTERNAL F
C
M=2
Q=0.5D0
X=0.
YX(1)= COS(Q)+1.D0
YX(2)=-SIN(Q)+1.D0
H=1.0D0
XK=X+H
YT(1)=COS(Q*(2.D0*XK-1.D0))+1.D0
YT(2)=SIN(Q*(2.D0*XK-1.D0))+1.D0
IMAX=16
JSTART=0
INIAPR=1
MEXACT=0
C
CALL DE06D(F,M,K,INIAPR,IMAX,JSTART,YX,X,AU,AJK,H,XP,MEXACT,
1 1,1.D0,1.D0,M,1,1,1,1, 1.D0,BUL,R,IERR)
C
DELY(1)=YT(1)-YX(1)
DELY(2)=YT(2)-YX(2)
PRINT 47
PRINT 51,XP,X
PRINT 35, H,K,INIAPR,IMAX,JSTART
PRINT 60,YX,YT,DELY
NX=1
KP2=K+2
KP1=K+1
PRINT 110,NX
DO 20 L=1,M
PRINT 47
PRINT 150,L
PRINT 85,(J-1,AU(L,J),AJK(L,J),J=1,KP1), KP1,AU(L,KP2)
20 CONTINUE
PRINT 180
C ***** BESSEL FUNCTIONS ******************
CALL SF33D(Q,KP2,KP2,JK,EXACT,RAB,IERR)
C **************
V1=1.D0
V0=Q/2.D0
W1=1.D0
W2=W1/1.D0
C
DO 8 J=1,KP2
JK(J)=JK(J)*V1*W2
V1=V1*V0
W2=W2/DBLE(J)
8 CONTINUE
C **************
EXACT(1)=2.D0*JK(1)+2.D0
C
DO 3 J=2,KP2,2
3 EXACT(J)=0.D0
C
DO 4 J=3,KP2,4
4 EXACT(J)=-2.D0*JK(J)
C
DO 5 J=5,KP2,4
5 EXACT(J)= 2.D0*JK(J)
C
DO 6 J=1,KP2
6 ERROR(J)=EXACT(J)-AU(1,J)
PRINT 160
PRINT 161,(J-1,AU(1,J),EXACT(J),ERROR(J),J=1,KP2)
PRINT 48
C ***** BESSEL FUNCTIONS ******************
CALL SF33D(Q,KP2,KP2,JK,EXACT,RAB,IERR)
C *****
V1=1.D0
V0=Q/2.D0
W1=1.D0
W2=W1/1.D0
C
DO 800 J=1,KP2
JK(J)=JK(J)*V1*W2
V1=V1*V0
W2=W2/DBLE(J)
800 CONTINUE
C **************
EXACT(1)=2.D0
C
DO 300 J=3,KP2,2
300 EXACT(J)=0.D0
C
DO 400 J=4,KP2,4
400 EXACT(J)=-2.D0*JK(J)
C
DO 500 J=2,KP2,4
500 EXACT(J)= 2.D0*JK(J)
C
DO 600 J=1,KP2
600 ERROR(J)=EXACT(J)-AU(2,J)
PRINT 170
PRINT 161,(J-1,AU(2,J),EXACT(J),ERROR(J),J=1,KP2)
STOP
47 FORMAT(1X,64(1H-)/)
48 FORMAT(1X,89(1H-)/)
50 FORMAT(1X,'NX=',I3,4X,'XS='/(3D25.16))
51 FORMAT(1X,'XP=',D25.16,5X,'X=',D25.16)
35 FORMAT(1X,'H=',D25.16,3X,'K=',I2,3X,'INIAPR=',I1,3X,'IMAX=',I2,
1 3X,'JSTART=',I1)
60 FORMAT(1X,'YX=',2D25.16/1X,'YT=',2D25.16/1X,'DELTA Y =',2D25.16)
85 FORMAT(5X,I2,5X,2D25.16)
110 FORMAT(/1X,'COEFFICIENTS AU AND AJK ON ',I3,' SEGMENT ')
C
150 FORMAT(' Number of',7X,'Chebyshev coefficients'
* ' for ',I1,' component'/
1 14X,51(1H-)/
2 ' coefficient',11X,'for Y',21X,'for Y'''/
3 1X,64(1H-)/)
160 FORMAT(//' Number of',19X,'Chebyshev coefficients for 1',
* ' component'/
1 14X,76(1H-)/
2 ' coefficient',6X,'approximate',18X,'exact',16X,
3 'absolute error'/
4 1X,86(1H-)/)
161 FORMAT(5X,I2,5X,3D25.16)
170 FORMAT(' Number of',19X,'Chebyshev coefficients for 2',
* ' component'/
1 14X,76(1H-)/
2 ' coefficient',6X,'approximate',18X,'exact',16X,
3 'absolute error'/
4 1X,86(1H-)/)
180 FORMAT(/1X,86(1H*)/)
END
SUBROUTINE F(X,Y,Z,M)
DOUBLE PRECISION X,Y(2),Z(2),R,Q
Q=0.5D0
R =2.D0*X-1.D0
Z(1)=-2.D0*Q*(Y(2)-1.D0)+(1.D0-DEXP(1.D0-Y(1)+DCOS(Q*R)))/(1.D0+X)
Z(2)= 2.D0*Q*(Y(1)-1.D0)+(1.D0-DEXP(1.D0-Y(2)+DSIN(Q*R)))/(1.D0+X)
RETURN
END
Результаты:
----------------------------------------------------------------
XP= .0000000000000000D+00 X= .1000000000000000D+01
H= .1000000000000000D+01 K=11 INIAPR=1 IMAX=16 JSTART=1
YX= .1877582561890373D+01 .1479425538604203D+01
YT= .1877582561890373D+01 .1479425538604203D+01
DELTA Y = .2220446049250313D-15 .0000000000000000D+00
COEFFICIENTS AU AND AJK ON 1 SEGMENT
----------------------------------------------------------------
Number of Chebyshev coefficients for 1 component
---------------------------------------------------
coefficient for Y for Y'
----------------------------------------------------------------
0 .3876939614481626D+01 -.8381458775717082D-15
1 -.4585666969845331D-16 -.4845369153497486D+00
2 -.6120804691736534D-01 -.6547191987778950D-15
3 -.2775557561562891D-16 .5127459989174128D-02
4 .3214729527285557D-03 -.3216522913903480D-15
5 -.9412738329658138D-17 -.1610725448276222D-04
6 -.6721369257266663D-06 -.1333975247971853D-15
7 -.4137513937923756D-17 .2403173467777737D-07
8 .7516446311980316D-09 -.1754713453532009D-16
9 .2682270999638617D-18 -.2089352055964174D-10
10 -.5226345925825359D-12 -.2720331013401911D-16
11 -.6182570485004343D-18 .1186314365969263D-13
12 .2471488262435964D-15
----------------------------------------------------------------
Number of Chebyshev coefficients for 2 component
---------------------------------------------------
coefficient for Y for Y'
----------------------------------------------------------------
0 .2000000000000000D+01 .1876939614481626D+01
1 .4845369153497479D+00 -.3436311023057026D-15
2 -.8447459582967137D-17 -.6120804691736582D-01
3 -.5127459989174504D-02 -.2760514256419655D-15
4 -.1930811603266178D-16 .3214729527282392D-03
5 .1610725448270591D-04 .3287843088062292D-16
6 .0000000000000000D+00 -.6721369258790109D-06
7 -.2403173465988319D-07 .3287843088062292D-16
8 .9050547041412199D-18 .7516445977185249D-09
9 .2089353515759578D-10 .3916680348103885D-17
10 .6757967666373710D-17 -.5226679549233393D-12
11 -.1187881715734862D-13 -.2664020263068445D-15
12 -.5550042214725927D-17
**************************************************************************************
Number of Chebyshev coefficients for 1 component
----------------------------------------------------------------------------
coefficient approximate exact absolute error
--------------------------------------------------------------------------------------
0 .3876939614481626D+01 .3876939614481626D+01 .0000000000000000D+00
1 -.4585666969845331D-16 .0000000000000000D+00 .4585666969845331D-16
2 -.6120804691736534D-01 -.6120804691736528D-01 .6245004513516506D-16
3 -.2775557561562891D-16 .0000000000000000D+00 .2775557561562891D-16
4 .3214729527285557D-03 .3214729527285751D-03 .1946142899611480D-16
5 -.9412738329658138D-17 .0000000000000000D+00 .9412738329658138D-17
6 -.6721369257266663D-06 -.6721369257237697D-06 .2896535042254487D-17
7 -.4137513937923756D-17 .0000000000000000D+00 .4137513937923756D-17
8 .7516446311980316D-09 .7516446309595219D-09 -.2385096714138463D-18
9 .2682270999638617D-18 .0000000000000000D+00 -.2682270999638617D-18
10 -.5226345925825359D-12 -.5226354721645607D-12 -.8795820247308316D-18
11 -.6182570485004343D-18 .0000000000000000D+00 .6182570485004343D-18
12 .2471488262435964D-15 .2476765118959866D-15 .5276856523901907D-18
-----------------------------------------------------------------------------------------
Number of Chebyshev coefficients for 2 component
----------------------------------------------------------------------------
coefficient approximate exact absolute error
--------------------------------------------------------------------------------------
0 .2000000000000000D+01 .2000000000000000D+01 .0000000000000000D+00
1 .4845369153497479D+00 .4845369153497477D+00 -.1665334536937735D-15
2 -.8447459582967137D-17 .0000000000000000D+00 .8447459582967137D-17
3 -.5127459989174504D-02 -.5127459989174487D-02 .1734723475976807D-16
4 -.1930811603266178D-16 .0000000000000000D+00 .1930811603266178D-16
5 .1610725448270591D-04 .1610725448271494D-04 .9032759349519859D-17
6 .0000000000000000D+00 .0000000000000000D+00 .0000000000000000D+00
7 -.2403173465988319D-07 -.2403173465552604D-07 .4357150715565922D-17
8 .9050547041412199D-18 .0000000000000000D+00 -.9050547041412199D-18
9 .2089353515759578D-10 .2089353517865796D-10 .2106217358883670D-19
10 .6757967666373710D-17 .0000000000000000D+00 -.6757967666373710D-17
11 -.1187881715734862D-13 -.1188370792446492D-13 -.4890767116303012D-17
12 -.5550042214725927D-17 .0000000000000000D+00 .5550042214725927D-17
2) Начальное условие задачи Коши для системы (7)
XN = - 0,5, y1(XN) = 1 + cos q (2 * XN - 1), y2(XN) = 1 + sin q (2 * XN - 1).
Рассматриваются два элементарных сегмента: [- 0,5, 0] и [0, 1]. Выполняются два обращения к подпрограмме DE06D без контроля точности. Первое обращение осуществляется из начальной точки X = - 0,5 с шагом H = 0,5. Второе обращение осуществляется из точки X = 0 с шагом H = 1. Приводятся вызывающая программа и результаты счета. В качестве выходных данных, полученных после каждого обращения к подпрограмме DE06D, приводятся: значения параметров XP, X (которые последовательно представляют границы каждого из двух элементарных сегментов); значения параметров H, K, INIAPR, IMAX, при которых было вычислено приближенное значение решения YX в конце каждого элементарного сегмента; приближенное значение решения YX в конце каждого элементарного сегмента; точное значение решения YT в конце каждого элементарного сегмента и абсолютная погрешность DELY приближенного значения YX. После этого для каждого элементарного сегмента приводятся приближенные значения коэффициентов Чебышёва AU для компонент решения y1(x), y2(x) и приближенные значения коэффициентов Чебышёва AJK их производных y'1(x), y'2(x). Даются также точные значения коэффициентов Чебышёва решения в (8), (9) и абсолютные погрешности приближенных значений AU, вычисленных за два обращения к подпрограмме DE06D.
DOUBLE PRECISION RAB,YT,DELY(2),
1 H,XP,XK,Q,JK(32),EXACT(32),
2 ERROR(42),V0,V1,W1,W2,
3 YX,X,AU,AJK,R
PARAMETER (K=11)
DIMENSION RAB(13),YT(2),
1 YX(2),AU(2,K+2),AJK(2,K+1),
2 R(2*K*K+7*K+3*K*2+5*2)
C__________________________________________________
C
DOUBLE PRECISION WC1,WC2,WC3,HD2,HD4,HOLD, WC30,WC12,WC22,WC32
C
COMMON/COM70D/WC1,WC2,WC3,LASN,HD2,HD4,HOLD
COMMON/COM75D/LASN0,LASN2,WC30,WC12,WC22,WC32
C__________________________________________________
C
COMMON/COM06D/I(22)
C__________________________________________________
C
LOGICAL BUL
EXTERNAL F
C
M=2
Q=0.5D0
X=-0.5D0
YX(1)=COS(Q*(2.D0*X-1.D0))+1.D0
YX(2)=SIN(Q*(2.D0*X-1.D0))+1.D0
H=0.5D0
XK=X+H
YT(1)=COS(Q*(2.D0*XK-1.D0))+1.D0
YT(2)=SIN(Q*(2.D0*XK-1.D0))+1.D0
IMAX=14
JSTART=0
INIAPR=1
MEXACT=0
C
CALL DE06D(F,M,K,INIAPR,IMAX,JSTART,YX,X,AU,AJK,H,XP,MEXACT,
1 1,1.D0,1.D0,M,1,1,1,1, 1.D0,BUL,R,IERR)
C
DELY(1)=YT(1)-YX(1)
DELY(2)=YT(2)-YX(2)
PRINT 44
PRINT 51,XP,X
PRINT 35, H,K,INIAPR,IMAX,JSTART
PRINT 60,YX,YT,DELY
PRINT 47
NX=1
GO TO 15
C_______________________________________
C
1 H=1.D0
XK=X+H
YT(1)=COS(Q*(2.D0*XK-1.D0))+1.D0
YT(2)=SIN(Q*(2.D0*XK-1.D0))+1.D0
IMAX=16
C
CALL DE06D(F,M,K,INIAPR,IMAX,JSTART,YX,X,AU,AJK,H,XP,MEXACT,
1 1,1.D0,1.D0,M,1,1,1,1, 1.D0,BUL,R,IERR)
C
DELY(1)=YT(1)-YX(1)
DELY(2)=YT(2)-YX(2)
PRINT 47
PRINT 45
PRINT 51,XP,X
PRINT 35, H,K,INIAPR,IMAX,JSTART
PRINT 60,YX,YT,DELY
PRINT 47
NX=2
15 KP2=K+2
KP1=K+1
PRINT 110,NX
DO 20 L=1,M
PRINT 47
PRINT 150,L
PRINT 85,(J-1,AU(L,J),AJK(L,J),J=1,KP1), KP1,AU(L,KP2)
20 CONTINUE
GO TO (1,2),NX
2 PRINT 180
C ***** BESSEL FUNCTIONS ******************
CALL SF33D(Q,KP2,KP2,JK,EXACT,RAB,IERR)
C **************
V1=1.D0
V0=Q/2.D0
W1=1.D0
W2=W1/1.D0
C
DO 8 J=1,KP2
JK(J)=JK(J)*V1*W2
V1=V1*V0
W2=W2/DBLE(J)
8 CONTINUE
C **************
EXACT(1)=2.D0*JK(1)+2.D0
C
DO 3 J=2,KP2,2
3 EXACT(J)=0.D0
C
DO 4 J=3,KP2,4
4 EXACT(J)=-2.D0*JK(J)
C
DO 5 J=5,KP2,4
5 EXACT(J)= 2.D0*JK(J)
C
DO 6 J=1,KP2
6 ERROR(J)=EXACT(J)-AU(1,J)
PRINT 160
PRINT 161,(J-1,AU(1,J),EXACT(J),ERROR(J),J=1,KP2)
PRINT 48
C ***** BESSEL FUNCTIONS ******************
CALL SF33D(Q,KP2,KP2,JK,EXACT,RAB,IERR)
C *****
V1=1.D0
V0=Q/2.D0
W1=1.D0
W2=W1/1.D0
C
DO 800 J=1,KP2
JK(J)=JK(J)*V1*W2
V1=V1*V0
W2=W2/DBLE(J)
800 CONTINUE
C **************
EXACT(1)=2.D0
C
DO 300 J=3,KP2,2
300 EXACT(J)=0.D0
C
DO 400 J=4,KP2,4
400 EXACT(J)=-2.D0*JK(J)
C
DO 500 J=2,KP2,4
500 EXACT(J)= 2.D0*JK(J)
C
DO 600 J=1,KP2
600 ERROR(J)=EXACT(J)-AU(2,J)
PRINT 170
PRINT 161,(J-1,AU(2,J),EXACT(J),ERROR(J),J=1,KP2)
STOP
44 FORMAT(/1X,'RESULTS AFTER THE FIRST CALL SUBROUTINE DE06D'//)
45 FORMAT(/1X,'RESULTS AFTER THE SECOND CALL SUBROUTINE DE06D'//)
47 FORMAT(1X,64(1H-)/)
48 FORMAT(1X,89(1H-)/)
50 FORMAT(1X,'NX=',I3,4X,'XS='/(3D25.16))
51 FORMAT(1X,'XP=',D25.16,5X,'X=',D25.16)
35 FORMAT(1X,'H=',D25.16,3X,'K=',I2,3X,'INIAPR=',I1,3X,'IMAX=',I2,
1 3X,'JSTART=',I1)
60 FORMAT(1X,'YX=',2D25.16/1X,'YT=',2D25.16/1X,'DELTA Y =',2D25.16)
85 FORMAT(5X,I2,5X,2D25.16)
110 FORMAT(/1X,'COEFFICIENTS AU AND AJK ON ',I3,' SEGMENT ')
C
150 FORMAT(' Number of',7X,'Chebyshev coefficients'
* ' for ',I1,' component'/
1 14X,51(1H-)/
2 ' coefficient',11X,'for Y',21X,'for Y'''/
3 1X,64(1H-)/)
160 FORMAT(//' Number of',19X,'Chebyshev coefficients for 1',
* ' component'/
1 14X,76(1H-)/
2 ' coefficient',6X,'approximate',18X,'exact',16X,
3 'absolute error'/
4 1X,86(1H-)/)
161 FORMAT(5X,I2,5X,3D25.16)
170 FORMAT(' Number of',19X,'Chebyshev coefficients for 2',
* ' component'/
1 14X,76(1H-)/
2 ' coefficient',6X,'approximate',18X,'exact',16X,
3 'absolute error'/
4 1X,86(1H-)/)
180 FORMAT(/1X,86(1H*)/)
END
SUBROUTINE F(X,Y,Z,M)
DOUBLE PRECISION X,Y(2),Z(2),R,Q
Q=0.5D0
R =2.D0*X-1.D0
Z(1)=-2.D0*Q*(Y(2)-1.D0)+(1.D0-DEXP(1.D0-Y(1)+DCOS(Q*R)))/(1.D0+X)
Z(2)= 2.D0*Q*(Y(1)-1.D0)+(1.D0-DEXP(1.D0-Y(2)+DSIN(Q*R)))/(1.D0+X)
RETURN
END
Результаты:
RESULTS AFTER THE FIRST CALL SUBROUTINE DE06D
XP= -.5000000000000000D+00 X= .0000000000000000D+00
H= .5000000000000000D+00 K=11 INIAPR=1 IMAX=14 JSTART=1
YX= .1877582561890373D+01 .5205744613957970D+00
YT= .1877582561890373D+01 .5205744613957970D+00
DELTA Y = .0000000000000000D+00 .0000000000000000D+00
----------------------------------------------------------------
COEFFICIENTS AU AND AJK ON 1 SEGMENT
----------------------------------------------------------------
Number of Chebyshev coefficients for 1 component
---------------------------------------------------
coefficient for Y for Y'
----------------------------------------------------------------
0 .3440601623170462D+01 .1342059372335287D+01
1 .1690818267858915D+00 -.1814968541164727D+00
2 -.1137320976131721D-01 -.1059524195184482D-01
3 -.4420444436200430D-03 .4745020646026972D-03
4 .1483978914782698D-04 .1382469503621041D-04
5 .3457975458437643D-06 -.3711881277662207D-06
6 -.7735964763302993D-08 -.7206797540160214D-08
7 -.1287287390210027D-09 .1381808723230371D-09
8 .2159545790970223D-11 .2011845015939279D-11
9 .2794812154524764D-13 -.3005829905724372D-13
10 -.3753036123793196D-15 -.4197353185506070D-15
11 -.4769719528984170D-17 -.3401006689815467D-16
12 -.3542715301891111D-18
----------------------------------------------------------------
Number of Chebyshev coefficients for 2 component
---------------------------------------------------
coefficient for Y for Y'
----------------------------------------------------------------
0 .6579406276647125D+00 .1440601623170462D+01
1 .1814968541164724D+00 .1690818267858918D+00
2 .1059524195184448D-01 -.1137320976131710D-01
3 -.4745020646027070D-03 -.4420444436198796D-03
4 -.1382469503643224D-04 .1483978914786932D-04
5 .3711881278164416D-06 .3457975459519126D-06
6 .7206797388126834D-08 -.7735964788347499D-08
7 -.1381807922053367D-09 -.1287286781754800D-09
8 -.2011821717662790D-11 .2159575151355637D-11
9 .3000002286107092D-13 .2791175493857545D-13
10 .3489962089936113D-15 -.4264946414696963D-15
11 -.4846530016701094D-17 -.7941780913456320D-17
12 -.8272688451516999D-19
----------------------------------------------------------------
RESULTS AFTER THE SECOND CALL SUBROUTINE DE06D
XP= .0000000000000000D+00 X= .1000000000000000D+01
H= .1000000000000000D+01 K=11 INIAPR=1 IMAX=16 JSTART=1
YX= .1877582561890373D+01 .1479425538604203D+01
YT= .1877582561890373D+01 .1479425538604203D+01
DELTA Y = .2220446049250313D-15 .0000000000000000D+00
----------------------------------------------------------------
COEFFICIENTS AU AND AJK ON 2 SEGMENT
----------------------------------------------------------------
Number of Chebyshev coefficients for 1 component
---------------------------------------------------
coefficient for Y for Y'
----------------------------------------------------------------
0 .3876939614481626D+01 -.8381458775717082D-15
1 -.4585666969845331D-16 -.4845369153497486D+00
2 -.6120804691736534D-01 -.6547191987778950D-15
3 -.2775557561562891D-16 .5127459989174128D-02
4 .3214729527285557D-03 -.3216522913903480D-15
5 -.9412738329658138D-17 -.1610725448276222D-04
6 -.6721369257266663D-06 -.1333975247971853D-15
7 -.4137513937923756D-17 .2403173467777737D-07
8 .7516446311980316D-09 -.1754713453532009D-16
9 .2682270999638617D-18 -.2089352055964174D-10
10 -.5226345925825359D-12 -.2720331013401911D-16
11 -.6182570485004343D-18 .1186314365969263D-13
12 .2471488262435964D-15
----------------------------------------------------------------
Number of Chebyshev coefficients for 2 component
---------------------------------------------------
coefficient for Y for Y'
----------------------------------------------------------------
0 .2000000000000000D+01 .1876939614481626D+01
1 .4845369153497479D+00 -.3436311023057026D-15
2 -.8447459582967137D-17 -.6120804691736582D-01
3 -.5127459989174504D-02 -.2760514256419655D-15
4 -.1930811603266178D-16 .3214729527282392D-03
5 .1610725448270591D-04 .3287843088062292D-16
6 .0000000000000000D+00 -.6721369258790109D-06
7 -.2403173465988319D-07 .3287843088062292D-16
8 .9050547041412199D-18 .7516445977185249D-09
9 .2089353515759578D-10 .3916680348103885D-17
10 .6757967666373710D-17 -.5226679549233393D-12
11 -.1187881715734862D-13 -.2664020263068445D-15
12 -.5550042214725927D-17
**************************************************************************************
Number of Chebyshev coefficients for 1 component
----------------------------------------------------------------------------
coefficient approximate exact absolute error
--------------------------------------------------------------------------------------
0 .3876939614481626D+01 .3876939614481626D+01 .0000000000000000D+00
1 -.4585666969845331D-16 .0000000000000000D+00 .4585666969845331D-16
2 -.6120804691736534D-01 -.6120804691736528D-01 .6245004513516506D-16
3 -.2775557561562891D-16 .0000000000000000D+00 .2775557561562891D-16
4 .3214729527285557D-03 .3214729527285751D-03 .1946142899611480D-16
5 -.9412738329658138D-17 .0000000000000000D+00 .9412738329658138D-17
6 -.6721369257266663D-06 -.6721369257237697D-06 .2896535042254487D-17
7 -.4137513937923756D-17 .0000000000000000D+00 .4137513937923756D-17
8 .7516446311980316D-09 .7516446309595219D-09 -.2385096714138463D-18
9 .2682270999638617D-18 .0000000000000000D+00 -.2682270999638617D-18
10 -.5226345925825359D-12 -.5226354721645607D-12 -.8795820247308316D-18
11 -.6182570485004343D-18 .0000000000000000D+00 .6182570485004343D-18
12 .2471488262435964D-15 .2476765118959866D-15 .5276856523901907D-18
-----------------------------------------------------------------------------------------
Number of Chebyshev coefficients for 2 component
----------------------------------------------------------------------------
coefficient approximate exact absolute error
--------------------------------------------------------------------------------------
0 .2000000000000000D+01 .2000000000000000D+01 .0000000000000000D+00
1 .4845369153497479D+00 .4845369153497477D+00 -.1665334536937735D-15
2 -.8447459582967137D-17 .0000000000000000D+00 .8447459582967137D-17
3 -.5127459989174504D-02 -.5127459989174487D-02 .1734723475976807D-16
4 -.1930811603266178D-16 .0000000000000000D+00 .1930811603266178D-16
5 .1610725448270591D-04 .1610725448271494D-04 .9032759349519859D-17
6 .0000000000000000D+00 .0000000000000000D+00 .0000000000000000D+00
7 -.2403173465988319D-07 -.2403173465552604D-07 .4357150715565922D-17
8 .9050547041412199D-18 .0000000000000000D+00 -.9050547041412199D-18
9 .2089353515759578D-10 .2089353517865796D-10 .2106217358883670D-19
10 .6757967666373710D-17 .0000000000000000D+00 -.6757967666373710D-17
11 -.1187881715734862D-13 -.1188370792446492D-13 -.4890767116303012D-17
12 -.5550042214725927D-17 .0000000000000000D+00 .5550042214725927D-17
3) В третьем примере решается задача Коши
y' = qy, y(0) = eq, q = 4, 0 ≤ x ≤ 7
с контролем точности. Точное решение задачи Коши y(x) = eq(1+x). Используется пять сегментов разбиения промежутка интегрирования и на всех пяти сегментах решение приближается частичными суммами ряда Чебышёва одного и того же порядка K + 1 при K = 18; погрешность этого решения вычисляется с помощью оценивающего решения, которое представляется частичной суммой тоже одного и того же порядка K2 + 1 при K2 = 25. Пример содержит пять обращений к подпрограмме с контролем точности. Приводятся вызывающая программа, подпрограмма вычисления правой части дифференциального уравнения, а также результаты счета на произвольном элементарном сегменте после каждого обращения к подпрогрaмме. Результаты включают границы элементарного сегмента [XP, X], приближенное значение Y решения в конце X сегмента, точное значение YT решения в X, относительную погрешность (YT - Y) / Y приближенного решения, рекомендуемую длину следующего элементарного сегмента H, значения параметров K, IMAX, K2, IMAX2, INIAPR, IERR, BUL, значения параметра JSTART на входе в подпрограмму и на выходе из нее, значения переменных NACCEP, NREJEC из общего блока и число обращений к правой части N. Приводятся также вычисленные на каждом элементарном сегменте коэффициенты Чебышёва для решения и его производной.
DOUBLE PRECISION X,EPS,THRESH,H,HMIN,XP,AU,AJK,R
C
C
LOGICAL BUL
C
PARAMETER (K=18,K2=25,M=1)
DIMENSION AU (M,K+2), AJK (M,K+1),
1 R(2*K*K+5*K+3*K*M+K2*K2+5*K2+4*K2*M+15*M)
C__________________________________________________________________
C
DOUBLE PRECISION WC1,WC2,WC3,HD2,HD4,HOLD,WC30,WC12,WC22,WC32
C
COMMON/COM70D/WC1,WC2,WC3,LASN,HD2,HD4,HOLD
COMMON/COM75D/LASN0,LASN2,WC30,WC12,WC22,WC32
C___________________________________________________________________
C
COMMON/COM06D/I(22)
C
COMMON/STAT76/NACCEP,NREJEC
C
C___________________________________________________________________
C
COMMON/BLOCK/N,Q
DOUBLE PRECISION Q,Y,YT,XK
C
EXTERNAL F
C
M0=M
Q=4.0D0
N=0
NREJEC=0
NACCEP=0
C
X=0.D0
Y=EXP(Q)
XK=7.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
C
JSTART=0
BUL=.FALSE.
DO 20 NX=1,5
PRINT 54,JSTART
C
CALL DE06D(F,M0,K0,INIAPR,IMAX,JSTART,Y,X,AU,AJK,H,XP,MEXACT,
1 IU,EPS,THRESH,M0,MESTER,K20,IMAX2,NATTEM,HMIN,BUL,R,
2 IERR)
C
YT=EXP(Q*(X+1.D0))
PRINT 44,NX
PRINT 50,IERR,N,NACCEP,NREJEC,BUL
PRINT 51,XP,X
PRINT 35, H,K0,INIAPR,IMAX,JSTART,K20,IMAX2
PRINT 36, H
PRINT 52,Y,YT,(YT-Y)/Y
C
PRINT 110,NX
DO 10 L=1,M
PRINT 47
PRINT 150,L
PRINT 85,(J-1,AU(L,J),AJK(L,J),J=1,K0+1), K0+1,AU(L,K0+2)
10 CONTINUE
C
PRINT 55
IF (NX .NE. 4) GO TO 20
C
H=XK-X
BUL=.TRUE.
20 CONTINUE
C
STOP
C
35 FORMAT(1X,'H =',D25.16,3X,'K =',I2,3X,'INIAPR=',I1,3X,'IMAX =',I2,
1 3X,'JSTART (exit) =',I1/1X,31X,'K2=',I2,13X,' IMAX2=',I2)
36 FORMAT(1X,'It is advisable to take the next segment H = ',D22.16/)
44 FORMAT(/1X,'RESULTS AFTER ',I1,'-nd CALL SUBROUTINE DE06D'//)
47 FORMAT(1X,64(1H-)/)
50 FORMAT(1X,'IERR=',I2,3X,'N=',I5,3X,'NACCEP=',I3,3X,'NREJEC=',I3,
1 3X,'BUL=',L3)
51 FORMAT(1X,'XP=',D25.16,5X,'X=',D25.16)
52 FORMAT(1X,'Y =',D25.16/1X,'YT=',D25.16/1X,'(YT-Y)/Y=',D25.16/)
53 FORMAT(1X,'It is advisable to take the next elementary segment '
1 'equal to '/1X,' H=',D25.16)
54 FORMAT(1X,'JSTART (input) =',I1)
55 FORMAT(1X,64(1H*)/)
85 FORMAT(5X,I2,5X,2D25.16)
110 FORMAT(/12X,'COEFFICIENTS AU AND AJK ON ',I3,' SEGMENT ')
C
150 FORMAT(' Number of',7X,'Chebyshev coefficients'
* ' for ',I1,' component'/
1 14X,51(1H-)/
2 ' coefficient',11X,'for Y',21X,'for Y'''/
3 1X,64(1H-)/)
C
END
C
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
Результаты:
JSTART (input) =0
RESULTS AFTER 1-nd CALL SUBROUTINE DE06D
IERR= 0 N= 666 NACCEP= 1 NREJEC= 0 BUL= F
XP= .0000000000000000D+00 X= .1000000000000000D+01
H = .1462072113407607D+01 K =18 INIAPR=1 IMAX =28 JSTART (exit) =1
K2=25 IMAX2= 3
It is advisable to take the next segment H = .1462072113407607D+01
Y = .2980957987041719D+04
YT= .2980957987041728D+04
(YT-Y)/Y= .3051014827201587D-14
COEFFICIENTS AU AND AJK ON 1 SEGMENT
----------------------------------------------------------------
Number of Chebyshev coefficients for 1 component
---------------------------------------------------
coefficient for Y for Y'
----------------------------------------------------------------
0 .1839300696370417D+04 .7357202785481671D+04
1 .1283417414302829D+04 .5133669657211319D+04
2 .5558832820675875D+03 .2223533128270354D+04
3 .1716508501676545D+03 .6866034006706186D+03
4 .4093073156462493D+02 .1637229262584993D+03
5 .7927923909155056D+01 .3171169563661981D+02
6 .1291112018849544D+01 .5164448075398151D+01
7 .1812517960576923D+00 .7250071842307492D+00
8 .2234944644573490D-01 .8939778578276779D-01
9 .2456224491789299D-02 .9824897967232359D-02
10 .2434260195698533D-03 .9737040783530104D-03
11 .2196429608173141D-04 .8785718443822619D-04
12 .1818762681781754D-05 .7275050756828591D-05
13 .1391438984679396D-06 .5565757127019844D-06
14 .9891938164326202D-08 .3956803649573376D-07
15 .6567404212323482D-09 .2627175499717027D-08
16 .4089880504035658D-10 .1636112217928698D-09
17 .2396438554643678D-11 .9651977134206291D-11
18 .1307305681532330D-12 .6534000770996862D-12
19 .8022023286874694D-14
****************************************************************
JSTART (input) =1
RESULTS AFTER 2-nd CALL SUBROUTINE DE06D
IERR= 0 N= 1332 NACCEP= 2 NREJEC= 0 BUL= F
XP= .1000000000000000D+01 X= .2462072113407607D+01
H = .1532966265247129D+01 K =18 INIAPR=1 IMAX =28 JSTART (exit) =1
K2=25 IMAX2= 3
It is advisable to take the next segment H = .1532966265247129D+01
Y = .1033321008678100D+07
YT= .1033321008678101D+07
(YT-Y)/Y= .5633066629307278D-15
COEFFICIENTS AU AND AJK ON 2 SEGMENT
----------------------------------------------------------------
Number of Chebyshev coefficients for 1 component
---------------------------------------------------
coefficient for Y for Y'
----------------------------------------------------------------
0 .5095960823292119D+06 .2038384329316864D+07
1 .4098258344442519D+06 .1639303337777021D+07
2 .2292912802021268D+06 .9171651208085136D+06
3 .9617334343582617D+05 .3846933737433065D+06
4 .3195489189291315D+05 .1278195675716525D+06
5 .8749770823017981D+04 .3499908329207163D+05
6 .2032391002635105D+04 .8129564010540043D+04
7 .4093162695215231D+03 .1637265078085889D+04
8 .7269704474073458D+02 .2907881789627752D+03
9 .1154084340356675D+02 .4616337361410507D+02
10 .1655685234835414D+01 .6622740939284641D+01
11 .2166055656541846D+00 .8664222626000115D+00
12 .2603837925257364D-01 .1041535171011070D+00
13 .2894799815643235D-02 .1157919934606169D-01
14 .2992947964282581D-03 .1197179289141259D-02
15 .2891713907189341D-04 .1156686171161492D-03
16 .2621971404995999D-05 .1048792939031817D-04
17 .2239451148232430D-06 .8957763096562044D-06
18 .1807690881716652D-07 .7239128632652303D-07
19 .1382744594359486D-08
****************************************************************
JSTART (input) =1
RESULTS AFTER 3-nd CALL SUBROUTINE DE06D
IERR= 0 N= 1998 NACCEP= 3 NREJEC= 0 BUL= F
XP= .2462072113407607D+01 X= .3995038378654736D+01
H = .1508421188765385D+01 K =18 INIAPR=1 IMAX =28 JSTART (exit) =1
K2=25 IMAX2= 3
It is advisable to take the next segment H = .1508421188765385D+01
Y = .4756312916275247D+09
YT= .4756312916275291D+09
(YT-Y)/Y= .9273451497032786D-14
COEFFICIENTS AU AND AJK ON 3 SEGMENT
----------------------------------------------------------------
Number of Chebyshev coefficients for 1 component
---------------------------------------------------
coefficient for Y for Y'
----------------------------------------------------------------
0 .2283153386407679D+09 .9132613545630736D+09
1 .1860185017542536D+09 .7440740070170157D+09
2 .1069698752921242D+09 .4278795011684965D+09
3 .4645916804004535D+08 .1858366721601802D+09
4 .1604973746530335D+08 .6419894986121231D+08
5 .4580262211113702D+07 .1832104884445409D+08
6 .1110523488616919D+07 .4442093954467108D+07
7 .2336949821032878D+06 .9347799284127976D+06
8 .4339963089823036D+05 .1735985235926468D+06
9 .7207905996657899D+04 .2883162398648638D+05
10 .1082226112071841D+04 .4328904448188316D+04
11 .1482195794845022D+03 .5928783178757106D+03
12 .1865712734438479D+02 .7462850933277682D+02
13 .2172315938322924D+01 .8689263730999869D+01
14 .2352560605988909D+00 .9410242565271432D+00
15 .2381147000857250D-01 .9524591586284714D-01
16 .2261992349068948D-02 .9047974398043834D-02
17 .2024301729855873D-03 .8097148698169576D-03
18 .1712302581414227D-04 .6848667213432691D-04
19 .1372970277624434D-05
****************************************************************
JSTART (input) =1
RESULTS AFTER 4-nd CALL SUBROUTINE DE06D
IERR= 0 N= 2664 NACCEP= 4 NREJEC= 0 BUL= F
XP= .3995038378654736D+01 X= .5503459567420121D+01
H = .1512727617257633D+01 K =18 INIAPR=1 IMAX =28 JSTART (exit) =1
K2=25 IMAX2= 3
It is advisable to take the next segment H = .1512727617257633D+01
Y = .1984569961813829D+12
YT= .1984569961813862D+12
(YT-Y)/Y= .1645384603317564D-13
COEFFICIENTS AU AND AJK ON 4 SEGMENT
----------------------------------------------------------------
Number of Chebyshev coefficients for 1 component
---------------------------------------------------
coefficient for Y for Y'
----------------------------------------------------------------
0 .9614306841453001D+11 .3845722736581171D+12
1 .7799344472085522D+11 .3119737788834184D+12
2 .4443771893936655D+11 .1777508757574646D+12
3 .1907393434757767D+11 .7629573739030995D+11
4 .6502821532118846D+10 .2601128612847512D+11
5 .1829887179448613D+10 .7319548717794297D+10
6 .4372504798242567D+09 .1749001919296878D+10
7 .9064955902453510D+08 .3625982360980564D+09
8 .1658089635656771D+08 .6632358542618788D+08
9 .2711805404723104D+07 .1084722161884197D+08
10 .4009004613238497D+06 .1603601845251839D+07
11 .5405659881262979D+05 .2162263952339132D+06
12 .6698502761906634D+04 .2679401102719428D+05
13 .7677470299048989D+03 .3070988132609295D+04
14 .8184192932842436D+02 .3273677233231283D+03
15 .8153477041859235D+01 .3261391954020655D+02
16 .7623494593312734D+00 .3049405498408305D+01
17 .6714759872954737D-01 .2685999699242529D+00
18 .5589907782479134D-02 .2237515196611639D-01
19 .4412645753107231D-03
****************************************************************
JSTART (input) =1
RESULTS AFTER 5-nd CALL SUBROUTINE DE06D
IERR= 0 N= 3330 NACCEP= 5 NREJEC= 0 BUL= T
XP= .5503459567420121D+01 X= .7000000000000000D+01
H = .1519889233714398D+01 K =18 INIAPR=1 IMAX =28 JSTART (exit) =1
K2=25 IMAX2= 3
It is advisable to take the next segment H = .1519889233714398D+01
Y = .7896296018267691D+14
YT= .7896296018268069D+14
(YT-Y)/Y= .4788637598251465D-13
COEFFICIENTS AU AND AJK ON 5 SEGMENT
----------------------------------------------------------------
Number of Chebyshev coefficients for 1 component
---------------------------------------------------
coefficient for Y for Y'
----------------------------------------------------------------
0 .3842661031210711D+14 .1537064412484283D+15
1 .3110528300639395D+14 .1244211320255758D+15
2 .1764181737953541D+14 .7056726951814186D+14
3 .7528482813802057D+13 .3011393125520830D+14
4 .2550044414557033D+13 .1020017765822807D+14
5 .7126444734746635D+12 .2850577893898591D+13
6 .1690714118893398D+12 .6762856475573130D+12
7 .3479544965157729D+11 .1391817986062768D+12
8 .6317274240572766D+10 .2526909696227131D+11
9 .1025433937757140D+10 .4101735751016442D+10
10 .1504475789562973D+09 .6017903158224826D+09
11 .2013147026984537D+08 .8052588107817611D+08
12 .2475517428784340D+07 .9902069711953914D+07
13 .2815494125856304D+06 .1126197651777321D+07
14 .2978172755044162D+05 .1191269088557214D+06
15 .2944053381332099D+04 .1177621249174327D+05
16 .2731356948213740D+03 .1092541701730341D+04
17 .2387094743611971D+02 .9548266446217895D+02
18 .1971719625379019D+01 .7890468580648303D+01
19 .1543599161706787D+00
Из приведенных результатов наглядно видно, что на всех выбранных подпрограммой DE06D сегментах коэффициенты Чебышёва, как для решения Y, так и для его производной Y', быстро убывают.
4) Решается задача Коши из примера 3 с контролем точности. Используется разбиение отрезка интегрирования на семь сегментов, при этом порядок частичных сумм ряда Чебышёва, представляющих приближенное решение, или порядок частичных сумм для оценивающего решения изменяется от сегмента к сегменту. Поэтому каждое обращение к подпрограмме DE06D (кроме последнего, седьмого) осуществляется с нулевым значением параметра JSTART. Приводится та же информация, что и в примере 3.
DOUBLE PRECISION X,EPS,THRESH,H,HMIN,XP,AU,AJK,R
C
C
LOGICAL BUL
C
C take into consideraitoin: all arraies dependent on K
C must be described for maximal value K
C
PARAMETER (K=18,K2=27,M=1)
C
DIMENSION AU (M,K+2), AJK (M,K+1),
1 R(2*K*K+5*K+3*K*M+K2*K2+5*K2+4*K2*M+15*M)
C__________________________________________________________________
C
DOUBLE PRECISION WC1,WC2,WC3,HD2,HD4,HOLD,WC30,WC12,WC22,WC32
C
COMMON/COM70D/WC1,WC2,WC3,LASN,HD2,HD4,HOLD
COMMON/COM75D/LASN0,LASN2,WC30,WC12,WC22,WC32
C___________________________________________________________________
C
COMMON/COM06D/I(22)
C
COMMON/STAT76/NACCEP,NREJEC
C
C___________________________________________________________________
C
COMMON/BLOCK/N,Q
DOUBLE PRECISION Q,Y,YT,XK,H1
C
EXTERNAL F
C
Q=4.0D0
M0=M
C
N=0
NREJEC=0
NACCEP=0
C
X=0.
Y=EXP(Q)
XK=7.D0
K0=12
INIAPR=1
IMAX=23
IU=2
EPS=0.5D-11
THRESH=1.D0
MEXACT=1
MESTER=1
H=1.D0
HMIN=1.D-3
K20=25
IMAX2=3
NATTEM=3
C
C take into consideraitoin: parameter JSTART must be equal to zero
C
JSTART=0
BUL=.FALSE.
C
DO 20 NX=1,7
PRINT 54,JSTART
C
CALL DE06D(F,M0,K0,INIAPR,IMAX,JSTART,Y,X,AU,AJK,H,XP,MEXACT,
1 IU,EPS,THRESH,M0,MESTER,K20,IMAX2,NATTEM,HMIN,BUL,R,
2 IERR)
C
C
YT=EXP(Q*(X+1.D0))
PRINT 44,NX
PRINT 50,IERR,N,NACCEP,NREJEC,BUL
PRINT 51,XP,X
PRINT 35, H,K0,INIAPR,IMAX,JSTART,K20,IMAX2
PRINT 36, H
PRINT 52,Y,YT,(YT-Y)/Y
C
PRINT 110,NX
DO 10 L=1,M
PRINT 47
PRINT 150,L
PRINT 85,(J-1,AU(L,J),AJK(L,J),J=1,K0+1), K0+1,AU(L,K0+2)
10 CONTINUE
PRINT 55
C
GO TO (1,2,3,4,5,6,20),NX
C
C take into consideraitoin: parameter JSTART must be equal to zero
C
1 JSTART=0
C
K0 =16
IMAX=25
GO TO 20
C
C take into consideraitoin: parameter JSTART must be equal to zero
C
2 JSTART=0
C
K0= 17
IMAX=24
GO TO 20
C
C take into consideraitoin: parameter JSTART must be equal to zero
C
3 JSTART=0
C
K0= 18
IMAX=25
GO TO 20
C
C take into consideraitoin: parameter JSTART must be equal to zero
C
4 JSTART=0
C
K20 = 26
IMAX2= 3
GO TO 20
C
C take into consideraitoin: parameter JSTART must be equal to zero
C
5 JSTART=0
C
K20=27
GO TO 20
C
6 H1=XK-X
IF (H1 .GT. H) GO TO 20
H=H1
BUL=.TRUE.
20 CONTINUE
C
STOP
C
C
35 FORMAT(1X,'H =',D25.16,3X,'K =',I2,3X,'INIAPR=',I1,3X,'IMAX =',I2,
1 3X,'JSTART (exit) =',I1/1X,31X,'K2=',I2,13X,' IMAX2=',I2)
36 FORMAT(1X,'It is advisable to take the next segment H = ',D22.16/)
44 FORMAT(/1X,'RESULTS AFTER ',I1,'-nd CALL SUBROUTINE DE75D'//)
47 FORMAT(1X,64(1H-)/)
50 FORMAT(1X,'IERR=',I2,3X,'N=',I5,3X,'NACCEP=',I3,3X,'NREJEC=',I3,
1 3X,'BUL=',L3)
51 FORMAT(1X,'XP=',D25.16,5X,'X=',D25.16)
52 FORMAT(1X,'Y =',D25.16/1X,'YT=',D25.16/1X,'(YT-Y)/Y=',D25.16/)
53 FORMAT(1X,'It is advisable to take the next elementary segment '
1 'equal to '/1X,' H=',D25.16)
54 FORMAT(1X,'JSTART (input) =',I1)
55 FORMAT(1X,64(1H*)/)
85 FORMAT(5X,I2,5X,2D25.16)
110 FORMAT(/12X,'COEFFICIENTS AU AND AJK ON ',I3,' SEGMENT ')
C
150 FORMAT(' Number of',7X,'Chebyshev coefficients'
* ' for ',I1,' component'/
1 14X,51(1H-)/
2 ' coefficient',11X,'for Y',21X,'for Y'''/
3 1X,64(1H-)/)
C
C
END
C
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
Результаты:
JSTART (input) =0
RESULTS AFTER 1-nd CALL SUBROUTINE DE75D
IERR= 0 N= 426 NACCEP= 1 NREJEC= 0 BUL= F
XP= .0000000000000000D+00 X= .1000000000000000D+01
H = .9339989440424415D+00 K =12 INIAPR=1 IMAX =23 JSTART (exit) =1
K2=25 IMAX2= 3
It is advisable to take the next segment H = .9339989440424415D+00
Y = .2980957987051175D+04
YT= .2980957987041728D+04
(YT-Y)/Y= -.3169089101004236D-11
COEFFICIENTS AU AND AJK ON 1 SEGMENT
----------------------------------------------------------------
Number of Chebyshev coefficients for 1 component
---------------------------------------------------
coefficient for Y for Y'
----------------------------------------------------------------
0 .1839300696375054D+04 .7357202785508012D+04
1 .1283417414306673D+04 .5133669657232480D+04
2 .5558832820697687D+03 .2223533128281319D+04
3 .1716508501684951D+03 .6866034006743304D+03
4 .4093073156484638D+02 .1637229262593788D+03
5 .7927923909197766D+01 .3171169563678836D+02
6 .1291112018856089D+01 .5164448075423461D+01
7 .1812517960606050D+00 .7250071842422237D+00
8 .2234944644816022D-01 .8939778572651956D-01
9 .2456224486281866D-02 .9824897901096714D-02
10 .2434260144623886D-03 .9737042203723919D-03
11 .2196430237098456D-04 .8785732260117168D-04
12 .1818768226197809D-05 .7274916049071489D-05
13 .1391402218433298D-06
****************************************************************
JSTART (input) =0
RESULTS AFTER 2-nd CALL SUBROUTINE DE75D
IERR= 0 N= 984 NACCEP= 2 NREJEC= 0 BUL= F
XP= .1000000000000000D+01 X= .1933998944042441D+01
H = .1463018769583101D+01 K =16 INIAPR=1 IMAX =25 JSTART (exit) =1
K2=25 IMAX2= 3
It is advisable to take the next segment H = .1463018769583101D+01
Y = .1249908453213111D+06
YT= .1249908453209149D+06
(YT-Y)/Y= -.3170104968364770D-11
COEFFICIENTS AU AND AJK ON 2 SEGMENT
----------------------------------------------------------------
Number of Chebyshev coefficients for 1 component
---------------------------------------------------
coefficient for Y for Y'
----------------------------------------------------------------
0 .8037958824974234D+05 .3215183529989694D+06
1 .5424502552893991D+05 .2169801021157596D+06
2 .2230133679673537D+05 .8920534718694149D+05
3 .6490503023351323D+04 .2596201209340530D+05
4 .1453873109271962D+04 .5815492437087844D+04
5 .2640586850763683D+03 .1056234740305480D+04
6 .4028111989831159D+02 .1611244795932513D+03
7 .5293168336248575D+01 .2117267334499627D+02
8 .6106485449183694D+00 .2442594179657737D+01
9 .6276803386794723D-01 .2510721354558518D+00
10 .5816699639627018D-02 .2326679855159702D-01
11 .4906653858040516D-03 .1962661543525729D-02
12 .3797871070547109D-04 .1519148431233242D-03
13 .2715658130067093D-05 .1086264108524082D-04
14 .1804262604877445D-06 .7217141971826990D-06
15 .1119425955930262D-07 .4478026027765480D-07
16 .6515943505424323D-09 .2596067721905015D-08
17 .3561687178441543D-10
****************************************************************
JSTART (input) =0
RESULTS AFTER 3-nd CALL SUBROUTINE DE75D
IERR= 0 N= 2119 NACCEP= 3 NREJEC= 1 BUL= F
XP= .1933998944042441D+01 X= .3026545196739182D+01
H = .1146370222494962D+01 K =17 INIAPR=1 IMAX =24 JSTART (exit) =1
K2=25 IMAX2= 3
It is advisable to take the next segment H = .1146370222494962D+01
Y = .9881558369177736D+07
YT= .9881558369146470D+07
(YT-Y)/Y= -.3164112410904332D-11
COEFFICIENTS AU AND AJK ON 3 SEGMENT
----------------------------------------------------------------
Number of Chebyshev coefficients for 1 component
---------------------------------------------------
coefficient for Y for Y'
----------------------------------------------------------------
0 .5780818042052528D+07 .2312327216821008D+08
1 .4196551137329945D+07 .1678620454931976D+08
2 .1939743920960109D+07 .7758975683840424D+07
3 .6456828492877231D+06 .2582731397150884D+07
4 .1667768332217746D+06 .6671073328870935D+06
5 .3508413894436391D+05 .1403365557774531D+06
6 .6215489215671377D+04 .2486195686268363D+05
7 .9501741790476867D+03 .3800696716190284D+04
8 .1276744097381627D+03 .5106976389517394D+03
9 .1529817230266042D+02 .6119268921043457D+02
10 .1653611639886406D+01 .6614446559295304D+01
11 .1627797659830340D+00 .6511190643731926D+00
12 .1470855313245595D-01 .5883421283884438D-01
13 .1228127053072204D-02 .4912509167649404D-02
14 .9530279603193810D-04 .3812122431057219D-03
15 .6907324619409443D-05 .2762986426496639D-04
16 .4696239783306994D-06 .1878666937571438D-05
17 .3007292939514618D-07 .1198759793652471D-06
18 .1808565583073439D-08
****************************************************************
JSTART (input) =0
RESULTS AFTER 4-nd CALL SUBROUTINE DE75D
IERR= 0 N= 2731 NACCEP= 4 NREJEC= 1 BUL= F
XP= .3026545196739182D+01 X= .4172915419234144D+01
H = .1235733327394152D+01 K =18 INIAPR=1 IMAX =25 JSTART (exit) =1
K2=25 IMAX2= 3
It is advisable to take the next segment H = .1235733327394152D+01
Y = .9688900316055118D+09
YT= .9688900316024432D+09
(YT-Y)/Y= -.3167094533155484D-11
COEFFICIENTS AU AND AJK ON 4 SEGMENT
----------------------------------------------------------------
Number of Chebyshev coefficients for 1 component
---------------------------------------------------
coefficient for Y for Y'
----------------------------------------------------------------
0 .5507695694472854D+09 .2203078277789139D+10
1 .4078135968999383D+09 .1631254387599751D+10
2 .1950262075756098D+09 .7801048303024373D+09
3 .6756364317878412D+08 .2702545727151356D+09
4 .1821515163661657D+08 .7286060654646595D+08
5 .4005985175514870D+07 .1602394070205921D+08
6 .7426759175759568D+06 .2970703670303678D+07
7 .1188853381188680D+06 .4755413524753515D+06
8 .1673472466805787D+05 .6693889867211709D+05
9 .2101253259175835D+04 .8405013036627621D+04
10 .2380650710946103D+03 .9522602843055625D+03
11 .2456750422042004D+02 .9827001685569394D+02
12 .2327487241500224D+01 .9309948989056466D+01
13 .2037808989550966D+00 .8151236634540737D+00
14 .1658309884890176D-01 .6633245250463915D-01
15 .1260490432937310D-02 .5041966184080593D-02
16 .8988051643308878D-04 .3595019759359275D-03
17 .6035210244051902D-05 .2408196235137439D-04
18 .3824385829688405D-06 .1507421870883263D-05
19 .2275429323436656D-07
****************************************************************
JSTART (input) =0
RESULTS AFTER 5-nd CALL SUBROUTINE DE75D
IERR= 0 N= 3348 NACCEP= 5 NREJEC= 1 BUL= F
XP= .4172915419234144D+01 X= .5408648746628296D+01
H = .1222344503872057D+01 K =18 INIAPR=1 IMAX =25 JSTART (exit) =1
K2=26 IMAX2= 3
It is advisable to take the next segment H = .1222344503872057D+01
Y = .1358198193438054D+12
YT= .1358198193433751D+12
(YT-Y)/Y= -.3168039796049561D-11
COEFFICIENTS AU AND AJK ON 5 SEGMENT
----------------------------------------------------------------
Number of Chebyshev coefficients for 1 component
---------------------------------------------------
coefficient for Y for Y'
----------------------------------------------------------------
0 .7385214498624342D+11 .2954085799449741D+12
1 .5626476588406478D+11 .2250590635362594D+12
2 .2832066611716351D+11 .1132826644686539D+12
3 .1042855593594509D+11 .4171422374378013D+11
4 .3003174782840780D+10 .1201269913136295D+11
5 .7074463202524471D+09 .2829785281009670D+10
6 .1407193299953187D+09 .5628773199811665D+09
7 .2419536209274748D+08 .9678144837092318D+08
8 .3661009324781738D+07 .1464403729909931D+08
9 .4943952668126404D+06 .1977581067238737D+07
10 .6026691322536863D+05 .2410676528998309D+06
11 .6693657657416524D+04 .2677463062813714D+05
12 .6826707356569832D+03 .2730682942569212D+04
13 .6435613469215659D+02 .2574245407728013D+03
14 .5639749370523604D+01 .2255900530428335D+02
15 .4616941102191801D+00 .1846773508630577D+01
16 .3546015191904983D-01 .1418332502071280D+00
17 .2564693420986640D-02 .1025289976678323D-01
18 .1750003986023030D-03 .7031627028482035D-03
19 .1131587558177000D-04
****************************************************************
JSTART (input) =0
RESULTS AFTER 6-nd CALL SUBROUTINE DE75D
IERR= 0 N= 3970 NACCEP= 6 NREJEC= 1 BUL= F
XP= .5408648746628296D+01 X= .6630993250500353D+01
H = .1224525460641845D+01 K =18 INIAPR=1 IMAX =25 JSTART (exit) =1
K2=27 IMAX2= 3
It is advisable to take the next segment H = .1224525460641845D+01
Y = .1804650227235185D+14
YT= .1804650227229473D+14
(YT-Y)/Y= -.3165216887347327D-11
COEFFICIENTS AU AND AJK ON 6 SEGMENT
----------------------------------------------------------------
Number of Chebyshev coefficients for 1 component
---------------------------------------------------
coefficient for Y for Y'
----------------------------------------------------------------
0 .9875973615366320D+13 .3950389446146527D+14
1 .7494320213510355D+13 .2997728085404143D+14
2 .3744870485462777D+13 .1497948194185112D+14
3 .1366963361820775D+13 .5467853447283091D+13
4 .3899324361061464D+12 .1559729744424568D+13
5 .9094850718381935D+11 .3637940287352643D+12
6 .1790758183748640D+11 .7163032734993637D+11
7 .3047354370874895D+10 .1218941748349303D+11
8 .4562982357247394D+09 .1825192942895151D+10
9 .6097379302984758D+08 .2438951721168850D+09
10 .7354312361964374D+07 .2941724944607217D+08
11 .8081659198219120D+06 .3232663677802565D+07
12 .8154671371608855D+05 .3261868532254882D+06
13 .7605552804944042D+04 .3042220989152230D+05
14 .6593806950923336D+03 .2637522734543309D+04
15 .5340203935980719D+02 .2136076454278082D+03
16 .4057527906895036D+01 .1623033175617456D+02
17 .2903243740263771D+00 .1161984471604228D+01
18 .1960724484544482D-01 .7935519330203533D-01
19 .1255700875618072D-02
****************************************************************
JSTART (input) =1
RESULTS AFTER 7-nd CALL SUBROUTINE DE75D
IERR= 0 N= 4592 NACCEP= 7 NREJEC= 1 BUL= T
XP= .6630993250500353D+01 X= .7000000000000000D+01
H = .3321060745496823D+01 K =18 INIAPR=1 IMAX =25 JSTART (exit) =1
K2=27 IMAX2= 3
It is advisable to take the next segment H = .3321060745496823D+01
Y = .7896296018293066D+14
YT= .7896296018268069D+14
(YT-Y)/Y= -.3165645632090114D-11
COEFFICIENTS AU AND AJK ON 7 SEGMENT
----------------------------------------------------------------
Number of Chebyshev coefficients for 1 component
---------------------------------------------------
coefficient for Y for Y'
----------------------------------------------------------------
0 .8613410585753402D+14 .3445364234301361D+15
1 .2979974189365701D+14 .1191989675746280D+15
2 .5377474890583328D+13 .2150989956233333D+14
3 .6540696376622693D+12 .2616278550649085D+13
4 .5993282489630859D+11 .2397312995852262D+12
5 .4403202264276050D+10 .1761280905712864D+11
6 .2699288890354053D+09 .1079715556123627D+10
7 .1419491845594438D+08 .5677967383432746D+08
8 .6535185241938498D+06 .2614074090104342D+07
9 .2675432161266559D+05 .1070172973954678D+06
10 .9860332245442534D+03 .3944114536285400D+04
11 .3304344415690883D+02 .1321820547580719D+03
12 .1015205359386096D+01 .4046922683715820D+01
13 .2873995011261581D-01 .1252365112304687D+00
14 .7750034867670807D-03 -.3077507019042969D-02
15 .6639264667270388D-04 .7622957229614258D-02
16 -.9827424780178726D-05 -.1387286186218262D-01
17 -.3699995146846387D-04 .9327411651611328D-02
18 -.8735488693492202D-05 -.7054567337036133D-02
19 -.3862121288260887D-04
Из приведенных результатов четвертого примера наглядно видно, что на всех выбранных подпрограммой DE06D сегментах коэффициенты Чебышёва, как для решения Y, так и для его производной Y', быстро убывают. Отметим, однако, что в значениях последних коэффициентов Чебышёва (как для решения, так и для его производной), вычисленных на седьмом сегменте, по мере уменьшения их абсолютных величин начинает заметнее проявлять себя вычислительная погрешность.