Scientific journal
Modern problems of science and education
ISSN 2070-7428
"Перечень" ВАК
ИФ РИНЦ = 0,737

CONSTRUCTION ОNE-STEP NINE POINTS BLOCK METHOD FOR SOLVING STIFF SYSTEM OF ORDINARY DIFFERENTIAL EQUATIONS

Tursunov D.A. 1 Semenov M.E. 2
1 Osh State University
2 National Research Tomsk Polytechnic University
Исследование посвящено развитию теории численных методов в плане построения линейных неявных n-шаговых k-точечных блочных методов для решения жестких систем обыкновенных дифференциальных уравнений (ОДУ). Приводится пример линейного неявного одношагового девятиточечного блочного метода, записанного в виде формул дифференцирования назад. С помощью метода коллокаций определены коэффициенты блочного метода. Проведена проверка условия согласованности полученных коэффициентов метода, построена область устойчивости, вычислены константы погрешности, установлена сходимость метода и определен порядок точности. Проведены численные эксперименты решения ОДУ в системе MatLAB. Предлагаемый метод является самостартующим и может быть применен для численного решения задачи Коши для уравнений и систем ОДУ первого порядка, в том числе и для жестких систем ОДУ. Вычисление приближенного значения искомой функции в каждой k-ой точке внутри блока независимо друг от друга и может рассматриваться как отдельная задача.
Research is devoted to the development of the theory of numerical methods in terms of constructing a linear implicit n-step k-points block methods for solving stiff systems of ordinary differential equations (ODEs). The example of the one-step linear implicit nine points block method in the form of backward differentiation formulas was written. Coefficients of the presented method were defined with the collocation technique. The conditions for consistency coefficients, the region of stability, the error constants, the convergence and the order accuracy of the method were defined. Numerical experiments of solving ODEs had been carried out with a MatLAB program. The proposed method is the self-starting method and it can be applied for the numerical solution of the Cauchy problem for the first order ODEs as well as system, including for stiff ones. The calculation of the approximate value of the unknown function for each k-th point in a block not depends from each other and it can be considered as a separate task.
ordinary differential equations
order of approximation
stability
linear block methods

Введение

Многие научно-технические задачи описываются с использованием аппарата дифференциальных уравнений. При этом аналитическое решение задачи неизвестно или его не удается выразить в элементарных функциях. Научным сообществом проведено множество исследований, посвященных разработке методов, с помощью которых можно получить приближенное решение дифференциальных уравнений и систем. Особого внимания заслуживают исследования, направленные на разработку численных методов интегрирования жестких систем обыкновенных дифференциальных уравнений (ОДУ), а также систем, полученных в результате дискретизации уравнений в частных производных [1, 6, 8]. Жесткие задачи очень распространены во многих областях прикладных исследований: управление, биология, химико-технологические процессы, электротехника, динамика жидкостей, акустика, пластическая деформация и т.д.

Данное исследование посвящено развитию одношаговых многоточечных блочных методов. Предлагаемый метод является самостартующим и может быть применен для численного решения задачи Коши для уравнений и систем ОДУ первого порядка. Для создания блочных методов необходимо решить задачи: 1) определить коэффициенты уравнений для блоков, включающих k точек; 2) провести проверку коэффициентов метода на согласованность; 3) определить порядок точности (аппроксимации) и область устойчивости метода; 4) установить сходимость метода; 5) провести численные эксперименты.

Теоретический обзор

Для численного решения уравнений и жестких систем ОДУ можно использовать различные методы: формулы дифференцирования назад (ФДН или методы Гира) [4, 7, 8, 10], методы Розенброка, неявные методы типа Рунге-Кутты [2, 3], а также явные разностные схемы [1]. Данные методы устойчивы, шаг интегрирования можно выбирать произвольно, руководствуясь требованиями точности, а не устойчивости. Классические неявные многошаговые методы не являются самостартующими и для их использования требуется «участок разгона». В свою очередь явные методы являются менее трудоемкими при их реализации, но им свойственны определенные недостатки, в частности, условная устойчивость, что существенно ограничивает область их применения.

Блочным будем называть метод, с помощью которого в блоке из k равноотстоящих расчетных точек xn+1, ..., xn+i,…, xn+k вычисляются одновременно новые k значений функции yn+1, ..., yn+i,…, yn+k. Будем называть блочный метод одношаговым, если для расчета значений в новом блоке используется только последняя точка предшествующего блока. Если используются несколько точек предшествующего блока, то будем говорить о многошаговых блочных методах. Впервые идеи создания блочных методов решения уравнений и систем ОДУ были предложены в 1953 году У. Милном. Дальнейшее развитие этих идей принадлежат Дж. Россеру, который предложил оригинальные идеи использования методов Рунге-Кутты. Развитие линейных многошаговых методов для нахождения приближенного решения уравнений и систем ОДУ могут быть получены с использованием различных подходов: метод коллокаций, метод интерполяции степенных рядов, разложение в ряд Тейлора [7, 9, 10].

