Анализ устойчивости вычислительных схем

Инженерно-технологическая академия

Ростов-на-Дону  Таганрог

Издательство Южного федерального университета

2018

МИНИСТЕРСТВО НАУКИ И ВЫСШЕГО ОБРАЗОВАНИЯ

РОССИЙСКОЙ ФЕДЕРАЦИИ

Федеральное государственное автономное образовательное

учреждение высшего образования

«ЮЖНЫЙ ФЕДЕРАЛЬНЫЙ УНИВЕРСИТЕТ»

А. Н. ЦЕЛЫХ

В. С. ВАСИЛЬЕВ

Э. М. КОТОВ

АНАЛИЗ УСТОЙЧИВОСТИ 
ВЫЧИСЛИТЕЛЬНЫХ СХЕМ 

Учебное пособие

по курсу

ЧИСЛЕННЫЕ МЕТОДЫ

Содержание

4

УДК 519.62(075.8)+004.42(075.8)
ББК  22.193я73+32973-018я73

Ц349

Печатается по решению кафедры информационно-аналитических систем
безопасности Института компьютерных технологий и информационной 

безопасности Южного федерального университета 

(протокол № 6 от 22 февраля 2018 г.)

Рецензенты:

доктор технических наук, профессор, заведующий кафедрой 
информатики Таганрогского института имени А. П. Чехова 

(филиал) РГЭУ (РИНХ) Я. Е. Ромм

доктор технических наук, профессор, профессор кафедры

информационно-аналитических систем безопасности 
Южного федерального университета А. В. Боженюк

Целых, А. Н.

Ц349 
Анализ устойчивости вычислительных схем : учебное пособие 

по курсу «Численные методы» / А. Н. Целых, В. С. Васильев, 
Э. М. Котов ; Южный федеральный университет.  Ростов-на-Дону ; 
Таганрог : Издательство Южного федерального университета, 2018. 
 146 с.

ISBN 978-5-9275-2912-4

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

Пособие предназначено для студентов высших учебных заведений, 

обучающихся по направлению 10.03.01 «Информационная безопасность» 
(профиль «Информационно-аналитические системы безопасности») по 
курсу «Численные методы».

УДК 519.62(075.8)+004.42(075.8)

ББК  22.193я73+32973-018я73

ISBN 978-5-9275-2912-4

© Южный федеральный университет, 2018
© Целых А. Н., Васильев В. С., Котов Э. М., 2018
© Оформление. Макет. Издательство

Южного федерального университета, 2018

СОДЕРЖАНИЕ

СОДЕРЖАНИЕ ………………………………………………………..
3

ВВЕДЕНИЕ ……………………………………………………………..
5

ГЛАВА 1. АНАЛИЗ УСТОЙЧИВОСТИ МЕТОДАМИ СПЕКТРАЛЬНОГО АНАЛИЗА ……………………………………………..
7

1.1. Теоремы устойчивости …………………………………………
7

1.2. Спектральный анализ устойчивости …………………………..
12

1.3. Энергетические соотношения …………………………………
17

1.4. Сеточный аналог теоремы Ляпунова ………………………….
19

1.5. Проверка условий теоремы Самарского ………………………
20

1.6. Нелинейные модели ……………………………………………
22

1.7. Анализ симметричной нелинейной модели …………………..
28

ГЛАВА 2. ПРОГРАММНАЯ РЕАЛИЗАЦИЯ РЕШЕНИЯ ЗАДАЧ 
КОШИ ДЛЯ ОБЫКНОВЕННЫХ ДИФФЕРЕНЦИАЛЬНЫХ 
УРАВНЕНИЙ …………………………………………………………..
41

2.1. Постановка задачи ……………………………………………...
41

2.2. Анализ системы уравнений исходной модели ………………..
42

2.3. Свойства исходной модели …………………………………….
45

2.4. Программная реализация ………………………………………
51

ГЛАВА 3. ИСПОЛЬЗОВАНИЕ ФУНКЦИЙ БИБЛИОТЕКИ GNU 
SCIENTIFIC LIBRARY (GSL) …………………………………………
56

