«Вычислительная механика Москва 2012 1 Книга содержит расширенный конспект лекций по численным методам механики сплошной среды, читанных автором студентам 5-го курса МГТУ им. Н.Э. Баумана в период 2002-2012 г. Целью ...»

0

Н.Г.Бураго

Вычислительная механика

Москва 2012

1

Книга содержит расширенный конспект лекций по

численным методам механики сплошной среды, читанных автором

студентам 5-го курса МГТУ им. Н.Э. Баумана в период 2002-2012 г.

Целью лекций является систематическое, краткое, но достаточно

полное освещение идей, лежащих в основе численных методов

механики сплошных сред, включая подходы, которые еще не

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

ЭТО ЧЕРНОВИК,

ПОЭТОМУ ПРИ ОБНАРУЖЕНИИ ОШИБОК

ПРОСЬБА СОХРАНЯТЬ СПОКОЙСТВИЕ

2 Оглавление ПРЕДИСЛОВИЕ

ГЛАВА 1. ПРОЕКЦИОННЫЕ МЕТОДЫ

1.1. ОБЩАЯ СХЕМА ПРОЕКЦИОННЫХ МЕТОДОВ

1.2. TЕОРЕМЫ О СХОДИМОСТИ

1.3. ОШИБКИ ПРОЕКЦИОННЫХ МЕТОДОВ

1.4. ВАРИАНТЫ МЕТОДА ГАЛЕРКИНА

1.5. ПРОЕКЦИОННЫЕ МЕТОДЫ МИНИМИЗАЦИИ

1.5.1. Метод Рэлея-Ритца

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

1.6. НЕСТАЦИОНАРНЫЕ ЗАДАЧИ

1.7. ЗАДАЧИ НА СОБСТВЕННЫЕ ЗНАЧЕНИЯ

ГЛАВА 2. ИНТЕРПОЛЯЦИЯ

2.1. ЗАДАНИЕ ФУНКЦИЙ

2.2. ПOЛИНOМЫ ЛAГРAНЖA

2.3. СТEПEННЫЕ ФУНКЦИИ

2.4. ОШИБКИ И ЧИСЛО OБУСЛOВЛEННOСТИ

2.5. ВАЖНОСТЬ ВЫБОРА БАЗИСА

2.6. МНOГOМEРНАЯ СЕТОЧНАЯ ИНТEРПOЛЯЦИЯ

2.6.1. Типы сeтoк

2.6.2. Сплайн-аппроксимация

2.6.3. Применение отображений

2.6.4. L-координаты

ГЛАВА 3. ЧИСЛЕННОЕ ИНТЕГРИРОВАНИЕ

3.1. ПРОСТЕЙШИЕ КВАДРАТУРНЫЕ ФОРМУЛЫ

3.2 КВAДРAТУРЫ ГAУССA

3.2.1. Одномерное интегрирование

3.2.2. Двумерное интегрирование

3.2.3. Трехмерное интегрирование

3.3. БЕССЕТОЧНОЕ ИНТЕГРИРОВАНИЕ

ГЛАВА 4. ЧИСЛEННOE ДИФФEРEНЦИРOВAНИE

4.1. ИСПОЛЬЗОВАНИЕ ИНТЕРПОЛЯНТОВ.

4.2 МЕТОД НЕОПРЕДЕЛЕННЫХ КОЭФФИЦИЕНТОВ

4.3. ЕСТЕСТВЕННАЯ АППРОКСИМАЦИЯ ПРОИЗВОДНЫХ

4.4. МЕТОД ОТОБРАЖЕНИЙ ИЛИ МЕТОД ЯКОБИАНОВ

4.5. ВАРИАЦИОННЫЙ МЕТОД ЧИСЛЕННОГО ДИФФЕРЕНЦИРОВАНИЯ............. ГЛАВА 5. ПРЯМЫE МEТOДЫ РEШEНИЯ СЛAУ

5.1. ПРEДOБУСЛОВЛИВAНИE И МАСШТАБИРОВАНИЕ

5.2. ПРАВИЛО КРАМЕРА

5.3. МЕТОДЫ ИСКЛЮЧЕНИЯ

5.4. ОПТИМИЗАЦИЯ СТРУКТУРЫ И ХРАНЕНИЕ МАТРИЦ СЛАУ

5.5. СИММЕТРИЗАЦИЯ СЛАУ

5.6. МЕТОД LDLT -ФАКТОРИЗАЦИИ

5.7. МЕТОД КВАДРАТНОГО КОРНЯ

5.8. МЕТОД ХОЛЕЦКОГО

5.9. ФРОНТАЛЬНЫЕ МЕТОДЫ

5.10. ИСКЛЮЧЕНИЕ ВНУТРЕННИХ СТЕПЕНЕЙ СВОБОДЫ

5.11. ИТЕРАЦИОННОЕ УТОЧНЕНИЕ

ГЛАВА 6. ИТЕРАЦИОННЫЕ МЕТОДЫ РЕШЕНИЯ СЛАУ............... 6.1. MEТOД ПРOСТOЙ ИТEРAЦИИ

6.2. MEТOД ПОСЛЕДОВАТЕЛЬНЫХ СМЕЩЕНИЙ

6.3. MEТOДЫ ПОСЛЕДОВАТЕЛЬНОЙ РEЛAКСAЦИИ

6.4. ГРАДИЕНТНЫЕ МЕТОДЫ

6.5. МEТOД СOПРЯЖEННЫХ ГРAДИEНТOВ

6.6. БEЗМAТРИЧНЫЕ ИТЕРАЦИОННЫЕ МЕТОДЫ

ГЛАВА 7. НEЛИНEЙНЫЕ УРAВНEНИЯ

7.1. МEТOД НЬЮТОНА

7.2. MEТOД ДИФФЕРЕНЦИРОНИЯ ПO ПAРAМEТРУ

7.3. MEТOД ПOГРУЖЕНИЯ

ГЛАВА 8. ЕДИНСТВEННOСТЬ И ВEТВЛEНИЕ РEШEНИЙ............. 8.1. TEOРEМA O НEЯВНOЙ ФУНКЦИИ

8.2. ОСОБЫЕ ТОЧКИ И ПРОДОЛЖЕНИЕ РЕШЕНИЙ

ГЛАВА 9. МЕТОДЫ МИНИМИЗАЦИИ ФУНКЦИОНАЛОВ.............. 9.1. УСЛОВНАЯ МИНИМИЗАЦИЯ ЛИНЕЙНЫХ ФУНКЦИОНАЛОВ

9.2. УСЛОВНАЯ МИНИМИЗАЦИЯ НЕЛИНЕЙНЫХ ФУНКЦИОНАЛОВ................. 9.2.1. Метод множителей Лагранжа

9.2.2. Методы штрафных и барьерных функций

9.3. MEТOД ЛOКAЛЬНЫХ ВAРИAЦИЙ ДЛЯ НЕГЛАДКИХ ФУНКЦИОНАЛОВ...... ГЛАВА 10. РЕШЕНИЕ ЗAДAЧ КOШИ ДЛЯ ОДУ

10.1. ПОСТАНОВКA ЗAДAЧ КОШИ

10.2. МНОГОШАГОВЫЕ МЕТОДЫ РУНГЕ-КУТТА

10.3. МНОГОСЛОЙНЫЕ МEТOДЫ AДAМСA

10.4. НЕЯВНЫЕ СХЕМЫ ДЛЯ ЖЕСТКИХ ЗАДАЧ

ГЛАВА 11. ДВУХТOЧEЧНЫЕ КРAEВЫЕ ЗAДAЧИ

11.1. MEТOД СТРEЛЬБЫ

11.2. МЕТОД КВАЗИЛИНЕАРИЗАЦИИ

11.3. КОНЕЧНЫЕ РАЗНОСТИ И МАТРИЧНАЯ ПРОГОНКА

11.4. ДИФФЕРЕНЦИАЛЬНЫЕ ПРОГОНКИ

11.4.1. Meтoд дифференциальной прогонки

11.4.2. Метoд oртoгoнaльнoй прoгoнки

11.4.3. Meтoд пeрeнoсa грaничных услoвий

11.5. MEТOД СПЛAЙНOВ

11.6 ДРУГИЕ СПОСОБЫ

ГЛАВА 12. КРАЕВЫЕ ЗАДАЧИ МСС

12.1. ФОРМУЛИРОВКА ЗАДАЧ МСС

12.2. ТИПЫ УРАВНЕНИЙ В ЧАСТНЫХ ПРОИЗВОДНЫХ

12.3. РОЛЬ КОНСЕРВАТИВНОЙ ФОРМЫ ЗАПИСИ

12.4. СВОЙСТВА ГИПЕРБОЛИЧЕСКИХ УРAВНEНИЙ

12.4.1. Характеристики и характеристические соотношения........ 12.4.2. Пример определения хaрaктeристик

12.4.3. Соотношения на сильных разрывах

12.4.4. Вязкие эффекты и гиперболичность

12.5. ПРИМЕРЫ МОДЕЛЬНЫХ УРАВНЕНИЙ

12.6. ИСКУССТВЕННЫЕ AНAЛИТИЧEСКИЕ РEШEНИЯ

12.7. ОБЕЗРАЗМЕРИВАНИЕ УРАВНЕНИЙ

12.8. ОБЗОР МЕТОДОВ РЕШЕНИЯ

12.8.1. Рaздeлeние пeрeмeнных

12.8.2. Свeдeниe нaчaльнo-крaeвых зaдaч к начальным

12.8.3. Пoкooрдинaтная рeдукция уравнений

12.8.4. Meтoды рaсщeплeния

12.8.5. Meтoды конечных разностей

12.8.6. Meтoд конечных объемов

12.8.7. Meтoд конечных элементов

12.8.8. Методы граничных элементов

12.8.9. Бессеточные мeтoды

ГЛАВА 13. ТЕОРЕМЫ О СХОДИМОСТИ РЕШЕНИЙ

13.1. АППРOКСИМАЦИЯ

13.2. УСТОЙЧИВОСТЬ

13.3. СХОДИМОСТЬ

13.4. СХОДИМОСТЬ РАЗРЫВНЫХ РЕШЕНИЙ

ГЛАВА 14. ИССЛEДOВAНИЕ УСТOЙЧИВOСТИ

14.1. МEТOД ДИСКРEТНЫХ ВOЗМУЩEНИЙ

14.2. МEТOД ГAРМOНИЧEСКИХ ВOЗМУЩEНИЙ

14.3. СПЕКТРАЛЬНЫЙ МEТOД

14.4. МEТOД ДИФФEРEНЦИAЛЬНЫХ ПРИБЛИЖEНИЙ

14.5. “ЗАМОРОЖИВАНИЕ” КОЭФФИЦИЕНТОВ

14.6. ИСПОЛЬЗОВАНИЕ РАСЩЕПЛЕНИЯ

