Многошаговые дискретные схемы и явные, последовательные сплайн-методы решения обыкновенных дифференциальных уравнений на нерегулярном шаблоне Текст научной статьи по специальности «Математика»

Вычислительные технологии

Том 1, № 2, 1996

МНОГОШАГОВЫЕ ДИСКРЕТНЫЕ СХЕМЫ И ЯВНЫЕ, ПОСЛЕДОВАТЕЛЬНЫЕ СПЛАЙН-МЕТОДЫ РЕШЕНИЯ ОБЫКНОВЕННЫХ ДИФФЕРЕНЦИАЛЬНЫХ

УРАВНЕНИЙ НА НЕРЕГУЛЯРНОМ ШАБЛОНЕ*

В. И. Киреев

Московский государственный авиационный институт, Россия

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

Выдающийся математик и механик, академик Н. Н. Яненко в своей многогранной научной деятельности определил наиболее важные направления развития вычислительной механики и математики. В частности, Н. Н. Яненко придавал большое значение развитию схем повышенной точности. Так, в одной из последних статей [1] он относил эту проблему к числу весьма важных, требующих решения. В соответствии с этим в данной работе в дополнение к традиционным способам конструирования численных схем решения обыкновенных дифференциальных уравнений (ОДУ) и в развитие [2, 3] предложен новый подход, основанный на адаптации параметрических соотношений для интегродифференциальных сплайнов (ИД-сплайнов) и их следствий для конструирования новых дискретных схем и дискретно-непрерывных сплайн-схем второго и третьего порядков.

1. Дискретные схемы второго и третьего порядков, построенные на нерегулярном шаблоне

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

аУ

— = ^ (х,у), У (а) = Уо, х е [а,Ь]. ах

*© В. И. Киреев, 1996.

Полученные в работе дискретные схемы второго и третьего порядков следуют из параметрических соотношений, относящихся к дифференциальным и интегродифференциальным сплайнам, и следствий из этих соотношений — аппроксимационных формул для производных, а также из обобщенной на нерегулярный шаблон квадратурной формулы парабол, экстраполяционной формулы и др. [4]. Так, из квадратных дифференциальных сплайнов и условий непрерывности производных следуют параметрические соотношения функционально-дифференциального типа:

т. + т. _ 2А/г тг _ тг_1 _ А/г+1 _ А/ (!)

кг ' кг+1 кг к2+1

(Здесь и ниже / обозначает аппроксимируемую функцию в задаче теории приближения,

Хг+1

т г _ & тг _ /í, 1г+1 _ I / (x)dx, кг+1 _ хг+1 _ хг, А/г _ /г _ Л_Ъ Аг+1 _ —

Хг

параметр нерегулярности сетки.)

Из соотношений (1) в работе [2] получены аппроксимационные формулы для производных на нерегулярном шаблоне (хг_1,хг,хг+1) (к _ уаг):

1

т ¿_1 _

Н+1

Н2+1, +(Н+)!* __!_{

к /г_1 + к к /г А /г+1

кг кг кг+1 Аг+1

те г

г+1

Н

т г+1

тг _

А , + Нг+1Акг+1 , + 1 , ■

_Аг+1/г_1 +--ГГ-/г + Т-/г+1

кг кг+1 Аг+1

г+1

6

М3.г

1

Н+1 2 "

(Нг+1)2 н2(г+1)

Аг+1 /г_ 1--ТГ-/г +--Т-Уг+1

к к

Н +1

гкг +1

кг / ^ кг кг+1

к

+1

1

/г + 7 Уг+1 к +1

(2)

(3)

(4)

здесь Нр(г+1) _ ккг + ркг+1, Ыъ,

вир |/'"(х)|, в скобках рядом с формулами

[Хг-1,Хг+1]

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

Рассмотрение нескольких типов квадратных ИД-сплайнов приводит к обобщенной на к _ уаг двухинтервальной квадратурной формуле парабол

Г+1 + — Г

к2+1 г + к2 г_1 3

—/ 1^1 I 1

к к к +1

/г + 7 /г+1 к +1

(5)

к интегродифференциальной параметрической связи

Г

+1

1г 1

, _ т (кг+1?Пг+1 + 2Щ+1тг + кгтг_1)

к +1 к 6

квадратурным формулам

к?

_1

6Н +1

н?г+1). + Н+1н?(г+1) . .

"/г_1 + 1.2^ /г I. /г+1

к2

к2к

+1

к

+1

1