3.1. Система уравнений ……………………………………………..
56

3.2. Простые вычислительные схемы ……………………………...
57

3.3. Симметричные явно-неявные схемы ………………………….
64

ГЛАВА 4. ЧИСЛЕННЫЕ МЕТОДЫ РЕШЕНИЯ СИСТЕМ ЛИНЕЙНЫХ АЛГЕБРАИЧЕСКИХ УРАВНЕНИЙ …………………..
70

4.1. Численное интегрирование …………………………………….
70

4.2. Численное решение обыкновенных дифференциальных урав
нений …………………………………………………………………..…
71

4.3. Простейшее уравнение переноса ………………………………
72

4.4. Матричные разложения ………………………………………...
75

Содержание

4

4.5. Стационарные итерационные методы решения систем линей
ных алгебраических уравнений ………………………………………...
80

4.6. Нестационарные итерационные методы решения систем ли
нейных алгебраических уравнений …………………………………….
82

ГЛАВА 5. УСТОЙЧИВОСТЬ ЧИСЛЕННЫХ МЕТОДОВ ……….
87

5.1. Постановка задачи ……………………………………………...
87

5.2. Простейшая вычислительная схема …………………………...
88

5.3. Устойчивость простейшей вычислительной схемы …………..
89

5.4. Развитая вычислительная схема ……………………………….
89

5.5. Устойчивость развитой вычислительной схемы ……………...
92

5.6. Устойчивость численных методов решения систем обыкно
венных дифференциальных уравнений ………………………………..
98

ГЛАВА 6. ИНТЕРПОЛЯЦИЯ И ПРИБЛИЖЕНИЕ ФУНКЦИЙ ...
107

6.1. Постановка задачи интерполяции ……………………………..
107

6.2. Локальная интерполяция ……………………………………….
108

6.3. Глобальная интерполяция ……………………………………...
111

6.4. Полином Лагранжа ……………………………………………..
112

6.5. Метод наименьших квадратов …………………………………
113

ГЛАВА 7. МЕТОДЫ РЕШЕНИЯ ОБЫКНОВЕННЫХ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ …………………………………….
116

7.1. Постановка задачи ……………………………………………...
116

7.2. Аналитические методы ………………………………………...
117

7.3. Численные методы решения. Правило Рунге …………………
120

КОНТРОЛЬНЫЕ ВОПРОСЫ ………………………………………..
126

ЗАДАНИЯ ДЛЯ ЗАКРЕПЛЕНИЯ МАТЕРИАЛА ………………….
129

ПРАКТИЧЕСКАЯ РАБОТА. Решение задачи Коши ……………….
129

ЛАБОРАТОРНАЯ РАБОТА. Численные методы решения нелинейных уравнений …………………………………………………………...
131

ЗАКЛЮЧЕНИЕ ………………………………………………………...
143

СПИСОК ЛИТЕРАТУРЫ …………………………………………….
144

1.1. Теоремы устойчивости

5

ВВЕДЕНИЕ

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

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

При выборе численного метода необходимо учитывать две группы 

требований:

—
дискретная модель должна быть адекватной математической 
модели;

—
численный метод должен быть корректным и реализуемым.

Для обеспечения адекватности дискретная модель должна обладать 

свойствами сходимости численного метода, выполнения дискретных аналогов сохранения и качественно правильного поведения решения [1].

Сходимость численного метода, например, означает, что при 

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

Корректность численного метода означает, что дискретная задача 

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

Задача или алгоритм решения задачи называются вычислительно 

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

Введение

6

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

Устойчивость результата к малым возмущениям входных данных 

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

Из-за наличия ошибок округления при вводе входных данных в 

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

1.1. Теоремы устойчивости

7

АНАЛИЗ УСТОЙЧИВОСТИ МЕТОДАМИ 
СПЕКТРАЛЬНОГО АНАЛИЗА

1.1. Теоремы устойчивости

Теорема Ляпунова [1].
Пусть в Ω существует непрерывная вместе с частными производ
ными первого порядка положительно определённая функция V(𝐱) такая, 
что функция