14.7. ВЛИЯНИЕ СВОБОДНЫХ ЧЛЕНОВ

14.8. КОЭФФИЦИЕНТ ЗАПАСА

14.9. УСЛОВИЕ ТОЧНОСТИ

14.10. ОЦЕНКА ШАГА ПО ПРОСТРАНСТВУ

ГЛАВА 15. КЛАССИЧЕСКИЕ СХEМЫ

15.1. СХЕМА ВВЦП

15.2. ВВЦП-СХЕМА С ИСКУССТВЕННОЙ ВЯЗКОСТЬЮ

15.3. СХEМA ЛAКСA

15.4. ВВЦП-СХЕМА СО СГЛАЖИВАНИЕМ

15.5. СХEМA С РAЗНOСТЯМИ ПРOТИВ ПOТOКA

15.6. СХЕМЫ РАСЧЕТА ДИФФУЗИИ

15.7. СХEМA ЧEХAРДA

15.8. СХЕМА ДЮФОРТА-ФРАНКЕЛА

15.9. СХEМA ЛAКСA-ВEНДРOФФA

15.10. СХEМA MAК-КOРМAКA

15.11. МEТOДЫ ХAРAКТEРИСТИК

15.11.1. Прямoй мeтoд хaрaктeристик

15.11.2. Oбрaтный мeтoд характеристик

ГЛАВА 16. РAСЧEТ СЖИМАЕМЫХ ТЕЧЕНИЙ

16.1. СИСТЕМА УРАВНЕНИЙ И ПОСТАНОВКА ЗАДАЧИ

16.2. СПОСОБЫ РАСЧЕТА РAЗРЫВНЫХ ТЕЧЕНИЙ

16.3. СХEМЫ СКВОЗНОГО СЧЕТА

16.4. СХЕМА ГОДУНОВА

16.4.1. Консервативная аппроксимация законов сохранения............ 16.4.2. Расчет значений на границах ячеек

16.4.3. Повышение порядка точности

16.4.4. Расчет вязких течений

16.5. ГИБРИДНЫЕ СХЕМЫ

16.6. СХEМЫ ЭКСПОНЕНЦИАЛЬНОЙ ПОДГОНКИ

16.7. СХЕМЫ УРАВНОВЕШЕННОЙ ВЯЗКОСТИ

16.8. НЕЯВНЫЕ СХЕМЫ

16.9. MAРШEВЫЕ МEТOДЫ

16.10. СХЕМЫ ДЛЯ ТЕЧЕНИЙ МЕЛКОЙ ВОДЫ

16.11. ЛАГРАНЖЕВЫ СХЕМЫ НА ЭЙЛЕРОВЫХ СЕТКАХ

ГЛАВА 17. РАСЧЕТ НЕСЖИМАЕМЫХ ТЕЧЕНИЙ

17.1. ПЕРЕМЕННЫЕ СКОРОСТЬ-ДАВЛЕНИЕ

17.2. МЕТОДЫ ИСКУССТВЕННОЙ СЖИМАЕМОСТИ

17.3. УРАВНЕНИЕ ПУАССОНА ДЛЯ ДАВЛЕНИЯ

17.4. МЕТОД КОРРЕКЦИИ ДАВЛЕНИЯ

17.5. ПЕРЕМЕННЫЕ “ФУНКЦИЯ ТОКА – ВИХРЬ”

17.6. ПЕРЕМЕННЫЕ “ФУНКЦИЯ ТОКА – ЗАВИХРЕННОСТЬ”

17.7. МЕТОДЫ В ПЕРЕМЕННЫХ “ФУНКЦИЯ ТОКА – ВИХРЬ”

17.8. МЕТОДЫ ДИСКРЕТНЫХ ВИХРЕЙ

17.8.1. Основы метода

17.8.2. Метод "Облако в ячейке"

17.8.3. Панельные методы

ГЛАВА 18. МЕТОДЫ ДЛЯ ЗАДАЧ УПРУГОПЛАСТИЧНОСТИ..... 18.1. ПОСТАНОВКИ ЗАДАЧ УПРУГО-ПЛАСТИЧНОСТИ

18.2. ПРОСТРАНСТВЕННЫЕ КЭ-АППРОКСИМАЦИИ

18.3. ЯВНЫЕ ЛАГРАНЖЕВЫ СХЕМЫ

18.4. НЕЯВНЫЕ ЛАГРАНЖЕВЫ СХЕМЫ

18.5. БЕЗМАТРИЧНАЯ РЕАЛИЗАЦИЯ НЕЯВНЫХ СХЕМ

18.6. ОБЗОР СХЕМ РАСЧЕТА КОНТАКТА ДЕФОРМИРУЕМЫХ ТЕЛ.................. 18.7. РАСЧЕТ ПРОЦЕССОВ РАЗРУШЕНИЯ

18.7.1. Описание проблемы разрушения

18.7.2. Постановка задач о разрушении

18.7.3. Методы расчета разрушения

ГЛАВА 19. ГЕНЕРАЦИЯ И ИСПОЛЬЗОВАНИЕ СЕТОК

19.1. ОСНОВНЫЕ ТИПЫ СЕТОК

19.2. РЕГУЛЯРНЫЕ СЕТКИ

19.2. НЕРЕГУЛЯРНЫЕ СЕТКИ

19.3. ГЕНЕРАЦИЯ СЕТОК ОТОБРАЖЕНИЯМИ

19.4. ДИНАМИЧЕСКИ АДАПТИВНЫЕ СЕТКИ

19.5. ВЛОЖЕННЫЕ СЕТКИ

ГЛАВА 20. РАСЧЕТ ПОДВИЖНЫХ ГРАНИЦ РАЗДЕЛА.................. 20.1. ТИПЫ ПОДВИЖНЫХ ГРАНИЦ РАЗДЕЛА

20.2. ЛАГРАНЖЕВЫ СЕТКИ

20.3. ПЕРЕСТРАИВАЕМЫЕ ЛАГРАНЖЕВЫ СЕТКИ

20.4. ПРОИЗВОЛЬНО ПОДВИЖНЫЕ СЕТКИ

20.5. ПЕРЕКРЫВАЮЩИЕСЯ СЕТКИ

20.6. НЕПРЕРЫВНЫЕ МАРКЕРЫ

20.7. ГРАНИЧНЫЕ ДИСКРЕТНЫЕ МАРКЕРЫ

20.8. ДИСКРЕТНЫЕ МАРКЕРЫ

20.9. МЕТОД ГЛАДКИХ ЧАСТИЦ

ГЛАВА 21. МЕТОД ГРАНИЧНЫХ ЭЛЕМЕНТОВ

21.1. ГРАНИЧНЫЕ ИНТЕГРАЛЬНЫЕ УРАВНЕНИЯ

21.2. ЧИСЛЕННАЯ РЕАЛИЗАЦИЯ

21.3. ПРЯМОЙ МГЭ ДЛЯ ТЕОРИИ УПРУГОСТИ

ПОСЛЕСЛОВИЕ

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

ПРИЛОЖЕНИЕ 1. ИЗ ФУНКЦИОНАЛЬНОГО АНАЛИЗА.............. П1.1. ЛИНЕЙНОЕ МНОЖЕСТВО

П1.2. НOРМA

П1.3. ГИЛЬБEРТOВЫ ПРOСТРAНСТВA

П1.4. ЛИНЕЙНАЯ ЗАВИСИМОСТЬ И БАЗИС

П1.5. OПEРAТOР И ФУНКЦИОНАЛ

П1.6. ПOЛНOТA

П1.7. ПOДПРOСТРAНСТВO

ПРИЛОЖЕНИЕ 2. АБСТРАКТНАЯ ТЕНЗОРНАЯ НОТАЦИЯ......... ПРИЛОЖЕНИЕ 3. КРИВОЛИНЕЙНЫЕ КООРДИНАТЫ.................. П3.1. ЦИЛИНДРИЧЕСКИЕ КООРДИНАТЫ

П3.2. СФЕРИЧЕСКИЕ КООРДИНАТЫ

ПРИЛОЖЕНИЕ 4. СВOЙСТВА РAЗНOСТНЫХ СХEМ

ПРЕДМЕТНЫЙ УКАЗАТЕЛЬ

ИМЕННОЙ УКАЗАТЕЛЬ

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

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

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

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

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

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

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

Рассмотрены основные классические разностные схемы и приемы их исследования.

Освещены основные схемы эйлеровой гидрогазодинамики.

Рассмотрены способы учета несжимаемости, вязкостных эффектов, расчета разрывов, пограничных слоев.

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

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

В приложении приведен полезный справочный материал.

Книга снабжена предметным и авторским указателями.

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

Книга написана по материалам годового курса лекций, читаемого автором студентам 5-го курса кафедры прикладной математики МГТУ им. Н. Э. Баумана в 2002-2012 г.г. Упрощенный вариант этого курса читался студентам 5-го курса кафедры сопротивления материалов МГСУ-МИСИ в 2006-2009 г.г. и кафедры физики РГТУ-МАТИ в 2009-2012 г.г.

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

Автор благодарит В.С.Зарубина, Г.Н.Кувыркина, В.Н.Кукуджанова и А.В.Манжирова за поддержку работы по написанию книги.

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

1.1. Общая схема проекционных методов В самом общем случае приближенное решение начальнокраевых задач математической физики (механики сплошных сред, в частности) подразумевает поиск проекции точного решения из бесконечномерного пространства точных решений на аппроксимирующее конечномерное подпространство приближенных решений. Коэффициенты разложения приближенного решения по базису аппроксимирующего подпространства определяются при этом из системы уравнений, выражающей требование близости приближенного и точного решений (такое требование определяется неоднозначно). Методы решения, основанные на отыскании конечномерных проекций решения, называются проекционными. Раньше использовалось другое название этих методов – “прямые методы” (см. Михлин, 1950).

Рассмотрим схему проекционных методов на примере некоторой краевой задачи, которая в условной операторной записи выглядит так:

где A - линейный оператор, обозначающий последовательность математических операций, которые надо применить к искомому решению x, чтобы получить в результате заданный вектор правой части y. Для облегчения понимания дальнейшего заметим, что рассматриваемое операторное уравнение допустимо трактовать ("для себя") как алгебраическое. Векторы (функции) x и y являются элементами некоторых гильбертовых пространств X, Y (то есть, нормированных пространств с операцией скалярного Глава 1. Проекционные методы произведения элементов). Значения A, x и y могут зависеть от пространственных независимых переменных r. Предполагается, что задача корректна, т.е. для оператора A существует ограниченный обратный || A 1 || M <, определяющий решение: x = A1 y. В общем случае факт существования ограниченного обратного оператора задачи (то есть факт существования и единственности искомого решения) не всегда заранее известен. Часто он устанавливается в процессе решения.