Л3 ( 1 Н^+1 Н^+1 Н 2(г+^

р+1 = Лг+1 1 г + Нг Н3г г + Н3г г | (7)

/г = 6Н+ Л/г-1 + /г + ~Л2+Г/т ь (7)

/-1 = Л (/^-1 + 2/г) - |т Мз,г) ,

= Нг/г-1 + (тт + 2тпг-1) ^Мз,г н+1 = /Н|Г1) , - 3(Н-+1)2 + 2Лг+1Нг-1 , +

1г 6Н+1 ^ Лг-1 /г-2 Лг-1Лг /г-1 +

+ 6ЛгНг-1 + Лг+1(3Лг-1 + 6Нг + 2Лг+1)Д (§)

ЛгЛг+1 /

лево- и правосторонним одноинтервальным интегрофункциональным аппроксимацион-ным формулам для тг-1, тг+1

- 6/г-1 2(/ +2/ ) _ 2 ( + ) 6/г-1 (Л2 М N

"- = "ЛГ - £(/г + 2/г-1), ^ = £(/г-1 + 2/г) - Мз^ ,

трехинтервальным интегральным формулам для тг, ?_г+1 2

m г = a

= a

h2 - h2+1 _i+2 , 3fci+iHi+2 + (fe?+2 - _ - НЦ2H2++1),i ' h 1 i+1 + h 1 h 1 i-1 hi+1 hi+1 hi

Hi+1Hf(h2 - h2+2) + 3hi+! Hi+1 , h2+! - h

Ti+2 __ (h2 h2+2) + 3hi+1 Hi+ H+1 + h2+1 h2+2 li Ji+1 z, Ji + г, Ji-1

hi+2 i+1 hi+1 i hi

(9)

(10)

где A = h2+1(H2+1 + 2hi+2) + hi+1(h2 + h2+2) + hihi+2(H3(i+1) + hi+2). Формулы (9), (10) при h = const становятся двухинтервальными:

^i = ^ (If1 - Ii_1>, mn i+1 = 1 (If? - 1+1). 2. Явные двухшаговые схемы второго порядка

Схема 2Я2А (2 — "шаговость", Я — "явная", 2 — порядок точности, А — модификация) получается аппроксимацией (dy/dx)Xn по формуле (2):

Уп+1 = Уп - ¿ПцДУп + Hn+l8n+1Fn (h = var), Уп+1 = Уп_1 + 2hFn (h = const).

При h = const (регулярный шаблон) схема 2Я2А есть двухшаговая схема Эйлера. Из оценки порядка аппроксимации (2) следует, что порядок аппроксимации схемы при hn+1 < hn повышается на единицу без изменения количества точек шаблона.

2

Из аппроксимации второй производной (4) следует функционально-дифференциальная схема 2Я2Б:

Уп+i = -¿п+1Уп-1 + —h— Уп + H,n+1F/(xra,yra) (h = var),

hn 2

Уп+i = -y/n-1 + 2y„ + h2Fn (h = const).

Из второго параметрического соотношения из (1) следует схема 2Я2В:

yn+i = yn + ¿П+I A/ + 'П+I ( - ) (h = var)'

hn+i hn

Уп+i = Уп + АУп + hAF„ (h = const).

Путем преобразования двух соотношений (1) легко получается схема Адамса—Бэшфорта (2Я2Г):

h2 ( Hra+i 1 \

Уп+i = Уп + "^т1 7—г2— Fn - — Fn-И (h = var). 2 \h™h™+i h2 /

3. Неявные одно- и двухшаговая схемы второго порядка

Первая схема — это классическая одношаговая схема Эйлера—Коши (1НЯ2) (НЯ - "неявная"), следующая из первого параметрического соотношения из (1), а вторая двухшаговая (2НЯ2) — из аппроксимационной формулы (3):

Уп+i = Уп + "у1 [Fn + F (xn+i, Уп+i)],

hn+1

Уп+i =

tr2(n+i) Hn

(Hn+i)2

h h Уп - ¿n+i/n-i + H2+iF(Xn+i ,Уn+l)

i tn 10'

',nhn+1

Последняя схема при h = const преобразуется к виду

1 4 2

Уп+i = - 3У/n-i + 3Уп + 3hF(х2+1,УП+1).

4. Явная трехшаговая схема третьего порядка (3Я3)

Эта схема вытекает из квадратурной формулы (8) и имеет вид

h2

/ = / + ' п+1 уп+" = Уп + 6Hn+1

Яэ2^+1) F 3(H2+")2 + 2hn+iH2-i _ +

-h-Fn-2--h-h-Fn-1 +

hn—1 hn-1hn

6'пН22-1 + hn+1 (Нз(22- i) + 2hn+1)

+-h"h-

hnhn+1

При h = const 3Я3 преобразуется к известной схеме, легко получающейся интегроинтер-поляционным методом [4]:

h

Уп+1 = Уп + —(5Fn-2-s16Fn-i + 23Fn).

5. Неявные двух- и одношаговые схемы третьего порядка

Из обобщенной формулы парабол (5) следует первая двухшаговая схема (2НЯ3А)

Л2

г2 Л - | ""п+1

Уп+1 = Уп - ¿п+1ДУп + -у-

1 2H n+i 1

^п-1 + -Г-^- К + т-F (Xn+1,yn+l)

которая при т = const преобразуется к известной:

-

Уп+1 = Уп + 3 [Fn-1 + 4Fn + F(хп+1,Уп+1)] •

Одноинтервальная, трехточечная квадратурная формула (7) определяет вторую двухша-говую схему 2НЯ3Б:

yn+1 = +

-3

тп+1

бдп+1

1 н "+1^ п+1 н 2(п+1)

— ^F"-1 +--- -2 " F" +--72-F (xn+1,yn+1)

тп тп+1

п+1 п+1

которая при т = const упрощается и принимает следующий вид [4]:

-

Уп+1 = Уп + 12 [—Fn-1 + 8F" + 5F(Х"+1, Уп+1)] •

Две неявные одношаговые схемы 1НЯ3А и 1НЯ3Б выведены из квадратурных формул (8), а двухшаговая схема 2НЯ3Д получена из параметрической связи (6) —

1НЯ3А : 1НЯ3Б : 1НЯ3Д :

- -2

Уп+1 = Уп + ^+1 [Fn + 2F (Ж"+1,Уп+1)]--"+1 F'(xn+1,yn+1),

3

б

Уп+1 = Уп + -n+1F" + [2F" + F'(xn+1, Уп+1)],

б

-

Уп+1 = Уп + ^п+1АУп + -п+i + 2Н£+^п + -п+1^(жп+1, Уп+1)].

б1

Все эти схемы имеют повышенно-дифференциальный тип, так как в их правые части входят вторые производные от решения У (ж).

В работе выполнен анализ устойчивости некоторых двухшаговых схем на нерегулярном шаблоне путем проверки условия корней однородных характеристических уравнений. В табл. 1 приведены значения корней для шести двухшаговых схем.

Таблица 1

Схема 2я2а 2я2б 2я2в 2ня2 2ня3а 2ня3д

qi h = const 0<£п+1<2

91 ¿2 0п+1 ¿п+1 0п+1 1 1 1 1

q2 -1 1 1 1/3 0 ^ 4/5 ¿2 -°п+1 ¿п+1

Из анализа q1 и q2 можно сделать следующие выводы.

Для обеспечения устойчивости схем 2Я2А, 2НЯ3А на нерегулярном (в общем случае) шаблоне необходимо выбрать параметр неравномерности сетки ^п+1<1.

Чтобы обеспечить устойчивость схем 2Я2Б, 2Я2В, 2НЯ3Д, для которых при - = const условие корней не выполняется, необходимо формировать неравномерную сетку с $п+1 <1.

Схема 2НЯ2 при формировании сетки с 0<$п+1<2 для обеспечения устойчивости не требует выполнения условия ^п+1<1.

6. Явные последовательные сплайн-методы второго и третьего порядков

В основу данных дискретно-непрерывных методов кладутся дискретные решения второго и третьего порядка (шаг "С-1") и соответствующие по порядку сплайн-функции — квадратные и кубические (шаг "С-2").

В качестве этих сплайн-функций могут быть приняты

52(ж) = Уп + т„(ж - Жп) + т— ( ^Уга+1 - т-п) (ж - ж„)2 (г = 0,Ж - 1), (11)

ьп+1 \ Ьп+1 /

ь ь

с / ч - ( \ , (3Ауп+1 3тп 2 .

5з(ж) = Уп + тп(ж - Жп) + —Т2--Т---Т- (ж - Жп) +

V Ьп+1 ьп+1 ьп+1 /

+12^ (-АДУп+1 + 2?тп + Агйп+И (ж - Жп)3 (г = 0, N - 1). (12)

ьп+1 V Ьп+1 /

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

Пусть п = 0, N. На каждом отрезке [жп,жп+1]:

а) На шаге С-1 рассчитываются У/п+1, Уп+1 = тп+1 по одной наиболее приемлемой явной (несоставной) или явно-неявной (составной) схеме, скомпонованной из совокупности схем — явной и неявной. Решение, полученное на этом шаге, может быть названо опорным. Порядок точности схемы должен соответствовать порядку сходимости сплайн-функции.

б) На шаге С-2 на каждом отрезке [жп, жп+1] определяется сплайн-функция (11) или (12). Для этого необходимо рассчитать коэффициенты многочленов 52(а^,ж) (г = 0, 2) и

£3(аг,ж)

(г = 0, 3). Коэффициенты а0 определяются значением Уо и значениями Уп, вычисленными на предыдущем шаге, а коэффициенты а1 для 52(а^,ж) — дифференцированием 52(а^,ж) на предыдущем отрезке (т0 = ^(ж0,У0)). Коэффициенты а1 для 53(а^,ж) вычисляются по формуле тп = ^(жп,Уп) (п = 0, N). Коэффициенты а2 для 52(а^, ж)|[Хп'Хп+1] и а2, аз для 53 (а^, ж) |[Хп' Хп+1] (п = 1, N - 1) вычисляются с помощью подстановки в их формулы значений Уп+1, тп для ж) и Уп+1, тп, ттп+1 для 53(а*,ж).

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

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

^у/^ж = ж + У, у(0) = 1, ж € [0, 0.4065].

Задача решалась на неравномерной сетке с $п+1 = 0.9 (ж = 0, 0.15, 0.285, 0.4065) по составным (явным) схемам типа "предиктор-корректор", образующих следующие сочетания: П1К2, П2К2, П2К3, П3К3. (Здесь П и К указывают на шаги "предиктор"и "корректор", а следующие за ними цифры — на порядок точности схем.) Из сочетаний схем следует, что для схемы К2 принимались предикторы первого и второго порядков, а для К3 — второго и третьего.

Результаты расчетов (п = 1, 2, 3) по указанным составным схемам приведены в табл. 2, где указаны явные и неявные схемы, составляющие шаги "предиктор"и "коррек-тор"соответственно. В последнем столбце этой таблицы приведено точное решение, а в скобках — расхождение численных и точных решений (в процентах).

Таблица 2

Схема хп п1 к2 п2 к2 п2 кз пз кз Точное решение

1я1 1ня2 2я2а 1ня2 2я2а 2ня3а зяз 2ня3а

0.150 1.172711 (0.08) 1.174097 (0.036) 1.173700 (0.0027) 1.173700 (0.0027) 1.173668

0.285 1.372442 (0.15) 1.375370 (0.061) 1.374445 (0.0057) 1.374467 (0.00415) 1.374524

0.4065 1.593437 (0.20) 1.597867 (0.079) 1.596537 (0.0043) 1.596576 (0.0019) 1.596606

Анализ приведенных в табл. 2 результатов указывает на повышение точности расчетов при повышении порядков схем как на шаге "П", так и на шаге "К". Минимальную погрешность имеют результаты, полученные по схеме П3К3, однако схема П2К3 тоже дает хорошие результаты.

Таблица 3

Схема жп П1К2 П2К2 П2К3 Точное значение

0.150 1.302812 (1.57) 1.321290 (0.18) 1.323700 (0.0024) 1.323668

0.285 1.656175 (0.47) 1.660534 (0.06) 1.659445 (0.0048) 1.659524

0.4065 1.981588 (1.07) 2.001969 (0.06) 2.003037 (0.0035) 2.003106

Получено также непрерывно-дискретное решение, в котором на шаге С-1 принимались схемы П1К2, П2К2 — для сплайна (11) и П2К3 — для сплайна (12) .

Непрерывные сплайн-решения в силу их коллокационных свойств можно использовать для анализа поведения производных. Результаты численных значений производных, вычисленных дифференцированием 52 (ж) и 53 (ж) , приведены в табл. 3. В скобках ниже численных значений производных указаны проценты их отличия от точных.

Как видно из табл. 3, точность вычисления производных, полученная с использованием 52(ж) и опорного решения П2К2, существенно превышает точность, полученную по 52(ж) и П1К2, а точность для 5з(ж) и П3К3 примерно на порядок лучше точности для 52(ж) и П2К2.

Список литературы

[1] ЯНЕНКО Н. Н. Проблемы вычислительной механики. В "Научные основы прогр. техники и технол.", Машиностроение, М., 1986.

[2] Киреев В. И., Патрикеева Т. К. Интегродифференциальные консервативные сплайны и их применение в интерполяции, численном дифференцировании и интегрировании. В "Вычисл. технологии", 4, №10, ИВТ СО РАН, Новосибирск, 1995, 233-244.

[3] КирЕЕв В. И. Интегральный метод приближения функций алгебраическими многочленами и биквадратными сплайнами. Вестник МАИ, 1, №1, 1994, 48-57.

[4] Самарский А. А., Гулин А. В. Численные методы. Наука, М., 1989.

Поступила в редакцию 30 июня 1996 г.