Коэффициенты блочного метода

Запишем задачу Коши для ОДУ первого порядка в виде

y’=f(x, y), (1)

где 0 < x ≤ b и задано начальное условие y(0) = y0. Для записи n-шагового k-точечного блочного метода в виде ФДН будем использовать: вычисленные значения искомой функции Ym=(yn+1, ..., yn+i, …, yn+k)T, Ym-1=(yn−k+1, yn−k+2, ..., yn)T и производные Fm=(fn+1, ..., fn+i, …, fn+k)T в блоке равноотстоящих расчетных точек xn+1, ..., xn+i,…, xn+k. Здесь m = 1, 2, … номер текущего блока, k – количество точек в блоке, n + i = k∙(m – 1) + i – номер текущей точки в блоке, xn+1 – начало блока, xn+i, i = 2, 3, …, k – 1 – точки внутри блока, xn+k – конец блока. Обозначим через h =xn+i – xn+i–1 размер шага интегрирования, тогда размер блока равен (k –1)∙h. В общем случае n-шаговый k-точечный блочный метод можно записать [5]:

. (2)

Здесь k×k матрицы коэффициентов A(i) и B(i), i = 0, 1, …, n, явный вид которых необходимо определить. Применение формулы (2) генерирует новый блок, включающий k новых значений Ym = (yn+1, ..., yn+i, …, yn+k)T в блоке из k равноотстоящих точек.

Мы предлагаем использовать метод коллокаций для построения одношагового девятиточечного (k = 9) блочного метода в виде ФДН. Предположим, что точное решение у(х) локально представлено на интервале [x0, x0 + k·h] непрерывным решением Y(х) вида

, (3)

где bj – неизвестные коэффициенты, которые необходимо определить и φj(x) – некоторые полиномы степени j=0, 1, ... , k. Мы предлагаем конструировать одношаговый девятиточечный блочный метод в виде ФДН, выбирая φj(x) = xj и используя следующие коллокационные соотношения

Y(xn+j) = yn+j, j = 0, 1, ..., k – 1 и , (4)

где yn+j является приближением для точного решения y(xn+j), fn+k = f(xn+k, yn+k). Использование соотношений (4) приводит к системе линейных алгебраических уравнений, которую необходимо решить для получения коэффициентов bj, j=0, 1, ... , k. Далее полученные коэффициенты bj, j=0, 1, ... , k подставляются в выражение (3) и после некоторых алгебраических вычислений при k = 9 получаем следующее выражение для Y(х):

, (5)

где αj(x)и β9(x) – коэффициенты метода. Далее формула(5) используется для записи первого уравнения одношагового девятиточечного метода в виде ФДН:

yn+9 = (280yn – 2835yn+1 + 12960yn+2 – 35280yn+3 + 63504yn+4 – 79380yn+5 +
+ 70560yn+6 – 45360yn+7 + 22680yn+8 + 2520∙h∙fn+9)/7129. (6)

Для записи остальных восьми выражений конструируемого метода воспользуемся опять выражением (5), применив к нему операцию дифференцирования в точках x=xn+j, j=1, 2, …, 8 соответственно. При этом будем использовать следующее выражение дифференциального исчисления.Приведем в явном виде остальные коэффициенты одношагового девятиточечного метода в форме ФДН:

fn+1 = (–796yn– 427253yn+1/35 + 28336yn+2 – 98336yn+3/3 + 97160yn+4/3–23849yn+5–

–23849yn+5 + 184912yn+6/15 – 12368yn+7/3 + 4924yn+8/7 – 35∙h∙fn+9)/(7129h), (7)

fn+2 = (801yn/8 – 3587yn+1/2 – 154791yn+2/20 + 49483yn+3/3 – 48895yn+4/4 +
+48013yn+5/6 – 46543yn+6/12 + 6229yn+7/5 – 4969yn+8/24 + 10∙h∙fn+9)/(7129h),

fn+3 = ( –2423yn/84 + 10851yn+1/28 – 3081yn+2 – 259573yn+3/60 + 21135yn+4/2 –
–20757yn+5/4 + 6709yn+6/3 – 18867yn+7/28 + 15087yn+8/140 – 5∙h∙fn+9)/(7129h),

fn+4 = (817yn/56 – 3659yn+1/21 + 1039yn+2 – 14426yn+3/3 – 1325yn+4 + 7003yn+5 –
–6793yn+6/3 +12746yn+7/21 – 5113yn+8/56 + 4∙h∙fn+9)/(7129h),