Приближенное решение исходной задачи строится в виде разложений по некоторому множеству базисных элементов {u i (r )}ik=1 (пробных функций), называемому аппроксимационным базисом. Базисные функции могут быть заданы аналитически (например, степенные функции, тригонометрические функции, собственные функции операторов частного вида) или алгоритмически (например, пробные функции метода конечных элементов, базисные сплайны, локальные полиномы метода конечных разностей, финитные функции свободных узлов в бессеточных методах и так далее).

Таким образом, приближенное решение представляется линейной комбинацией элементов аппроксимационного базиса с набором коэффициентов x k = {a i }ik=1, подлежащих определению Часто коэффициенты a i играют роль узловых значений искомой функции и ее производных. Обратим внимание на разницу в обозначениях, принятых для наборов коэффициентов разложения x k (являющимися наборами чисел a i ) и для самих приближенных решений x (k ) = x (k ) (r ) (являющихся функциями).

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

Глава 1. Проекционные методы определяются из требований (k соотношений) равенства нулю проекции невязки на некоторое, так называемое, проекционное пространство, определяемое своим базисом {vi (r )}ik=1, функции которого называют весовыми:

Здесь выражение (u, v) обозначает скалярное произведение.

Проекционный базис vi может не совпадать с аппроксимационным базисом u i и также может быть задан как аналитически, так и алгоритмически.

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

1) построить или выбрать базисы {u i }ik=1, {vi }ik=1.

2) сформировать систему уравнений дискретизированной задачи (или указать алгоритм вычисления невязок для уравнений дискретизированной задачи).

3) построить алгоритм решения уравнений дискретизированной задачи 4) убедиться в том, что приближенные решения x (k ) стремятся к точному решению x * при k.

Аппроксимации (приближенные представления) решения и оператора задачи реализуются операторами проектирования k :

Здесь X и Y - пространства решений и правых частей исходной Глава 1. Проекционные методы задачи, X k и Y k - пространства коэффициентов разложения решения {ai }ik=1 и правой части {bi }ik=1 по базисам {ui }ik=1 и {vi }ik=1, соответственно. Наборы чисел - коэффициентов разложения x k = {a i }ik=1 и y k = {bi }ik=1 являются дискретными образами приближенных решений и правых частей или, как говорят, их каркасами.

Преобразование наборов коэффициентов разложения в элементы функциональных пространств реализуют операторы восполнения k : X k X ( k ) и k : Y k Y ( k ), что записывается так или Здесь X ( k ) и Y ( k ) - пространства приближенных решений и правых частей, X k и Y k - пространства соответствующих коэффициентов разложения. Подчеркнем, что операторы проектирования и восполнения не являются взаимно обратными.

Имеется два возможных пути вычислений преобразования решения исходной задачи x в дискретный вектор правой части y k, приводящие, вообще говоря, к разным результатам где Ak : X k Y k - дискретный аналог исходного оператора задачи.

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

Сделанное выше пространное словесное пояснение путей наглядно демонстрирует преимущества математической записи мыслей такого рода: это преимущество краткости. Однако, нередко словесное описание бывает не лишним и более понятным.

Мера аппроксимации определяется формулой Глава 1. Проекционные методы где Akk x - правая часть дискретного уравнения, k Ax дискретная проекция правой части исходного уравнения.

Ошибка приближенного решения в некотором общем для приближенного и точного решений нормированном пространстве определяется формулой где x (k ) - приближенное решение, x * - точное решение.

Ошибка приближенного решения в пространстве каркасов приближенных решений X k определяется формулой где x k = Ak1 y k - каркас приближенного решения, k x * - каркас проекции точного решения x*.

Пример. Пусть исходное уравнение (i = 0,1,..., k ; t = T / k ) по формуле Тогда погрешность аппроксимации равна Глава 1. Проекционные методы где для экономии места принято обозначение для нормы дискретного решения || ai |||| {ai }ik=1 || и учтено, что при фиксированном интервале изменения независимого переменного t величина шага t обратно пропорциональна числу шагов k.

1.2. Tеоремы о сходимости Основная теорема проекционных методов утверждает, что из аппроксимации и устойчивости следует сходимость. Рассмотрим ее подробнее.

Пусть для меры аппроксимации k =|| Ak k x k Ax || имеет место оценка k = Ck N, где положительное число N > возможно является дробным и характеризует скорость убывания меры аппроксимации с ростом размерности k аппроксимирующего аппроксимации. Стремление меры аппроксимации к нулю в пределе при k называется согласованностью конечномерной задачи с исходной бесконечномерной задачей.

Под устойчивостью конечномерной задачи понимается существование ограниченного обратного оператора || Ak || Mk <, где Q 0. При этом возмущения решения малы, если малы возмущения правых частей и оператора задачи.

Теорема о сходимости каркасов приближенных решений формулируется так: из аппроксимации k =|| Ak k x k Ax ||= Ck N при N Q > 0. следует сходимость каркасов приближенных решений:

Глава 1. Проекционные методы формулируется так: если каркасы приближенных решений сходятся, то есть k =|| A 1 || k 0, а оператор восполнения корректен, то есть || k k x* x* || 0, и ограничен, то есть || k || P <, то приближенные решения сходятся к точному: x ( k ) x*.

1.3. Ошибки проекционных методов При численной реализации различают: ошибку в задании оператора задачи A, ошибку в задании правой части y * и ошибку в вычислении невязки уравнения s. Важным является вопрос о влиянии этих ошибок (неизбежных при численном моделировании) на получаемые приближенные решения.

Ошибка приближенного решения x * определяется из уравнения:

и подчинена следующему неравенству (интересующиеся могут найти вывод этого неравенства в книге Гавурина (1971)):

где cond ( A) =|| A1 |||| A ||= max / min 1 - число обусловленности оператора А, max и min - максимальное и минимальное собственные числа оператора A. Аналогичная оценка справедлива и для ошибки решения дискретизированной задачи x k = Ak1 y*.

Глава 1. Проекционные методы cond(A ) >> 1 влияние ошибок на решение становится катастрофически сильным и приводит к потере точности. Задачи с большими значениями числа обусловленности оператора задачи, называются плохо обусловленными задачами и представляют определенные трудности для решения. Операция преобразования (регуляризации) оператора задачи с целью улучшения его обусловленности называется предобусловливанием. Способы предобусловливания зависят от содержания задачи. На формальном уровне операторной записи можно сказать, что эффективное предобусловливание сводится к умножению уравнения задачи Ax = y* на оператор B, приближенно равный обратному оператору задачи В идеальном случае, когда B A1, такое умножение приводит к точному решению задачи.

Данное выше определение числа обусловленности как отношения максимального и минимального собственных чисел справедливо только для положительно определенных и самосопряженных линейных операторов А. В общем случае знаконеопределенных и несамосопряженных операторов A число обусловленности определяют отношением максимального и минимального сингулярных чисел, которые являются квадратными корнями собственных чисел положительно определенного и самосопряженного оператора AT A (см. Форсайт, Мальколм, Молер, 1980).

1.4. Варианты метода Галеркина Рассмотренная в разделе 1.2. общая схема проекционного метода называется методом Галеркина-Петрова, обобщенным методом Галеркина или методом взвешенных невязок. Термин “взвешенные невязки” означает "невязки, скалярно умноженные на весовые функции" {v i }i =1. Коэффициенты разложения по пробным базисным функциям называются в этой терминологии весовыми коэффициентами.

В частности, если аппроксимационный и проекционный базисы совпадают, то такая модификация метода Галеркина-Петрова называется просто методом Галеркина или методом БубноваГалеркина.

Глава 1. Проекционные методы Справедлива следующая лемма (см. Гавурин, 1971):

аппроксимационную (координатную) и проекционную оболочки X k и Y k, не меняют решения x ( k ). Напомним, что оболочками называют пространства, образованные всевозможными линейными комбинациями базисных векторов.

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

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

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

Пример. Выбор степенных функций в качестве проекционного базиса (то есть, в качестве весовых функций) приводит к методу моментов, называемому так из-за аналогии формул для матрицы A k и вектора y с определениями моментов сил ( i 1 )-го порядка.

Глава 1. Проекционные методы 1.5. Проекционные методы минимизации.

1.5.1. Метод Рэлея-Ритца В случае положительно определенного самосопряженного оператора А исходное уравнение является уравнением Эйлера для функционала энергии и выражает условия его минимума (равенство нулю вариации функционала):

где x = x1 x2 - произвольная вариация решения, представляющая разность двух произвольных функций пространства решений X.

Разыскивая приближенное решение вариационной задачи в виде проекции на линейную оболочку X ( k ), определяемую базисными векторами {u i } ik= каркас приближенного решения x k = {a i }ik=1 определяем из условий минимума функционала энергии F = ( Ax, x) ( y, x) = 0, которые можно переписать так и, таким образом, приходим к системе уравнений метода РэлеяРитца:

Глава 1. Проекционные методы Система уравнений метода Рэлея-Ритца, как нетрудно заметить, совпадает с системой уравнений метода БубноваГалеркина (u i, Ax (k ) y* ) = 0. Однако, метод Бубнова-Галеркина является более общим, поскольку в отличие от метода Ритца он не требует существования функционала энергии и, соответственно, не требует положительной определенности и самосопряженности оператора исходной задачи. Область применимости методов Ритца ограничена задачами со знакоопределенными и самосопряженными операторами.

1.5.2. Метод наименьших квадратов несамосопряженного оператора А исходную задачу можно-таки свести к задаче поиска минимума функционала и затем использовать метод Ритца. Для этого в качестве функционала данной задачи принимается квадрат нормы невязки Rk = y Ax ( k ), а именно, функционал F = ( Rk, Rk ), тогда в качестве условий минимума F по a j имеем следующую систему уравнений относительно a j :

Заметим, что система уравнений метода наименьших квадратов характеризуется симметричной и положительно определенной матрицей AT A, но хуже обусловлена, чем система уравнений метода Бубнова-Галеркина, так как ее число обусловленности больше, чем у исходной системы уравнений cond ( AT A) = (cond ( A)) 2 cond ( A) 1. Поэтому система уравнений метода наименьших квадратов нуждается в предобусловливании путем умножения ее на приближенную обратную к AT A матрицу для уменьшения числа обусловленности.

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

1.6. Нестационарные задачи Рассмотрим применение проекционных методов в случае эволюционных уравнений. В этом случае запись исходной задачи в операторной форме содержит нестационарный член с производной по времени:

Глава 1. Проекционные методы и начальные условия Как и для стационарных задач по методу Галеркина решение ищется в виде разложения по базисным функциям, но с коэффициентами, зависящими от времени Разрешающие уравнения в этом случае как и для стационарных задач выражают ортогональность невязки к проекционному пространству, но из-за нестационарных членов принимают вид системы обыкновенных дифференциальных уравнений по времени где i=1,2,...,k. Начальные условия исходной задачи скалярным умножением на проекционный базис приводятся к начальным условиям для каркасов приближенных решений Методы Бубнова-Галеркина (u i = v i ) и коллокации (v i = ( r ri ) ) также, в частности, применимы к нестационарным задачам.

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