𝑊(𝑡, 𝒙) = (𝑔𝑟𝑎𝑑𝑉, 𝒇(𝑡, 𝒙))

удовлетворяет неравенству

𝑊(𝑦, 𝒙) ≤ 0 для 𝑡 > 0, 𝒙 ∈ 𝛺.

Тогда тривиальное решение системы

𝑑𝒙
𝑑𝑡 = 𝒇(𝑡, 𝒙)

устойчиво.

Теорема Самарского [2].
Если 𝐀 = 𝐀∗ – самосопряжённый положительный оператор и суще
ствует оператор 𝐁−1, то для устойчивости схемы (двухслойная итерационная схема в канонической форме)

𝑩

𝒚𝑘+1−𝒚𝑘

𝜏
+ 𝑨𝒚𝑘 = 𝟎, 𝑘 = 0;  1;  2; …, 𝒚0 = 𝒖0

в 𝐻𝐴: ‖𝒚𝑘+1‖𝐴 ≤ ‖𝒚𝑘‖𝐴 необходимо и достаточно, чтобы выполнялось 
неравенство

(𝑩𝒚, 𝒚) −

1

2 𝜏(𝑨𝒚, 𝒚) ≥ 0 для всех 𝒚 ∈ 𝐻𝐴 (или 𝑩 ≥

1

2 𝜏𝑨).

Интерпретация теоремы Ляпунова. 

Пусть

1. Анализ устойчивости методами спектрального анализа

8