fn+5 = ( –831yn/70 + 1861yn+1/14 – 2114yn+2/3 + 7339yn+3/3 – 7255yn+4 +
+15833y5/10 + 13838y6/3 – 6499y7/7 + 5239y8/42 – 5∙h∙fn+9)/(7129h),

fn+6 = (2563yn/168 – 11481yn+1/70 + 3261yn+2/4 – 7549yn+3/3 + 22395yn+4/4 –
–22017yn+5/2 +280573yn+6/60 +20127yn+7/7 – 16347yn+8/56 + 10∙h∙fn+9)/(7129h),

fn+7 = ( –901y n/28 + 4037yn+1/12 – 8029yn+2/5 + 55783yn+3/12 – 55195yn+4/6 –
–54313yn+5/4 –52843yn+6/3 +1178937yn+7/140 + 5869yn+8/4 – 35∙h∙fn+9)/(7129h),

fn+8 = (1041y n/8 – 9334yn+1/7 + 18578yn+2/3 – 258412yn+3/15 + 64015yn+4/2 –
–126266yn+5/3 +123326yn+6/3–33556yn+7+4134649yn+8/280 + 280hfn+9)/(7129h).

Система (6)-(7) определяют одношаговый девятиточечный метод в виде системы ФДН. Вычислительная сложность метода коллокаций (3)-(5) для нахождения коэффициентов метода определена сложностью нахождения обратной матрицы и имеет кубический порядок O(k3), зависит от размера входа (в нашем случае это количество точек в блоке k). Исследование роста сложности при увеличении размера входа с k = 2 до k = 9, свидетельствуют, что асимптотически сложность не изменяется с ростом порядка блочного метода.

Согласованность коэффициентов и порядок точности

Для проверки согласованности коэффициентов блочного метода (6)-(7) воспользуемся подходом, предложенным в монографии [5]. Для этого вычислим константы погрешности многошаговых методов Сp+1 и сравним их с нулем:

для p=0, 1,…, 8 и C10≠0. (8)

Выполнение условий (8) позволяет установить порядок аппроксимации. Последовательно подставив в выражение (8) коэффициенты ai и bi, i=1, 2,…, 9 из выражения (6) и варьируя p, можно показать, что Сp+1 = 0, p=0, 1, 2,…, 8, но при этом С10 = −252/7129. В этом случае можно заметить, что выражение (6) определяет метод девятого порядка точности. Аналогичные вычисления можно провести для системы (7). При этом Сp+1=0, p=0, 1, …, 8 для всех уравнений из системы (7) и С10=(3722/320805, −7489/2566440, 7549/5988360, −7633/8982540, 7759/8982540, −7969/5988360, 8389/2566440, −9649/641610) соответственно.

Область устойчивости и сходимость

В качестве основы определения области устойчивости блочных методов будем использовать метод характеристического полинома [5]. Для этого рассмотрим задачу

y'= ly, . (9)

Запишем выражение (2) для одношагового (n = 1) блочного метода:

, (10)

явный вид k×k матриц коэффициентов A(1), B(0) и B(1) определен в системе (6)-(7). Применим выражение (10) к тестовой задаче (9), будем иметь

Ym+1=(I – z B)-1A Ym, (11)

где , zC, I – единичная матрица. Явный вид k×k матриц A, В может быть получен нормализацией A(1), B(0) и B(1), определенных в (6)-(7). Доминирующее собственное значение матрицы (I – z B)-1A будет определять рациональную функцию R(z), определяющую область устойчивости метода (6)-(7) на комплексной плоскости

R(z)=R1(z)/R2(z), (12)

здесь R1(z) = 15120 – 60480z + 114660z2 + 136080z3 + 112245z4 + 6728z5 + 2953z6 + 9132z7 + 1680z8, R2(z) = 15120 – 75600z + 182700z2 – 283500z3 + 316365z4 – 269325z5 + 180920z6 –97725z7 + 42774z8 – 15120z9.

На рис. 1 приведен график области устойчивости и полюсы функции устойчивости предложенного девятиточечного метода (6)-(7) в сравнении с областью устойчивости восьмиточечного блочного метода [4], ближайшего по точности.

Надпись:  
Рис. 1. Области устойчивости и полюсы одношаговых блочных методов: k= 8 (пунктиная линия) [11] и k= 9 (сплошная линия)