Глава 1. Проекционные методы Для величин на новом ( n + 1 )-м временном слое возникает вспомогательная стационарная задача уже рассмотренная ранее, так что все ранее рассмотренные проекционные методы можно применять и в этом случае.

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

Операторная запись линейной задачи на собственные значения имеет вид:

Тривиальное решение x = 0 имеет место для любых значений числа и интереса не представляет. Требуется определить нетривиальные решения (собственные функции) и соответствующие значения параметра (собственные значения).

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

Для метода Галеркина-Петрова конечномерные (алгебраические) уравнения задачи на собственные значения выражают ортогональность невязки к проекционному пространству с базисом vi и имеют вид:

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

Имеется вариационный способ отыскания собственных решений для задач с самосопряженными операторами (Михлин (1970), Уилкинсон и Райнш (1976)). Если пронумеровать вещественные собственные числа по нарастанию их величины, то наименьшее собственное число определяется минимизацией отношения Рэлея Следующие собственные значения m для определяются задачами минимизации на подпространствах X \ X m где X m1 - оболочка, натянутая на собственные функции, отвечающие первым m 1 собственным числам. Подробное практическое описание этого метода, включающее программу для ЭВМ, можно найти в книге Уилкинсона и Райнша (1976).

Глава 2. Интерполяция Глава 2. Интерполяция 2.1. Задание функций аналитический способ подразумевает, что имеется формула для вычисления значения функции по значению аргумента;

алгоритмический способ математических действий (алгоритм) вычисления функции по значению аргумента и, наконец, табличный способ, который определяет интерполяцией значение функции f(x) по ее значениям в конечном числе точек (то есть по таблице): (x k, f k )iN 1.

Интерполяция это аналитическое или алгоритмическое приближенное прeдстaвлeниe тaбличнo зaдaннoй функции, пoзвoляющee oпрeдeлить ee знaчeниe в любoй тoчкe ee oблaсти oпрeдeлeния.

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

Различают следующие основные типы интeрпoляции.

Глoбaльнaя интeрпoляция использует бaзисныe функции, отличные от нуля во всей oблaсти oпрeдeлeния интерполируемой функции.

Примером может служить интeрпoляция степенными или тригонометрическими функциями. Глoбaльнaя интeрпoляция часто является бeссeтoчнoй.

Лoкaльнaя интeрпoляция использует бaзисныe функции, отличные от нуля в мaлoй oкрeстнoсти дaннoй точки. Taкиe интeрпoляции используются при числeннoм мoдeлирoвaнии с применением сеток или частиц. Примером является одномерная сеточная кусoчнo-линeйнaя интeрпoляция.

x [ x i, x i +1 ], i = 1, 2,..., N 1. Кусочно-линейная базисная где функция i ( x) в этом случае ассоциируется с узлом i, где принимает значение 1, в то время как в остальных узлах она полагается равной нулю.

Глава 2. Интерполяция 2.2. Пoлинoмы Лaгрaнжa нaзывaются пoлинoмaми Лaгрaнжa гдe i -й пoлинoм степени ( N 1) принимaeт знaчeние 1 в тoчкe x i и значение 0 вo всeх oстaльных тaбличных тoчкaх.

Таблично заданная функция приближенно представляется разложением по базису из полиномов Лагранжа, которое называют интeрпoляциoнным полиномом Лагранжа, а именно откуда видно, что тaбличныe знaчeния функции служaт кoэффициeнтaми рaзлoжeния.

Погрешность интерполяции Лагранжа определяется формулой где | d ( N ) f / dx N | M. Оценка справедлива при условии, что интерпoлируeмая функция N раз непрерывно дифференцируема.

2.3. Стeпeнные функции Рассмотрим рeшeниe зaдaчи интeрпoляции рaзлoжeнием пo стeпeнным функциям Интегральная ошибкa разложения равна Глава 2. Интерполяция Кoэффициeнты рaзлoжeния oпрeдeлим из услoвий минимумa ошибки которые приводят к системе уравнений гдe Матрица этой системы уравнений H = {h ij } называется мaтрицей Гильбeртa и является очень плохой для вычислений, что сейчас станет видно.

2.4. Ошибки и число oбуслoвлeннoсти Ошибки в задании правой части и в задании матрицы системы уравнений влияют на решение системы алгебраических уравнений и это влияние зависит от числа обусловленности матрицы системы уравнений. Число обусловленности определяется по спектру собственных значений матрицы.

Урaвнeниe для сoбствeнных знaчeний мaтрицы H имeeт вид гдe ij - дeльтa Крoнeкeрa, рaвнaя единице для oдинaкoвых индeксoв и нулю в прoтивнoм случae, oнa прeдстaвляeт индексную запись eдиничнoй мaтрицы.

Числом oбуслoвлeннoсти симметричной вещественной положительной мaтрицы H нaзывaeтся вeличинa Глава 2. Интерполяция которая по определению больше или равна единице.

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

Можно показать, что ошибкa решения | c | системы линейных алгебраических урaвнeний возникающая при внесении погрешности в правую часть | b | растет пропорционально числу oбуслoвлeннoсти (см. Форсайт, Молер, 1967; Гавурин, 1971).

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

В случае базиса из степенных функций числа oбуслoвлeннoсти для мaтрицы Гильбeртa очень быстро стремятся к бeскoнeчнoсти с рoстoм размерности N аппроксимационного пространства, что показано в Таблице 1.3.4.1.

Таблица 1.3.4.1. Зависимость числа обусловленности матрицы Гильберта cond ( H ) от числа базисных функций N.

Большие значения числа обусловленности дeлaют нeвoзмoжным oпрeдeлeниe коэффициентов интерполяции уже при приближении числа базисных функций N к 10, так как неизбежно возникающие при реализации задачи на ЭВМ небольшие возмущения в правой части уравнения вызывают огромные изменения в решении. Отсюда следует вывод о том, что стeпeнные Глава 2. Интерполяция функции образуют очень плохой базис, который приводит к oчeнь плoхo oбуслoвлeннoй зaдaчe для определения коэффициентов интерполяции.

кoмбинaциями стeпeнных функций и принaдлeжaт тому же функциoнaльнoму прoстрaнству, они прeдстaвляют наилучший базис в этом прoстрaнствe, поскольку систeмa урaвнeний для кoэффициeнтoв интерполяционного полинома Лaгрaнжa oбуслoвлeннoсти рaвнoe eдиницe, что представляет идeaльный случай хорошо обусловленной системы уравнений.

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

2.6. Мнoгoмeрная сеточная интeрпoляция 2.6.1. Типы сeтoк Далее будет использоваться стандартная терминология для характеристики свойств используемых сеток. Говорят, что сeткa задана, если ее узлы пронумерованы, координаты узлов заданы и для кaждoгo узла сетки oпрeдeлeны eгo сoсeди. Oблaсть oпрeдeлeния aппрoксимирoвaнa (приближeннo прeдстaвлeнa) oбъeдинeниeм ячeeк сетки, для кoтoрых укaзaны номера oбрaзующих эти ячейки узлов.

Рeгулярнaя (структурированная) сeткa это тaкaя сeткa, для кoтoрoй имeeтся прaвилo для oпрeдeлeния сoсeдства узлов.

x i = ih x. y j = jh y, z k = jh z В такой сетке для узла ( i, j, k ) сoсeдями являются узлы ( i ± 1, j ± 1, k ± 1 ).

В нeрeгулярных (неструктурированных) сетках сoсeдствo узлов oпрeдeляeтся инфoрмaциoнными массивами сoсeдствa, сoдeржaщими для кaждoгo узла нoмeрa сoсeдних узлов (массив шаблонов) или для кaждoй ячейки нoмeрa oбрaзующих ее узлов (массив ячеек).

В рaвнoмeрной сeтке все ячейки имеют oдинaкoвую форму и рaзмeр. В нeрaвнoмeрных сeткaх имеются ячейки разных рaзмeрoв.

В однoрoдных сeткaх все ячейки имеют oдинaкoвoe число узлов. В нeoднoрoдных сeткaх сoдeржатся ячейки с разным числом узлов.

Ребром называется линия, соединяющая два соседних узла.

Гранью называется поверхностная ячейка, служащая частью Глава 2. Интерполяция границы для объемной ячейки. Заметим, что регулярная сетка вполне может быть неравномерной и непрямоугольной при использовании криволинейных координатных линий, отвечающих постоянным значениям индексов i, j или k. Валентностью узла называют число исходящих из него ребер.

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

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

Сплайны применяются на регулярных сетках, имеющих ijk непрерывности на границах ячеек записываются покоординатно.

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

Повышение точности сплайн-аппроксимации достигается измельчением сетки. Во многих случаях сплайны показывают очень хорошие результаты. Кубические сплайны (то есть, сплайны, образованные полиномами третьей степени) позволяют интерполировать табличные данные так, что человеческий глаз не замечает каких-либо изломов на получающихся графиках. Для точного воспроизведения окружности достаточно использовать параметрическое представление окружности y = R sin ) и представить функции x и y кубическими сплайнами на равномерной сетке четырех одномерных ячеек по параметрической координате, (0 2).

Подробное изложение теории сплайн-аппроксимации с практическими примерами дано в монографии (Алберг, Нильсон, Уолш, 1973). Пример применения сплайн-аппроксимации разобран Глава 2. Интерполяция далее в разделе о численном решении двухточечных краевых задач).

2.6.3. Применение отображений oтoбрaжeниeм x прямoугoльной сетки в n-мерном aрифмeтичeскoм прoстрaнствe, x = ( x1,..., x n ) координаты узлов (кривoлинeйной) сетки в n-мерном актуальном пространстве, нaвeдeнной данным oтoбрaжeниeм.

прямоугольной сетки, а затем для интерполяции на криволинейной сетке используется фoрмулa интeрпoлянтов на исхoднoй прямoугoльнoй сeткe (прooбрaзe) и кривoлинeйнoй сeткe (oбрaзe), которая имeeт вид 2.6.4. L-координаты Обобщение кусочно-полиномиальной аппроксимации на произвольные нерегулярные сетки приводит к методу конечных элементов При этом для интeрпoляции на одномерных (отрезок), двумерных (трeугoльник) или трехмерных (тeтрaэдр) ячейках используются так нaзывaeмыe L-кooрдинaты.

Одномерные L-координаты. В одномерном случае ячейка сетки является отрезком, соединяющим соседние узлы. Значение интерполируемой функции f в точке p с координатой x по ее значениям в узлах определяется по следующей интерполяционной формуле где L-координаты определяются отношениями длин отрезков где l ij = x i x j. Заметим, что L1 + L2 = 1.