{𝜉′ = 𝑓1(𝑡, 𝜉, 𝜂),
𝜉(0) = 𝛿1,

𝜂′ = 𝑓2(𝑡, 𝜉, 𝜂),
𝜂(0) = 𝛿2

и функция V(ξ, η) существует. Умножим левую и правую части первого 
уравнения на V′ξ, а второго – на V′η, и сложим полученные уравнения:

𝑉′𝜉𝜉′ + 𝑉′𝜂𝜂′ = 𝑉′𝜉𝑓1(𝑡, 𝜉, 𝜂) + 𝑉′𝜂𝑓2(𝑡, 𝜉, 𝜂).

Правая часть результирующего уравнения и является функцией 

W(t, ξ, η) = V′ξf1(t, ξ, η) + V′ηf2(t, ξ, η), а левая представляет собой полную 
производную функции V(ξ, η) по t

𝜕𝑉

𝜕𝜉

𝑑𝜉

𝑑𝑡 +

𝜕𝑉

𝜕𝜂

𝑑𝜂

𝑑𝑡 =

𝑑𝑉

𝑑𝑡.

Таким образом, если W(t, ξ, η) ≤ 0, то 

dV

dt = W(t, ξ, η) ≤ 0, т. е. 

функция V(ξ, η) с ростом t > 0 не возрастает. С другой стороны, как положительно определённая функция V(ξ, η) ≥ 0 ограничена снизу значением 0. Поэтому возмущения в решении просто не могут нарастать.

Схема доказательства теоремы Самарского. 

Представим 𝐲k в виде 𝐲k =

1

2 (𝐲k + 𝐲k+1) −

1

2 (𝐲k+1 − 𝐲k) и умножим 

скалярно левую и правую части двухслойной итерационной схемы на 
(𝐲k+1 − 𝐲k)

(𝑩 𝒚𝑘+1 − 𝒚𝑘

𝜏
, 𝒚𝑘+1 − 𝒚𝑘) + (𝑨 (1

2 (𝒚𝑘 + 𝒚𝑘+1) − 1

2 (𝒚𝑘+1 − 𝒚𝑘)) , 𝒚𝑘+1 −  𝒚𝑘) =

= (𝟎, 𝒚𝑘+1 − 𝒚𝑘) = 0.

Перегруппируем слагаемые:

(𝑩(𝒚𝑘+1 − 𝒚𝑘), 𝒚𝑘+1 − 𝒚𝑘) −

1

2 𝜏(𝑨(𝒚𝑘+1 − 𝒚𝑘), 𝒚𝑘+1 − 𝒚𝑘) +

+

1

2 𝜏(𝑨(𝒚𝑘 + + 𝒚𝑘+1), 𝒚𝑘+1 − 𝒚𝑘) = 0.

Раскроем скобки в последнем слагаемом

(𝑨(𝒚𝑘 + 𝒚𝑘+1), 𝒚𝑘+1 − 𝒚𝑘) = (𝑨𝒚𝑘, 𝒚𝑘+1) − (𝑨𝒚𝑘, 𝒚𝑘) + (𝑨𝒚𝑘+1, 𝒚𝑘+1) −

−(𝑨𝒚𝑘+1, 𝒚𝑘).

Если оператор 𝐀 = 𝐀∗ – самосопряжённый, то

1.1. Теоремы устойчивости

9

(𝑨𝒚𝑘, 𝒚𝑘+1) = (𝒚𝑘, 𝑨𝒚𝑘+1) = (𝑨𝒚𝑘+1, 𝒚𝑘).

Поэтому

(𝑨(𝒚𝑘 + 𝒚𝑘+1), 𝒚𝑘+1 − 𝒚𝑘) = (𝑨𝒚𝑘+1, 𝒚𝑘+1) − (𝑨𝒚𝑘, 𝒚𝑘).

А если оператор 𝐀 – положительно определённый, пространство 

может быть наделено нормой ‖𝐳‖A

2 = (𝐀𝐳, 𝐳)

(𝑨(𝒚𝑘 + 𝒚𝑘+1), 𝒚𝑘+1 − 𝒚𝑘) = ‖𝒚𝑘+1‖𝐴

2 − ‖𝒚𝑘‖𝐴

2.

Таким образом,

(𝑩(𝒚𝑘+1 − 𝒚𝑘), 𝒚𝑘+1 − 𝒚𝑘) −

1

2 𝜏(𝑨(𝒚𝑘+1 − 𝒚𝑘), 𝒚𝑘+1 − 𝒚𝑘) +

+

1

2 𝜏‖𝒚𝑘+1‖𝐴

2 − − 

1

2 𝜏‖𝒚𝑘‖𝐴

2 = 0.

Наконец, если для любого 𝐳 выполняется (𝐁𝐳, 𝐳) −

1

2 τ(𝐀𝐳, 𝐳) ≥ 0, то 

1

2 𝜏‖𝒚𝑘+1‖𝐴

2 −

1

2 𝜏‖𝒚𝑘‖𝐴

2 ≤ 0  

или 

‖𝒚𝑘+1‖𝐴

2 ≤ ‖𝒚𝑘‖𝐴

2, 

что применительно к доказательству устойчивости означает ненарастание 
возмущений с ростом k.

Схема доказательства альтернативной теоремы Самарского. 
Системы обыкновенных дифференциальных уравнений, в которых 

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

𝐵(𝑦𝑘+1 − 𝑦𝑘) + 𝜏𝐴𝑦𝑘 = 𝐵(𝑦𝑘+1 − 𝑦𝑘) −

1

2 𝜏𝐴𝑦𝑘+1 +

1

2 𝜏𝐴𝑦𝑘+1 + +

1

2 𝜏𝐴𝑦𝑘 +

+ 

1

2 𝜏𝐴𝑦𝑘 = 0  

или, окончательно, в виде

(𝑩 −

1

2 𝜏𝑨) (𝒚𝑘+1 − 𝒚𝑘) +

1

2 𝜏𝑨(𝒚𝑘+1 + 𝒚𝑘) = 𝟎. 

1. Анализ устойчивости методами спектрального анализа

10

Скалярно умножим обе части равенства на (𝐲k + 𝐲k+1):

((𝑩 −

1

2 𝜏𝑨) (𝒚𝑘+1 − 𝒚𝑘), 𝒚𝑘+1 + 𝒚𝑘) +

1

2 𝜏(𝑨(𝒚𝑘+1 + 𝒚𝑘), 𝒚𝑘+1 + 𝒚𝑘) = 0. 

Если теперь матрица 𝐀 – положительно определённая (или хотя бы 

полуопределённая), то

(𝑨(𝒚𝑘+1 + 𝒚𝑘), 𝒚𝑘+1 + 𝒚𝑘) ≥ 0, 

((𝑩 −

1

2 𝜏𝑨) (𝒚𝑘+1 − 𝒚𝑘), 𝒚𝑘+1 + 𝒚𝑘) ≤ 0. 

А если при этом матрица (𝐁 −

1

2 τ𝐀) – самосопряжённая, то

((𝑩 −

1

2 𝜏𝑨) (𝒚𝑘+1 − 𝒚𝑘), 𝒚𝑘+1 + 𝒚𝑘) = ((𝑩 −

1

2 𝜏𝑨) 𝒚𝑘+1, 𝒚𝑘+1) −

− ((𝑩 −

1

2 𝜏𝑨) 𝒚𝑘, 𝒚𝑘) ≤ 0. 

Наконец, если матрица 𝐁 −

1

2 τ𝐀 – положительно определённая 

(строго), то в норме, порождаемой положительно определённой самосо
пряжённой матрицей 𝐁 −

1

2 τ𝐀, выполняется

‖𝒚𝑘+1‖𝑩−1

2𝜏𝑨

2
− ‖𝒚𝑘‖𝑩−1

2𝜏𝑨

2
≤ 0  

или

‖𝒚𝑘+1‖𝑩−1

2𝜏𝑨 ≤ ‖𝒚𝑘‖𝑩−1

2𝜏𝑨, 

что составит основу получения условий устойчивости.

Таким образом, совместное требование самосопряжённости и стро
гой положительной определённости матрицы 𝐀 вместе с требованием не
строгой положительной определённости матрицы 𝐁 −

1

2 τ𝐀 заменилось 

совместным требованием нестрогой положительной определённости матрицы 𝐀 и требованием самосопряжённости и строгой положительной 

определённости матрицы 𝐁 −

1

2 τ𝐀. Но условие, найденное академиком 

А.А. Самарским:

(𝐵𝑧, 𝑧) > 1

2 𝜏(𝐴𝑧, 𝑧) 

1.1. Теоремы устойчивости

11

для двухслойной итерационной схемы оказывается необходимым в любом случае. Заметим, что это и предшествующие условия вместе с положительной определённостью (строгой или нестрогой) матрицы 𝐀 требуют 
положительной определённости (строгой) матрицы 𝐁.

Пример линейной модели. 
Пусть сообщающиеся сосуды представляют собой П-образную 

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

1

2 (𝜌𝑆(𝐿 + 2𝐻))(𝑦′)2 + (𝜌𝑆𝑦)𝑔𝑦 = 𝑐𝑜𝑛𝑠𝑡,  
(1.1) 

где ρ – плотность жидкости; S – площадь сечения трубки; L – длина горизонтального участка; H – высота невозмущенных столбов жидкости в 
двух вертикальных участках; y – отклонение свободных поверхностей от 
положения равновесия; g – ускорение свободного падения. После сокращения получим

1

2 (𝐿 + 2𝐻)(𝑦′)2 + 𝑔𝑦2 = 𝑐𝑜𝑛𝑠𝑡. 

Продифференцируем обе части равенства по времени

(𝐿 + 2𝐻)𝑦′𝑦′′ + 2𝑔𝑦𝑦′ = 0. 

Поскольку величины L ≥ 0, H > 0, g > 0 неотрицательны, то можно 

обойтись одним коэффициентом

(𝑦′′ + 𝜔2𝑦)𝑦′ = 0, 

где ω = √2g (L + 2H)
⁄
.

Считая y′ ≠ 0, окончательно получим

𝑦′′ + 𝜔2𝑦 = 0. 

Уравнение высокого порядка может быть записано в виде системы 

обыкновенных дифференциальных уравнений в нормальной форме:

{𝑦′ = 𝑣,
𝑦(0) = 𝑦0,

𝑣′ = −𝜔2𝑦,
𝑣(0) = 𝑣0. 

Для разрешения системы можно использовать чисто явную схему