Области устойчивостей указанных методов лежат вне границ кривой S = ={R(z) ≤ 1: zC} [4]. Области устойчивости методов практически совпадают, при этом предложенный метод (6)-(7) имеет более высокий порядок аппроксимации. Для функции устойчивости (12) определены полюсы с координатами (0.767, 0.000), (–0.454, ±1.463), (0.219, ±1.047), (0.549, ±0.686) и (0.716, ±0.341) соответственно. Наличие указанных полюсов накладывает ограничение на область устойчивости метода, при этом предлагаемый метод (6)-(7) является А(α)-стабильным, где α=72.760 (рис. 1). Блочный метод (2) обладает нуль-устойчивостью, если корни его характеристического полинома

(13)

лежат внутри единичной окружности |Rj| ≤ 1, а корни, принадлежащие единичной окружности |Rj| = 1 – простые, j=1, 2, …, k [5]. Рассмотрим блочный метод (11) при h→0, в этом случае будем иметь и характеристический полином (13):

ρ(R) = det(R∙I– A(1)) =2520∙R8(1 – R)/7129 = 0.

Таким образом, выполнение условий согласованности коэффициентов (9) и нуль-устойчивости (13) позволяет сделать вывод, что метод (11) является сходящимся.

Результаты численных расчетов

Для проведения тестирования предложенного блочного метода (6)-(7) разработан алгоритм численного решения уравнений и систем ОДУ. Данный алгоритм программно реализован в системе MatLAB. Кратко приведем полученные результаты численного эксперимента (табл. 1-2) на задачах различной степени жесткости, имеющих аналитические решения [4, 9, 10]. Для сравнения точности полученных численных решений использована характеристика [4] , которую будем называть максимальной абсолютной погрешность численного решения. Данная погрешность выбирается среди погрешностей для всех расчетных точек интервала интегрирования.

Задача1. Автономное линейное ОДУ: y'= –9y, y(0)=e. Аналитическое решение у=e1-9x.

Задача 2. Автономное нелинейное ОДУ: y'=50/y–50y, y(0)=.

Аналитическое решение y=(1+e–100x)1/2.

Таблица 1

h

Задача 1

Задача 2

MaxE

MaxE [9]

MaxE

MaxE [9]

10–2

1.6291e-11

9.33734e-02

6.0156e-04

1.83156e-02

10–3

3.9879e-13

2.08215e-02

2.5320e-11

3.45336e-02

10–4

2.2906e-12

2.25130e-03

2.0606e-13

8.42901e-03

10–5

1.3794e-11

2.26891e-04

7.0144e-13

9.19344e-04

10–6

3.1240e-10

2.27068e-05

3.2572e-13

9.27340e-05

Задача 3. Жесткая автономная система нелинейных ОДУ

Аналитическое решение y1=e–2x, y2=e–x.

Задача 4. Жесткая система нелинейных неоднородных ОДУ

Аналитическое решение y1=1/(1+х), y2=1+х.

Таблица 2

h

Задача 3

h

Задача 4

MaxE

MaxE [10]

MaxE

MaxE [10]

10–2

1.5364e-12

5.86134e-05

10–3

2.9382e-12

5.33017e-06

10–4

1.1761e-11

3.03964e-06

10–5

5.7333e-12

5.39953e-10

10–6

9.6801e-12

1.12154e-07

10–7

1.0836e-13

1.35196e-11

Из результатов расчетов (табл. 1-2) следует, что с помощью предлагаемого одношагового девятиточечного блочного метода (6)-(7) при величине шага h=10-6…10-2 можно получить численное решение тестовых задач 1-4 с более низкой погрешностью MaxE по сравнению с существующими методами.

Заключение

В работе использован коллокационный подход для получения формул линейных неявных n-шаговых k-точечных блочных методов. На примере одношагового девятиточечного метода в виде формул дифференцирования назад определены коэффициенты блочного метода, проведена их проверка на согласованность, определен порядок аппроксимации, установлена сходимость и область устойчивости предложенного метода. Предложенный блочный метод определяет метод девятого порядка точности, является А(α)-стабильным, α=72.760. Вычисление приближенного значения искомой функции в каждой k-ой точке внутри блока независимо друг от друга и может рассматриваться как отдельная (параллельная) задача. Данные методы можно использовать для решения уравнений и систем ОДУ, в том числе жестких систем.

Работа выполнена при поддержке гранта РФФИ по проекту № 13-01-90903 мол_ин_нр, государственного задания «Наука» № 1.604.2011 (Россия).

Рецензенты:

Шумилов Б.М., д.ф.-м.н., профессор кафедры прикладной математики, ГОУ ВПО Томский государственный архитектурно-строительный университет, г. Томск.

Гришин А.М., д.ф.-м.н., профессор, зав. кафедры физической и вычислительной механики. ФГБОУ ВПО НИ Томский государственный университет, г. Томск.