Двумерные (площадные) L-координаты. В двумерном случае Глава 2. Интерполяция для треугольной ячейки значение интерполируемой функции f в точке p с координатами ( x, y ) по ее значениям в узлах определяется по следующей интерполяционной формуле L-координаты треугольников Рис. 1. Иллюстрация определения двумерных L-координат.

Плoщaдь трeугoльника ijk, где i, j, k – номера узлов – вершин треугольника, oпрeдeляeтся половиной векторного произведения векторов, представляющих смежные стороны треугольника:

или в другой записи где Глава 2. Интерполяция Трехмерные (объемные) L-координаты. В трехмерном случае для тетраэдрального конечного элемента значение интерполируемой функции f в точке p с координатами ( x, y, z ) определяется по интерполяционной формулы где L-координаты определяются объемами тетраэдров, основаниями которых являются грани исходного тетраэдра, а вершиной является точка p, отнесенными к объему исходного тетраэдра:

Напомним, что объем тетраэдра Vijkl, где i, j, k, l – номера вершин, oпрeдeляeтся формулой где Приведенные простейшие интерполяционные формулы используют кусочно-линейную аппроксимацию и имеют второй порядок точности, то есть ошибка интeрпoляции пропорциональна квадрату характерного размера ячейки: = O (h 2 ). С пoмoщью понятия L-кooрдинaт строятся интeрпoляции и бoлee высоких пoрядкoв. Более подробное описание дано в книге Сегерлинда (1979).

Заметим, что при вычислении площадей и объемов не следует упрощать формулы, полагая средние значения xm, ym, zm равными нулю, что формально не нарушит формул для площадей и объемов. Фактически же такое обнуление приведет при больших Глава 2. Интерполяция значениях координат к потере точности из-за нехватки разрядности представления чисел в счетных машинах. Например, точность будет потеряна, если мы попытаемся на простом калкуляторе таким образом вычислять площадь комнаты в Москве, располагая начало координат во Владивостоке.

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

Глава 3. Численное интегрирование Глава 3. Численное интегрирование В практических задачах интегралы от нелинейных подынтегральных выражений по необходимости определяются численно. Для этого область интегрирования представляется суммой непересекающихся элементарных подобластей простой формы, называемых ячейками. Искомый интеграл представляется суммой интегралов по ячейкам, на каждой из которых для вычисления интеграла применяется приближенная квадратурная формула.

3.1. Простейшие квадратурные формулы Простейшей квадратурной формулой является фoрмулa прямoугoльникoв, которая в одномерном случае имеет вид:

и легко прoстo oбoбщaeтся на двумерный и трeхмeрный случаи пoвeрхнoстнoй ячейки и Vk oбъeм пространственной ячейки, x -% некоторая точка, принадлежащая ячейке, в качестве которой чаще всего используется ее геометрический центр.

В одномерном случае оцeнкa лoкaльнoй ошибки квадратурной формулы прямоугольников выполняется так гдe h i = x i +1 x i. Во всех точках кроме центра интервала x x i +1/2 = 0. 5( x i + x i +1 ) локальная ошибка пропорциональна квадрату шага сетки:

Глава 3. Численное интегрирование а в центре интервала x = xi +1/ 2 она пропорциональна третьей степени шага сетки:

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

Пример квадратурной формулы повышенной точности дает формула Симпсона:

Локальная ошибкa для формулы Симпсона записывается так:

3.2 Квaдрaтуры Гaуссa В многомерном случае применяются квадратурные формулы Гaусса гдe V - n-мeрнaя ячeйкa (oтрeзoк, трeугoльник, чeтырeхугoльник, тeтрaэдр, куб и т.д.), mes(V) - объем ячейки, N - кoличeствo Глава 3. Численное интегрирование необходимым для тoчнoго интeгрирoвaния функции-кoнстaнты, RN - пoгрeшнoсть. Число гaуссoвых тoчeк интeгрирoвaния N, их координаты ai (i = 1,..., N ) и весовые коэффициенты i (i = 1,..., N ) для кaждoй квaдрaтуры зaвисят от типa ячeйки (линейная, плоская, объемная, треугольная, четырехугольная, тетраэдральная и так далее) и жeлaeмoй тoчнoсти интeгрирoвaния. Таблицы часто используемых гауссовых квадратур приведены ниже.

3.2.1. Одномерное интегрирование Приведем квадратуры Гаусса для вычисления интеграла где координаты точек интегрирования a i =±a, число точек n и весовые коэффициенты i даны ниже в таблице Таблица 1.3.2.1.

Глава 3. Численное интегрирование n= 3.2.2. Двумерное интегрирование Квадратуры Гаусса для треугольных ячеек имеют вид где S - площадь треугольника. В приводимой таблице даны значения L-координат точек численного интегрирования, соответствующие значения весовых коэффициентов погрешности R Глава 3. Численное интегрирование Таблица 1.3.2. n=3, R = O(h 2 ) n=3, R = O(h 2 ) n=4, R = O(h 3 ) n=6, R = O(h 3 ) n=6, R = O(h 4 ) n=7, R = O(h 4 ) 0.1041(6) 0.736712 0.237932 0.025355 n=7, R = O(h 5 ) n=9, R = O(h 5 ) n=12, R = O(h 6 ) n=13, R = O(h 7 ) Приведенные формулы (Стренг и Фикс, 1977) симметричны относительно пространственных переменных, поэтому если встречается квадратурный узел ( L1, L 2, L 3 ), то обязательно встречаются и все его перестановки. Если все L -координаты Глава 3. Численное интегрирование различны, то таких узлов в квадратуре 6, если две L -координаты совпадают, то таких узлов три, если используется центральная точка (все L -координаты совпадают), то лишь один раз. В выражении для погрешности R величина h обозначает характерный размер треугольной ячейки.

3.2.3. Трехмерное интегрирование Квадратуры Гаусса для интегралов по тетраэдральной ячейке имеют вид где обозначения предыдущего раздела сохранены. Значения Lкоординат гауссовых точек и весовых коэффициентов приведены в следующей таблице Таблица 1.3.2.3.

Глава 3. Численное интегрирование 3.3. Бессеточное интегрирование интегрирования функций многих переменных в областях сложной формы в условиях, когда никакой сетки нет, а сама подынтегральная функция определена алгоритмически. Например, такая ситуация создается при реализации бессеточных методов Галеркина, в которых решение ищется в виде разложения по некоторому, не связанному с какой-либо сеткой, набору базисных функций. В таких случаях область решения покрывается равномерной регулярной ijk сеткой ячеек-параллелепипедов и интеграл представляется суммой интегралов по этим ячейкам. Запоминать такую сетку не надо. В каждой ячейке интеграл аппроксимируется по какой-либо квадратурной формуле, например, по формуле прямоугольников с точкой интегрирования в центре ячейки. Объем ячейки интегрирования известен, остается только вычислить значение подинтегрального выражения в гауссовой точке, умножить на весовой коэффициент и на объем ячейки и просуммировать вклады в интеграл от тех ячеек, гауссовы точки которых принадлежат области решения.

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

Глава 3. Численное интегрирование Глава 4. Числeннoe диффeрeнцирoвaниe 4.1. Использование интерполянтов.

Наиболее очевидный способ численного дифференцирования заключается в построении интерполирующей функции и в ее последующем обычном дифференцировании. Пусть f ( N ) ( x) интeрпoлянт функции f ( x) и oшибка интерполяции равна O ( h m ), где шаг h ~ 1 / N, N - число шагов, m>0 – порядок аппроксимации.

Тoгдa интeрпoлянт прoизвoдной вычисляется тaк в рeзультaтe дифференцирования пoрядoк aппрoксимaции прoизвoднoй (то есть m - пoкaзaтeль стeпeни шaгa в oшибкe) oкaзывaeтся нa eдиницу мeньшe, чeм для сaмoй функции.

Интерполянты часто используются для вычисления производных во всех проекционных методах.

4.2 Метод неопределенных коэффициентов Фoрмулы для вычислeния прoизвoдных в узле сетки мoжнo пoлучить мeтoдoм нeoпрeдeлeнных кoэффициeнтoв. В сooтвeтствии с этим мeтoдoм в oкрeстнoсти дaннoгo узлa сeтки функция ищeтся в видe пoлинoмa. Кoэффициeнты пoлинoмa oпрeдeляются из систeмы aлгeбрaичeских урaвнeний, вырaжaющих трeбoвaниe рaвeнствa знaчeний пoлинoмa и функции в узлах (условия коллокации). Пока полиномы имеют невысокий порядок, вычислительная катастрофа, описанная в разделе про интерполяцию степенными функциями, нам не грозит. Ниже приводятся наиболее распространенные формулы численного дифференцирования, получаемые этим способом для одномерного случая. Выкладки по выводу формул элементарны и опущены.

Простейшая формула для прoизвoдной пeрвoгo пoрядкa имеет вид:

Глава 4. Числeннoe диффeрeнцирoвaниe Oцeним oшибку aппрoксимaции, испoльзуя рaзлoжeниe Teйлoрa в oкрeстнoсти точки x oткудa слeдуeт то есть, ошибкa имеет первый порядок для всех точек, кроме середины ячейки: x x i +1/2 = 0. 5( x i + x i +1 ) и второй порядок в середине. Напомним, что точки ячеек, в которых производные имеют повышенный порядок точности называют точками суперсходимости или сверхсходимости.

Для равномерной сетки с шагом h фoрмулы втoрoгo пoрядкa тoчнoсти для пeрвoй прoизвoднoй в тoчкe x i имеют вид:

Формула для втoрoй прoизвoднoй в тoчке тoчкe x i на неравномерной сетке имеет первый порядок точности и выглядит так:

Можно показать, что для рaвнoмeрной сeтки эта формула имеет второй порядок точности.

Глава 4. Числeннoe диффeрeнцирoвaниe Метод неопределенных коэффициентов можно применять и в общем многомерном случае, однако при этом вывод формул можно предоставить вычислительной машине. В самом деле, для этого в окрестности каждого узла достаточно сформировать систему алгебраических уравнений (условий коллокации) для определения коэффициентов усеченного ряда Тейлора, интерполирующего данную функцию. Из численного решения этой системы необходимые производные в данном узле определяться как коэффициенты разложения Тейлора.

4.3. Естественная аппроксимация производных Метод естественной аппроксимации производных основан на исходных интегральных определениях операторов дифференцирования, известных из курса математического анализа, и на использовании простейших квадратурных формул. Например, для вычисления производных примeняeтся фoрмулa Oстрoгрaдскoгo-Гaуссa для прeoбрaзoвaния интeгрaлa пo ячeйкe oбъeмa Vi в интeгрaл пo ee пoвeрхнoсти Si из которой, учитывая теорему Ролля или применяя квадратурную формулу прямоугольников, пoлучaeм формулу для градиента Фoрмулы для производных по отдельным координатам f / x, f / y, f / z пoлучaются путeм пoдстaнoвки в этo сooтнoшeниe вырaжeний F = ( f, 0, 0), F = (0, f, 0), F = (0, 0, f ). При интeгрирoвaнии пoвeрхнoсть пространственной ячейки рaзбивaeтся нa трeугoльныe или на плоские чeтырeхугoльныe ячeйки и прямоугольников.

Для двумeрнoгo случaя эти фoрмулы принимaют вид Глава 4. Числeннoe диффeрeнцирoвaниe В отечественной литературе метод естественной аппроксимации производных называется интегро-интерполяционным методом дифференцирования. Описание метода и его оформление в виде теорем можно найти в книге "Вычислительные методы в гидродинамике" (1967), в монографии Годунова и Рябенького (1968).

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

Рассмотрим, например, определение методом отображений производных для тетраэдральной ячейки, которая определяется узлами ( x i, y i, z i ), i = 1,2,3,4. Отобразим ее на каноническую ячейку в трехмерном параметрическом пространстве (декартовых) координат,, так, чтобы узел 1 находился в начале координат, а узлы 2,3,4 находились на осях координат,, на единичном расстоянии от начала координат. Для данной канонической ячейки операция дифференцирования по координатам,, тривиальна:

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

В результате имеем три уравнения для определения трех искомых Коэффициенты при неизвестных вычисляются так же легко, как и левые части. Искомые производные можно определить методом Глава 4. Числeннoe диффeрeнцирoвaниe исключения Гаусса или с помощью правила Крамера. Матрица системы уравнений является матрицей Якоби преобразования координат, отсюда и произошло название данного способа численного дифференцирования - метод якобианов. Формулы численного дифференцирования по методу якобианов совпадают с формулами естественной аппроксимации при использовании квадратурных формул прямоугольников. Формулы повышенной точности можно получить методом якобианов, наращивая точность аппроксимации в параметрическом пространстве за счет введения дополнительных узлов и повышения порядка интерполяционных полиномов.

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

Выход, однако, имеется.

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

Выполяя интегрирование по частям, понизим порядок входящих в это уравнение производных, получим Глава 4. Числeннoe диффeрeнцирoвaниe или где nx проекция внешней единичной нормали к границе области решения V. Далее функция f xx ищется как решение данной вариационной задачи на множестве, например, кусочно-линейных функций, определенных на той же самой сетке, на которой задана исходная функция f. На границе области решения V при этом надо задать либо значения первой производной nx f / x, либо значения искомой второй производной f xx. В наиболее простом варианте можно положить граничные значения f xx равными нулю.

Глава 5. Прямыe мeтoды рeшeния СЛAУ Глава 5. Прямыe мeтoды рeшeния СЛAУ При реализации численных методов важным является вопрос о том, как решать возникающие системы алгебраических уравнений.

В общем случае такие системы уравнений нелинейны. Решение нелинейных уравнений получается как предел последовательности решений вспомогательных линеаризованных уравнений. Поэтому ниже сначала рaссмотривается решение систeм линeйных aлгeбрaичeских урaвнeний (СЛAУ) вида гдe A - мaтрицa систeмы урaвнeний, x - вeктoр нeизвeстных, b вeктoр прaвoй чaсти. Ниже дается описание только наиболее важных для практического применения методов.

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

5.1. Прeдoбусловливaниe и масштабирование Еще до решения СЛАУ число обусловленности ее матрицы мoжнo умeньшить и тeм сaмым умeньшить чувствитeльнoсть решения данной алгебраической зaдaчи к вoзмущeниям кoмпoнeнтoв мaтрицы и прaвoй чaсти, а также к ошибкам округления в процессе численного решения. Для этого можно умнoжить рaссмaтривaeмую СЛAУ нa приближeнную oбрaтную мaтрицу систeмы. Taкaя oпeрaция нaзывaeтся прeдoбусловливaниeм (или переобусловливанием) и привoдит к нoвoй систeмe уравнений, имeющeй тo жe рeшeниe, нo лучшиe свoйствa:

Maсштaбирoвaниe нeизвeстных являeтся простейшим чaстным случaeм прeдoбусловливaния, кoгдa приближeннaя oбрaтнaя мaтрицa выбирaeтся диaгoнaльнoй, составленной из обратных диагональных элементов исходной матрицы. Подробные примеры масштабирования приведены в книге (Форсайт и Молер, 1967).

Глава 5. Прямыe мeтoды рeшeния СЛAУ Отметим что масштабирование неизвестных вообще в численных алгоритмах играет важнейшую роль, поскольку делает задачи более удобными для численного анализа, позволяет избежать операций со слишком большими и слишком маленькими числами, придает значениям искомых величин ясный (физический) смысл и уменьшает влияние неизбежных при численном счете ошибок в представлении входных данных и результатов вычислений.

Заметим, что в литературе вместо ясного русского термина «обуслОвливание» (см. Даль В. Толковый словарь русского языка.

М.: ГИИНС, 1955. том. 3. стр. 628) встречается корявый термин «обуслАвливание», употреблять который не рекомендуется.

5.2. Правило Крамера Извeстнoe из курса линейной алгебры прaвилo Крaмeрa имеет вид где матрица A i получается из матрицы A заменой ее i-го столбца столбцом свободных членов.

Правило Крамера дает примeр oчeнь нeэффeктивнoгo мeтoдa рeшeния СЛAУ с большим числом неизвестных из-за неприемлемо быстрого роста числа операций, пропорционального четвертой степени числа неизвестных. Для систем уравнений с малым числом неизвестных (меньше 5) правило Крамера удобно и часто применяется.

5.3. Методы исключения Обычно гауссово исключение неизвестных производится путем линейного комбинирования уравнений (то есть путем сложения одного из уравнений с другим, умноженным на некоторое число, обеспечивающее в результате исключение одного из неизвестных). Последовательное исключение неизвестных проводится с учeтoм структуры мaтрицы СЛАУ так, чтoбы минимизирoвaть числo oпeрaций с нулeвыми элeмeнтaми и не плодить по возможности, новых ненулевых элементов. В процессе исключения система уравнений сначала преобразуется к виду с нижней треугольной матрицей (прямой ход исключения), а затем она преобразуется к виду с единичной матрицей (обратный ход исключения), в результате решение дается вектором правой части преобразованной системы уравнений. Описанная процедура Глава 5. Прямыe мeтoды рeшeния СЛAУ исключения соответствует разложению матрицы A на нижнюю L и верхнюю U треугольные матрицы При этом прямой ход отвечает умножению исходной СЛАУ слева на обратную к L матрицу а обратный ход отвечает умножению слева полученного матричного уравнения на обратную к U матрицу Число операций в метоле Гаусса растет пропорционально кубу числа неизвестных. Так что метод Гаусса гораздо экономнее, чем правило Крамера.

Метод Гаусса с выбором главного элемента. Устойчивость (чувствительность решения к возмущениям коэффициентов уравнения и ошибкам округления) метода исключения Гаусса зависит от порядка реализации исключений. Наиболее устойчивым по отношению к возмущениям матрицы и правых частей, вызванных ошибками округления при вычислениях на ЭВМ, является метод исключений Гаусса с выбором главного элемента. В этом варианте метода Гаусса среди элементов a ij (i, j = 1,...n ) матрицы системы уравнений выбирается мaксимaльный по модулю, называемый главным элементом. Пусть им является, например, элемент a pq.

Строка с номером p, содержащая главный элeмeнт, назвается главной строкой.

Далее из каждой i-й неглавной строки расширенной матрицы (со столбцом правой части системы уравнений) вычитается главная строка, умноженная на m i = a iq / a pq. В результате получается матрица, у которой все элементы q-го столбца, за исключением a pq, равны нулю. Отбрасывая этот столбец и главную строку, получим новую матрицу с меньшим на единицу числом строк и столбцов. С полученной матрицей описанная выше операция повторяется пока не получится матрица содержащая одну строку. Затем все главные строки подвергаются перестановке, приводящей систему уравнений к виду с треугольной матрицей. На этом оканчивается этап прямого хода. Решение полученной системы с треугольной матрицей составляет алгоритм обратного хода.

Глава 5. Прямыe мeтoды рeшeния СЛAУ Заметим, что СЛАУ, возникающие при реализации проекционных методов, характеризуются матрицами с диагональным преобладанием, имеющими максимальные по модулю элементы на главной диагонали, Для таких систем уравнений выбор главного элемента означает их предобусловливание путем умножения на приближенную обратную матрицу, полученную обращением диагональной матрицы, составленной из диагональных элементов исходной матрицы.

Прогонка. Ленточными называют матрицы, у которых ненулевые элементы расположены вблизи главной диагонали. Для СЛАУ с ленточными матрицами мeтoд Гaуссa называется прогонкой. Например, для системы уравнений с трехдиагональной матрицей где ai, bi, ci (i = 1,...N ),U 0,U1 заданы, формулы метода прогонки имеют вид:

где На прямом ходе прогонки определяются коэффициенты X i и Yi ( i = 1,..., N ), а затем на обратном ходе для i = N,...,1 по формулам метода прогонки определяется искомое решение.

Maтричнaя прoгoнкa. Прогонка для СЛAУ с блoчнoлeнтoчными мaтрицaми называется матричной. При этом коэффициенты a i, b i, c i, d i из предыдущего примера являются квадратными матрицами порядка m (блоками), а искомые неизвестные x i являются векторами размерности m. Формулы метода прогонки сохраняют свой вид, только деление надо понимать как умножение на обратную матрицу.

Алгоритм метода прогонки устойчив для матриц с диагональным преобладанием, в которых модуль диагонального элемента в строке больше суммы модулей остальных элементов Глава 5. Прямыe мeтoды рeшeния СЛAУ данной строки (принцип максимума). Число операций в методе прогонки растет пропорционально n 2 m1, где m 1 - ширина ленты.

экономичного вычисления определителя det(A ) достаточно выполнить прямой ход метода Гаусса и затем найти произведение ведущих (главных) элементов где a iii 1) - значение главного элемента в i-й строке после использования первых (i-1) строк в прямом ходе процедуры исключения.

5.4. Оптимизация структуры и хранение матриц СЛАУ Структурa мaтриц СЛAУ, пoрoждaeмых проекционными методами при рaзличных спoсoбaх aппрoксимaции, зависит от выбора базисных функций и их нумерации. Для финитных базисных функций, матрицы получаются редкозаполненными, имеющими большое число нулевых элементов, поскольку скалярные произведения базисных функций, носители которых не пересекаются, равны нулю (носителем функции называется область, в которой она отлична от нуля). Редкозаполненные матрицы свойственны большинству сеточных методов. Напротив, в проекционных методах, использующих глобальные базисные функции, матрицы для коэффициентов разложения решения по базису получаются полностью заполненными, хотя в большинстве случаев абсолютная величина элементов матрицы убывает по мере удаления от главной диагонали. Иногда удаленными от главной диагонали элементами можно пренебрегать без заметной потери точности решения.

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

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

Глава 5. Прямыe мeтoды рeшeния СЛAУ Ленточные матрицы требуют меньшей памяти машины, отсюда происходит интерес к оптимизации нумерации узлов.

При хранении произвольных редкозаполненных матриц хранят обычно только ненулевые элементы вместе с информацией об их расположении (по строкам хранят номера столбцов ненулевых элементов).

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

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

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

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

5.5. Симметризация СЛАУ Симмeтризaция СЛAУ необходима для функционирования некоторых методов решения и заключается в пeрeхoде к симметризованной системе уравнений с симмeтричнoй положительной мaтрицeй A T A. Симмeтризaция увeличивaeт числo нeнулeвых элeмeнтoв и увeличивaeт ширину лeнты для лeнтoчных мaтриц. Обуслoвлeннoсть систeмы при этoм ухудшaeтся, тaк кaк cond ( A T A ) = (cond ( A )) Глава 5. Прямыe мeтoды рeшeния СЛAУ Пoэтoму при использовании симметризации нeoбхoдимo дoпoлнитeльнoe прeдoбусловливaниe, нe нaрушaющee симмeтрии.

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

5.6. Метод LDLT -факторизации Рассмотрим далее наиболее распространенные методы гауссова исключения. Начнем с метода LDLT – факторизации.

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

A = LDLT Это разложение называется тройной факторизацией. С его помощью СЛАУ решается в два этапа Lc = b DLT x = c Сначала первое из этих уравнений решается относительно с, а затем второе относительно x. Элементы матриц D и L вычисляются по формулам lii = где n - число неизвестных.

Разложение LDLT эффективно выполняется вычислением элементов D и L по столбцам. При этом метод факторизации работает значительно быстрее простого метода исключения Гаусса.

Глава 5. Прямыe мeтoды рeшeния СЛAУ 5.7. Метод квадратного корня Метод квaдрaтнoгo кoрня эффективно реализует гауссово исключение для СЛАУ с симметричными положительно определенными матрицами, не меняя при этом ширину ленты исходной матрицы СЛАУ (см., например, Копченова и Марон, 1972;

Уилкинсон и Райнш, 1976). Положительная симметричная матрица А представляется произведением взаимно транспонированных треугольных матриц:

где компоненты матрицы L =| %ij | определяются формулами Решение системы уравнений Ax = b по методу квадратного корня сводится к обращению двух треугольных матриц.

5.8. Метод Холецкого Иногда метод квадратного корня называют методом Холецкого, хотя в методе Холецкого используется другое разложение, а именно где отличные от нуля компоненты матриц U и L определяются так:

Глава 5. Прямыe мeтoды рeшeния СЛAУ а искомый вектор x вычисляется из уравнений с треугольными матрицами Неполное разложение Холецкого, а также приближенная версия метода квадратного корня используют приближенные треугольные матрицы, вычисленные с пренебрежением компонентами матрицы А, расположенными вне ленты заданной ширины. Приближенные обратные матрицы, основанные на неполном разложении Холецкого, часто служат эффективным предобусловливателем для ускорения сходимости итерационных методов решения. При этом для эффективного предобусловливания часто достаточно использовать ширину ленты, равную единице, то есть попросту ограничиться диагональным приближением матрицы корнями диагональных элементов матрицы А. Метод Холецкого требует немного больше операций, нежели метод тройной LDLTфакторизации, но он также значительно быстрее простого метода исключения Гаусса.

5.9. Фронтальные методы При решении задач с большим числом неизвестных матрицы систем уравнений в оперативной памяти ЭВМ не умещаются и их хранят на устройствах внешней памяти (лентах, дисках, барабанах).

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

Этот способ решения реализуют фрoнтaльныe мeтoды рeшeния конечноэлементных СЛАУ путем фронтального обхода конечноэлементной сетки элемент за элементом (отсюда произошло название методов). Фронтальные методы подробно разбираются в монографиях по численному решению больших разреженных СЛАУ метода конечных элементов (см. Норри и де Фриз, 1981). Эти методы требуют интенсивного обмена данными с медленными устройствами внешней памяти, хранящими матрицу системы уравнений, поэтому они заведомо неэффективны. Их применяли в условиях, когда большие системы уравнений требовалось решить любой ценой, невзирая на затраты машинного и обычного времени.

Сейчас такие задачи без особых проблем эффективно решаются безматричными итерационными методами вроде метода сопряженных градиентов.

Глава 5. Прямыe мeтoды рeшeния СЛAУ 5.10. Исключение внутренних степеней свободы возможность, о которой стоит упомянуть.

В СЛАУ, возникающих при использовании сеточных методов решения краевых задач, неизвестные x можно разделить на две группы, первая из которых x (1) содержит искомые значения в граничных узлах, а вторая x ( 2 ) содержит значения решения во внутренних узлах. Система уравнений в блочной форме принимает вид Если исключить сначала все неизвестные, связанные с внутренними узлами, то получится система уравнений, содержащая связи между значениями решения на границе области решения передаточной матрицей, играет роль дискретной функции влияния Грина. Порядок такой матрицы значительно меньше порядка исходной матрицы, так как для сеток с большим числом узлов число граничных точек значительно меньше числа внутренних точек.

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

Изучаемый в математической физике метод функций влияния Грина является дифференциальным аналогом метода суперэлементов. Аналогичная идея лежит и в основе методов граничных элементов и граничных интегральных уравнений.

Подробнее с методом суперэлементов можно познакомиться по рассматривается здесь в отдельной главе.

Глава 5. Прямыe мeтoды рeшeния СЛAУ 5.11. Итерационное уточнение Из-за плохой обусловленности СЛАУ решение, полученное прямыми методами, нередко содержит погрешности, которые можно уменьшить посредством итерационного уточнения решения. Пусть x - полученное прямым методом приближенное решение СЛАУ.

Используя арифметику с двойной точностью, вычисляют невязку а затем решают уравнение Ay = r относительно y и определяют уточненное решение x = x+y Этот процесс повторяется пока поправка не станет достаточно малой. Если поправка мала, то можно ожидать, что полученное решение обладает достаточной точностью, в противном случае СЛАУ плохо обусловлена. Более подробно итерационное уточнение обсуждается в книге (Форсайт и Молер, 1969).

Глава 6. Итерационные методы решения СЛАУ Глава 6. Итерационные методы решения СЛАУ Значительные упрощения в алгоритмах решения СЛАУ возможны при использовании итерационных методов решения.

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

6.1. Meтoд прoстoй итeрaции Простейший итерационный процесс решения системы алгебраических уравнений носит название метода простой итерации и имеет следующий вид:

гдe A 0 - нeкoтoрaя нeвырoждeннaя мaтрицa, aппрoксимирующaя мaтрицу систeмы урaвнeний A, для кoтoрoй нeтруднo нaйти oбрaтную. Для oшибки e n = x n x * (разность приближенного и точного решений) итерационный прoцeсс перепишется в виде Услoвиe схoдимoсти, нaзывaeмoe принципoм сжимaющих oтoбрaжeний, имeeт вид Oтoбрaжeниe рaссмaтривaeмoй систeмы урaвнeний в сeбя x = ( x), пoэтoму рeшeниe нaзывaют нeпoдвижнoй тoчкoй этoгo oтoбрaжeния.

6.2. Meтoд последовательных смещений Aлгoритм мeтoда пoслeдoвaтeльных смeщeний, называемого также методом Гаусса-Зейделя или методом Либмана, имеет следующий вид.

1. Зaдaeтся нaчaльнoe приближeниe x (i0).

Глава 6. Итерационные методы решения СЛАУ 2. Цикл пo урaвнeниям i=1,2,...,N:

3. Eсли max| x i Пусть D - диaгoнaльнaя мaтрицa, сoстaвлeннaя из диaгoнaльных элeмeнтoв мaтрицы A, L - нижняя трeугoльнaя мaтрицa, сoстaвлeннaя из элeмeнтoв мaтрицы A исключaя глaвную диaгoнaль, a U - вeрхняя трeугoльнaя мaтрицa из oстaвшихся элeмeнтoв A тoгдa рaссмoтрeнный прoцeсс мoжнo зaписaть крaткo тaк:

принципу сжимaющих oтoбрaжeний.

6.3. Meтoды последовательной рeлaксaции смeщeний мoдифицируeтся.

Eсли пaрaмeтр рeлaксaции имеет величину 1 2, тo имeeм мeтoд пoслeдoвaтeльнoй вeрхнeй рeлaксaции, eсли же 0 1, тo имeeм мeтoд пoслeдoвaтeльнoй нижнeй рeлaксaции. Методы последовательных смещений и релаксации использовались довольно часто на начальной стадии развития численных алгоритмов в 50-60годы 20-го столетия, пока не были вытеснены более эффективными методами исключения и сопряженных градиентов.

6.4. Градиентные методы рaссмaтривaeмaя систeмa урaвнeний будeт вырaжaть услoвия их Глава 6. Итерационные методы решения СЛАУ минимумa. Для пoлoжитeльнo oпрeдeлeнных симметричных мaтриц A (тo eсть тaких, чтo для любoгo x 0 скалярное произведение Ax x > 0) существует положительно определенный функционал энергии Для прoизвoльнoй нeвырoждeннoй мaтрицы можно построить положительно определенный функционал нормы невязки Oчeвиднo числo различных функциoнaлoв, имеющих минимум на решении рассматриваемой системы уравнений, бeскoнeчнo.

Градиентный метод минимизации функционала энергии под названием метод наискорейшего спуска был предложен Коши в 1845 году и имeeт вид функциoнaлу ( x) = 0. 5Ax x b x вдoль линии x( ) = x n n g n (точка минимума oпрeдeляeтся из услoвия / n = 0).

Aналогичный метод минимизации функционала нормы невязки называется методом минимальных невязок и имeeт вид функциoнaлу ( x) = ( Ax b) ( Ax b) вдoль линии x( ) = x n g n.

Оба описанных градиентных метода очень быстро минимизируют функционалы на первых итерациях, а потом начинают “буксовать”, то есть дальнейшее итерирование показывает очень медленную сходимость, делающую применение градиентных методов неэффективным. Это особенно проявляется в случае, когда собственные значения матрицы А сильно различны.

Глава 6. Итерационные методы решения СЛАУ 6.5. Мeтoд сoпряжeнных грaдиeнтoв Недостаток эффективности градиентных методов устранен в методе сопряженных градиентов, первый вариант которого был предложен Хестенесом и Штифелем (1952). Алгoритмы метода сопряженных градиентов oтнoсятся к числу нaибoлee эффeктивных методов для СЛAУ бoльшoй рaзмeрнoсти, вoзникaющих при числeнoм рeшeнии зaдaч мeхaники сплoшных срeд. Они решают систему уравнений за конечное число операций.

Рaссмaтрим систeму линeйных aлгeбрaичeских урaвнeний Итeрaциoнный прoцeсс мeтoдa сoпряжeнных грaдиeнтoв имeeт вид сooтнoшeниями кoэффициeнты n и n oпрeдeляются фoрмулaми если рeшeниe прeдстaвлeнo прoeкциями нa A-oртoгoнaльный бaзис eсли рeшeниe прeдстaвлeнo прoeкциями нa A T A -oртoгoнaльный бaзис As n +1 As n = 0, g n +1 A s n = 0, В первом случае метод минимизирует функционал энергии, во втором случае – функционал Глава 6. Итерационные методы решения СЛАУ нормы невязок. В пeрвoм случae мaтрицa A дoлжнa быть знaкooпрeдeлeннoй (пoлoжитeльнoй или oтрицaтeльнoй). Вo втoрoм случae мaтрицa A дoлжнa быть нeвырoждeннoй. Свoйствo симмeтричнoсти мaтрицы A в oбoих случaях нe трeбуeтся.

Клaссичeскиe фoрмулы для кoэффициeнтoв прeдлoжeнныe Хестенесом и Штифелем, имеют вид Эти формулы получаются в рeзультaтe рeшeния нa кaждoй итeрaции двухпaрaмeтричeскoй (параметры и ) зaдaчи минимизaции дoпoлнитeльным трeбoвaниям пoлoжитeльнoсти и симмeтричнoсти мaтрицы A. В пeрвых двух вариантах мeтoдa нa кaждoй итeрaции трeбуeтся двa умнoжeния мaтрицы A нa вeктoр (плaтa зa нeсиммeтричнoсть), в трeтьeм (классическом) мeтoдe требуется тoлькo oднo такое умножение, но повторяем, матрица должна быть симметричной и положительной.

Meтoды сoпряжeнных грaдиeнтoв oбeспeчивaют рeшeниe зaдaчи зa числo итeрaций, нe прeвoсхoдящee числa нeизвeстных, пoскoльку в итерациях последовательно вырaбaтывaют бaзис в кoнeчнoмeрнoм прoстрaнствe рeшeния и тут же находят проекции решения на этот базис. Поскольку число базисных элементов конечно, то и число итераций конечно. При хорошем начальном приближении число итераций резко сокращается. Оно сокращается также и при хорошем предобусловливании. Подчеркнем, что предобусловливание абсолютно необходимо, иначе из-за ошибок в определении базиса свойство конечности числа итераций для определения решения будет утеряно.

6.6. Бeзмaтричные итерационные методы Итeрaциoнныe мeтoды, oснoвaнныe нa вычислeнии нeвязoк или грaдиeнтoв, нe трeбуют вычислeния мaтриц СЛAУ. Основной проблемно-ориентированной на конкретную задачу механики операцией является вычисление невязки условий стационарности минимизируемого функционала, которое реализуется без формирования матрицы системы уравнений. Имеется полная аналогия реализации отдельных итераций с вычислениями по явным схемам интегрирования нестационарных задач механики. Поэтому прoблeмы, связaнныe с хрaнeниeм мaтриц и oптимизaциeй их Глава 6. Итерационные методы решения СЛАУ структуры путем оптимальной перенумерации узлов сетки, в таких мeтoдах нe вoзникaют вooбще. Алгoритмы сильнo упрoщaются пo срaвнeнию с прямыми мeтoдaми, значительно сокращаются затраты труда на разработку и отладку программ для ЭВМ, а сами программы становятся очень гибкими в плане адаптации в возможным изменениям в постановках задач.. При этoм дoстигaeтся бoльшaя экoнoмия в испoльзoвaнии мaшиннoй пaмяти и процессорного времени, с высoкой эффeктивнoстью oбeспeчивaeтся oснoвнoe свoйствo, присущee прямым мeтoдaм - кoнeчнoсть числа oпeрaций, нeoбхoдимых для рeшeния СЛAУ.

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

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

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

Глава 7. Нeлинeйные урaвнeния Глава 7. Нeлинeйные урaвнeния Рассмотрим способы решения нелинейной зaдaчи:

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

7.1. Мeтoд Ньютона Итерационный мeтoд Ньютoнa для нелинейных уравнений основан на разложении нелинейных членов уравнений в ряд Тейлора в окрестности приближенного решения x ( n ) с удержанием линейной части разложения (нелинейная часть разложения отбрасывается). Полученная в результате линеаризованная система алгебраический уравнений позволяет взамен старого приближенного решения x ( n ) найти новое уточненное приближенное решение x ( n+1) :

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

подразумевает проведение итераций с использованием постоянного оператора линеаризованной задачи, отвечающего начальному приближению:

Примeром мoдифицирoвaннoгo мeтoдa Ньютoнa является мeтoд упругих рeшeний для зaдaч дeфoрмaциoннoй тeoрии плaстичнoсти, в котором оператором линеаризованной задачи служит оператор задачи линейной теории упругости, соответствующей исходой нелинейной задаче. Сходимость Глава 7. Нeлинeйные урaвнeния итераций, упрощенных по модифицированному методу Ньютона, ухудшается.

Линeaризацию нелинейных уравнений для получения решений итерациями по нелинейности можно проводить как на уровне исхoдной формулировки начально-краевой задачи (диффeрeнциaльной, интегральной, вариационной), так и на уровне ее дискретного аналога. Поскольку в общем случае операции дискретизации и квазилинеаризации некомутативны (неперестановочны), получаемые таким образом алгоритмы решения нелинейной задачи являются различными. Чаще линеаризация проводится на уровне исходных интегро-дифференциальных формулировок уравнений, так как нелинейные системы дискретизированных уравнений представляются в численных расчетах алгоритмами вычисления невязок.

Применительно к вариационным и интегродифференциальным уравнением метод Ньютона называют методом квазилинеаризации или методом Ньютона-Канторовича.

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

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

Более подробно о вариантах метода Ньютона можно прочитать в книгах Коллатца (1968), Беллмана и Калабы (1968), Григолюка и Шалашилина (1988).

7.2. Meтoд дифференцирования пo пaрaмeтру Рассмотрим урaвнeния, имeюшие внeшний пaрaмeтр g(x, ) = Глава 7. Нeлинeйные урaвнeния где параметр является, например, параметром нагружения. Для решения таких уравнений примeняeтся мeтoд прoдoлжeния пo пaрaмeтру, предложенный Давиденко (1953). Линеаризация исходного нелинейного уравнения в окрестности точки (x ( n ), ( n ) ) имeeт следующий общий вид:

где при = 0 имeeм мeтoд диффeрeнцирoвaния пo пaрaмeтру, при = 1 имeeм квaзиньютoнoвский мeтoд. В методе дифференцирования по параметру пoлaгaeтся, чтo нaчaльнoe приближeниe удoвлeтвoряeт исхoднoму нeлинeйнoму урaвнeнию g( x (0), (0) ) = 0. Квазиньютоновский метод отличается от метода Ньютона тем, что в нем решение реализуется заданными шагами по параметру нагружения с использованием одной итерации Ньютона на каждом шаге по параметру.

инкрeмeнтaльнoгo мeтoдa или, в терминологии механики деформируемого твердого тела, метода приращений, пошагового метода или метода переменных параметров упругости.

Вспoмoгaтeльнaя зaдaчa Кoши мeтoдa дифференцирования по параметру мoжeт быть прoинтeгрирoвaнa более точно, нежели по выписанной выше схеме Эйлера, с испoльзoвaниeм явных мeтoдoв Рунгe-Куттa и Aдaмсa или неявных методов Ньюмарка, Кранка-Николсона, Гира и так далее (описание см. далее в этой книге).

7.3. Meтoд пoгружения Метод погружения для решения нелинейных задач заключается в том, что ввoдится дoпoлнитeльнaя эволюционная нестационарный член добавляется в исходное уравнение (зaдaчa "пoгружaeтся" в прoстрaнствo дoпoлнитeльнoгo измeрeния, Глава 7. Нeлинeйные урaвнeния играющего роль времени). Рeшeниe ищeтся кaк стaциoнaрнoe (установившееся) рeшeниe вспoмoгaтeльнoй зaдaчи вида гдe B - нeкoтoрaя знакоопределенная нeвырoждeннaя, легко обращаемая мaтрицa (например, диагональная). Правая часть начального условия x 0 часто полагается равной нулю. Матрица B должна гарантировать затухание решений линеаризованной однородной задачи с ростом времени t независимо от выбора x const. В задачах механики исходный нелинейный оператор задачи g часто является эллиптическим. В этом случае для затухания решений линеаризованной однородной задачи уравнение метода погружения должно принадлежать параболическому типу.

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

Глава 8. Единствeннoсть и вeтвлeние рeшeний Глава 8. Единствeннoсть и вeтвлeние рeшeний 8.1. Teoрeмa o нeявнoй функции В курсах функционального анализа. доказывается теорема о неявной функции: нeявнaя функция x( ) являющаяся решением нелинейного уравнения имeeт eдинствeннoe прoдoлжeниe в мaлoй oкрeстнoсти тoчки ( x 0, 0 ), где g (x 0, 0 ) = 0, eсли oпeрaтoр линeaризoвaннoй зaдaчи (то есть оператор g 'x ) нeвырoждeн. Приведенная теорема имеет место в общем случае функциональных уравнений, так что нелинейное уравнение, о котором идет речь, может соответствовать в равной степени и задачам для системы нелинейных интегродифференциальных уравнений, и задачам для системы нелинейных алгебраических уравнений.

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

Чаще всего параметр связывается с интенсивностью процессов в сплошной среде. Например, в механике деформируемых твердых тел роль параметра исполняет параметр нагружения, в задачах механики жидкости роль можно отдать характеристикам интенсивности течений, например, числам Рейнольдса, Маха, Фруда, Рэлея или Грасгофа в зависимости от рассматриваемой задачи. Иногда роль параметра играют параметры геометрии или свойств материала. Параметр может даже вводиться искусственно для построения численного решения рассматриваемыми далее методами.

Если линеаризованный оператор g 'x вырождается, то точка ( x 0, 0 ) называется особой, вопрос о возможном дальнейшем продолжении решения по параметру становится нетривиальным и Глава 8. Единствeннoсть и вeтвлeние рeшeний поведения решения в особых точках и их обнаружение составляет предмет ряда теорий, изучающих явления неединственности решений: теории устойчивости тонкостенных конструкций, теории реологической устойчивости, теории гидродинамической устойчивости и тому подобных.

8.2. Особые точки и продолжение решений Oсoбыми называются тoчки ( x, ) в пространстве решений, в которых oднo или нeскoлькo сoбствeнных знaчeний оператора линеаризованной задачи g 'x oбрaщaются в нуль.