Федеральное агентство по образованию
Государственное образовательное учреждение высшего
профессионального образования
«МАТИ» - Российский государственный технологический
университет им, К.Э. Циолковского
Кукуджанов В.Н.
Численные методы в механике
сплошных сред
Курс лекций
Учебное пособие
Москва 2006

УДК 519.6(075)
ББК 22.193
Рецензенты:
У.Г. Пирумов -д.ф.-м.н., чл.-корр. РАН, зав. лабораторией
Вычислительной математики и программирования Московского
авиационного института им. С. Орджоникидзе
В.И. Кондауров - д.ф.-м.н., проф. Московского физико-технического
института / Государственного университета
Кукуджанов В.Н.
Численные методы в механике сплошных сред; Учебное пособие. - М.:
«МАТИ» - РГТУ, 2006. - 158 с.
Данное пособие основано на курсе лекций, который в течении ряда лет
автор читает студентам «МАТИ» - РГТУ им. К.Э. Циолковского и студентам
Московского физико-технического института/Государственного университета,
специализирующихся в области физики и механики сплошных сред. В ней
изложены основы вычислительных методов и их применение для решения
задач термомеханики сплошных сред.
УДК 519.6(075)
ББК 22.193
© В.Н. Кукуджанов, 2006
© МАТИ, 2006

Оглавление
Предисловие 5
Обозначения 8
1 Понятия о разностных схемах 9
§1 Конечно-разностная аппроксимация дифференциального
оператора 9
§2 Устойчивость и сходимость решения разностной схемы . . 17
§3 Численное интегрирование задачи Коши для системы
дифференциальных уравнений первого порядка 21
§4 Решение задачи Коши для жестких систем ОДУ 29
§5 Примеры составления разностных схем для уравнений в
частных производных 37
§6 Исследование устойчивости разностных схем 44
§7 Упражнения 58
2 Методы решения систем алгебраических уравнений 63
§8 Норма и число обусловленности матрицы 63
§9 Прямые методы решения систем линейных уравнений . . 68
§10 Итерационные методы решения линейных алгебраических
систем уравнений 73
§11 Решение нелинейных уравнений 86
§12 Метод Ньютона и его модификации 92
§13 Методы минимизации функций (методы спуска) 100
§14 Упражнения 106
3

4 Оглавление
3 Методы решения краевых задач для систем ОДУ 111
§15 Численное решение двухточечных краевых задач. Устой^
чивые и неустойчивые алгоритмы 111
§16 Краевая задача для системы линейных уравнений .... 115
§17 Решение краевых задач для нелинейных уравнений .... 117
§ 18 Прогонка для обыкновенных дифференциальных уравнений 119
§19 Решение краевых задач для эллиптических уравнений . . 126
§20 Жесткие краевые задачи 139
§21 Упражнения 147
Литература 153
Предметный указатель 155

Предисловие
5
Предисловие
Книга написана на основе курсов лекций, которые автор читает в
течении ряда лет студентам Московского физико-технического
института /государственного университета/ (http://www.mfti.ru/) и студентам
"МАТИ" Российского государственного технологического
университета имени К.Э.Циолковского (http://www.mati.ru/), специализирующихся
в области физики и механики сплошных сред. В ней изложены основы
вычислительных методов и их применение для решения задач
термомеханики сплошных сред.
Вычислительная механика (computational mechanics) включает в
себя механику сплошных сред (МСС) и методы решения
вычислительной математики; предмет ее исследований — моделирование различных
термомеханических процессов, происходящих в сплошных средах. Она
отличается от классической механики сплошных сред, оперирующей
аналитическими методами, сложностью постановок задач и широтой охвата
решаемых проблем, что стало возможным благодаря огромным
возможностям вычислительных методов и компьютерной техники. В то же
время вычислительная механика наука менее строгая и законченная, чем
каждая из ее составляющих (вычислительная математика и МСС).
Большинство методов решения актуальных нелинейных задач современной
вычислительной механики пока далеки от математически строгого
обоснования, и являются в значительной степени, скорее
экспериментальными, чем теоретически обоснованными методами. Вместе с тем для
большого класса задач существует вполне законченная линейная теория,
которая при умелом использовании оказывается весьма полезной и для
решения сложных нелинейных задач. Поэтому вполне естественно, что
в настоящем курсе изучению основ линейной теории вычислительных
методов уделяется серьезное внимание.
Стиль изложения, принятый в книге, рассчитан больше на
инженеров-физиков, нежели на математиков-вычислителей. Мы не
стремились при изложении методов к безупречной строгости доказательств.

6
Предисловие
Доказательства, если и приводятся, то даны мелким шрифтом, или
указывается ссылка на литературу. Строгому изложению теорем с
детальной формулировкой всех необходимых условий предпочитается
изложение идей и качественных соображений, которые иллюстрируются на
конкретных примерах, связанных с важными приложениями механики.
Вместе с тем, достаточно подробно оговариваются те условия и ограничения,
которые существенны для практического применения излагаемого
метода, недостаточное внимание к которым может привести к ошибочным
результатам. Такой стиль непедантичного изложения вполне оправдан
тем, что учебник предназначен для физиков, которых больше интересует
содержательная сторона рассматриваемых методов, нежели формально
строгий подход их обоснования.
В пособии излагается часть курса, относящаяся к теории
разностных методов и их приложению к механике сплошных сред. Другая часть
курса, относящаяся к вариационным методам и методу конечных
элементов в настоящее пособие не вошла. Она изложена в [9]. Кроме того,
основные приложения вычислительных методов к задачам МСС
изложены в учебном пособии автора "Численное моделирование задач механики
деформируемых тел".
Часть курса, относящаяся к разностным схемам, по содержанию
близка к материалу традиционных курсов по численным методам, но
акцент сделан на методах эффективных для решения задач МСС. К
нетрадиционным методам решения задач МСС относятся методы
решения жестких и сингулярно-возмущенных краевых задач, нелинейных
волновых нестационарных задач, исследование устойчивости с помощью
дифференциальных приближений.
В книге основное внимание уделяется механике твердого
деформируемого тела. Это связано с одной стороны с научными интересами
автора, а с другой, с тем обстоятельством, что на русском языке
учебников и монографий, посвященных этой области механики значительно
меньше, чем вычислительной механике жидкости и газа.
Текст книги разбит на главы. Параграфы имеют сквозную нуме-

Предисловие
7
рацию. Нумерация формул и рисунков двойная — номер параграфа и
номер формулы в параграфе.
В основном тексте мелким шрифтом отмечена дополнительная ин^
формация (доказательства, нюансы применения методов и т.д.), которая
может быть пропущена при первом прочтении книги.
Приводится список литературы, на которую имеются ссылки в
тексте книги. Одни ссылки указывают на дополнительную учебную
литературу и приведены для читателя, более глубоко интересующегося
изложенным материалом. Другие ссылки приведены информации, кем были
получены те или другие результаты.
Книга снабжена предметным указателем, содержащим названия
методов и терминов вычислительной математики и механики.
Книга рассчитана на студентов четвертого-пятого курсов механико-
математических и физико-технических факультетов университетов,
знакомых с основами механики сплошной среды в рамках семестрового
курса и с понятиями вычислительной математики в объеме общего курса
математического анализа.
Автор благодарен своим коллегам и ученикам, принимавшим уча^
стие в обсуждении материала. Большую помощь оказал А.Л.Левитин.

Обозначения
Обозначения
В книге используются обозначения, общепринятые в математике]
вычислительной математике, в физике, в механике. Приведем некоторые
наиболее употребительные обозначения.
Полужирный - А или о: матрица
или тензор соответственно с
элементами Aij или Gij.
Стрелка - и или х: вектор с
элементами Щ ИЛИ Xi.
Крышка й: Фурье-образ скаляра
или вектора.
Нижний индекс
— xf элемент вектора х\
— aij или Dijki: элемент
матрицы/тензора <т, D;
— fx или f'x: частная производная
функции / по переменной ж;
— fg или /1: частная производная
функции / по вектору х.
Верхний индекс
— х2 или хп: степень;
— Рп(х): полином п-ой степени;
— йт или (г*1, щ)т или Ат:
операция транспонирования вектора
или матрицы;
— и* или А*: комплексное
сопряжение вектора или матрицы.
Фигурные скобки {х}: множество.
Квадратные скобки
— [а, Ь]: интервал;
— [и] = и+ — и~: разрыв
величины.
Круглые скобки - (и, /): скалярное
произведение векторов/матриц.
Точка - и • f или Н • и:
произведение векторов/матриц.
Точка сверху — и или (uv)':
материальная производная по времени.
Обозначения в разностных схемах
/г, Ах или Дг/: шаг разностной схемы по пространству;
т или At: шаг разностной схемы по времени;
Ар: шаг по параметру нагрузки;
щ: значение и в узле к разностной сетки;
Uij или Uij: значение и в узле (i,j) разностной сетки;
ип: значение и на гг-ом временном слое.

Глава 1
Понятия о разностных схемах
§1. Конечно-разностная аппроксимация
дифференциального оператора
Конечно-разностная аппроксимация (10). Оценка погрешности аппроксимации (12).
Экстраполяционная формула Ричардсона (15).
Для решения задач механики деформируемых тел вычислительны^
ми методами необходимо произвести дискретизацию тела на элементы,
с тем чтобы свести ее к решению системы алгебраических уравнений.
Исторически с самого зарождения механики был заложен
континуальный подход к решению ее задач. Тела рассматривались как непрерывный
континуум материальных точек и задачи формулировались как задачи
анализа непрерывных функций. Основным инструментом их
исследования было дифференциальное и интегральное исчисление. За несколько
столетий был создан мощнейший аппарат математического
исследования задач физики и механики методами решения дифференциальных
уравнений. Дискретный же анализ до появления вычислительных
машин практически развивался очень слабо. Положение изменилось после
появления ЭВМ и в настоящее время он бурно развивается]
Можно различить два направления этого развития: 1) прямое
физическое моделирование; 2) математическое моделирование. В первом
подходе сплошные тела рассматриваются как дискретная совокупность
материальных частиц, к которым непосредственно применяются законы
9

10
Гл.1 Понятия о разностных схемах
физики и выводятся дискретные уравнения, минуя формулировку в
виде математической задачи для функций непрерывного аргумента.
Однако в современной вычислительной механике основным является второй
подход к дискретизации. При этом подходе производится
дискретизация уже сформулированной континуальной математической задачи. Это
позволяет наиболее полно использовать те достижения
математического анализа, которые были получены в течение столетий. Именно этот
подход является основным и используется в дальнейшем.
Конечно-разностная аппроксимация. Для того чтобы
аппроксимировать задачу, сформулированную в виде операторного уравнения
функции непрерывного аргумента и
L(u) = f
дискретными уравнениями, необходимо
1) заменить область определения
непрерывного изменения аргумента и\
на дискретное множество значений и^.
Например, если оператор определен в
области 6J, показанной на рис. 1.1, то
ее можно заменить множеством
узлов квадратной сетки о;^,
покрывающей эту область. Граница области Г
заменяется ломаной Fh-
2) Ввести в рассмотрение
функции дискретного аргумента или
сеточные функции Uh, определенные на мно^
Рис. 1.1. жестве uh.
3) Заменить дифференциальный оператор L(u) на
разностный Lh(uh), определенный на дискретном множестве Uh
L(u) -> Lh(uh)

§1. Разностная аппроксимация диффер. оператора
11
После этого задача сводится к алгебраической системе уравнений
относительно значений функций щ в точках дискретного множества си^\
Эта общая схема должна быть математически строго
формализована. Для этого вводится понятие сетки и сеточной функции. Сеткой будем
называть множество и^ = {х^, £ со} (г = 1,..., 7V), а функции
непрерывного аргумента U(x) поставим в соответствие с помощью некоторого
оператора P/j дискретное множество значений функции щ = Рн(и), которое
и будет сеточной функцией.
Функции и(х) непрерывного изменения аргумента являются
элементами некоторого функционального пространства Н. Множество
сеточных функций образует некоторое векторное пространство Н^
размерность его совпадает с числом множества точек сетки, а компонентами
вектора являются значения Uh(xi) в этих точках.
В пространстве Н^ вводится норма сеточных функций Цгл^Ця^,
которая является аналогом нормы \\и\\н в пространстве Ну так чтобы
выполнялось условие их согласования
];™|К||яЛ= \\и\\н
Например: 1) норме в пространстве "С" непрерывных функций
\\и\\с = maxxGu; \и(х)\ отвечает норма \\uh\\Ch
\\uh\\Ch = max\u(xi)\
2) норме в L<i соответствует норма L2h
If V/2 г-1 Vl
\\u\\L = I / u2 dx\ -► \\uh\\Lh = i^2v%hi\
3) норме в W% — норма Wih
I x \ 1/2 /N-l k \ Х/2|
\\u\\w = \ dx u2dx\ -► \\uh\\Wh = l^2hk^2u2hi\

12
Гл.1 Понятия о разностных схемах
и т.д. Эти нормы индуцированные соответствующими скалярными
произведениями в пространствах L2j W2 и их сеточными аналогами в
векторных пространствах L2h, W2h-
Оценка погрешности аппроксимации. Основной задачей теории
разностных схем является оценка близости решения конечно-разностной
и дифференциальной задач. Однако эти решения определены в разных
пространствах Н и Hh имеют разные нормы, а оценить разность этих
решений можно только в какой-то общей для них норме. Чтобы обойти эту
сложность, возможны два различных пути. Первый состоит в том, чтобы
решение заданное на множестве Uh в виде сеточной функции щ
доопределить на всем множестве и до функции непрерывного аргумента и(х)
с помощью некоторого оператора интерполирования Я,[щ] —► й(х)\
Задача заключается в том, чтобы по значениям сеточной функ^
ции Uh(xi) восстановить непрерывную функцию й(х). Конечно, эта
функция будет отличаться от функции и(х) и надо оценить норму их разности
\\и(х) — й(х)\\н при х Е и. Такое доопределение функции неоднозначно
и оператор R(uh), с помощью которого осуществляется такое
доопределение является оператором интерполирования. Теория
интерполирования это классическая математическая теория, которая детально
разработана и продолжает развиваться, в частности, благодаря новым
задачам, связанным с применением метода конечных элементов.
Доопределение сеточной функции может быть сделано различны^
ми путями, например, путем линейной или квадратичной интерполяции
и т.д. Тогда нужные оценки можно получить в пространстве Н. Именно
такой подход используется, например, при аппроксимации в методе
конечного элемента (МКЭ). В методе КЭ щ определена во всей области
решения uj в виде кусочнонепрерывных функций. Это позволяет
использовать технику доказательств сходимости, устойчивости и т.д.,
основанную на рассмотрении функций из континуальных множеств.
Второй подход заключается в том, что оператор Ръ,(и) = Щ проек^
тирует и(х) на сетку и получается сеточная функция щ(х{). Заметим,

§1. Разностная аппроксимация диффер. оператора
13
что обратный оператор доопределения Я(щ) = и не может быть
построен по Pfc(it), т.к. они определены в разных пространствах. Рн{и)
действует из Н —► Hh, a R{uh) Hh —► Н.
В численном анализе применяются оба этих метода. Первый применяется в
методе конечных элементов (МКЭ). В МКЭ мы оперируем с функциями из Н,
узловые значения Uh доопределяются до и(х), т.е. строится оператор, который ставит
в соответствие векторному пространству Hh пространство Н. Оценки и
доказательства сходимости проводятся в пространстве Н. Оператору Rh{u) соответствует
построение функций формы фг{х) на выбранном шаблоне, т.е. по узловым значениям
восстанавливаются функции на КЭ. Эти функции образуют базис в Н и к ним
применяются непрерывные операторы интегрирования, дифференцирования и т.д. Эти
вопросы детально будут рассмотрены во второй части книги, где рассматривается
применение МКЭ к решению задач механики сплошных сред.
В теории конечно-разностных уравнений поступают обратным
способом: не доопределяют щ до и(х), а проектируют и(х) на и^ с помощью
оператора Рн{и) —► щ и рассматривают все в пространстве Н^\
Оператор проектирования Ph(u) в простейшем случае, когда мно^
жество точек Xi сетки Uh G а;, есть Ръ,{и) = u(xi).
Оператор проектирования Ph(u) может быть и более
сложным, например, оператором усреднения с весом по окружа- /+1|
ющим узлам, как показано на рис. 1.2, где х* — точка в центре ^^^KNc
правильной шестиугольной сетки внутри области определения *+^|\ k\^l*'
функции и(х) \. \\Х
Ф*) = ^-^ / \
где Si — площадь равностороннего треугольника с вершинами ^^
в точках г, г + 1, ж*.
Теперь можно поставить вопрос об определе- Рис. 1.2.
нии операции проектирования для
дифференциального оператора Ph(L) = Lh(uh), т.е. о его замене разностным. Это может
быть сделано бесконечным числом способов. Например, даже в
простейшем случае первой производной на 3-х точечном шаблоне можно полу-

14
Гл.1 Понятия о разностных схемах
чить семейство разностных операторов, зависящее от параметра
dv
Ph{L) = Lh(vh) = -Ш- = vx разность впере,
Ph(L) = Lh(vh) = г~1 = vx - разность назад
Lha) = avx + (1 - a)vx
Lha) — семейство разностных операторов, зависящее от параметра а
(О ^ а ^ 1), h — шаг сетки.
Также просто можно получить аппроксимационнные формулы для
производных второго
ph(L ) = Lh(vh) =
и более высокого порядков (см. упражнения к главе 2), а следовательно
для любого дифференциального оператора Lv.
Возникает вопрос о погрешности при такой замене, т.е. вопрос о
погрешности аппроксимации.
Погрешностью аппроксимации оператора L разностным
оператором Lh будем называть норму ||^||яЛ сеточной функции ф^
Фн = Ыщ)-Рн(Ь(и))н (1.1)|
где щ = Ph{u)i u(x) — функция непрерывного, a Uh(xi) —
дискретного аргумента. Эта характеристика погрешности аппроксимации во всей
области определения сеточного оператора Ьь(щ).
Если \\ф\\нк = О(hk), то будем говорить, что Lh, аппроксимирует L\
с порядком к. Таким образом, глобальная аппроксимация связана с
понятием нормы, а следовательно, с областью и с ее разбиением, поэтому
она отличается от локальной аппроксимации в окрестности точки.

§1. Разностная аппроксимация диффер. оператора
15
Локальную погрешность аппроксимации ф(х{) в точке Х{ легко
оценить , разлагая vi±\ = v(xi± h) в ряд Тейлора. Например для разности
вперед получим
= v'(xi) + ОЩ
^h{xi) = vx- v\xi) = Lh(uh) - Ph(L(u)) = 0(h)\
Аналогично это делается для любого разностного оператора.
Важно подчеркнуть разницу между локальной погрешностью аппроксимации
и глобальной погрешностью аппроксимации на всей сеточной области.
Покажем, что выбор нормы Ц^Ця^ имеет весьма существенное значение и что в|
разных нормах порядок погрешности аппроксимации может оказаться разным в
зависимости от того, в каком пространстве она оценивается. Это особенно существенно
при рассмотрении нерегулярных сеток [18].
Порядок аппроксимации зависит от шаблона — совокупности точек
сетки, используемого для аппроксимации дифференциального оператора
разностями. Увеличивая число точек шаблона можно повышать порядок
аппроксимации. Это позволяет уменьшить число рассчитываемых точек
при сохранении точности расчета. Однако, такие действия не всегда
оказываются выгодными, т.к. увеличение числа точек шаблона приводит к
усложнению аппроксимационных формул и к увеличению времени
расчета для одной точки сетки. Существует и другая возможность повысить
порядок аппроксимации, предложенная Ричардсоном.
Экстраполяционная формула Ричардсона. Для получения
решения с высоким порядком аппроксимации можно не увеличивать число
используемых точек шаблона, а применить расчет на вложенных сетках,
т.е. последовательно произвести расчет на сетках с шагом h, h/2, h/3
и т.д. С помощью вычисленных решений Uh, Щ/2, Щ/%,... можно
построить экстраполяционную формулу, которая позволяет найти решение щ
более высокого порядка аппроксимации, чем вычисленные решения.
Например, имея решение для двух сеток щ и i%/2, вычисленных по одной
у(х{) + v'(xi) h + v"(xi)— + 0(/i3) - v{xi)

16
Гл.1 Понятия о разностных схемах
и той же симметричной схеме второго порядка точности
uh(x) = и(х) + h2v{x) + 0(/i4)
uh/2(x) = и(х) + (-) v(x) + 0(hA)
составляя линейную комбинацию
1 4
Uh = ~^uh + -uh/2 = u(x) + 0(/i4)
получаем решение с четвертым порядком точности. При использовании
трех сеток для схемы первого порядка аппроксимации можно получить
решение с точностью до 0(h3) по формуле
Uh = auh + buh/2 + cuh/3 = u(x) + 0(h3) (1.2)|
где весовые коэффициенты определяются из системы уравнений]
а+Ъ+с=1
а+¥+У=о
Л2. /1ч2
а+\2 Г+ U С=°
Откуда находим
1 9
а = -, о = —4, с = -
2' ' 2
Во многих случаях такой способ позволяет повысить точность без
существенных затрат, т.к. расчет на двух-трех вложенных сетках, как
правило, выполняется для проверки сходимости используемого метода.
Следует, однако, отметить, что при увеличении числа членов п в
формуле (1.2) весовые коэффициенты быстро растут с уменьшением шага
сетки /i/n, а так как они входят в формулу с переменными знаками, это
приводит к влиянию ошибок округления на окончательный результат.
Чтобы избежать этого, дробление сетки производят в отношении /i/2n,
тогда коэффициенты растут медленнее, но каждый расчет требует
большего времени [12].

§2. Устойчивость и сходимость решения разностной схемы
17
§2. Устойчивость и сходимость решения разностной
схемы
Устойчивость (17). Теорема сходимости Лакса (17). Пример неустойчивой
разностной схемы (18).
Устойчивость. Возникает вопрос: следует ли из условия
аппроксимации (1.1) | l^ull = 0(hk), что решение разностного уравнения также
отличается от точного решения на 0(hk)! Ответ на него отрицательный.
Условие аппроксимации дифференциального оператора
оказывается недостаточным, а лишь необходимым для того, чтобы решение
разностного уравнения Lh(uh) = 0 сходилось бы к решению
соответствующего дифференциального уравнение L(u) = 0 при h —► 0. Для
сходимости необходимо выполнение еще одного условия, а именно, чтобы малые
погрешности, вносимые аппроксимацией в разностное уравнение, не
приводили бы к большим отклонениям в его решении.
Это важное свойство разностной схемы гарантирует устойчивость
разностного уравнения
Lh(uh) = fh
Оно тесно связано с непрерывной зависимостью решения от правых
частей уравнения. Малое возмущение в правой части уравнения Sfh
приводит к малому возмущению решения этого разностного уравнения 5щ.
Определение устойчивости. Разностная схема устойчива, если вы-\
полняется условие
\\Ы\ < C\\5fh\\ (2.1)
для любых Sfh £ Hh, здесь Suh = и — щ, Sfh = f — fh-
Теорема сходимости Лакса. Имеет место следующая теорема,
принадлежащая П.Лаксу: Если некоторый разностный оператор
Lh{uh) = fh аппроксимирует оператор L(u) = f и полученная с его
помощью разностная схема устойчива, то решение щ сходится к и\

18
Гл.1 Понятия о разностных схемах
Дано: 1)
\\Lh(Ph(u))-Ph(L(u))\\Hh = 0(hk), fc>0 (2.2)
Аппроксимация определяется для сеточной функции щ, Ph(L(u)) —
проекция дифференциального оператора на сеточное пространство Н^ L^ —
разностный оператор.
2) \\Suh\\ ^ с||£Д||. Решение щ разностного уравнения Lh(uh) = fh\
устойчиво по правой части.
Требуется доказать, что \\щ — Ph(u)\\ = 0(hk) при h —► 0, где и —\
решение L(u) = f.
Доказательство: Из условия аппроксимации (2.2) получаем
\\Lh (Ph(u)) - Ph (L(u)) || = \\Lh(uh) - Ph(f)\\ =
= \\fh-Ph(f)\\ = \\Sfh\\ = 0(hk) (2.3]
Обозначая ошибку решения разностного уравнения и^ через 8и^
5uh = Ph(u) - uh
Из условия устойчивости (2.1) и соотношения (2.3), находим
\\Suh\\^C\\5fh\\=C.O(hk)
\\uh-Ph(u)\\ = 0(hk)
что и требовалось доказать.
Другими словами ошибка вносимая аппроксимацией в правую часть
разностного уравнения малая порядка 0(hk), тогда ошибка в решении разностного уравнения
в силу условия устойчивости (2.1) будет того же порядка.
Пример неустойчивой разностной схемы. Рассмотрим
устойчивость разностной схемы на простейшем примере обыкновенного
дифференциального уравнения первого порядка
у' + ау = 0 при у = у0 (2.4)
Решение его элементарно у = уо е~ах.
Аппроксимируем (2.4) на трехточечном равномерном шаблоне с шагом h

§2. Устойчивость и сходимость решения разностной схемы
19
следующим семейством разностных операторов, зависящих от
параметра к (0 ^ к ^ 1), являющимся линейной комбинацией конечных
разностей вперед и назад.
, dy yn+1-yn yn-yn_i
2х- 1 - ah х- 1 ^ .\
Уп+1 Уп Н Уп-i = 0 (2.5)
х х
Решение этого разностного уравнения с постоянными коэффициентами
ищем в виде уп = Ci\n. Для определения Л получаем
2 1 — 2х + ah Л 1-х
А
Л' + Л = 0 (2.6)
х х
1,2 = — \-(1-2я + аК)± (l + 2ah(l-2K) + a2H2)l,2\
2х L ч J
Общее решение уравнения (2.6) запишется через две произвольные
константы С\ и С%
Уп = С1\п1+С2\1 (2.7)
Исследуем характер полученных решений при h —► О
Ai = 1 - ah + 0(a2h2) A2 = ^—^(1 + ah) + 0[a2h2
x
Учитывая, что nh = хп и lim^0(l — ah)llh = e~a, находим
Уп = Сг {е~ах- + 0(a2h2)) + С2 (е^ + 0{a2h2)) (^—^ "'* (2.8)1
Для определения Ci и Сг мы имеем только одно условие (2.4). Это
говорит о том, что одно из двух решений (2.7) появляется вследствие
принятой аппроксимации. Это "паразитное" решение.
Оно связано с тем обстоятельством, что формальный порядок
разностного уравнения (2.5) второй. Он определяется числом произвольных
констант в общем решении (2.7). В то же время порядок
дифференциального уравнения (2.4) первый и его общее решение зависит только от
одной произвольной константы, для определения которой имеется одно

20
Гл.1 Понятия о разностных схемах
начальное условие (2.4). Поэтому необходимо найти второе условие, из
которого можно было бы определить вторую произвольную константу в
решении (2.7).
Возможно несколько способов определения этой константы. Сде^
лать это можно с разным порядком аппроксимации. Прежде всего
очевидно, что это условие должно быть задано в соседней с х = 0 точке
сетки х = h1 у = у\(К) и должно мало отличаться от значения у0-
Если положить у\ = ус? то ошибка будет 0(h). Ее можно определить по
двухточечной схеме при к = 1 в (2.5), откуда получим, что
уг = 2/0(1 - ah) + 0(h2)
Ясно, что недостающее начальное условие необходимо аппроксимировать
согласовано с порядком аппроксимации используемой разностной
схемы, иначе порядок точности решения может оказаться сниженным. Эта
ситуация типична для схем, формальный порядок которых выше, чем
порядок аппроксимируемого диференциального уравнения. Такие схемы
часто используются для повышения порядка аппроксимации решения
задачи и здесь необходимо следить за согласованным порядком
аппроксимации для дополнительных начальных условий.
Итак, второе условие для определения С{ следует получить с
помощью двухточечной схемы (2.5) при к = 1
при х = 0, у = у0; при х = h, у1 = (1 - ah) yQ
Подставляя эти условия в (2.7), находим
а=уо + 0(a2h2) C2 = 0{o?h2)
Из (2.8) видно, что решение при \ ^ к ^ 1 при h —► 0 отличается от
точного решения уравнения (2.4) на 0(а2/г2), т.е. схема (2.5) устойчива и
аппроксимирует решение (2.4) со вторым порядком. Но если 0 < к < |,
то второе решение по абсолютной величине катастрофически растет с
уменьшением h при фиксированном хп, что свидетельствует о
неустойчивости схемы (2.5).
Действительно из (2.8) видно, что малая ошибка вследствие ап^
проксимации приводит к малому отклонению величины С^ от нуля при

§3. Численное интегрирование задачи Коши для ОДУ
21
1 > к ^ | и ошибка стремится к нулю при h —>• 0, а при х < | она
катастрофически растет.
Таким образом, хотя условие аппроксимации (2.4) выполняется при
к < |, но условие устойчивости нарушается и решение разностного
уравнения (2.5) не сходится к решению (2.4).
§3. Численное интегрирование задачи Коней для
системы дифференциальных уравнений первого
порядка
Схемы Эйлера (21). Схема Адамса (23). Построение схемы высокого порядка
методом разложения в ряд (25). Схема Рунге-Кутты (26).
Систему к обыкновенных дифференциальных уравнений (ОДУ)
первого порядка запишем в виде одного векторного уравнения
du ? и={иъ...ик) I
—- = туиЛ), где -> (3.1)
dt Д / = (/ь...Л) 1
Требуется определить функцию u(t) по начальным данным, заданным
при t = О
гГ(0) = щ O^t^T (3.2)|
щ — постоянный заданный вектор (размерности к).
К задаче (3.1)-(3.2) сводится широкий класс задач механики
твердого тела. Например, движение небесных тел, искусственных спутников,
ракет, а также задачи динамики механических систем материальных
точек, к которым приводятся задачи динамики сплошной среды и многие
другие задачи.
Рассмотрим некоторые из методов решения этой задачи, начиная
с самых простых, на равномерной сетке с постоянным шагом г = Т/п
(Т — общее время; п — число шагов по времени).
Схемы Эйлера. Дифференциальный оператор заменим на введенной
сетке разностным, одним из рассмотренных в §1 способов. Например,

22
Гл.1 Понятия о разностных схемах
заменим производную односторонней разностью вперед, а правую часть
в (3.1) возьмем в г-ой нижней точке. Тогда получим явную схему Эйлера
^^ = №,и) (з.з j
Здесь и далее нижний индекс обозначает номер шага по времени; щ —
уже вычисленные (известные) величины в предыдущих точках г,
соответствующих времени t^\ щ+\ — неизвестная величина.
Если правую часть взять в (г + 1)-ой точке, получим неявную схему\
Эйлера
^±l^l = f(ut+1,tl+1) (3.4J
При замене производной центральной разностью, получим явную или
неявную схему Эйлера с центральной разностью
— = f[uh U); — = /(г/г+ь ti+i) (3.5)|
Возможно и бесчисленное множество других схем, но пока ограничимся
приведенными выше.
Наиболее просто определяются значения щ+\ по явной схеме Эйле^
ра на всей сетке, последовательно находим
щ+1 = щ + т$(щ,и) г = 0,1... га (3.6)
в правой части (3.6) все величины вычислены на предыдущем шаге г.
При использовании неявной схемы получим
йг+1 = йг + т/(щ+ъ ti+i) г = 0,1... га (3.7)
в этом случае уравнение (3.7) для определения щ+\ является
нелинейным и его можно решить одним из итерационных методов, например,
методом простой итерации или методом Ньютона. Важной особенностью,
упрощающей решение уравнения (3.7) является то обстоятельство, что
для него всегда известно хорошее начальное приближение Щ^ = щ,
мало отличающееся от точного решения задачи щ+\.

§3. Численное интегрирование задачи Коши для ОДУ
23
В схеме (3.5) мы снова, как и в предыдущем параграфе,
встречаемся с ситуацией, когда формальный порядок разностного уравнения выше,
чем порядок исходного дифференциального уравнения (3.1) и
необходимо задавать дополнительное граничное условие.
Оценим эффективность приведенных схем с точки зрения
необходимого объема вычислений для получения решения, оставляя в стороне
другие характерные свойства этих схем. Очевидно, что основной
объем вычислений связан с вычислением вектор-функции / в правой части
уравнения и наиболее эффективной будет схема, в которой эта функция
вычисляется минимальное число раз для получения щ+\. В явной
схеме на каждом шаге это вычисление производится только один раз, в то
время как в неявной схеме это число равно числу итераций,
необходимых для получения решения с нужной точностью. В этом плане явная
схема экономичнее. Однако, точность определения решения зависит еще
и от порядка аппроксимации. Схемы высокого порядка аппроксимации
могут оказаться более эффективными, чем низкого порядка, хотя бы и
минимального объема вычислений.
В дальнейшем для простоты изложения будем рассматривать
интегрирование одного уравнения (3.1). Обобщение на случай системы
уравнений осуществляется формально.
Схема Адамса. Рассмотрим экономичную численную схему
высокого порядка аппроксимации. Экономичной считается схема, для которой
повышение порядка аппроксимации на единицу требует только одного
вычисления правой части уравнения (3.1).
Неизвестное щ+i ищется не по одному предыдущему значению щ\
а по (к + 1) значениям щ.. .Щ-k уже вычисленным на к предыдущих
шагах. Построим по этим найденным к значениям интерполяционный
многочлен
к
Ьк(щ ... щ-k) = ^2 lk(t) Ui~P (3>8)|
p=0

24
Гл.1 Понятия о разностных схемах
lpk(t) — базисные интерполяционные полиномы.
Базисные интерполяционные многочлены Лагранжа lpk(t) вычисляв
ются в интервале U^t^- U-k по формуле
pft\ = (t-ti)'-(t- tp-i)(t - tp+i) • • • (t + tj-k) ^ J
{tp — и) • • • (tp — tp-i)(t — tp+i) • • • (tp — ti-k)
Формулу (3.8) можно использовать для экстраполяции функции и\
на интервале [£j,^+i].
и(и + т)=щ+ J ^ dt = щ + [ f(u, О d£ (ЗЛО)
и и
Заменяя подинтегральную функцию f(u,r) ее интерполяционным
многочленом по формуле (3.8), получим
к и+т
U(U +т)=Щ + Т^2 fi-p / 1Ш) <% = Щ + ^(«о/г Н CLkfi-k) (3.11)
где fp = f(u(tp),tp); константы ai,...ai числа, независящие от
шага интегрирования т, которые вычисляются интегрированием
базисных интерполяционных полиномов ^(t), удовлетворяющих условиям
l*k(tm) = 5mpi m,p = i — к,.. .г, 5mp — символ Кронекера (см.
упражнение 7 на стр. 59).
При такой аппроксимации правой части погрешность будет равна]
\\f(u(t))-Lk(t)\\=0(rk)
тогда
u(t^ + T)-u(t^) = + а^_к + I
т
т.е. порядок аппроксимации будет равен числу точек, использованных
при интерполяции.
Эта схема требует запоминания к предыдущих значений правых
частей уравнения. Осложнение состоит в том, что при начале счета мы

§3. Численное интегрирование задачи Коши для ОДУ
25
не имеем этих значений. Чтобы схема начала работать первые к
значений должны быть вычислены нестандартным образом по какой-то
другой схеме; либо по схеме, порядок которой последовательно возрастает
от 1 до к при выполнении первых к шагов; либо методом Рунге-Кутты,
который будет изложен ниже.
Рассмотренное семейство схем носит название схемы Адамса. Они
являются схемами произвольного порядка аппроксимации (число точек,
по которым строится интерполяционный многочлен может быть
произвольным) и формально описывается разностными уравнениями, порядок
которых соответствует порядку аппроксимации и следовательно
содержит лишние паразитные решения, о которых говорилось выше в §2J
Рассмотрим этот вопрос на примере схемы Адамса второго порядка:
^^ = lf(ut)-1-f(ul.1) (3.13J
Если положить f(uj) = ащ, то получим разностную схему второго
порядка, как и в схеме (2.5) и следовательно для ее решения требуются
задавать дополнительные начальные условия при t = т согласованно с
порядком аппроксимации уравнения (3.13) равным 0(т2), как это было
сделано при решении уравнения (2.5).
Построение схемы высокого порядка методом разложения в
ряд. Разностную схему высокого порядка точности 0(тк) для
уравнения (3.1) можно получать разложением щ+\ = u(U + г) в ряд Тейлора
в точке ti
2 k—1
ul+l = щ + u'f + uf- + --- + uj*"1)^-^ + 0{тк) (3.14
и вычислением производных и\, и", ... vSk~^ в точке U путем
последовательного дифференцирования уравнения (3.1)
du d2u ,
-77 = ДМ), -Тл = f^u + ft
dt ,3 dt (3.15)
-^ = fuu{u'f + fuu" + ftt

26
Гл.1 Понятия о разностных схемах
Подставляя (3.15) в (3.14) и ограничиваясь первыми тремя членами
разложения, получим разностую схему третьего порядка аппроксимации
для уравнения (3.1).
^^ = /(«i,*0-(/»/ +Л) J+
т 21
+ {fuuf + к (/„/ + ft) + /«) ^ + 0(Т3) (3.16J
Откуда видно, что в формуле (3.16) число членов в
коэффициентах ряда очень быстро растет и вычисления становятся не
эффективными, т.к. для повышения порядка на единицу требуется многократное
вычисление правой части и ее производных. Этого можно избежать, если
вычислять производные не в дополнительных точках интервала [£^, ^+i]-
Схема Рунге-Кутты. Схема 3.16 принадлежит к разностным схемам
вида
^^ = P(U,Ul) (3.17J
где правая часть Р{Ь{,щ) = Pi[f{ui),U], строится зависящей
определенным образом от правой части /(щ) исходного дифференциального
уравнения (3.1) и аппроксимирует его с некоторым порядком 0(тк)\
К такому классу разностных схем принадлежит, например схема]
Эйлера предиктор-корректор. Решение вычисляется по двухшаговому
алгоритму:
1) Предиктор: Вычисляем г^+1/2 на полушаге по явной схеме; правая
часть в формуле (3.6) берется в виде |r/(^,t«);
^г+1/2 = Щ + \rf(Ui,ti)
2) Корректор: Вычисляется щ+\ в центральной точке (г + 1/2) при
правой части т f(ui+1/2, ti+1/2), где иг+1/2 = щ + \rj[u%, t%)\
Окончательно правая часть в уравнении (3.17) будет
щ+1 -щ ( т т\
( т
Р(щ, U) = f \Ui + -/(щ, U), U + ■

§3. Численное интегрирование задачи Коши для ОДУ
27
Нетрудно проверить, что эта схема имеет второй порядок
аппроксимации. Разложим функцию P(ui,ti) в ряд Тейлора как функцию двух
переменных Щ И t{.
P(uh U) = f [Ui + -f(iii, U), U + -J =
= [/(^, U) + ^fu(uu ti)f(ui, U) + ^ft(uu U) + 0(т2)]
с другой стороны имеем
г2
Щ+1 = щ + тщ(и) + -juttiU) + • • • + 0(тк) = щ + т Р(щ, U) (3.18)|
Производные по t при t = U легко определить из исходного
дифференциального уравнения (3.1)
щ(и) = f(ui,ti)
щ(и) = 1и(щ, и)щ + ft(ui, U) = /и(щ, и)/(щ, U) + ft(uh и)\
Откуда подставляя найденные решения в (3.18), получим оценку для
остаточного члена, определяющего порядок аппроксимации схемы
Эйлера предиктор-корректор
Гг+1 = ~(щ+1 ~ Щ) ~ Р(Щ, U) = 0{Т2)
т
Этот метод является наиболее простым методом, принадлежащим
к семейству схем Рунге-Кутты высокого порядка точности. Идея
метода состоит в замене многократного дифференцирования правой части
ее вычислением в к промежуточных точках интервала [щ^щ+г].
Комбинация этих значений может быть выбрана так, что будет эквивалентна
с точностью до остаточного члена выбранного порядка 0{тк)
соответствующему отрезку ряда Тейлора для щ+\ в (3.18). Этот способ требует
лишь вычисления функции f{u^ti) в нескольких дополнительных
точках, число этих точек равно к и соответствует выбранному А;-ому порядку

28
Гл.1 Понятия о разностных схемах
аппроксимации. Сначала вычисляются вспомогательные величины)
Ащ = кг = /(щ,и)т
Аи2 = к2 = f(ui + 0.5fci, U + 0.5т)т
Л^з = к3 = f(Ui + 0.5/с2, U + 0.5т)т
Ащ = к^ = /(щ + 0.5&3, ^ + ^)^
Окончательное значение г^+i с четвертым порядком точности
вычисляется по формуле
Щ+1 = Щ + -(к1 + 2к2 + 2/с3 + к4) (3.19)
Коэффициенты при ki подбираются из условия, чтобы правая часть
в (3.19) совпадала с разложением в ряд Тейлора с точностью до 0{тъ)^
что обеспечивает четвертый порядок аппроксимации для разностной
схемы уравнения (3.1). Очевидно, что для вычисления функции Р(и)
требуется, чтобы функция в правой части f(u,t) была бы трижды
дифференцируемой по своим аргументам.
Таком образом, схемы Рунге-Kymmu обладают такой же
экономичностью как и схемы Адамса. Каждый шаг в этой схеме повышает
порядок аппроксимации на единицу и требует при этом только одного
вычисления правой части уравнения (3.1) в точке интервала.
Устойчивость для схем Рунге-Кутты следует из следующей
теоремы.
Теорема Система разностных уравнений (3.17) устойчива, если\
1) Функция Р(и) удовлетворяет условию Липшица
\\Р(х)-Р(у)\\^С\\х-у\\
2) шаг интегрирования т достаточно мал и удовлетворяет
условию Ст <С 1. Доказательство этой теоремы приведено в курсах по
вычислительной математике, см., например, [2, 20].

§4. Решение задачи Коши для жестких систем ОДУ
29
§4. Решение задачи Коши для жестких систем ОДУ
Численное решение (31). Исследование устойчивости (32). Сингулярно-возмущенные
системы (33). Растяжение стержня из нелинейного вязкопластического
материала (34).
Среди систем ОДУ следует выделить класс, так называемых,
жестких систем, численное интегрирование которых вызывает особые
затруднения, в следствии того, что скорость изменения решения этих
уравнений крайне неравномерна. Поле направления интегральных
кривых для таких систем почти мгновенно меняет свое направление при
подходе к определенным траекториям. Решение задачи Коши или
краевой задачи включает в себя области очень быстрого изменения на малом
отрезке, которые сменяются областями медленной эволюции. Области
быстрого изменения в механике называют погранслоями или
внутренними слоями, а области медленного изменения — квазистационарным
режимом.
Рассмотрим систему уравнений
Систему (4.1) будем называть жесткой в окрестности решения
и = щ, если для матрицы Якоби G = -Л\ = /^ число
и\й0 °
тах|А^(£)| н , J
N(t) = . ';,( » 1 1 < г ^ п 4.2
где А^ — собственные числа матрицы G.
Это означает, что отдельные составляющие решения сильно разли^
чаются масштабами изменения по t. Спектр собственных значений
матрицы G можно условно разделить на 2 части: жесткую, для которой
выполняется условия
Re\i(u) ^ -L |1тЛг(г?)| < \Ке\г(и)\
и мягкую, для которой
\Хг{и)\ < /<L

30
Гл.1 Понятия о разностных схемах
Очевидно, что \{(й) зависит от решения и и нелинейная система (4.1)
может в разных частях фазового пространства обладать различной
жесткостью.
Отношение п = L/1 называют показателем жесткости систе]
мы. В прикладных реальных задачах п может достигать очень
больших значений от 107 до 1015. Тогда интегрирование с
обычной точностью стандартным методом на отрезке времени Т
потребует шага т, при котором выполняется условие тЦ/яЦ <С 1. Учитывая
||/й|| ~ тах|Л^| « L, получим т <С 1/L и необходимое число шагов
будет равно m = Т/т ^> TL ~ 1015. Это совершенно неприемлемо, если мы
интересуемся квазистационарным режимом для которого Т ~ 0(1), а не
структурой погранслоев. Для определения квазистационарного режима
вполне достаточно, чтобы m ~ Ю3. Поэтому нашей задачей будет
построение алгоритма, который позволял бы вести расчет с таким большим
шагом т.
Т.к. жесткость системы определяется матрицей щ , т.е. при за-
данном и = щ линейной частью /#, то в первом приближении
достаточно исследовать линеаризованную задачу и рассматривать вместо
системы (4.1) систему
— = Ш^)-й (4.3)|
Исследование устойчивости проводится на модельном уравнении
du Л ,Л \
— = \и (4.4)
dt v 1
где и — скалярная функция, Л — комплексное число, т.к. собственные
числа матрицы Якоби могут быть комплексными. При этом
рассматривается множество всех точек комплексной плоскости /х = т • А, для
которых разностная схема уравнения (4.1) будет устойчивой. Например,
явный метод Эйлера будет устойчив лишь в круге единичного радиуса
с центром в точке (—1,0), поэтому он непригоден для интегрирования с
большим шагом. Если метод устойчив на всей полуплоскость Re /x < 0, то
он называется А-устойчивым или абсолютно устойчивым. Т.к. решение

§4. Решение задачи Коши для жестких систем ОДУ
31
уравнения (4.4) устойчиво при Re Л < 0, то А-устойчивость ее
разностной схемы означает, что она устойчива при любом т > О, поскольку для
разностной схемы устойчивость определяется произведением Хт = /х\
Численное решение. Для численного решения жестких задач
необходимо использовать А-устойчивые или А{а)-устойчивые алгоритмы.
А(а)-устойчивой называется разностная схема, для которой область
устойчивости в комплексной плоскости /i ограничена углом а (рис. 4.1)
|arg(-/i)| <a.
К таким алгоритмам относятся неявные схе- . I A ti
мы Гира высокого порядка аппроксимации. Мно- \
гошаговые неявные схемы Гира высокого
порядка аппроксимации строятся аналогично схеме
Адамса с помощью интерполяционного
многочлена для функции f(u,t) в правой части уравне- ^—
ния (4.1), с той только разницей, что в число
узлов, по которым ведется интерполяция
включается и узел (п + 1), в котором ищется реше- /
ние. Интерполяционный полином (га + 1)-ой сте- /
пени аппроксимируется по узловым значениям
/(un+l)f„+l))/(u„)in)) >/(u„+l-mt„+l-m)_ Тоща рис 41
в отличии от явных схем Адамса, рассмотренных устойчивости
\
\
СМ
О
тц
R^H
Области
ком-
в §3, получим неявные схемы (га + 1)-ого порядка плекснои плоскости /х
и
п+1
- гг
Т
^2аг!^
fn+l—% у.п+1—%
(4.5
=0
Например, при га = 3 получим неявную схему четвертого порядка
точности
и
п+1
■и10
1
(9/n+1 + 19Г - bfn~l + /■
п-2\
т 24
Явная схема аппроксимируется по узловым точкам f(un, t
(4.6
f(un-l,tn-l),...f(un-m,tn-m). Сумма в (4.5)
начинается

32
Гл.1 Понятия о разностных схемах
При т = 3 получим схему третьего порядка точности
7/п+1 _ ?/п 1
= - (23Г " 16Г"1 + 5Г"2) (4.7|
Для решения системы неявных уравнений (4.6) можно
использовать схему предиктор-корректор. Предиктор вычисляется по явной
схеме (4.7).
~П+1 _ П 1
= j~2 (23Г " 16Г"1 + 5Г"2)
После чего решение уточняется по неявной схеме (4.6)
П+1 _ П 1
= ^ (9Г+1(й"+1,а:"+1) + 19Г - 5Г"1 + Г"2}
Исследование устойчивости. Рассмотрим исследование
устойчивости на примере неявной схемы Гира второго порядка точности.
ип+2 - \un+l + \ип = \rf{un+\ tn+2) (4.8)1
Исследование проведем для линеаризованного уравнения с правой
частью f(u) = A?i, где Л — комплексное число.
Решение уравнения (4.8) имеет вид
v,n = C1(r1)n + C2(r2)n
где 7*1 и Г2 — корни уравнения
2 ± л/1 + 2т А
Г1'2= 3-2гЛ ' ТА = ^
Как уже указывалось, решение уравнения (4.8) представляет ин^
терес для жестких систем при малых |/z| <^ 1 для определения
структуры погранслоев. Решение разностного уравнения (4.8) аппроксимирует
точное решение уравнения (4.4) и = ext. Эту область плоскости
комплексного переменного /i называют "областью точности". Легко видеть,

§4. Решение задачи Коши для жестких систем ОДУ
33
что первый корень уравнения (4.9) при малых /j, аппроксимирует точное
решение
г1(д) = е"(1 + 0(/.))
щ = Ci(n)" = С\еХтп (1 + 0{ii2)) , t = nr
Второй корень при \ji\ <С l
Г2 = \ + 0(р2), и2 = С2(г2)п^0
Откуда следует, что схема при малых /j, устойчива. Чем шире "область
точности" \/л\ < /iQ, в которой решение аппроксимируется с нужной
точностью, тем лучше схема.
Вторая область это область больших |/i| >> 1. Ее называют об-\
ластью устойчивости (квазистационарный режим). Легко видеть, что
при \fi\ ^> 1 корни (4.9) г\ ~ Г2 ~ (2|/х|)-1/2 С 1 и следовательно схема
устойчива. Однако, могут быть схемы устойчивые не для всех значений
ji в полуплоскости Re/x < 0. Например они устойчивы только в угловой
области |arg(—/i)| < а или в области Re/i < —а2. Такие схемы также
пригодны для получения медленно меняющегося решения
(квазистационарный режим). Необходимо, чтобы зона устойчивости содержала бы
достаточно широкую окрестность линии Im/i = 0, Re/x < 0. Более
полное изложение методов решения жестких систем можно найти в [15].
Сингулярно-возмущенные системы. Во многих физических
задачах, относящихся к классу жестких малый параметр г явно входит в
систему уравнений. Это упрощает интегрирование таких систем
уравнений. Пусть имеем систему двух уравнений такого вида
ей = f(u, v)
М } (4.10
v = (p(u,v)
В первое уравнение входит малый параметр е, как коэффициент при
производной, либо можно считать, что правая часть содержит большой

34
Гл.1 Понятия о разностных схемах
параметр L = е~1 ^> 1. Обе функции f(u,v) и ip(u,v) и их
производные — величины порядка 0(1). Спектр матрицы Якоби этого уравнения,
определяется из условия
det(Lfu~X Lfv )=0 (4.11)
Имеем жесткую часть, отвечающую функции L/u, а компонента
отвечающая функции (р будет мала. Квазистатический режим
определяется уравнением f(u, v) = 0. Оно разбивает плоскость и, v на две части.
Область f(u,v) > 0 является жесткой при Ai < 0. Теория таких систем
хорошо разработана [15].
Случай сингулярно-возмущенной системы с явно выделенным
большим параметром аналогичен общему случаю системы (4.1),
рассмотренных выше. Здесь большой параметр L играет ту же роль, что и | Re Amax|,
а качественная картина поведения решения в этом случае достаточно
прозрачна. Проведем это исследование на конкретном примере.
Растяжение стержня из нелинейного вязкопластического
материала. Рассмотрим решение сингулярно-возмущенной системы на
примере задачи о равномерном растяжении стержня из нелинейного
вязкопластического материала со стационарной диаграммой общего вида
а = as(s) [11], показанной на рис. 4.2.
Рис. 4.2. Стационарная и квазистационарная диаграммы материала]

§4. Решение задачи Коши для жестких систем ОДУ
35
Примем, что к концу стержня приложено напряжение,
изменяющееся по заданному закону а = ao(t).
Задача в безразмерных переменных может быть сведена к решению
системы, состоящей из двух уравнений с большим параметром 5 ^> 1
дбр дб дет
-т=т-т=Нщпа)ф^-а^)) ,412
да У ' А
= a0(t) t > 0 cr|t=0 = сг0
Первое уравнение системы (4.12) связывает напряжение а с
деформацией е в упруговязкопластическом материале, второе задает скорость
изменения приложенного напряжения. Функции в правых частях
уравнений Ф(сг, е)1 &'$(£) и их производные величины 0(1). Спектр вариационной
матрицы Якоби (4.11) определяется уравнением
det^o-A^=-(^-A)A = 0
Здесь параметр 6 = to/т в правой части (4.12), равный отношению
характерного времени Ц к времени релаксации т, для многих материалов
большая величина 5 ^> 1 и система уравнений (4.12) является
сингулярной системой вида (4.10).
Функция
' Ф(я) z > О
Ф(г) ,
Откуда следует, что скорость пластической деформации ёр = 0 при
|сг| < crs(e) и следовательно напряжение связано с деформацией законом
Гука в безразмерных переменных а = е.
Учитывая, что Ф£ = Фг dz/ds = —Фг das/ds: получим, что система]
жесткая при das/de > 0 и не является таковой при das/de < О, Фг > 0.
Поле траекторий интегральных кривых в фазовой плоскости а — s\
легко исследуется. Кривая а = as(s) разбивает плоскость на две части.
Слева от нее а — as(e) > 0, а справа sp = 0.

36
Гл.1 Понятия о разностных схемах
Вне малой окрестности 0{5~1) кривой а = as(e) поле направлений
интегральных кривых почти горизонтально, скорость изменения е очень
велика (порядка 0{8~1)) и направлена в сторону возрастания, т.е.
происходит быстрый рост пластической деформации. За малое время 0(6~1)
стержень из состояния Ai(cr°,£°) переходит по почти горизонтальной
прямой в точку С (а*, е*) в окрестности кривой а = crs(e). В этой
окрестности а = 0(1) и ё = 0(1), так как <&(z) = 0(£-1), и напряженно
деформируемое состояние изменяется вдоль возрастающей ветви CD
кривой а = as(s) до точки D, где das/ds = 0. Дальнейшее движение вдоль
кривой на падающем отрезке становится неустойчивым, система
теряет жесткость и происходит быстрое движение по горизонтали в точку
Е(а*,е*), рис. 4.2.
Чтобы описать изменение е на участках быстрого изменения (е0, £*)|
достаточно положить в первом уравнении (4.12) а = ао начальному
значению и проинтегрировать полученные уравнения
£
-^ = +5Ф (а0 - as(e)) i= —
dt J 5Ф(а0 -as{e))]
где учтено, что в начальный момент времени £q = ctq, поскольку, как
это следует из (4.12), мгновенное деформирование происходит по закону
Гука.
В точке (е = £*, о"о = crs (£*)), в зависимости от асимптотического
вида функции Ф(^) = aza при z —> 0 интеграл будет сходится при а < 1
и расходится при а ^ 1. Соответственно время, за которое е —> £* будет
конечным или бесконечным в масштабе 0(S~1). Однако и при а ^ 1
стремление е —► е* очень быстрое экспоненциального характера и
"эффективное время" такого перехода всегда конечная величина в масштабе
0(6-*).
Если переход совершается из точки неустойчивости, то имеем внут^
ренний погранслой. Начальной для него будет точка неустойчивости
а = а* и £ = е* при t = t* и время перехода на устойчивую ветвь будет

§5. Разностные схем уравнений в частных производных
37
вычисляться по формуле
s
de
U I 5Ф(а*
- <rs(e))
Таким образом, при медленном нагружении напряжением,
квазистационарная зависимость а = а(е) будет иметь вид кривой OA\CDE\
(рис. 4.2). Для нее характерно повышение предела текучести по
сравнению со стационарной зависимостью а = <JS(£) и появление площадки
идеального скольжения.
§5. Примеры составления разностных схем для
уравнений в частных производных
Волновое уравнение (37). Схема Лакса (38). Схема Лакса-Вендроффа (39).
Уравнение теплопроводности (42).
Рассмотрим простейшие разностные схемы для уравнений в
частных производных эволюционного типа, зависящих от одной
пространственной переменной и от времени. Начнем с уравнений акустики или
распространения упругих волн в стержне.
Волновое уравнение. Полная система уравнений распространения
продольных волн в стержне состоит из уравнения движения, уравнения
совместности поля скоростей и деформаций и закона Гука. При
одноосном растяжении/сжатии она имеет вид
dv да де dv ^ . \
рт=Ш' т = д~х> а = Ее (5Л1
Систему уравнений после исключения напряжения а можно пред^
ставить в виде системы двух волновых уравнений относительно скорости
v и деформации е
dv 9 де о Е
я7 = а я а = ~
I _ ajx * н
dt дх

38
Гл.1 Понятия о разностных схемах
Аппроксимация уравнений (5.2) строится на прямоугольной сетке
в плоскости xt. Шаг сетки в направлении оси х обозначим через /i, в
направлении t через т. Индекс г обозначает номера точек по оси ж,
индекс п — номера по оси t (рис. 5.1).
(г+1/2,и+1)
(М,и) (/+1/2,и) (Н-1,и)
0>+1)
ТТ"
' й+1
сЯ-— Н-3? и+щ
•-1/2J I ]г+1/2 "
1+4
(б)
Рис. 5.1. Шаблоны для схемы Лакса (а) и Лакса-Вендроффа (5)|
Определение Шаблоном разностной схемы будем называть точки,
задействованные в одном разностном уравнении для получения решения
в точке на (п + 1)-ом слое.
Схема Лакса. Шаблон может быть двухслойным или многослойным
по времени £, в зависимости от числа используемых слоев.
Аппроксимируем систему (5.2) на двухточечном шаблоне по х (рис. 5.1а), когда на
п-ом слое используются точки г + 1 и г — 1.
п+1 _ п п п
°г Ui+l/2 = ^ £i+l ~ £г-1
2h
£i ~ £i+l/2 V,
Щ+1/2 = \{vi+i + Vi-i) (5.3)
i+l ui-l
2h
Схема (5.3) называется схемой Лакса.
Определим порядок аппроксимации этой разностной схемы по а: и t.|

§5. Разностные схем уравнений в частных производных
39
Разлагая v™+l) v™+l ,2 и т.д. в ряд Тейлора в точке (г, гг), находим^
d2v\n
дх2/
= а2
'(-?)
[(§)>
[дх3
" Ь2
3!
+ о h4
Откуда следует, что локальная аппроксимация имеет порядок о(т + h2/r).
Член о(/г2/т) стремится к нулю, если только порядок h2 меньше или
равен порядку т при /г —» 0 и т —► 0. Такая аппроксимация называется
условной.
Схема Лакса-Вендроффа. Используем трехточечный шаблон по х,
состоящий из точек г — 1, г, г + 1, и получим другую разностную схему,
показанную на рис. 5.16
(5.4
<+1-«f o^-e?-!
т "а 2/г '
1)>@)>^!
(shs);h*-i
£г ~ £i _ Vi+l ~ Vi-1
т 2/г
[(1И£);м
:i)xsn^
Очевидно, что разложение в точке (г,гг) дает порядок аппроксимации
о(г + /г2).
Чтобы получить аппроксимацию одного и того же второго поряди
ка и по т и по /г поступим следующим образом: выразим остаточные
члены |^2^ и f-§2^, в левых частях уравнений (5.4) их через вторые
производные по ж, используя исходную систему уравнений (5.2), и
аппроксимируем полученные выражения вторых производных конечными
разностями с точностью до o(h2).
(si+i - 2si + Si-i) + o(h2]
(vi+i - 2vi + Vi-i) + o(h2)\
тд2е
2dt2 ~
rd2v
2dt2 ~
та2 д2е
2 дх2~
та2 d2v
2 дх2 "
та2
~ 2/г2
та2
" 2/г2

40
Гл.1 Понятия о разностных схемах
Окончательная разностная схема будет такой
■ v? т а2
^(<+1-2< + ^) = ^^
т
2 \£i+l Z£i + £i-l
s?+1-£? та2 -n+1
(5.5
т 2h2K г+1 г г"1У 2h
Нетрудно показать из (5.4), что порядок аппроксимации на
решениях системы (5.2) будет теперь о(т2 + /г2).
Здесь мы по-существу несколько сузили понятие аппроксимации и
использовали аппроксимацию на классе точных решений
дифференциальных уравнений (5.2). В этом случае вместо термина "условие
аппроксимации" используется термин "условие согласованности" (смотри [17]).
Оно так же, как и "условие аппроксимации", показывает, насколько
хорошо точное решение задачи удовлетворяет конечно-разностным
уравнениям.
Разностную схему (5.5) называют схемой Лакса-Вендроффа. Ту же
самую разностную схему можно получить иным способом, введя
промежуточный слой с индексом п + \ и используя двухшаговую схему,
(рис. 5.16"). На первом полушаге (предиктор) находим решение в точках
г + 1/2 и г — 1/2 по схеме Лакса (5.3), а затем находим окончательное
решение на слое п + 1 с помощью схемы "крест"
П+1 П РП+1/2 РП+1/2 П+1 П 7)П+1/2 11П+1/2
<+ ~ %П _ а2 £г+1/2 ~ gz-l/2 . gj ~ei _ ^г+1/2 ~ Vl/2
г 2h ' т 2h
Исключая величины в точках с полу целыми индексами, найдем^
а2 р+1-СЛ + т а2 /4+1 " W + v?_
2/i У 2 V h2
^г ^г_ — ^г+1 ^г-1 , га / £г+1 Zgz + £г-1
,+1 еп „п _.п о ^ о.п, .п X (5-6)|
г 2/i 2 \ /г2
Полученная двухшаговая схема (5.6) совпадает со схемой второго
порядка Лакса-Вендроффа (5.5).
Можно построить неявную схему второго порядка точности
0(т2 + Ь?) на четырехточечном шаблоне (рис 5.2а), когда на п-ш и

§5. Разностные схем уравнений в частных производных
41
п + 1-м слое используются по два узла и производная по х в (5.2)
аппроксимируется выражением
-О;
—О—
i /+1/2 /+1
(а)
-п+\
■W+I
/-1
I
(б)
/+1
Рис. 5.2. Шаблоны для неявных схем первой (а) и второй производной (б"]
иг+1
dv
дх " h
п+1 _ п
Ui+l/2 Ui+l/2
..n+l
+ (1 - К)
+ (!■
О ^ к ^ 1
-i+l '
n+l
ьг+1/2
^г+1/2
„п+1 .
„п+1
(5.7
+ (1
1
^г+1/2 = 2^+1+Vi)
При к = 1/2 порядок аппроксимации будет вторым, это легко доказать,
раскладывая все величины в точке (г + 1/2, п + 1/2). Схема (5.7)
существенно отличается от предыдущих (5.3)-(5.6) тем, что на верхнем
п + 1-м слое каждое уравнение содержит не одно, а два неизвестных,
и полная система уравнений для (п + 1)-го временного слоя не
распадается на систему рекуррентных уравнений. Такие схемы называются
неявными, в отличие от явных схем, для которых решение определяется
в каждой точке (п + 1)-го слоя независимо от других точек этого слоя, т.е
матрица системы уравнений для определения величин с индексом п+1
имеет диагональный вид.
Явные разностные схемы позволяют определять решение на

42
Гл.1 Понятия о разностных схемах
(п + 1)-ом временном слое, если решение известно во всех точках
предыдущего гг-го слоя, т.е. позволяют решать разностную задачу Коши,
либо решать задачу с периодическими условиями, заданными на отрезки
оси х. Решение начально-краевой задачи, когда кроме начальных
условий задаются еще и краевые условия на краях отрезка х = 0 и х = 1,
требует построение в этих точках специальных схем. Схема для внутренних
точек отрезка оказывается непригодной для определения решения в
точках, расположенных на границах отрезка, т.к. одна или несколько точек
шаблона для внутренней точки оказываются вне пределов отрезка]
В случае неявных схем число точек, используемых в разностной
схеме на (п + 1)-ом слое больше единицы и для их определения
получается система алгебраических уравнений, для замыкания которой
необходимы краевые условия. Решение можно определить только для всего
временного (п + 1)-го слоя одновременно решая систему алгебраических
уравнений. Это свойство неявной схемы находится в противоречии со
свойством волнового уравнения, решение которого в точке х, t зависит
от решения в тот же момент времени t в других точках стержня, т.к.
скорость распространения возмущений в упругом теле конечна. В
дальнейшем этот вопрос будет исследован более детально, а пока ограничимся
сделанным замечанием.
Уравнение теплопроводности. Рассмотрим семейство разностных
схем для уравнения теплопроводности
дТ _ д2Т
dt дх2
Введем безразмерные переменные x = j,i=y£,T = f-, где / и То
характерные длина и температура. Тогда получим
ж = ^ М
В дальнейшем черту над безразмерными переменными опускаем.
Вторую производную аппроксимируем на шеститочечном шаблоне
(рис 5.26"), в котором используются точки %—1, г, г+1 на гг-ом и (п + 1)-ом

§5. Разностные схем уравнений в частных производных
43
слоях по времени
грП+1 грп грп+1 г\грп+1 _|_ грп+1 грП слгрп _|_ грп
1{ -1{ =Hl_i±1 Щ +1i-i +{l_x)1i+i пг +A-i (5g)
т hz hz
Схема при 0 < к ^ 1 неявная и имеет порядок 0{Ь? + т), при к = 0 она
становится явной, а при я = \ имеет второй порядок аппроксимации
0(h2 + r2).
Нетрудно видеть, что если задача (5.8) решается на отрезке оси х £
[0,1] на сетке с узлами i = 1... N1 то разностных уравнений (5.9) при
О < х ^ 1 можно записать только (N — 2) для узлов i = 2 ... N — 1.
Чтобы система уравнений была бы замкнутой, к ним следует добавить два
уравнения, которые следуют из краевых условий, задаваемых на
концах отрезка х = 0 и х = 1, тогда число уравнений будет равно числу
неизвестных.
Для уравнения теплопроводности это свойство неявной схемы
вполне физично, т.к. тепловые возмущения распространяются с
бесконечной скоростью и все точки стержня в момент времени t влияют друг
на друга. Для уравнения теплопроводности наоборот является нефизич-
ной явная схема при к = 0. Однако это обстоятельство существенно
только в случае быстрого изменения решения или для высокочастотных
решений. Для гладких низкочастотных решений оно малосущественно и
для решения таких задач применимы как явные, так и неявные схемы.
Отметим, что увеличивая число точек шаблона, можно повышать
порядок аппроксимации, но при этом сильно усложняется система
разностных уравнений и усложняется ее исследование. Поэтому повышение
порядка аппроксимации не всегда оправдано. В дальнейшем при
аппроксимации эволюционных уравнений будем как правило, ограничиваться
двухслойными схемами первого и второго порядка точности.

44
Гл.1 Понятия о разностных схемах
§6. Исследование устойчивости разностных схем|
Устойчивость разностной схемы от 2х переменных (45). Условие устойчивости фон
Неймана (46). Устойчивость волнового уравнения. Условие устойчивости
Куранта (47). Устойчивость уравнения теплопроводности (50). Принцип замороженных
коэффициентов (51). Устойчивость при решении краевых задач (55). Выбор шага
при решении уравнения теплопроводности по неявной схеме (56). Выбор шага при
решении волнового уравнения (58).
Двухслойная система разностных уравнений в общем случае может
быть записана в виде
В^1 - В0гГ = О
или в более точной форме
N iVi
Y, bitP (^+1) - Y1 вот/5 №) = ° (6Л)|
/3=0 /3=0
где ТР(и) = u(x + f3h) есть оператор сдвига по оси Xj Bq , B^ матрицы
той же размерности, что и вектор неизвестных г?, элементы которых
постоянны, но могут зависеть от размера шагов г и /г,, /3 — целое число. Для
явной схемы множество точек N состоит из одной точки и матрица Bi
диагональная, для неявных из нескольких точек соседних с точкой Х{.
Если к уравнению (6.1) применить преобразование Фурье по х\
оо
й(к) = — / u{x)e~ikxdx
2тт J
(где крышка над переменной обозначает ее изображение в плоскости
комплексного переменного к) и учесть, что оператору сдвига соответствует
изображение
Т0{и)
ОО
-3- [ u{x + f3h)e-lkxdx = elkphU(k)
2тт J
то получим
Н1йп+1(к)-П0йп(к) = 0 (6.2)|

§6. Исследование устойчивости разностных схем
45
где
N N
Hi = ^ Вг еХР(^ Phk)^ Н0 = Yl B0 eXP(^ Phk\
/3=0 /3=0
А: параметр преобразования Фурье.
Разрешая (6.2) относительно йп+1, получим систему рекуррентных
уравнений для изображений йп+1
un+\k) = G(r,h1k)un(k) (6.3)|
где G = Н_1Но матрица перехода на п + 1 слой в пространстве
изображений. Применив п раз оператор G к й°(к), получим решение на
(п + 1)-ом слое в виде произведения
йп+1(к) = Gn(r,h,k)u°(k) (б.4)|
Устойчивость разностной схемы от 2х переменных. Используем
понятие устойчивости, по-существу эквивалентное ранее приведенному
в §2, но в иной, более удобной для дальнейшего, формулировке\
Разностную схему (6.1) будем называть устойчивой, если при\
некотором Т\ > 0 иТ > 0 бесконечное множество операторов перехода
Gn(r, k) 0 < г < п и О^пт^Т
будет равномерно ограниченным ||Gn(r, к)\\ < С, где С не зависит от т
и к. Это условие является необходимым и достаточным для
устойчивости.
Уравнение (6.2) является аналогом (6.1) в пространстве
изображений и G(t, к) является матрицей, зависящей от параметра
преобразования Фурье к. Условие устойчивости требует, чтобы операторы,
отвечающие матрицам Gn(r, к), были бы равномерно ограничены на конечном
интервале t для всех к.
Норма матрицы ||А(&)|| определяется формулой
||A(A)||=maxJ^§^ (6.5)1
УфО \v\

46
Гл.1 Понятия о разностных схемах
где \v\ = (Иг;2)1/2 — длина вектора.
Спектральным радиусом матрицы А называют R = max |Л^|; где\
\ — собственные числа матрицы А. Очевидно, что R ^ ||А||.
Спектральный радиус Rn матрицы п будет Rn = Rn. Кроме того,
||д2|| |А(АгЛ| |А(АгТ)| |Аг7| lAtfl |АгЯ
А2 =тах' ;^, л=тах' , \ _ л l—-1 < maxj M '
v^o \v\ ^о |Аг7| |гТ| ^б И |г?|
т.к. пространство векторов гТ шире, чем Av.
Следовательно, ||Gn(r, k)\\ < ||G||n i?n < ||Gn|| < ||G|H
Условие устойчивости фон Неймана. Необходимым условием
устойчивости является существование константы С, ограничивающей
спектральный радиус матрицы Gn(r, k)
Rn(r, k)^C
Д(т, к) < С1/п 0 < п < -
г
в частности, должно выполняться
R < Ст/Т
На конечном интервале 0 < г < т\ показательная функция в правой
части от г может быть ограничена линейной функцией
R < Ст/Т < 1 + Ci r
Откуда следует, что необходимое условие устойчивости разностной
схемы будет выполняться, если все собственные числа матрицы перехода G
удовлетворяют условию
Я = шах|Аг| ^ 1 + 0(т) (б.б)|
г
Это условие было получено фон Нейманом и носит его имя.
Если матрица G нормальная, т.е. перестановочная со своей
комплексно-сопряженной и транспонированной G G* = G*G, то
спектральный радиус равен норме матрицы и условие Неймана будет не
только необходимым, но и достаточным.

§6. Исследование устойчивости разностных схем
47
Заметим, что если решение исходного разностного уравнения (6.1)
искать в виде
v^p = ifm(k)XneXp(i0kh)
то непосредственно получим уравнение для определения
характеристических чисел Л матрицы G(r, к) из исходной системы разностных
уравнений
[ЛЕ - G(fc)] и°т = 0 det[AE - G(fc)] = 0
Этот прием удобен при практическом исследовании устойчивости
конкретных схем.
Исследуем устойчивость разностных схем, приведенных в §5 урав^
нений акустики и теплопроводности, для которых исследовалась
аппроксимация.
Устойчивость волнового уравнения. Условие устойчивости
Куранта. Рассмотрим схему Лакса. Подставляя в (5.3)
iC^ = \n+l^meMikh(3)
получаем
£0(А — cos kh) — vq— sin kh = 0
ri N
vo(X — cos kh) — Sq a2—- sin kh = 0
Из условия равенства нулю определителя системы (6.7) находим
а2т2
(А — cos kh)2 H—— sin2 kh = 0
Л = cos kh ± i —- sin kh (6-8)
CLT
|Л| ^ 1, если— ^ 1
т.е. схема условно устойчива.
Условие ат/h ^ 1 называется условием Куранта и является
необходимым условием устойчивости. Оно устанавливает зависимость

48
Гл.1 Понятия о разностных схемах
между шагом по пространству и шагом по времени и справедливо для
любого гиперболического уравнения. Смысл его в том, что шаг по
времени т должен быть таким, чтобы область зависимости решения в
точке х на п + 1 слое разностного уравнения была бы не меньше, чем
область зависимости дифференциального уравнения^ которая
определяется наклоном характеристик, выпущенных из точки х до пересечения с
п-ш слоем рис. 6.1.
Можно дать простое физическое объяснение этому условию. Если
оно нарушается и характеристики уравнений (5.2) проходят, как
показано на рис. 6.1 пунктиром, то решение разностного уравнения можно
сделать как угодно сильно отличающимся от решения
дифференциального уравнения, приложив на отрезках А г и А\{г + 1), которые не входят
в область определения разностного уравнения, а следовательно никак
не влияют на решение в точке О, достаточно большое возмущение. Это
приведет к как угодно большому отклонению разностного и
дифференциальных уравнений, т.е. к неустойчивости решения.
A i h i+l/2 i+1 A\
Рис. 6.1.
Докажем, что в случае волнового уравнения (6.7) матрица G нор^
мальная, и следовательно, условие Куранта (6.8) будет не только
необходимым, но и достаточным для сходимости. Действительно, легко про-

§6. Исследование устойчивости разностных схем
49
верить, что матрица
/ — coskh —j^smkh
I— g^jLsinkh —coskh
заменой переменных vq = v\a^ sq = S\ симметризуется
-coskh — ^ sin kh
^JA sin kh — cos kh
Gl ,__
h
и очевидно, что Gi GJ = G^Gi
Исследуем теперь устойчивость разностного уравнения (5.4),
которые отличаются от (5.3) только тем, что в разностной производной по t
используется значение в средней точке г на п-м слое, а не полусумма в
точках г + 1 и г — 1, как в (5.3), тогда
т ъ
£о(\ — 1) — Щ-г smkh = О
г;о(Л — 1) — 6q a1— sin kh = 0
aV (6'9)
det G = (Л - l)2 + —— sin2 kh = 0
hl
CLT
A = 1 ± i — sin /c/i
Откуда видно, что условие Неймана (6.6) нарушается, |Л| > 1 и схема
неустойчива.
Схема Лакса-Вендроффа. Рассмотрим устойчивость схемы (5.5). В
эти уравнения входит разностное представление второй производной.
Поскольку оно в дальнейшем часто встречается, введем для него отдельное
обозначение
А -> ит+\ — 2um + um-i
Aum =
Фурье-изображение для него будет
* Л eikh_2 + e-ikh 2(l-coskh)^ 4Л . 7kh ,
AUm = Um — = — Um = —Г2и™ Sin — (6.10)

50
Гл.1 Понятия о разностных схемах
Тогда для (5.5) находим
(А-1)
кг sin kh) vq +
к2 . 2kh
— sin -—
2 2
vq — (кг sin kh) £q = 0
X2 . 9 fc/i
-— Sir -—
2 2
(A-l)
Приравнивая нулю определитель, получаем
Ai,2 = 1 — х2(1 — cos а) ± г к sin a
£о = 0
Ai?2 = 1 - 2х sin^ — ± 2 г к sin
|A|2 = l-4x2(l-x2)sin4^l
• а (л • 2а\
1П — 1 — Sin —
2 V 27
kh
a\i/2
/i2
(6.11)
(6.12
При изменении 0 ^ a ^ 27Г в комплексной плоскости |А| описывает
эллипс, который лежит внутри единичной окружности |А| = 1, при к < 1,
где к — число Куранта. При к = 1 эллипс совпадает с единичной
окружностью, это свидетельствует о том, что схема при к = 1 недиссипатив-
на, т.е. амплитуда каждой Фурье-компоненты сохраняется точно. Здесь
также матрица G G* = G*G и условие Куранта достаточно для
устойчивости.
Устойчивость уравнения теплопроводности. Перейдем к
исследованию устойчивости разностной схемы (5.9) для уравнения
теплопроводности (5.8) в безразмерных переменных
-*- rn -L- гп
£ЛТ£ + (1 _ £)ЛТ,
П+1
(6.13
Подставляя в это уравнение решение в виде
та = 0П+1ехр(г^/3)

§6. Исследование устойчивости разностных схем
51
и используя формулу (6.10), получим
4т
А
4т г . 9 а . ^ . 9 ал
2т _ 1 £ 2 2 . Г • 9 а
Л[1 + (1 - С)Р ] = 1 - СР где р2 = -^ sin2 - (б.14]|
1-&>2+р2-р2 = р2
1 + (1-ф2 1 + (1"Ф2
Из (6.14) следует, что условие Неймана будет выполнено, если
0 * Tuh)? *2 Ч
Неравенства (6.15) должны выполняться при любом р, чтобы не было
ограничения на шаг по времени т. Левое неравенство выполняется при
любом 0 < £ < 1, а правое дает ограничение на величину (1 — £)|
2р2(1 - О + 2 ^ р2 (6.16)|
Откуда следует, что должно выполняться условие £ ^ |.|
Таким образом, схема безусловно или абсолютно устойчива при
£ ^ | и лишь условно устойчива, когда £ > |. Например, для явной
схемы при £ = 1
л 1 4г • 2 а
X = 1~¥Sm 2
Откуда следует, что для выполнения условия Неймана необходимо,
чтобы в размерных переменных т ^ |-. Т.е. на шаг по t следует
существенное ограничение, которое приводит к слишком малому шагу и явная
разностная схема для уравнения теплопроводности становится
неэффективной. С другой стороны, хотя неявная схема (6.13) при £ ^ 1/2 не дает
никаких ограничений на шаг по времени, она приводит к необходимости
решения на каждом шаге по времени системы алгебраических уравне-
Принцип замороженных коэффициентов. Приведенный
спектральный метод был развит для исследования линейных уравнений с

52
Гл.1 Понятия о разностных схемах
постоянными коэффициентами. Однако он оказывается полезным и при
исследовании устойчивости гораздо более широкого класса задач для
линейных и нелинейных уравнений. В нелинейных уравнениях или
линейных уравнениях с переменными коэффициентами для исследования
устойчивости разностного уравнения следует действовать по
следующему правилу: все входящие в уравнение коэффициенты, зависящие omx,t
и функции и, полагают постоянными или, как говорят,
"замороженными". Тогда уравнение становится линейным с постоянными
коэффициентами и к нему применяют спектральный признак устойчивости.
При этом условие устойчивости будет зависеть от замороженных
коэффициентов и следовательно от x,t и и. Шаг по времени г необходимо
выбрать таким, чтобы условие устойчивости соблюдалось при всех
значениях коэффициентов, которые используются в расчете.
Например, рассмотрим нелинейное уравнение теплопроводности с
коэффициентом теплопроводности, зависящим от координат и
температуры я = х(£, ж, Т)
дТ д ( f ^дТ\
Исследование проводится при х(£, х,Т) = x(^,xJ^,T^) так же как и в
предыдущем примере, считая, что к постоянно и зависит от
переменных £, х и Т как от параметров.
Исследуем на устойчивости явную схему. Функция q(t,x,u) в пра^
вой части уравнения не влияет на устойчивость, поэтому ее можно
отбросить. Тогда при фиксированном к получим
rpn+l _грп '
гг- h2 \ rn—1 L1m ' 1m—l
Условие устойчивости будет выполняться при ограничении
Т ^ 2x(t,x,T) ^ 2 maxn(tn,x,T) ^ ' *
(х,Т)
Это правило называется принципом "замороженных" коэффициентов.

§6. Исследование устойчивости разностных схем
53
Таким образом при расчете на каждом слое tn шаг т будет, вообще
говоря, переменным и определяется maxx(^,x, Т). Это может сильно
ограничивать т, если к(х, Т) принимает большое значение на небольшом
участке, в то время как во всем остальном стержне к мало и позволяет
вести расчет с существенно большим шагом, чем это диктуется
условием (6.17). Поэтому понятно желание получить явную, но безусловно
устойчивую схему. Такая схема была бы наиболее эффективной.
Возникает вопрос, можно ли построить такую разностную схему?
Аппроксимируем уравнение теплопроводности (5.8) следующей
трехслойной схемой
гПП — 1
грП+1
1m+l
1m+l
т
71+1 .
_ rpn — l _i_ rpn
2r - A* №\
Эта схема использует крестообразный шаблон, ее особенность состоит в
необычной аппроксимации второй производной в правой части
уравнения. Вместо обычного представления второй производной в виде|
1т+1 '
h2
взято выражение, в котором второй член 2TJ^ заменен на (Z^+1 + Z^_1)
Исследуем устойчивость схемы (6.18). Подставляя
T^Xn+1 ехр(г khm), находим
грП+1
1т+1
Л-
2т
А
+ е"
Л2(1 + q)- 2\q cosa + (1 - q) = О
2т
Ai
q cos
aiv/T^
qz sin a
1 + 9
Проанализируем полученное выражение
I q cos a±iq sin a
если q2 sin2 a > 1
если q2 sin2 a < 1
|A|<
|A|<
i + q
q cos a ± 1
1 + 9
l + «
i + q
<i
< i

54
Гл.1 Понятия о разностных схемах
Отсюда следует, что явная схема (6.18) абсолютно устойчива. Однако
оказывается, что несмотря на то, что устойчивость не налагает
никаких условий на шаг по времени, эти ограничения следуют из условия
аппроксимации. Действительно, проверим аппроксимацию правой части
уравнения (6.18)
Т(£, х + К) - T(t + т, х) - T(t - г, х) + Г(£, х + К)
г2 _2,
= Тхх - Ttt- + 0(к
Откуда ясно, что для аппроксимации необходимо, чтобы выполнялось
условие г ~ o(h), иначе при т ~ h (6.18) будет аппроксимировать не
уравнение теплопроводности, а телеграфное уравнение
гиперболического типа, содержащее вторую производную по времени
Tt + Ти^ - Тхх = 0(т2 + h2)
Таким образом, хотя разностная схема (6.18) абсолютно
устойчива, но она аппроксимирует уравнение (5.8) условно. Поэтому она так
же, как и явная условно устойчивая схема, неэффективна, т.к.
требуется очень малый шаг по времени т. Для уравнения теплопроводности не
удается построить эффективную явную схему и приходится
пользоваться неявными схемами. Наиболее часто используют простейшую неявную
схему (6.13) при £ = 0, она имеет первый порядок аппроксимации по
времени и второй по пространству. При £ = 1/2 схема использует
шеститочечный шаблон и, как легко видеть, имеет второй порядок по обеим
переменным.
Условие устойчивости схемы по спектральному признаку вообще]
конечно, не гарантирует устойчивости реального расчета, но оно
является необходимым условием, говорящим в пользу устойчивости.
Определенные неприятности здесь могут внести нелинейность при
"замороженных" коэффициентах, и особенно, аппроксимация краевых условий,
которые в спектральном методе не учитываются.

§6. Исследование устойчивости разностных схем
55
Спектральный метод был предназначен для исследования
устойчивости решения задачи Коши для уравнения в частных производных.
Однако большинство задач механики сплошных сред это краевые или
начально-краевые задачи. Поэтому возникает вопрос о влиянии
аппроксимации граничных условий на устойчивость решения. Исследование
этого вопроса в общем виде достаточно сложно. Здесь мы
ограничимся простейшим практическим способом определения устойчивости.
Устойчивость при решении краевых задач. Этот способ
заключается в том, что помимо стандартного исследования спектральной
устойчивости, исследуется устойчивость с помощью того же метода на
границах области интегрирования.
Рассмотрим, например, явную схему уравнения теплопроводности
с краевым условием
дТ
х = 0: —-/ЗТ = 0
дх
аппроксимация которого будет
1i — io
h
f3T™ = 0 (1 + h(3)T£ - T™ = 0 (6.19)
Решение (6.19) ищем в виде Г = T°\neika.
Подставляя получаем
eia = 1 + /3h а=- ln(l + (ЗК) = -i/3h + 0(h2)\
г
Вычислим для найденного а точку спектра Л(а). Для явной схемы
из (6.14) при £ = 1 получаем
4тг 9 а 4тг ff32h2
Л — 1 = — -пг sin — = — -гт
+ 0(h2)\ =т(р2 + ОЩ
К1 2 К2 V 4
\=l + rf52 + 0{h2) |Л| <1 + 0(т)
т.е. схема устойчива при левом краевом условии.

56
Гл.1 Понятия о разностных схемах
Аналогично рассматривается правое краевое условие. Его можно
проверять для левой полупрямой, полагая к = 0, — 1, — 2..., тогда
а = ij3h и Л остается прежним: Л = 1 + rf32 + 0(h2)\
Более полная неформальная теория устойчивости краевых задач
достаточно сложна и не будет здесь рассматриваться, см. [16].
Решение краевой задачи для системы уравнений (6.13) приводит
при фиксированном п к решению системы алгебраических уравнений
с трехдиагональной матрицей, которая эффективно решается методом
прогонки. Этот метод получил широкое, распространение и имеет очень
важное значение для вычислительной математики. Многие краевые
задачи, решаемые разностными методами, приводятся к алгебраическим
системам с матрицей, близкой к матрице с диагональной структурой,
которые эффективно решаются методом прогонки. Его основная идея
имеет различные обобщения. Известны скалярная, векторная и матричная
прогонки. Поэтому этот метод будет рассмотрен позже детально, после
того как будут разобраны общие методы решения разностных уравнений.
В заключение обсудим некоторые практически важные вопросы,
связанные со спектральным анализом устойчивости.
Выбор шага при решении уравнения теплопроводности по
неявной схеме. Для устойчивого счета по явной схеме условие Куранта
дает ограничение на шаг т ^ 0.5/г2. В то же время неявная схема не
накладывает никаких ограничений на шаг по времени т. Возникает
вопрос, какой же шаг следует использовать в реальном расчете по неявной
схеме? Здесь ограничение накладывает уже не условие устойчивости, а
условие точности расчета. Для того, чтобы ответить на поставленный
вопрос необходимо провести анализ точных решений дифференциального
и разностных уравнений, представленных в виде рядов Фурье.
Точное решение дифференциальной задачи (5.8) имеет следующий
вид
T = Y^ Cle~Xkt sin ктгх \к = к2тг2 (6.20)|
к

§6. Исследование устойчивости разностных схем
57
Решение разностного уравнения (6.13) при £ = 1
7£ = X)C?sin(bmft)
к
Подставляя его в (6.13), получим уравнение для определения С%\
СI - Q _ 4 2 knh +1
^ = C°(l + 4jsin2^y" (6.21j
Откуда видно, что С^ зависит от г и /г, а следовательно от t и ж, в то
время, как в точном решении С% зависит только от t и не зависит от х
С% = С°ке-к27тН
Из формулы (6.21), следует что выражение для С% не будет зависеть
от /г, а следовательно от ж, только при |/с7г/г| <С 1, тогда (6.21) можно
переписать так
1 + 4— sm2 —— = 1 + к2тт2т tte~k7Tt
hz 2
при к2тт2т <С 1 получаем выражение близкое к точному решению.
Оно показывает, что для реальной аппроксимации h следует вы^
бирать так, чтобы выполнялось условие \кттк\ <С 1, т.е. h <C К С другой
стороны из решения разностного уравнения следует, что шаг по
времени т должен удовлетворять условию тк2тт2 <С 1 для того чтобы решение
мало отличалось от точного. Тогда т « 0(h2). Например, если кп = 100,
h <С 10~2, т <С 10~4, т.е. соотношение такое же как из условия
устойчивости для явной схемы. Это условие является естественным для задачи
теплопроводности. Исключение составляют очень гладкие решения,
когда к « 1, тогда можно брать т & h. При такой гладкости целесообразно
пользоваться неявной схемой.

58
Гл.1 Понятия о разностных схемах
Выбор шага при решении волнового уравнения. Для того чтобы
определить каким должен быть шаг по времени для неявной абсолютно
устойчивой схемы при решении волновых уравнений нет необходимости
в сравнении точных решений дифференциального и разностного
уравнений. Достаточно вспомнить, что область зависимости точного
дифференциального уравнения определяется характеристиками и т = h/c
(рис. 6.1). Поэтому условие устойчивости Куранта для явной схемы
является одновременно ограничением на шаг г и по условию точности,
если точность оценивается в метрике Ch- Если же интерес представляют
только интегральные характеристики движения, то точность решения
оценивается в метрике L2h- Тогда шаг может быть взят больше т > h/c.
§7. Упражнения
1. Получить разностное представление третьей производной на 4хто-
чечном шаблоне равномерной сетки с номерами (г — 2, г — 1, г, г + 1).|
Определить порядок аппроксимации, используя
_ fu _ <l_ /<Ри\
dx3 dx \dx2 J
и разностную формулу для второй производной.
2. Получить разностное представление для четвертой производной на
бтиточечном шаблоне (г — 2, г — 1, г, г + 1, г + 2), используя полу^
ченное в задаче 1 разностное представление для третьей производи
ной.
3. Получить разностное представление для четвертой производной
используя представление
uiv = ^_ (#У>
х dx2 \dx2
и разностную формулу для второй производной на пятиточечном
(г — 2, г — 1, г, г + 1, г + 2) и семиточечном (г — 2, г — 2, г — 1, г, г + 1,|
г + 2, г + 3) шаблонах. Сравнить порядки аппроксимации.

§7. Упражнения
59
4. Записать разностный оператор, аппроксимирующий уравнение
Пуассона на прямоугольной сетке с шагом h по х и шагом Н по у\
Доказать, что порядок аппроксимации полученной схемы второй]
д и д и
L{u) = d^ + W = f{x'y)
Т2 (ui+l,j ~ %ui,j + Ui-l,j) + Тр (ui,j+l ~ %ui,j + ui,j-l) = fi,}\
L(u) + 0(h2)+0(H2) = f(x,y)
5. Скомбинировать расчет на вложенных сетках с шагом h и h/2 по
схеме первого порядка аппроксимации, чтобы повысить порядок ап^
проксимации до 0{Ь?) (формулы Ричардсона).
6. Определить устойчивость разностных схем
j ha
Уг+1 = У + "у [Уг + Уг+1)
■ + а у i = О
для уравнения
2h
у' + ау = 0 у(0) = уо
отыскивая точное решение разностных уравнений в виде yi = X]
(см. §2).
7. Получить формулу Адамса третьего порядка аппроксимации.
Использовать базисные интерполяционные полиномы Лагранжа
_ (t-T)(t-2T)...(t-kr))
к[) (t-ртМрт)
где uj(i) — выражение, стоящее в числителе после сокращения на
множитель (t — рт). 1р(гт) = 5{Р, г = 1... к, р < к, 5ip — символ Кро-|
некера.

60
Гл.1 Понятия о разностных схемах
8. Получить разностные уравнения метода Адамса третьего порядка
точности и формулы для вычисления первых трех значений функ^
ции щ (г = 0,1, 2), необходимых для начала счета по схеме третьего
порядка. Применить полученную схему для решения уравнения^
du
~dt
= u2 + t2
9. Получить схему Рунге-Кутты третьего порядка точности методом
удвоенного предиктор-корректора. Применить полученную схем^
для решения уравнения
du
~dt
t2 + их + и2
10. Привести телеграфное уравнение
д2и д2и
к системе трех уравнений первого порядка и исследовать для нее
устойчивость разностной схемы Лакса.
11. Исследовать устойчивость явной схемы для параболического
уравнения с постоянными коэффициентами
dt дх2 дх
12. Исследовать методом замороженных коэффициентах устойчивость
схемы Лакса-Вендроффа для нелинейного волнового уравнения]
ди 2 /ди\ д2и
dt2 \дх) дх2
предварительно перейдя к системе двух уравнений
13. Исследовать устойчивость разностной схемы крест телеграфного
уравнения.
д2и д2и

§7. Упражнения
61
14. Исследовать устойчивость явной разностной схемы методом
замороженных коэффициентов для уравнения
у' + {ку2)у = Ъх
где к и Ь — постоянные.
15. Определить порядок аппроксимации и исследовать устойчивость
разностной схемы для системы волновых уравнений
dt дх ) At Ax
де dv I el+l _ £v v'^ - vv
dt дх у At Ax
16. Получить систему двух уравнений первого порядка для стержня из
упруговязкого материала, описываемого уравнением
да де 1
т=Ет-та
Здесь Е — модуль упругости, т — время релаксации.
Выписать схему второго порядка Лакса-Вендроффа.
17. Получить для системы уравнений упруговязкого стержня (см.
упражнение 16) явную разностную схему на трехточечном шаблоне
"уголок". Определить порядок аппроксимации. Исследовать
устойчивость.
18. Для одномерного движения вязкой жидкости в канале с
определяющим уравнением
да 1
получить неявную схему на бтиточечном шаблоне (рис 5.26).
Исследовать устойчивость.

62 Гл.1 Понятия о разностных схемах

Глава 2
Методы решения систем
алгебраических уравнений
§8. Норма и число обусловленности матрицы
Относительная ошибка решения при возмущении правых частей. Число
обусловленности матрицы (63). Относительная ошибка решения при возмущении
коэффициентов матрицы (65). Пример (66).
Прежде чем рассматривать конкретные методы решения систем ал^
гебраических уравнений исследуем некоторые общие вопросы, связанные
с их решением.
При решении уравнения Ах = Ь возникает вопрос об устойчивости
решения по отношению к малым возмущениям коэффициентов матрицы
А + 8 А или правых частей Ъ + 5Ь.
Пусть А положительно определенная (Air, х) > О и
симметричная матрица А Ат = АТА.
Относительная ошибка решения при возмущении правых
частей. Число обусловленности матрицы. Рассмотрим возмущение
правой части Ъ + 5Ъ и определим какая ошибка возникает при этом в
решении х + 8х
А(х + 5х) = Ъ + 5Ъ А5х = 5Ь 5х = А~г5Ъ\
\\Sx\\ = \\A-1Sb\\^\\A-1\\\\Sb\\
63

64
Гл.2 Методы решения систем алгебраических уравнений
Если норма ||А-11| велика, тогда и \\5х\\ будет велико. Наибольшее
увеличение длины вектора будет в направлении наибольшего собственного
вектора A_1^i = Amaxxi, т.е. тогда когда вектор 5Ь направлен вдоль
наибольшего собственного вектора х\.
5Ъ = ах\, 5х = аА~гх\ = —
м
где 0 < Лх ^ Л2 ... ^ Лп собственные числа матрицы А, а Л^ = 1/Л^
собственные числа обратной матрицы А-1.
Если Л^ —► 0, т.е. матрица А близка к вырожденной, то 5х —> оо. Но
важны относительные, а не абсолютные ошибки, т.е. требуется оценить
относительную ошибку 4гЗ по сравнению ^-А Наихудший случай, когда
ошибка ||&с|| максимальна, а норма \\х\\ минимальна. Последнее будет,
когда х лежит в другом экстремальном направлении, соответствующему
минимальному Хп
х = А_16 Ъ = afn х = аА~
тогда
\\Sx\\ < Лтах \\6Ь\\
Ы\ " Лтт ||&||
т.е. ошибка тем больше, чем больше число С --
вается числом обусловленности матрицы А.
Если матрица несимметричная, то норма определяется как
максимальное увеличение длины вектора || Ах|| относительно его
первоначальной длины ||х|| и обозначается ||А||
ЦА||=тах»Ц (8.2J
Для несимметричной матрицы этот максимум может достигаться уже
не на максимальном собственном векторе хп и Лтах ^ ||А||, поэтому в
выражении для числа обусловленности Лтах следует заменить на ||А||,
a Amin на || А"11|
С=||А||||А-1|| (8.3)|
1- _ Л - _ "
?п — &Лп£п — д
-^maxl
(8.1)
= Amax/Amin, которое назы-

§8. Норма и число обусловленности матрицы
65
Для несимметричных матриц систему уравнений Ах = b можно
заменить симметризованной АтАж = АТЬ = Ь\ и найти число
обусловленности через собственные числа. Но следует отметить, что такая
симметризация растягивает спектр и увеличивает число обусловленности.
Из вида (8.1) следует, что чем более растянут спектр, т.е. чем
больше отношение Amax/Amin = С, тем менее устойчиво решение к ошибкам
округления и другим возмущениям.
Относительная ошибка решения при возмущении
коэффициентов матрицы. Рассмотрим возмущение коэффициентов матрицы
А + 5 А и определим относительную ошибку 4jrP- по сравнению с
величиной \xw
(А + 5 А) (х + 5х) = Ъ А5х = -5А(х + 6х)
5х = -А~г5А(х + 5х)
\\5х\\ <||A-1||||^A||||f+ 5x\\ <||А-1|
WSAW
ы\
с
\5А\\
(8.4
и в этом случае относительная ошибка пропорциональна числу
обусловленности матрицы А.
Приведем примеры хорошо и плохо обусловленных матриц\
1 1
1 1.0001
0.0001 1
1 1
Покажем, что А1 — плохо обусловленная, а А11 — хорошо
обусловленная матрицы.
Собственные числа матрицы определяются из уравнения
A2-2.0001A + 10"4 = 0
Aij2 = 1.00005 ± у^-ООООб)2 - Ю"4
Откуда С1
10-
2-104

66
Гл.2 Методы решения систем алгебраических уравнений
На примере решения системы убедимся в ее неустойчивости^
Атх ■
2
2.0001
для возмущенной правой части 5Ь = (0,10 4)т находим, что решение
изменится на величину 0(1)
АТХ :
2.0002
Собственные числа матрицы А11 удовлетворяют уравнению
Л2-1.0001 Л-0.9999 = 0
|1 + л/5|
То же самое возмущение правой части на 5Ь = (0,10 4)т дает
возмущение решения 5х = (0,10~4)т • С7/, т.е. того же порядка.
Пример. Интересно оценить число обусловленности для трехдиаго-
нальной матрицы, которая получается при решении простейшей 2-х
точечной краевой задачи для уравнения второго порядка вида
u"{x) = f(x\ u(0) = 0
u\l) = 0
При разбиении отрезка [0,1] на п частей матрица А разностной системы
будет трехдиагональной размерности п х п
I 2 -1 0 ••• 0 0 0 0 \
-1 2 -1 ••• 0 0 0 0
0 0 0
0 0 0
\ о о о
-12-10
0-12-1
0 0-12/

§8. Норма и число обусловленности матрицы
67
Правая части f(x) ~ 0(1), т.е. ||2>|| ~ 0(1).
Можно показать, что для матрицы А собственные числа будут
АШах = 4, Amin = тт2/п2. Тогда ошибке округления в правой части
равной дЬ = 10~9 при п = I2, относительная ошибка в решении на
основании (8.1) будет ||&е||/||х|| ~ 10~5, но уже при п = 104 получаем, что
||&г||/||:с|| « Ю-1, т.е. излишнее увеличение числа разбиений может
привести к потере точности за счет ошибок округления. Если же взять
матрицу А для уравнения uIV(x) = /(ж), то у нее Amin ~ ^4 и ошибки
округления скажутся уже при п = 102. Этот пример показывает, что число
разбиений брать с "запасом" далеко небезобидно, это может привести к
плохо обусловленной системе алгебраических уравнений.
Регуляризация плохообусловленной системы уравнений. При
решении плохообусловленных систем уравнений имеет смысл проводить
их регуляризацию, которая заключается в малом возмущении
коэффициентов матрицы Ai = А + аЕ, где Е — единичная матрица, А —
симметричная положительно определенная матрица и а > 0. Собственные
значения Л$ матрицы Ai равны Л^ + а, где Л^ собственные числа
матрицы А.
Решение регуляризованной системы уравнений можно представить
в виде разложения по собственным векторам ё\ матрицы AJ
E(b,ei) Ъе{
ciei, где а = —— = ——
. Ai Х{ + а
Если а выбрать таким, чтобы выполнялось условие Атах ^> а >> Атц
(такое неравенство для плохообусловленной матрицы А выполнимо), то
члены отличающиеся малым А^ уже не будут вносить существенного
изменения в решение при возмущении правых частей Ь + 5Ь, которые имеют
место для исходной матрицы А. В то же время добавление а в члены,
отвечающие большим А^ будут несущественны и решение может
оказаться вполне приемлемой точности. Оптимальный выбор а зависит от
конкретной задачи и должен осуществляться экспериментальным подбором
и сравнением результатов при разных а.

68 Гл.2 Методы решения систем алгебраических уравнений
§9. Прямые методы решения систем линейных
уравнений
Регуляризация плохообусловленной системы уравнений (67). Метод исключения
Гаусса. Факторизация матрицы (68). Метод Гаусса с выбором главного
элемента (70). Разложение Холецкого. Метод квадратного корня (71).
Метод исключения Гаусса. Факторизация матрицы.
Рассмотрим простейшие прямые методы решения системы алгебраических
уравнений
Ах = Ь (9.1)|
А — квадратная матрица п х п с вещественными коэффициентами,
clet А^О, Ъ — вещественный вектор.
Метод Гаусса заключается в последовательном исключении
неизвестных: сначала исключается х\ из (п — 1) уравнений, затем х2
из (п — 2) оставшихся уравнений и т.д. Исключение производится на
каждом следующем шаге из системы на единицу меньшего порядка, чем
на предыдущем. В результате система (9.1) преобразуется к системе,
матрица которой имеет верхнетреугольный вид
Ш=6,
Un U12
О U22
О О
и1п
U2n
U пп
Xi
(9.2
здесь вектор Ъ\ получен описанным преобразованием из вектора правых
частей Ь исходной системы 9.1.
Следующий шаг состоит в обратном порядке решения: из
последнего уравнения находится хП1 подставляется в (п — 1)-е уравнение, затем
находится хп-\ и т.д. Обращение верхнетреугольной матрицы
производится элементарно.
Весь алгоритм можно записать в терминах преобразований матриц.
Действительно, рассмотрим процесс обратный исключению при прямом

§9. Прямые методы решения линейных систем
ходе. Это позволит получить матрицу А умножением матрицы L на
матрицу U и получить алгоритм вычисления их элементов Uik и L^j.
Умножаем (п — 1)-е уравнение на Ln,n-i и складываем с n-м, затем
умножаем (п — 2)-е на Ьщп_<2 и на Ln_iin_2, и складываем с п-м и (п — 1)-м,
(п — 3)-е на ЬП;П_з, Ьп_1)П_з и Ьп_2,п-з и складываем соответственно с
n-м, (п — 1)-м и (п — 2)-м и т.д. Такой алгоритм совпадает с тем, что
получается при перемножении матриц L на U. В результате получаем
матрицу А
LU = A (9.3)|
/ 1 0 0 ••• 0 \
1/21 1 0 ••• О
^31 -^32 1 * * * О
\Аг1 ^гг2 ^пЗ ' ' "
«12 «13
О /722 и™ -" U,
2п
\ О О О
U пп/
О а22 «23
\ О О О
аы\\
«2п
Важным моментом является, что матрица L получилась
нижнетреугольной с единичной диагональю.
Рекуррентные формулы для определения матриц L и U
записываются в следующем виде
Un = ап
U\j = aij, Lji =
г-l
Uи = ац — у ^ LikUki
k=i
г-1
Uij = C^ij / J J^ikUkj
J = 2,3.
г = 2,3.
г = 2,3.
. .п
•Н
•Н

70
Гл.2 Методы решения систем алгебраических уравнений
Lji = — \aJ% - ^2 LjkUb j = {г + 1), 3 ... n\
Ult \ k=l J
Таким образом Ax = L\Jx = Lbi = b
Итак, процедура Гаусса — это расщепление (факторизация) мат]
рицы А на произведение нижнетреуголъной матрицы L на
верхнетреугольную U. Обе матрицы легко обращаются: после прямого хода
обращаем матрицу L и находим L_16 = &i, а при обратном ходе обращаем U
и находим х = U-1&i. Таким образом, исходная задача редуцируется на
две более простые: обращение верхне- и нижнетреугольной матриц U
и L.
Отметим, что часто при решении конкретных задач механики
сплошных сред меняются лишь правые части, а матрица А остается
неизменной — изменяется внешняя нагрузка, а уравнения и область, в
которой ищутся решения, остаются теми же. Тогда L и U обращаются
один раз и хранятся в памяти ЭВМ и решение отдельной задачи сводится
к перемножению матрицы на вектор.
Метод Гаусса с выбором главного элемента. Для того, чтобы
процесс Гаусса был бы осуществим, необходимо, чтобы все миноры j-ro
порядка в левом верхнем углу матрицы не обращались бы в нуль |Aj| ^ 0
j = 1... п. Процесс, например, нельзя провести для матрицы]
-CD
т.к. |Ai| = 0, хотя система с такой матрицей имеет очевидное решение,
но достаточно переставить строки и |Aj| ^ 0 {j = 1, 2). Более того, если
|Aj| = £, где е мало, то в процессе исключения приходится делить на
малую величину, в результате ошибок округления можно потерять
истинную информацию об исходных коэффициентах а^-, причем это может
быть вызвано не плохой обусловленностью самой матрицы А, а плохо

§9. Прямые методы решения линейных систем
71
выбранным алгоритмом. Например, для матрицы
вначале необходимо поменять местами строки, а затем проводить
процедуру Гаусса. Т.е. прежде чем производить исключение какого-то к-ого
неизвестного, необходимо провести сначала выбор главного элемента
среди всех элементов к-ого столбца, и то уравнение, в котором
коэффициент по модулю окажется наибольшим, поставить в первую строку.
Тогда Ujj > 0 будет максимальным и алгоритм Гаусса станет настолько
устойчивым, насколько это позволяет матрица А. Этот метод называ-
ютметодом Гаусса с выбором главного элемента.
Разложение Холецкого. Метод квадратного корня. Если
матрица А положительно определенная, как при решении задачи теории
упругость методом конечного элемента, то
det A = det L det U = Un • U22 •
т.к. все Ujj > 0, то процесс всегда осуществим. Кроме того, он не требует
перестановки строк, если А симметрична. Действительно, в этом случае
разложение А = LU запишем так А = LD(D_1U), где D
диагональная матрица с элементами С/ц ... Ujj. Тогда получаем
■ Ujj > 0,
D_1U =
/1 ии
0 1
0 0
\0 0
7"UQ OT\/T
/ 4dcM
LDU = U1
DLT
Uis '
Щ2 '
0 •
0 •
где
ищ \
U2n
Un—\^n
■ 1 /
U = D1
Т.К. А =
откуда по симметрии следует, что U = LT. Более того, в этом случае
матрицу А можно представить как
А = LDLT = LD^D1/2!/ = £ХТ (9.4)|

72
Гл.2 Методы решения систем алгебраических уравнений
Это разложение называется разложением Холецкого. Если число
обусловленности матрицы А есть С, то для прямого хода в методе Гаусса и
обратного оно одинаково и равно С1//2, т.к. £ = А1/2. Тогда собственные
1/2
числа li матрицы £ будут /^ = Л^ .
Пример разложения Холецкого.
£ 1\ - ( l °\ (£ °
1 Оу ~ \е-1 \) \0 г'1 J \0 1
LB^t1 °)(£l/2 ° W^ °
detLD1/2 = detX = l
Хотя сама матрица А требует выбора главного элемента, после
разложения Холецкого, решение получим без выбора главного элемента]
Расчетные формулы для элементов матриц D и L имеют вид
di = sqnon, In = 0п)1/2, hj
di = sqn I an - ^(kifdk j , 1ц
l%dl \ k=l J
f- i = 2,3...rJ
( г~г V/2
an - ^2(hi)2dk
\ k=i J
i = 2, 3 .. .n
j = i + 1... n
J Id-
Вычислительные затраты в методе Холецкого примерно в двое меньше,
чем в методе Гаусса без учета симметрии матрицы А.
Метод Гаусса часто приходится применять к неполностью
заполненным матрицам. При этом в разностных и вариационно-проекционных
методах матрица имеет квазидиагональную или ленточную структуру,
когда ненулевые элементы расположены вблизи диагонали. Такие
матрицы для своего обращения требуют меньшего числа операций. Поэтому
целесообразно преобразовать матрицу к виду, имеющему минимальную
ширину ленты. При малой ширине ленты, когда число ненулевых
элементов в строке к <С N, где N — размерность матрицы, выигрыш будет

§10. Итерационные методы решения линейных систем
73
весьма большим. Более того, когда лента имеет переменную ширину,
целесообразно запоминать весь профиль ленты, т.е. первый и последний
ненулевые элементы в каждой строке. Иногда, когда матрица редко
заполненная и не удается придать ей ленточный вид, запоминают
расположение каждого ненулевого элемента и производят исключение с учетом
их расположения [4]. Во многих разностных задачах матрица А имеет
совсем простую структуру трех-, пятидиагональную, либо блочнодиаго-
нальную. Для этих случаев разработан наиболее эффективный метод
решения метод прогонки.
Ниже мы рассмотрим некоторые задачи с применением этого
метода, очень важного для решения задач в конечных разностях, а пока
перейдем к другому способу решения систем алгебраических уравнений
итерационным методам.
§10. Итерационные методы решения линейных
алгебраических систем уравнений
Одношаговые итерационные процессы (74). Итерационные процессы Зейделя и
Якоби (76). Метод установления (77). Оптимизация скорости сходимости
установившегося процесса (80). Оптимизация нестационарного процесса (82).|
Для решения систем алгебраических уравнений, получающихся в
результате дискретизации уравнений механики деформируемых сред, во
многих случаях целесообразно применять итерационные методы.
Очень важным моментом при решении уравнений итерационными
методами является начальное приближение. От того насколько хорошим
будет оно часто зависит успех всего решения задачи. При решении
эволюционных задач механики с помощью неявных схем приходится решать
системы уравнений на каждом временном шаге At. В этих задачах
решения на предыдущем шаге служат прекрасным начальным приближением
и применение метода итераций здесь весьма эффективно. При решении
стационарных задач, не зависящих от времени, такого простого способа
нахождения начального приближения нет. Но и в этих задачах физиче-

74
Гл.2 Методы решения систем алгебраических уравнений
екая интуиция и ясное понимание существа задачи позволяют получать
такое приближение. В некоторых случаях бывает целесообразно
искусственно вводить в задачу время и решать ее как эволюционную до выхода
решения на установившийся режим (метод установления).
Важным преимуществом итерационных методов по сравнению с
прямыми методами является невысокие требования к объему
используемой памяти, и хотя эта проблема в связи с колоссальным прогрессом
вычислительной техники в последние годы потеряла свою остроту, тем
не менее она существенна и в настоящее время.
Одношаговые итерационные процессы. Наиболее существенное
значение итерационные методы имеют для решения нелинейнх задач,
здесь им нет альтернативы, по-существу это единственный метод их
практического решения. Но начинать изучение итерационных методов,
конечно, следует с решения линейных уравнений.
Пусть имеем систему линейных алгебраических уравнений
Ax = f (10.1|
где А положительно полуопределенная матрица А ^ 0. Матрицу А
можно рассматривать как линейный оператор в конечномерном
эвклидовом пространстве со скалярным произведением (х, у) = xiyi и
ассоциированной с ним нормой ||а?|| = (х, х)1/2, х — гг-мерный вектор.
Если для получения каждого следующего приближения хк+1
используется только предыдущее приближение хк, то итерационный
процесс будем называть одношаговым или двухслойным (в соответствии с
отвечающей ему разностной схемой). При использовании двух
предыдущих приближений хк и хк~г процесс называется двухшаговым или
трехслойным.
Каноническая форма записи двухслойного итерационного процесса]
имеет вид
В + Ахк = f к = 0,1... (10.2)
Тк+1

§10. Итерационные методы решения линейных систем
75
Матрица В и скаляр т^ — параметры итерационного процесса, которые
выбираются так, чтобы процесс был бы наиболее эффективным.
Форма (10.2) соответствует записи в виде явной разностной схемы. Таким
образом существует тесная связь между итерационными методами и
явными разностными схемами
Сходимость итерационного процесса будем оценивать в
энергетическом пространстве Не, порожденном симметричной положительно
определенной матрицей С, со скалярным произведением (а, Ь)с = (Са, Ь) и
соответствующей ему норме ||а||с = (Са, а)1/2.
Процесс будет сходящимся, если норма ||^||с ~~> 0 ПРИ & —> оо, где
^ = хк — х, хк — к-ое приближение. В силу того, что точное
решение х уравнения (10.1) неизвестно, оценку точности проводят не по
норме Holies а по норме невязки
w^k - /1 = на
которая легко вычисляется на каждой итерации. В качестве меры
сходимости примем относительную погрешность невязки е и оценку точности
можно проводить как
iiai^nai
Такой выбор оценки точности отвечает сходимости в энергетическом
пространстве с матрицей С = АТА.
||£*|| = (Affc - /, Ахк - /)1/2 = (Ахк - Ах, Ахк - А£)1/2 =|
= (Аг*, А#)1/2 = (АТА#,#)1/2 = ||г*||с|
Итерационный процесс, задаваемый формулой (10.2) можно опти^
мизировать за счет выбора матрицы ~Вк и параметра 7>. При изменении
этих параметров от итерации к итерации итерационный процесс
называют нестационарным, в отличии от стационарного процесса, при котором
они остаются неизменными.

76
Гл.2 Методы решения систем алгебраических уравнений
Итерационные процессы Зейделя и Якоби. Рассмотрим как
будут представляться в форме (10.2) классические итерационные
процессы Якоби и Зейделя. В итерационном процессе Якоби k-ое приближение
при решении уравнений (10.1) состоит в нахождении каждой компоненты
вектора xk+l из г-ого уравнения по формуле
г—1 п
Y ai5x) + alxxk+l + Y, ачх) = fi г = 1, 2 ...n (10.3)|
j=l j=i+l
в которой все компоненты кроме Х{ берутся из предыдущего
приближения.
В методе Зейделя в формуле (10.3) в первой сумме хк заменяется]
уже вычисленными значениями xk+1 (j = 1... г — 1).
Матрицу А представим в виде разложения на верхнетреугольную
матрицу L, диагональную D и нижнетреугольную U
A=L+D+U
Тогда итерационный процесс Якоби можно представить в форме (10.2)
со следующими параметрами В = D, г = 1, а процесс Зейделя с
параметрами В = L + D, т = 1.
Обобщение итерационных процессов Якоби и Зейделя на
нестационарный случай, в котором изменяется параметр rfe, получим в виде|
D - + Axk = f
/+1+lsk Ar = 0,l...rJ
(L + D) — + Axk = f
Tk+l
Еще одно обобщение такого итерационного процесса — метод верхней
релаксации записывается в виде
Подставляя А =
xk+l =
(Lrfc+i + D)-
I A rK — f
I J\X — J
Tk+l
L + D + U, явно выразим xk+l
L + —D
-l
'/-(U+(1-^)D)*'"

§10. Итерационные методы решения линейных систем
77
и обращать нужно только верхнетреугольную матрицу.
Эти процессы в стационарном случае будут сходящимися, если
матрица перехода G от А:-ого слоя к к + 1-му будет удовлетворять условию
Неймана.
xk+1 = Gxk + В-Ч+i/, где G = E-rfc+1B-1Aj
В методе Якоби
G = E-rD -lA
В методе Зейделя
G = E-r(L + D)"1A
В обоих случаях матрица G перестановочна со своей транспонированной
и условие Неймана не только необходимое, но и достаточное условие для
сходимости, если А = Ат — симметрично сопряженная матрица.
Стационарный итерационный процесс (10.2) соответствует методу
простой итерации. Он будет сходится, если по критерию Неймана все
собственные значения А^ матрицы G удовлетворяют условию |А$| < 1.
Метод установления. Наиболее просто анализируется сходимость
итерационного процесса при В = Е. Тогда уравнение (10.2) можно в
пределе при Th+i —> 0 заменить дифференциальным уравнением
Ж + АФ = / (10.4)
при начальном условии
t = 0, Ф = 0 (10.4а)
Докажем, что решение задачи (10.1) можно получить методом
установления, решая уравнение (10.4), а затем устремляя t —► оо.|
Доказательство. Решение задачи (10.4) будем искать в виде разложения по|
собственным векторам ип матрицы А
■^■tln = Лп11п
N N
71=1 П=1

78
Гл.2 Методы решения систем алгебраических уравнений
где an(t) — коэффициенты Фурье функции Ф(£). Подставляя разложения в (10.4) и
учитывая, что an(t) = (Ф • йп), получаем
E(^ + A„a„(t) -/„)«„ = 0
П—1 ^ '
В силу линейной независимости векторов базиса ип находим N дифференциальных
уравнений с постоянными коэффициентами для определения an(t)
Из начального
тогда решение
условия (
dan _
1, т лпап jn
10.4а) получаем
Оп(0)
= 0
п = 1
(10.4) принимает следующий вид
ап = Се-х^ + ^ где
In
ап = т-
- (1 - е"^*)
с =
ф =
In
N ,
...N
из
(1-е
аг
-Хт
,(0)
'«)й
При £ —>• оо, учитывая, что для положительно определенной матрицы А > 0
все собственные числа Лп > 0, получаем стационарное решение уравнения (10.4]
N fn
lim Ф = V" -^йп
*^°° n^i n
С другой стороны решение уравнения (10.1) тоже можно представить в виде
разложения по собственным векторам йп матрицы А. Разлагая х и / и подставляя
в (10.1), получим
х = ^2 Хп^п f = 5^ ^п
п=1
N
/ J %п-г*-У"п / J AnXnUn / j Jn^n •"
Хп
Таким образом, решение системы (10.4) в пределе при t —> оо
сходится к решению уравнения (10.1). Это позволяет утверждать, что
решение разностного уравнения, которое аппроксимирует дифференциальное

§10. Итерационные методы решения линейных систем
79
уравнение (10.4) также будет сходится к решению системы (10.1).
Уравнение (Ю.4) можно записать в разностном виде
+ АФП = / п = 1,2,... (10.5)
f — не зависит от £, поэтому правая часть в (10.5) без индекса.
Решение (10.5) запишем в виде рекуррентного соотношения
Ф"+1 = (Е-Атп)ФЧтв/ (Ю.6)|
где Е — единичная матрица, п — номер итерации или номер шага
интегрирования по параметру t. Формулу (10.6) можно рассматривать как
итерационную для решения уравнения (10.1). В (10.6) входит
неопределенный параметр тп, который должен быть выбран так, чтобы процесс
итераций сходился. Очевидно, что тп должно зависеть от свойств
матрицы А.
Перепишем уравнение (10.6) относительно новой переменной £п —
невязки решения уравнения (10.1) на n-й итерации
£п = АФП - /
тогда (10.6) можно представить в виде однородного уравнения
относительно £п
фп = A_1(<fn + f) Фп+1 - Фп = A_1(fn+1 - <fn)|
Подставляя в (10.5), получаем
Г+1 = (Е-Атп)Г
Разлагая ^п по собственным векторам матрицы А и приравнивая
коэффициенты при линейно независимых векторах щ, получим
N

80
Гл.2 Методы решения систем алгебраических уравнений
Откуда ясно, что для сходимости стационарного процесса при т^ = тп
необходимо и достаточно, чтобы
max|l-Afcr| < 1, к = 1... N
А/с
Если матрица А имеет границы спектра Ai ^ An(A) ^ Адг, то для
параметра т должно выполняться неравенство
г < A (10.7J
Другими словами, при таком выборе параметра т норма оператора
перехода ||Qn|| = ||Е — Атп|| ^ 1, и он является оператором сжатия (см.§11).
Однако остается широкий выбор для т и при выполнении условия (10.7).
Возникает вопрос, как выбрать т, чтобы процесс сходился бы с
наибольшей скоростью?
Оптимизация скорости сходимости установившегося процесса.
Скорость сходимости итерационного процесса будет определяться
наибольшим собственным значением и соответствующим ему вектором гГдг.
Можно выбрать т так, чтобы полностью подавить компоненту вектора
невязки £]у-, соответствующую максимальному собственному вектору
l-XNr = 01 т = -—
An
тогда остальные компоненты вектора невязки ££ будут затухать в
итерационном процессе как
и скорость сходимости будет определяться
m&xqk = 1- — = 1- —
к Лдг О
где С — число обусловленности матрицы А. Очевидно, что чем больше С
для матрицы А, тем медленнее будет сходимость выбранного процесса.

§10. Итерационные методы решения линейных систем
81
В общем случае скорость сходимости итерационного процесса
определяется нормой матрицы перехода Q
„пи IIQril llfn+1ll
q = max,» = ||Q|| = max^p = ™|x^
Если Q не меняется от итерации к итерации, то и q не зависит от п. Если
же шаг тп зависит от п, то и q зависит от п, тогда следует вычислять
среднее значение qcp
Pn(Qn) = (Е - пА)(Е - т2А)... (Е - гпА)
^ = lim \\Pn(Qn)\\1/n = lim || (E - пА)... (Е - ТпЩ1^
т.е. нужно вычислять асимптотическую скорость сходимости. Удобнее
использовать в качестве характеристики итерационного процесса
экспоненциальную скорость сходимости
S = -1п^Ср
При плохой обусловленности матрицы А средняя скорость
сходимости qcp мало отличается от единицы qcp « 1 — а, а — малая величина,
тогда S « а. Чем меньше S, тем медленнее сходится процесс.
Выше в формуле (10.7) г было выбрано не оптимальным образом]
выбор можно улучшить, найдя оптимальную скорость сходимости. В
самом деле, откуда следует, что нужно подавлять компоненту,
отвечающую Лдг: max |1 — г Л|? Скорость сходимости при одинаковом шаге тп — т
определяется qcp = lim ||Pn(A)||1/n = max |1 — т Л|, где Ai ^ Л ^ Адг. По-
п^оо Л
этому т нужно выбрать так, чтобы этот максимум при изменении А в
указанных пределах был бы минимален.
Это задача на нахождение минимакса Р(А). Нужно найти
minmaxP(A), где Р(Х) — линейная функция на отрезке А« ^ А ^ Адг
г Л
при Р(0) = 1. Геометрическое решение этой задачи очевидно (рис. 10.1).
Необходимо, чтобы прямая Р(А) проходила через срединную точку
интервала. Тогда максимальное значение |Р(А)|, достигаемое на концах
интервала, будет минимально.

82
Гл.2 Методы решения систем алгебраических уравнений
Рис. 10.1.
На рис. 10.1 видно, что если прямая пересекает ось Л левее средин^
ной точки (пунктирная линия), то значение Р(Х) на правом конце
интервала возрастает. Если правее, то возрастает значение на левом конце.
Поскольку тахР(А) находится всегда на концах интервала, то minmax
отвечает срединная точка пересечения.
Ajv + Ai
0
qk
2Хк
min max |P(A)|
Ai + Адг
2Ai
1
Xn + Ai
n/b i
Xn + Ai
max
к
Xn — Ai
Xn + Ai
1
i_
с
i + *
С) V C + ° \C2
Скорость сходимости при больших С будет S — ^
два раза больше, чем при первом выборе.
С + °\С2
2 + 0(С-2)
Оптимизация нестационарного процесса. До сих пор параметр т
от итерации к итерации был фиксирован. Возникает вопрос, нельзя ли
увеличить скорость сходимости за счет выбора переменного шага
итерации? Итерационный процесс с параметром тп, зависящим от номера

§10. Итерационные методы решения линейных систем
83
итерации, называется неустановившимся, в отличие от установившегося
при фиксированном т.
Если тп переменная величина, тогда
п
£п = Д(Е - тг А)|° = Р»(А)^
г=1
где Рп(А) — полином n-ой степени матрицы А. Поскольку Атип = \™иП)
где Хп и ип собственные значения и собственные векторы матрицы А, то
разлагая £п и £° по собственным векторам матрицы А, получим для
коэффициентов Фурье следующие формулы
d = Pn{^m)im
где
п
Рп^ш) = 1[(1-г1Хт) (10.8)|
г=1
Так же, как и раньше, чтобы получить оптимальную скорость
сходимости, нужно минимизировать максимум этого полинома
qcp = minmax |Pn(A)| Ai < Л < \n
Эта задача может быть легко сведена к известной задаче Чебышева
найти полином степени п на интервале [—1,1] и принимающий значение
Рп(0) = 1, максимальное значение которого было бы минимальным.
Чтобы свести нашу задачу к задаче Чебышева, достаточно провести
преобразование независимой переменной А
лп -\- \\ — 2л
Хп — Ai
которое переводит отрезок [Ai, An] в отрезок [+1,-1] и необходимо
нормировать решение, чтобы Рп(0) = 1, тогда
WAn + Ai-2AA
±п \ Xn-Xi )

84
Гл.2 Методы решения систем алгебраических уравнений
где Тп(у) — обозначает полином Чебышева степени п
тп{у) = —^со&{пксссоъу)
Корни полинома Чебышева п-й степени Тп(у) вычисляются по формуле
(2г - 1)тг
Уг = cos г = 1,2...п
У 2п
тогда корни Рп(Х) будут
1 "
2
л ,л л а (2г — 1)тт'
Адг + Ai — (Ajv — Aij cos
2n
Теперь Ti легко определяются через A^ с помощью представления (10 J
п ,
РП(А) = П(1-
г=1 ^
А\ 1
. , а следовательно, Ti = —
Хо ] X,
г = 1, 2 ... п\
Адг + Ai - (XN - Ai) cos (2г2п1)7Г
Таким образом поставленная задача решена — найдена формула, по
которой определяются шаги оптимального неустановившегося итерационного
процесса.
Оценим асимптотическую скорость сходимости такого
нестационарного итерационного процесса.
Наихудшая сходимость будет при значении |maxPn(A)| = q™ ,
которое полиномы Чебышева принимают на концах интервала [—1,-1-1]:
Tn(d=l) = (±1)п, поэтому на основании формулы (10.9) находим
максимальное значение Pn(Ai) = Tn(l)/Tn(r). Тогда
1 2 2
?С1
^
|Г„(г)| (г + VT^lf + (г - V^l)n " (г + V^l)n\
где
_ An + Ai
An — Ai

§10. Итерационные методы решения линейных систем
85
Выражая г через число обусловленности С, получим
1 +с Vc*2
с
Тогда асимптотическая оценка экспоненциальной скорости сходимости
при плохой обусловленности матрицы А будет такой
S = -lngcp = — ln2 + ln(r+\/r2 + l)
=-^п2+1п(1+§+у(1+§
1, , ( 2 2 \ /1
= —ln2 + ln [l + - + -= )+o -
lim S = In f 1 + -У = -?= + о (^] (10.10)1
V Vc/ Vc Vcy v 1
Сравнивая (10.10) с формулой для стационарного процесса
5 = | + 0(С-2),
получаем, что в нестационарном процессе скорость сходимости S в
квадрат раз больше.
Таким образом выбор т\ по формулам чебышевских полиномов
теоретически дает мощное ускорение итерационного процесса. На практике
оказывается, что реализация чебышевского итерационного процесса
может не дать полученного ускорения и даже иногда привести к
расхождению. Это связано с проблемой вычислительной устойчивости, которая
решается специальным порядком выбора параметров т\. Объяснение и
алгоритм выбора приводятся в [14].
Полученные результаты по ускорению итерационного метода, в
котором невязка вычисляется по формуле
Г+1 = (Е-г&А!)^

86
Гл.2 Методы решения систем алгебраических уравнений
применимы к любому методу простой итерации. Например, для метода
Якоби матрица Ai = D_1A, для метода Зейделя Ai = (L + D)_1A.|
Во всех рассмотренных нами примерах параметр т^ не зависел от
предыдущих приближений, а только от свойств матрицы А. Такие
итерационные процессы являются линейными, в отличие от нелинейных,
когда Т{ зависит от предыдущих приближений и корректируется на каждой
итерации в зависимости от получаемого решения. Нелинейные
итерационные процессы по скорости сходимости могут превосходить
оптимальный линейный процесс (см. [20]).
В заключение обсудим вопрос о точности, до которой следует про^
водить итерационный процесс. Поскольку предполагается, что система
алгебраических уравнений, которая решается, получена
аппроксимацией дифференциальных уравнений с некоторой точностью е = 0(hk), то и
итерации следует проводить с той же точностью. Практически это
означает, что итерационный процесс должен продолжаться, по крайней мере,
до такого п, пока не выполнится неравенство
НГ+1-П<£
где £ — невязка решения уравнения (10.1).
Дальнейшие итерации большого смысла не имеют, т.к. ошибка ап^
проксимации будет превалировать над точностью итерационного
процесса.
§11. Решение нелинейных уравнений
Нелинейные уравнения и итерационные методы (86). Сжимающие отображения,
принцип неподвижной точки (87). Метод простых итераций. Достаточное условие
сходимости (89).
Нелинейные уравнения и итерационные методы. Рассмотрим
систему нелинейных уравнений, которую запишем в векторном виде
f(xi,X2...xn) = 0 или /(f) = 0 (П.1)|

§11. Решение нелинейных уравнений
87
где х = {х\.. .хп) — вектор неизвестных, f(x) — известная вектор-
функция. В частном случае, когда / = Ах + Ъ система уравнений (11.1)
будет линейной.
Нелинейные уравнения решаются итерационными методами,
которые можно рассматривать как обобщение линейных итерционных
методов, рассмотренных в предыдущем параграфе. Формула (10.2)
двухслойного итерационного процесса для нелинейного уравнения примет вид
ztn+l _ ~гп ^
Вп+1 + /(f) = 0 п = 0,1... (11.2
или в разрешенном относительно хп+1 виде
fn+1 = Р(£")
где Р(х) = х — rB~j:1/(x) — нелинейный оператор, В — квадратная
матрица п х п, имеющая обратную. Для стационарного процесса В и г не
зависят от п.
Решение нелинейного уравнения (11.1) в общем случае не един^
ственно. Выделение области, в которой (11.1) имеет один единственный
корень представляет самостоятельную проблему, которую мы будем
считать решенной. Наша задача будет состоять в нахождении итерационным
методом в выделенной области этого корня.
Определим скорость сходимости итерационного процесса.
Будем говорить, что метод сходится с линейной скоростью, если
Xх1 — х* = О (ж™ — х*). Метод имеет квадратную сходимость, если
Xй — х* = О {{х4 — х*)2). Здесь х* - корень уравнения (11.1).]
Сжимающие отображения, принцип неподвижной точки. Итак
Р(х) — нелинейный оператор отображает n-мерное евклидово
пространство Еп в себя. Всякая точка х £ Еп} удовлетворяющая условию]
Р(х)=х (11.3)
называется неподвижной точкой оператора Р. Решение (11.3) есть
решение задачи о неподвижной точке. Решение любого нелинейного
уравнения f(x) = 0 можно рассматривать как решение задачи о неподвижной

88
Гл.2 Методы решения систем алгебраических уравнений
точке
X = X + f(x) = Р(х)
Для отыскания неподвижных точек применим метод
последовательных приближений. Пусть х° — пробное решение, х1 = Р(ж°), где х1
можно принять за уточненное значение х°. Процесс уточнений можно
продолжить
хп+1 = Р(хп), га = 1,2... (11.4)|
Вводя степени оператора Р, можно записать
£П+1 = pn+l(£0) ? pn^ =p.p... P(f0j
Определение Оператор Р, отображающий пространство Еп в себя,
называется сжимающим в замкнутом шаре
R(x0,r) = (£, ||x-fo|| ^ г)
если для любых х и у £ R(xq, г) выполняется условие Липшица\
\\Р(Х)-Р(у)\\^4х-У\\ (11-51
или при у = О, ||Р(х)|| ^ А||:г||; где А < 1 коэффициент сжатия\
Теорема 1) Пусть Р(х) оператор сжатия в R(x0, r) с коэффициентом
Л< 1.
2) Пусть хо удовлетворяет условию
\\Р(х0) - х0\\ < г0(1 - X) < г
Тогда последовательность {хГ} сходится к х* £ R, их* единственная
неподвижная точка в R(xo,r).
Условие (11.2) для нулевого приближения гарантирует, что все приближения!
не выходят из R(x0,r)
||F0+1-f0|| = ||Pr+1M-^o|| =
= ||Pr+1(f0) - РГШ + РГШ - ... + Р{х0) - foil ^
< \\Рг+1(х0) - РГШ\\ + • • • + \\Р(х0) - foil <
^ Лг(1 - Л)г0 + ... (1 - Л)г0 = (1 - Лг+1)г0 < г0,|

§11. Решение нелинейных уравнений
89
где каждая разность преобразована по следующей формуле
Рг+1(х0) - Рг(х0) = Рг(Р(х0) - х0) < Рг[(1 - Л)г0] < Хг(1 - A)r<J
Доказательство сходимости основано на признаке сходимости последова-|
тельность Коши {хп}. Для любого га и п > О
Wx™ - оГ^\\ = \\(хт- ж^1 + ... + (fm+n"1 - ^+"11 <
< Н^1 - х™\\ +... +11^+" - ж"1-"-"-1!! <
^ \т(\ - Х)г0 + ... + A^+^l - A)r0 <
£ Am(l - Хп)г0 < Хтг0
Для любого е > 0 существует такое N(e), что для любых т > N(e), и
п > 0 Цж™ — жт+п|| < е. Действительно, возьмем N таким, что XN < £/V0, тогда
\\хт — хт+п\\ < е. Откуда следует сходимость последовательности {х*1} к своему
пределу х*.
Доказательство единственности. Пусть решение неединственно, тогда имеем]
два решения х* и -г\
х* = Р(£*), Л = Р{Л)
||£* - ^|| = \\Р(х*) - P(*S)\\ ^ А||£* - Л\\ 0 < А < 1
откуда следует \\х* — -г\\\ = 0.
Пусть х* = lim Pm(rro), тогда будем говорить, что х* достижима]
из xq, а множество всех точек {xq}, откуда достижима точка ж*, будем
называть областью достижимости.
Из теоремы следует, что х* достигается из любой точки шара^
R(xo,r). Действительно,
\\Р(х) - Р(х0)\\ < А||£ - foil < Ar0, \\P(x)-X!\\ < Ard
т.е. Р отображает R(xq, т*о) на R^\xi, Аго), Р2 —► R^2\x2, A2ro) и т.д., т.е.
радиус шара R^ уменьшается с каждой итерацией и {х71} —► х*\
Метод простых итераций. Достаточное условие сходимости.
Покажем, что для выполнения условия Липшица достаточно, чтобы
выполнялось условие
ЦВД1КА/7К1,

90
Гл.2 Методы решения систем алгебраических уравнений
где п — размерность ж, т.е. норма матрицы Якоби должна быть
меньше, чем Л/п, (Л < 1). Здесь для матрицы Якоби функции Р(х) введено
обозначение
(dPi дРп\
дх\ дх\
Р^
dPi ^ дРп
\дхп дхп)
а норма понимается в смысле максимум-нормы для вектора и матрицы
||£|| = max \xi\, || A|| = max \atj\
г ij
Доказательство следует из разложения по формуле Тейлора функций в точке ж0,
входящих в условие Липшица.
Pi(Xi) = Рг(х°3) + J2 (Щ) ^ fe - *?) + 0(х3 - x°f
Pt(Vi) = Р,(*°) + J2 (Щ)^ (%■ - *?) + 0(У3 - х°)2
Откуда
\Pi(Xj) - р,(й)|< J2
9Рг(х°Л
дх.
\*з-У}\
Если записать это неравенство через максимум-норму, получим
'дР^
||Д(^)-р,Ы1К
дх.
) (Xj-Vj)
м- iiP
\\х — у\\ > max
9Pi($)
dxj
^X\\x-y\
Тогда ясно, что для выполнения условия Липшица достаточным
условием будет
max
дЩф
dxj
А
где $eR(xo,r0), А<1
При п = 1 должно выполняться условие

§11. Решение нелинейных уравнений
91
Геометрическая интерпретация метода простых итераций (11.4) по^
казана на рис. 11.1а и рис. 11.16. На рисунках видно, что при |i=>/(ir) | ^ 1
процесс сходится, кроме того сходимость знакопеременная монотонная,
когда Р'{х) < 0 (рис. 11.1а). Приближения ведут себя так
х0 — х* < 0, х\ — х* > 0, х2 — х* < 0, и т.д]
Ахп меняет знак после каждой итерации. При Р'{х) > 0 в
окрестности корня х\ на рис. 11.16 сходимость знакопостоянная Ахп > 0.
При ЦР^х)!! > 1 в окрестности большего корня х\ процесс
расходится (рис. 11.16).
Х0Х2Х* Х\ X Х*ХъХг Х\ХоХ*Хъ Х\
(а) (Ь)
Рис. 11.1.

92
Гл.2 Методы решения систем алгебраических уравнений
§12. Метод Ньютона и его модификации
Метод Ньютона-Рафсона (92). Модифицированный метод Ньютона-Рафсона (94).
Метод секущих (95). Двухступенчатые итерационные методы (95).
Нестационарный метод Ньютона. Выбор оптимального шага (97).
Метод Ньютона-Рафсона. Кроме метода простых итераций,
можно построить семейство более точных итерационных методов, которые
имеют квадратичную, а не линейную скорость сходимости.
В формуле (11.2) полагаем т = 1, а в качестве матрицы В берем
матрицу Якоби
в = j = д(Д ... fn)
д(х1...хп)
Тогда получим итерационный метод Ньютона-Рафсона
£»+! = £»_ [/'(гР)Г7(^) (12.13I
В одномерном случае
ox fx
Ахп
Уравнение (12.1) можно рассматривать как
х = Р{х) =х- J~\x) f{x)
и доказывать, что Р(х) — оператор сжатия. Справедлива следующая
теорема.
Теорема Если
1) матрица (J0)-1 невырожденная, т.е. det(J°)_1 Щ 0 в шаре R(xq,t)
IKJ0)"1!! ^B, гдеВ>0;
2)^-^ С, С>0> ij = l...n
OXiOXj

§12. Метод Ньютона и его модификации
93
3) ||^-^°||<Го, Г0>В
тогда процесс сходится.
Пусть х° — начальное приближение, а х — решение системы нели^
нейных уравнений f(x) = 0. Разлагая f(x) в ряд Тейлора в окрестности
х = гг°, получим
/(f) = f(x°) + J(f°)(f - f°) + 0((f - £°)2) = б
откуда при малом Ах = x — x° ограничиваясь линейным членом,
находим
х = x0-(J°)-1(x°)f(x0)
Последовательно применяя эту формулу к найденному хП1 получаем
рекуррентную формулу
хп+1 = хп- J"1^) /(£") (12.2)1
Процесс закончится, когда выполнится условие ||Ажгг+1|| ^ е\
Доказательство квадратичной сходимости метода Ньютона-Рафсона. Запи-|
шем (12.2) в виде
^+1 = Ф(^), где ф{х) = х- (И)"1 f(x)
<?{х) = {h)-1 fm (h)-1 /(х), <j?{5?) = 0
здесь использована формула дифференцирования обратной матрицы1
(А"1)* = А^А^А"1
Тогда
Af™+1 = f^1 - х* = £(£") - х* =
= ф{х*) + 0{х*) Ах" + ф"(х*) А(^)2 - £* =
= <р"(х*)А(хт)2 + 0((хп)3) (12.3)|
Откуда получаем, что скорость сходимости в методе Ньютона
квадратичная, если в малой окрестности х* удовлетворяются условия теоремы.
1 Запись /^ или f'{x) является матрицей, полученной при дифференцировании вектора / по
вектору х

94
Гл.2 Методы решения систем алгебраических уравнений
Модифицированный метод Ньютона-Рафсона. Недостатки
метода Ньютона заключаются в том, что на каждом шаге приходится
вычислять и обращать матрицу Якоби. Кроме того, для ее получения
необходимо дифференцировать исходную систему уравнений (11.1), которая
может иметь сложный вид. Модификация метода Ньютона направлена
на минимизацию вычислительной работы, а также на увеличение
окрестности корня, в которой можно задавать начальное приближение.
Простейшей модификацией метода Ньютона является случай,
когда А = J_1(f°) не зависит от п. В этом случае интерполяционный метод
стационарный и скорость его сходимости линейная, а не квадратичная,
как при А = J-1^).
На рис. 12.1 и 12.2 показана геометрическая интерпретация
условий сходимости метода Ньютона. На рис. 12.1 условия теоремы
выполнены, а на рис. 12.2а нарушено первое условие ff(x) ^ 0, и
последовательные приближения расходятся. На рис. 12.26 первое условие выполняется
f'(x) т^ 0, но нарушено второе условие теоремы: в окрестности корня х*
вторая производная f^ меняет знак и итерационный процесс расходится.
Л, * Х>2 J\fi Дл
Хо X
Рис. 12.1.
На рис. 12.1 видно, что модифицированный метод, когда^
Л_1(ж°) = [/(^°)]-1 фиксирована, процесс сходится значительно
медленнее (пунктирные линии), чем метод с переменной касательной.

§12. Метод Ньютона и его модификации
95
№
I
Х3^ —
J4
хГ~
/х^
}х0
J** »
X
(а)
(Ь)
Рис. 12.2.
Метод секущих. Другая модификация — метод секущих получается
заменой обратных производных в (12.1) разделенной разностью
хп+1 = хп- _ ^-^ /(£*) п = 0,1... (12.4)
Этот метод является двухшаговым итерационным методом. Для
получения приближения жп+1 необходимо использовать два предыдущих х71
и f"_1. Чтобы начать итерационный процесс на первом шаге необходимо
найти х1 с помощью одношагового итерационного процесса, например
модифицированного метода Ньютона.
На рис. 12.3 показана сходимость итерационного метода секущих.
Этот метод сходится медленнее, чем метод Ньютона, но быстрее чем
модифицированный метод.
Двухступенчатые итерационные методы. Метод Ньютона —
метод линеаризации нелинейной системы уравнений. Он сводит решение
нелинейного уравнения к многократному решению линейных уравнений,
т.е. обращению матрицы для линеаризованных уравнений. Линейную
задачу можно решать в свою очередь итерационными методами, используя
например методы ускоренной сходимости, рассмотренные выше в §10.
Тогда полный итерационный цикл решения нелинейного уравнения
становится двухступенчатым и состоит из внешнего цикла (заключается в
применении метода Ньютона к исходной систему уравнений, т.е. в
сведении решения нелинейной задачи к решению линейной) и внутреннего

96
Гл.2 Методы решения систем алгебраических уравнений
Рис. 12.3.
цикла решения линейной задачи, например неустановившемся методом
Зейделя с ускорением.
Для решения нелинейных уравнений можно непосредственно при^
менятв обобщение стандартнвгх линейнвгх итерационнв1х методов на
нелинейный случай. Например, нелинейный метод Зейделя будет
выглядеть следующим образом
Л(я?
п+1 „п+1
Li+li ' ' ' э хк)
о,
1,2... к
тогда
fl(xn1+\xn2,xl...,xnk)=0, Мхп1+1,хп2+1,х1...,хпк) = 0, ..J
На каждом шаге решается нелинейное уравнение с одним
неизвестным к = 1. На первом шаге относительно х7^1, на втором
относительно ^2+1. Эти п уравнений, каждое из которых содержит только по одному
неизвестному решаются методом Ньютона, сводя задачу к линейной. В
этом случае имеем тоже двухступенчатый метод, в котором внешние
итерации осуществляются нелинейным методом Зейделя, а внутренние —
методом Ньютона для системы функций, зависящих от одной
переменной.

§12. Метод Ньютона и его модификации
97
Ясно, что возможны и другие двухступенчатые методы с различи
ными сочетаниями внешних и внутренних итерационных процессов^
Нестационарный метод Ньютона. Выбор оптимального шага.
Пусть необходимо найти решение
/(f) = 0 (12.5)1
Рассмотрим семейство уравнений Ньютона с параметром г и проведем
оптимизацию процесса последовательных приближений по этому
параметру.
Решение уравнения (12.5) итерационным методом Ньютона с пара^
метром т можно трактовать как решение разностного уравнения^
Д£»+1 = -т/г1^) /(f), Х ~Х = -/^(f) /(f) (12.б|
Непрерывный аналог метода Ньютона заключается в том, что
вместо (12.2) берем уравнение (12.6) при т —> 0. Это матричное
дифференциальное уравнение легко интегрируется
% = -fs\^)m in/(f) = -*+c
f = f°, t = t°, /(f) = /(fV^0 (12.7J
f j .g/. . 9(/i.../fc)
df 9(xi...Xfe)
Легко видеть, что решение (12.7) при любом гг° сходится к решению (12.5)
при t —>• оо.
lim /(ж) = 0, lim ж(£) = х*
т.е. мы доказали сходимость процесса при т —>■ 0, но не при конечном г.
Хотя решение уравнения (12.5) теоретически достигается при t —► сю, но
убывание /(ж) экспоненциальное и /(ж) ~ 0 уже при конечном tj
Практически реализуется дискретный аналог уравнения (12.6)
-т71^1^) fix*1) = -т71^)

98
Гл.2 Методы решения систем алгебраических уравнений
хп+1 = хп- т^х*1) (12.8)1
Это нестационарный метод Ньютона. При тп = 1 получаем
классический стационарный метод Ньютона.
Доказано, что (12.8) сходится к (12.7) при тп —> 0 практически на
конечном интервале, если /# в окрестности \\х — х*\\ ^ \\х° — х*\\ имеет
обратную матрицу f^1: \f^l\ ^ 0.
Для практического вычисления необходимо (12.8) дополнить ал^
горитмом оптимального выбора тп на каждой итерации. Из принципа
уменьшения невязки в процессе итерации следует, что тп должно
уменьшаться пропорционально невязке и его следует брать удовлетворяющим
условию
тп = тп+1—— (12.9)
где 5п = 11/(^)11 = (/(ж71), J(я^))1/2 — норма невязки.
Из (12.9) следует, что шаг тп убывает с уменьшением невязки.
Для того, чтобы он не стал слишком малым, должно выполняться
условие, ограничивающее шаг снизу: 0 < G ^ тп ^ 1. Строгое
доказательство скорости сходимости алгоритма (12.8)-(12.9) отсутствует, но
численные расчеты показывают, что она быстрее, чем у стационарного метода
Ньютона. При уменьшении тп область сходимости расширяется, а
скорость падает.
Можно предложить приближенную формулу, уточняющую вы^
бор тп. Естественно считать оптимальным такой шаг тп, при котором
невязка вдоль выбранного направления минимальна. Направление
вектора приращения Ах™ на n-ом шаге определяется вектором а,
вычисляемого по формуле (12.8), квадрат невязки вдоль этого направления
равен
5n+V) = (/V+1)./V+1)) =
= (/(£" - тпа(х")), /(х" - тпа(хп))) (12.ю|
Строгое нахождение тп, минимизирующего 5п+11 сложная
задача. Но ее можно решить приближенно, заменяя действительную функ-

§12. Метод Ньютона и его модификации
99
цию 5п(т) параболой, разлагая в ряд в точке г = 0, полагая, что 5п(0)
отвечает невязке полученной на предыдущей итерации. Это следует из
формулы (12.10).
d2Sn
5п(т) = 5п(0)+т
d5n
+ т*
dr2
6n(0) = [f(^)}2 (12.11)
dSn
~d^
-2/(f - то") /*(£" - тсГ)сГ = -2/(f - тсГ) /(f)
~d^
(/(£"), /(f)) =-2(^(0) (12.12)1
Коэффициент при г2 определим по значению невязки при т = 1 по
формуле (12.11)
£п(1) = £п(0) - 2£п(0) + (<5°)"
Найденное значение (J0)" подставляем в (12.11)
6п{т) = 5п(0) - 25п(0)т + (5п(1) + 5п(0))т2
и находим г, при котором невязка 5п(т) минимальна
5fn{r) = -2^(0) + 2т{5п(1) + ^(0)) = 0
^(°) П2 13)|
при этом всегда выполняется 0 ^ тп ^ 1. Таким образом, на каждом шаге
необходимо по формуле (12.10) вычислять 5п для двух значений т = 0
и тп = 1. Эта формула обладает одним недостатком, можно построить
примеры, когда тп вдали от корня излишне мало.
Если тп мало, то скорость сходимости будет медленной вдали от
корня, поэтому его следует ограничивать некоторой величиной В
5п(0)
В,
6п(0) + 6п(1)\
Доказательство сходимости. В окрестности простого корня х = х* сходи^
мость при выборе шага тп+1 по формуле (12.13) так же, как в методе Ньютона
является квадратичной.

100
Гл.2 Методы решения систем алгебраических уравнений
Разлагая выражение (12.13) по степеням ёп(1)/6п(0), получим
6п(0) \5n(0)J
т.к. 5п(1) есть квадрат невязки на следующей итерации по отношении к 5п(0), в
стационарном методе Ньютона, для которого выше было доказано что 6п(1) <~ [£п(0)]2.
Тогда скорость сходимости неустановившегося процесса Ньютона (12.13) будет
квадратичной так же как и в установившемся процессе.
■2-+1 = g°- ^1^-1(^)^(2») = 2j+i - о(зГ - х») /тЧЛ/гИ
где х1^1 — значение полученное в (п + 1)-й итерации классическим методом
Ньютона.
§13. Методы минимизации функций (методы спуска)
Метод покоординатного спуска (100). Метод наискорейшего спуска (102).
Метод сопряженных градиентов (103). Итерационный метод с использованием
спектрально-эквивалентных операторов или с переобуславливанием (105).|
Метод покоординатного спуска. Выше были рассмотрены
итерационные процессы решения нелинейных уравнений. Перейдем теперь к
рассмотрению итерационных методов, в которых приближения строятся
в значительной степени с некоторым произволом. Это семейство методов
минимизации функций вида F(x) = p{x) f(x)^ к которым приводятся
многие задачи механики минимизации функционалов. Эти методы
объединяются под названием методов спуска, если для каждой
последующей итерации выполняется условие
F (^+1) < F (х11)
Здесь F(x) — скалярная функция, зависящая от векторного аргумента.
Направление спуска задается некоторым вектором Dn, а величина
поправки Ах™ — скаляром \п так, что
xn+l = xn + \nDn (13.1)|

§13. Методы минимизации функций (методы спуска)
101
Наиболее простой способ минимизации, когда Dn задается произвольно,
например, по одной из координат Xj
/5» = е? = { 0...0 1 0...0} r J
1... j ...Г
a \n определяется из приближенного условия максимального спуска в
точке Xй. Этот метод называется методом покоординатного спуска.]
Функция F(x) заменяется квадратичной функцией в соответствии
с разложением по Хп
F(£" + А»£>») = F(x») + А" (Ц) D» + \ (А«)2 Щ Т {^щ} D"
3 (13.3)
Если принять условие, что F[x) в окрестности Xй выпуклая функция, т.е.
функция удовлетворяющая неравенству
F(xt) = F(aX + (1 - a)Y) < aF(X) + (1 - a)F(Y) 2 ^ ° V?
X ^ x* ^ Y
то максимальный шаг будет определяться из условия
Зп = 0 (13.4)
dF(x» + \»D») (^\пп + )п\бп\т OF
д\п \дхп J L J dxi&Xj\
Откуда, учитывая что Dn определяется формулой (13.2), получаем j
компоненту (gradF)j
dF
дхп± , 3
и диагональный член матрицы Гессе Fjl-
9]=l^Dn
' 33
т яр
Dni
dxidxj
Xn
dF
_ ^ __ 771П
dxjdxj n
33
Ясно, что произвол в направлении Dn можно использовать для
того чтобы увеличить скорость спуска Лп. Для этого нужно взять j таким,

102
Гл.2 Методы решения систем алгебраических уравнений
чтобы д™ было максимальным среди 1 ^ j ^ г, где г — размерность
вектора х или пространства Ег
\\dF/dxj\\ = max |aF/&^|
з
Это наиболее простое использование произвола в методе
покоординатного спуска или методе покоординатной релаксации. Возможны и другие
модификации этого метода.
Метод наискорейшего спуска. Ускорить процесс спуска можно,
если направление спуска брать не по координатным направлениям, а в
направлении градиента OF'/dx11\ полагая
бп=-Л=& Н
а величину Лп определить, как и выше, из условия (13.4). Тогда получаем
r^nr^n BF \
Хп = 7f где Н = — матрица Гессе (13.6)
\gA ~HGn oxi0Xj\gn
Если F(x) квадратичная функция вида
F(x) = A + 2bYx + fTHf
где b — вектор, а Н — симметричная положительно определенная
матрица Н = Нт > 0, то отыскание ее минимума отвечает решению линейной
системы уравнений с матрицей Н. Тогда можно строго доказать
существование и единственность решения, а следовательно и сходимость
итерационного процесса.
На рис. 13.1 изображены последовательные приближения методе
наискорейшего спуска и линии уровня функции F{x). Итерационный
процесс (13.1), (13.5), (13.6) называют методом наискорейшего спуска
решения системы линейных уравнений Их = 6, когда F(x)
квадратичная функция.
Если F(x) не является квадратичной функцией, то и тогда процесс
может сходится при условии, что имеется начальная точка Xq настолько

§13. Методы минимизации функций (методы спуска)
103
Рис. 13.1.
близкая к ж*, что квадратичные члены существенно преобладают над
членами более высокого порядка, а матрица Гессе в окрестности х*
положительно определена.
Метод сопряженных градиентов. Предназначен для решения
линейных алгебраических уравнений с симметричной положительно
определенной матрицей. Этот метод позволяет минимизировать
квадратичную функцию п переменных за г итераций, где г ^ п. Запишем
квадратичную функцию в виде произведения сомножителей (х — я?*), где ж* —
точка, в которой F(x) принимает минимальное значение
F{x) = F0 + \{х - f*)T ' Н • (х - х*) (13.7)|
Градиент F(x) будет
G(x) = VF = Н • (х - х*) (13.8)|
Будем считать, что векторы /5°, D1,... Dn в формуле (13.1) линейно
независимы и образуют базис. Обозначим ак = XkDk1 тогда к-ую итера-

104
Гл.2 Методы решения систем алгебраических уравнений
цию можно записать через предыдущую
хк = xk~l + ak~l = ...= xk~q + dk~q + ... dk~l
Для n-ой итерации получим
га-1
Xй = хк+1 + ^ & п^к + 2
q=k+l
Gn = U-(xn-x,) = Шк+1 + Н • I ^ аЛ - ШJ
\q=k+l J
т.к. Нж* = 0 и Hxk+l = Gk+1, то n-ая итерация Gn может быть записана
в виде разложения по векторам базиса Dn
п-1
Gn = Gk+1+ Y1 Л"н^ fc = -l,0,...n-2
пХМ (13-9)|
Gn = G° + Yl A<?HZ^
Свертывая (13.9) с Dn', получим
т т n_1 т
[Gnl Gn=[(5fc+1] (5fc+1+^A4^j H№ (i3.io)
q=k+l
т.к. Dn — линейно независимые вектора, их можно выбирать так, чтобы
они были бы Н-сопряженными или Н-ортогональными, тогда|
\3q] UDk = 0 q^k (13.11)
Тогда из (13.10) получим
Gn] Dk = 0 (13.12)
Т.к. Z5°,... Dn образуют базис, то (13.12) означает, что
Gn = 0

§13. Методы минимизации функций (методы спуска)
105
а в силу (13.8) и Н > 0 получаем, что х71 = ж*, т.е. процесс за п итераций
приводит к точке, в которой F{x) имеет минимум.
Н-сопряженный базис, удовлетворяющий условиям (13.11) легко
можно получить из первоначального базиса D0,... Dn следующим
способом.
бо = _& = _VF(fo) I
3k+1 = Gk + pk&t [ \
(Зк выбирается так, чтобы удовлетворить условию Н-сопряженности (13.11).
Используя формулу (13.13), находим (Зк
[&]т
d
UDk
Dk\ UDk
/Г
Для вычисления хк шаг вычисляется по формуле (13.6).
Если F(x) неквадратичная функция, но начальное приближение
достаточно хорошее, то и в этом случае сходимость достигается, но
возможно не за п а за большее число итераций.
Итерационный метод с использованием
спектрально-эквивалентных операторов или с переобуславливанием. Как уже
отмечалось в §10, часто в итерационный процесс вводят в качестве
дополнительного параметра матрицу В и выбирают ее так, чтобы ускорить
сходимость процесса. Покажем это на примере метода простой итерации.
B^+1 = Вхп-т (Ах" - Ь\ (13.14)1
Матрица В должна быть легко обратимой. Тогда, умножая (13.14)
на В-1, получим
хп+1 = хп- тВ1 (Ах" -Ь\ = (Е- тВ-гА) Xй - тВ~1Ъ (13.15)1
где Е единичная матрица. Итерационный процесс (13.14) есть метод
простой итерации с матрицей (Е — тВ_1А) вместо матрицы А.

106
Гл.2 Методы решения систем алгебраических уравнений
Пусть отношение границ спектра /j, ^ Л^ ^ М или число
обусловленности М//2 ^ 1, Л^ — собственные значения положительно
определенной матрицы А > 0. В этом случае, как известно, итерационный
метод (13.15) при В = Е сходится медленно. Какой должна быть
матрица В > 0, чтобы процесс (13.15) сходился быстрее?
Обозначим
(Ах, х) . _ (Ах, х)
Mi = sup ; J _!, /ii = inf - ;
(Bx, xj x (Bx, xj
при В = E, M\ = M и iii = fi. Можно показать [2], что если В подобрана
так, что
Мг М J
— < — 13.16
Ml M
то итерационный процесс с переобусловленной матрицей В сходится
существенно быстрее. Множитель перехода от итерации к итерации q^
выражается отношением
?2-M^-i + ft/Mi (13Л71
и в силу неравенства (13.16) (ft < 0.ъ гДе Qi множитель перехода в
итерационном процессе при В = Е. Следовательно введение
переобуславливающей матрицы В ускоряет метод простой итерации.
Переобуславливание можно вводить и в другие рассмотренные вы^
ше итерационные процессы аналогично тому как это сделано в
отношении метода простой итерации.
Отметим, что в некоторых случаях выбор матрицы В оказывается
достаточно простым. Например, можно в качестве такой матрицы
использовать матрицу, состоящих из диагональных элементов матрицы А.
Этого часто оказывается достаточным для улучшения сходимости
итерационного процесса.
§14. Упражнения
1. Построить примеры плохо и хорошо обусловленных систем
линейных уравнений с двумя и тремя неизвестными. Исследовать их обу-

§14. Упражнения
107
словленность.
2. Доказать, что условие clet А = е<1, недостаточно для плохой
обусловленности матрицы. И обратно, условие clet A = 0(1)
недостаточно для хорошей обусловленности. Привести контрпримеры.
3. Доказать, что для матрицы А = ВТВ > 0 число обусловленности
(7(A) = С (В)2. Использовать определение нормы матрицы А
||А||2 = ?1^Г
4. Доказать, что если матрица А ^ 0, то ВТАВ ^ О
5. Исследовать сходимость метода Зейделя, когда для матрицы А
имеется диагональное преобладание по столбцам.
п
<l\aij\ > ^ ач> г = 1...П, q< 1
6. Решить систему уравнений Ах = Ъ используя разложение Холецко-
го
1 0.5 0 \ /2^
0.5 1 0.4 , Ъ =
О 0.4 2 у
7. Стержень из нелинейно упругого материала с диаграммой
а = Е(е + £о) растягивается напряжением <т* = 10Esq. Можно
ли определить деформацию растяжения 1) методом простой
итерации; 2) методом Ньютона?
8. Определить методом секущих модулей деформацию при
растяжении стержня напряжением а* = |, если диаграмма материала в
безразмерных переменных имеет вид
2-1/s е>1
Начальное приближение s^ = 1. Вычислить 5 итераций.

108
Гл.2 Методы решения систем алгебраических уравнений
9. Методом простой итерации
ип+1 = ип + т (А?1 - /)
найти решение системы двух уравнений Ах = f
хг + х2 = 1
/-. \ -1 ' е < 1
х\ + (1 + е)х2 = 1
Показать, что система плохообусловленна и определить параметр т
при котором скорость сходимости итерационного процесса опти^
мальна. Проверить это численным расчетом для начального при^
ближения х\ = 0.8, х2 = 0 при е = 0.1. Сравнить с точным
решением.
10. Найти решение системы нелинейных уравнений
х\ + х2 = 5
Х1+х2 + Ъ,Щх2)2 = Ъ
методом простой итерации. Определить значение параметра т при
котором процесс будет сходящимся. Определить число итераций для
получения решения с точностью до 10~4 при начальном приближе^
нии 40) = 0.9 и 40) = 0.
11. Найти решение системы двух нелинейных уравнений методом
секущих и модифицированным методом Ньютона в полуплоскости
xi > 0.
Ы2 + Ы2 = 1
х\ + х2 = 0.1
Вычислить 5 последовательных приближений при начальном при^
ближении х\ = —х2 = у/2/2. Оценить точность полученного
решения. Сравнить решения полученные по двух методам с точным
решением задачи.

§14. Упражнения
109
12. Получить решение задачи о релаксации напряжения неявным
методом Эйлера. Рассмотреть нелинейный упруговязкий стержне,
определяющее уравнение которого в безразмерных переменных
да де лп
при начальном условии a\t=o = 1
Полученную систему алгебраических уравнений решить методом
Ньютона для к = 3, 7, 5, к число разбиений отрезка t= [0,3],
а = 10, п = 5, tmax = 3.

110 Гл.2 Методы решения систем алгебраических уравнений

Глава 3
Методы решения краевых задач для
систем ОДУ
§15. Численное решение двухточечных краевых
задач. Устойчивые и неустойчивые алгоритмы]
Жесткая двухточечная краевая задача (111). Метод начальных параметров (113).
В главе 1 рассматривалось решение задач с начальными условиями
(задач Коши) для обыкновенных дифференциальных уравнений и для
уравнений в частных производных. Однако большинство задач МСС —
краевые или начально-краевые задачи.
Перейдем к решению краевых или граничных задач. Начнем с наи^
более простой задачи для обыкновенного дифференциального уравнения
второго порядка.
Жесткая двухточечная краевая задача. Численное решение
двухточечной краевой задачи для обыкновенного дифференциального
уравнения второго порядка связано с определенными трудностями, хотя
аналитическое решение не вызывает затруднений. Хотя исходная краевая
задача хорошо обусловлена, далеко не всякий алгоритм решения
обладает свойством хорошей обусловленности и устойчивости по отношению
к малым возмущениям. Это приводит к быстрому накоплению ошибок
округления и неадекватному решению задачи.
111

112
Гл.З Методы решения краевых задач для систем ОДУ
Рассмотрим сначала решение краевой задачи для
дифференциального уравнения с постоянными коэффициентами при наличии в
уравнении большого параметра а2 ^> 1. Пусть необходимо решить краевую
задачу на интервале х Е [0,1].
y"-a2y = 0 y(0) = Yo y(l) = 11 (15.1JI
Ее общее решение имеет вид
у = С\ ch ах + Ci sh ax
у' = а(С\ sh ах + С<2 ch ax)
Удовлетворяя краевым условиям (15.1), получаем систему двух
уравнений для определения произвольных констант Ci, Сч откуда находим
решение краевой задачи (15.1) в виде
( \ v л, i Yl ~ Y° ch a ъ
у(х) = Yq ch ах -\ sh ax
sha
Представим это решение в виде двух слагаемых, содержащих
множителями Уо и У\
, ч sh(l — x)a shax .. .^г „/ Чтл . ^ \
у(х) = Уо^-г—- + Y1—— = A(x)Y0 + B(x)Y1 15.2
sh a sh a
Исследуем его поведение при больших значениях параметра а ^> 1.
Найдем значения А(х) и В(х) при а —► оо
А(х) = lim shfl(1~a;) = о х е (0,1], А(0) = ll
a^oo sh a
с In /77»
B(x) = lim -— = 0 xe[0,1), Б(1) =
a^oo sha
т.е. Л (ж) и 5 (ж) разрывные функции на концах интервала в точках х = 0
и ж = 1 при a —► оо. Решение краевой задачи (15.1) при а ^> 1 показано
на рис. 15.1а. Функции А{х) и Б(х) при Yq и ^i B (15.2) ограничены при
a —► оо для любых ж £ [0,1], и решение краевой задачи (15.1) устойчиво.
Малое отклонение в У0 или ^i приводит к малому отклонению
решения у(х), ошибка не растет, хотя решение в окрестности концов
изменяется очень быстро и образует "пограничный слой" или "краевой эффект"
(рис. 1.1). Этот эффект имеет место и для более сложных систем
уравнений при наличии большого параметра.

§15. Численное решение двухточечных краевых задач
113
Yo
1 х
(а)
(Ь)
Рис. 15.1.
Метод начальных параметров. Для численного решения краевой
задачи (15.1) можно воспользоваться методом начальных параметров.
Он состоит в отыскании 2-х линейно независимых частных решений,
каждое из которых удовлетворяет двум условиям на левом конце
1) yi(0) = l, y'1(0) = 0. У1:
2) у2(0)=0, Уг(0) = 1. у2
chax
1
shax
(15.3)
а
Затем с помощью линейной комбинации двух этих решений
удовлетворяют граничным условиям (15.1) на правом и левом концах
^2,Ciyh
y(0) = Ciyi(0)+C2y2{0) = Y0
y{l) = Clyl(l)+C2y2{l) = Yl
(15.4
Из полученной системы уравнений находят произвольные константы С\
^olfe(l) - ^Jfe(O)
С2
!/i(0)ife(l)-i/i(l)jfe(0)'
Yoyi(l)-Yiyi(0)
lim d
lim C2
a—>oo
OO
CO
OO
OO
J/i(0)3te(l)-3/i(l)№(0)'
Так как 2/i(1), 2/2(1) ~^ °° при a —► 00 (рис 15.16), то в численном
решении только при абсолютно точных вычислениях можно
удовлетворить граничным условиям при а —> оо: lim C\ = Yq и lim C2 = Y\. А при
приближенных вычислениях сколь угодно малые отклонения Jyi(O) или
Sy2(0) могут привести к сколь угодно большой ошибке, т.к. с помощью

114
Гл.З Методы решения краевых задач для систем ОДУ
линейных комбинаций бесконечно больших чисел необходимо получить
конечные величины С\ ~ 0(1), С^ ~ 0(1).
Решение системы (15.4) находится через разности решений у\ и у^
при х = 1. Но 2/i и ?/2 — большие числа при а ^> 1. Следовательно С\ и Сч
будут определяться их "хвостами", которые зависят от ошибок
округления или других малых возмущений. Главные же части больших чисел
при вычитании исчезают и решение будет зависеть от ошибок
округления, а не от реальных физических условий задачи.
Ясно, что такой алгоритм получения С\ и Сг, неудовлетворителен.
Алгоритм метода начальных параметров опирается на получение двух
линейно независимых решений задачи Коши, а это приводит к плохо
обусловленной задаче при а ^> 1.
Действительно матрица А системы (15.4) имеет вид
2/1 (0) jfe(O)
l/i(l) 2/2(1)
ее собственные числа равны
> 2/1 (0) + 2/2(0)^ |7yi(0) + 2/2(0)
А1,2 = 7> ±
N |1/2
1/1 (0)2/2(1)- 2/2(0)yi(l)N'
Учитывая, что из (15.3) следует, что
l/i(0) ~ jto(0) ~ 0(1), 2/i(l) ~ 2/o(l) ~ 0(ea)j
получим
Ai,2 = 0(еа) ± [0{е2а) + 0(еа)]1/2 = 0(еа) ± [1 + |0(е"а)]1/2
Откуда следует Ai = 0(еа), Л2 = 0(1) и число обусловленности
задачи (15.4) С = 0(еа).
Таким образом, поскольку решение задач Коши для уравнения
второго порядка при а ^> 1 дает быстрорастущее решение порядка 0(еа),
то оно не должно входить в численный алгоритм решения краевой
задачи, т.к. это приводит к неудовлетворительному численному решению

§16. Краевая задача для системы линейных уравнений
115
для больших значений а. Конечно, при а ~ 0(1) метод начальных
параметров позволяет получать вполне корректное решение краевой задачи,
т.к. в этом случае матрица системы (15.4) будет хорошо обусловленной
С = 0(1).
Итак, для дифференциальных уравнений второго порядка вида
(15.1) при больших параметрах а получить хорошее численное
решение невозможно ни методом начальных параметров, ни каким либо
другим методом. Например, методом "пристрелки", основанным на
численном решении задачи Коши, и варьировании второго условия при х = О
в (15.3), до тех пор, пока не выполнится граничное условие на правом
конце х = 1.
Краевая задача (15.1) при больших значениях а2 ^> 1 относится к
так называемым "жестким" задачам, численное решение которых
связано с определенными трудностями, преодоление которых требует
специальных методов решения.
§16. Краевая задача для системы линейных
уравнений
Такая же ситуация возникает и в общем случае при решении кра^
евой задачи для линейной системы неоднородных дифференциальных
уравнений с переменными коэффициентами
du
— = А(х)и + а(х) 0<ж<1 (16.1)
ах
где и(х) и а(х) — n-мерные вектор-функции, А(х) — матрица
размерности п х п.
Краевые условия в общем виде могут быть записаны следующим
образом
Bu(0) + CU(l) = f (16.2)|
В и С матрицы с постоянными элементами пхп, f — n-мерный вектор.

116
Гл.З Методы решения краевых задач для систем ОДУ
Общее решение (16.1), как следует из общей теории линейных диф^
ференциальных уравнений, имеет вид
п
u(t) = и°(х) + ^ СА(х) (1б.3)|
г=1
где vP(х) частное решение неоднородной системы, щ(х) — система п-
линейно независимых решений однородной системы (16.1) при а(х) = О
и при любых однородных линейно-независимых краевых условиях
вида (16.2) при /= 0.
Задача (16.1)-(16.2) невырожденная и имеет единственное решение]
если определитель системы, полученной из (16.2) после подстановки в нее
решения (16.3), относительно неизвестных d не равен нулю.
Решение этой общей задачи методом начальных параметров, также
как и для уравнения второго порядка, рассмотренного выше, может быть
получено численно сведением к решению (п + 1)-ой задачи Коши.
Сначала одним из численных методов решения задачи Коши, на^
пример методом Эйлера, находится частное решение vP(x),
удовлетворяющее нулевым начальным условиям й(0) = 0. Затем находятся п
частных решений задач Коши щ однородной системы уравнений (16.1) при
а{х) = 0, удовлетворяющих условиям
^г(О) = k г = 1.. .п
где li — n линейно независимых единичных векторов n-мерного базиса.
Подставляя полученные решения в (16.3) и затем в (16.2), решаем
систему (16.2) относительно постоянных С{. Число обусловленности
получаемой матрицы будет большим, а задача плохо обусловленной в том случае,
когда среди решений щ(х) имеются быстрорастущие решения, т.е.
когда характеристические числа матрицы А(х) таковы, что среди них них
найдутся такие, для которых будут выполняться условия Re\ ^ С, где
С ^> 1, аналогично тому, как это было при решении уравнения второго
порядка. Если таких А^ нет, то задача хорошо обусловлена и изложенный
способ вполне приемлемый.

§17. Решение краевых задач для нелинейных уравнений
117
§17. Решение краевых задач для нелинейных
уравнений
Метод пристрелки (117). Метод квазилинеаризации (118).
Краевая задача для системы нелинейных уравнений
формулируется также как и для линейных уравнений, с той разницей, что линейные
функции в уравнении (16.1) заменяются нелинейными.
Пусть требуется найти решение системы уравнений вида
— = f(u,x) O^z^l (17.1}
при краевых условиях общего нелинейного вида
Ф(и(0),м(1)) = б (17.2)|
и(х), /(х,г?), Ф вектор-функции размерности п.
Метод пристрелки. Метод решения заключается в сведении этой
задачи к многократному решению задачи Коши. Начальные данные Коши
щ(0) = oil (г = 1.. .р) принимаются за неизвестную переменную а,
которую нужно определить из краевых условий (17.2). Определяем
решения уравнения (17.1) для последовательности значений ар, для которых
Ф (а, гГ(1, а)) ^ 0 меняет знак в интервале изменения параметра о\
ap-i ^ а ^ ар (1Т.3)|
Задача сводится к решению нелинейной системы п уравнений
относительно OL{
F(a) = $(a,u{l,a)) = 0 (17.4)|
заданной численно на интервале (17.3). Эта система уравнений может
быть решена одним из итерационных методов, например методом
Ньютона. Вычисление матрицы F^{a) также необходимо будет осуществить
численно. При этом придется находить решение вспомогательных
задач Коши для соответствующих приближенных значений нулей
функции (17.4) а = с№\ Вычисление каждого следующего приближения по

118
Гл.З Методы решения краевых задач для систем ОДУ
методу Ньютона
,#+1) = ,*<*) _ psl (<#>) F (<#>) (17.б|
связано с решением задачи Коши при й(0) = о№\ вычислением матрицы
Fa (а^) и ее обращением.
Естественно, что предлагаемый способ решения обладает, по край^
ней мере, теми же недостатками, которые были отмечены выше при
решении линейных уравнений и связаны со сведением решения краевой
задачи к решению задачи Коши, не считая дополнительных трудностей,
связанных с нелинейностью задачи.
Метод квазилинеаризации. Другой метод решения краевой задачи
для системы нелинейных уравнений основан на другой идее сведения к
краевой же задаче для системы линейных уравнений. Пусть требуется
решить для системы уравнений
!=/(**) Н
краевую задачу на интервале 0 ^ х ^ L при нелинейном условии общего
вида
F(y(0),y(L)) = 0 (17.7|
В качестве нулевого приближения берем некоторую функцию уо(х)\
по возможности близкую к действительному решению задачи (17.6)-
(17.7).
Линеаризуем уравнение (17.6) вблизи уо(х)
Ш(х) = у0(х) + 6у(х) I
j?0(x) + 8$ = f(y0(x),x) + fg [у0, х] 5у
Получаем для определения 5у систему линейных дифференциальных
уравнений с правой частью, fg - якобиева матрица функции f{y)\
Краевые условия тоже линеаризируются
F (Й(0) + ЩО), УоЩ + ВД) = О
F(y0(0),y0(L)) +4(0) (Й(ОШЬ)) $7(0)+ (17.9)|
+4w(yo(0),yo(i))^(L) = 0|

§18. Прогонка для ОДУ
119
Получаем линейную неоднородную систему краевых уравнений
относительно Syo(0) и 6yo(L) аналогичную (16.2). Задача (17.8)-(17.9)
решается также как и задача (1б.1)-(16.2).
После определения 8у(х) и соответственно у\{х) процесс
повторяется: ищем
у2(х) = уг(х) + 8ух{х) и т.д.
Таким образом, решение сводится к решению последовательности
неоднородных линейных краевых задач (17.8)-(17.9). Этот метод называется в
литературе методом квазилинеаризации. Сходимость метода
существенно зависит от того, насколько хорошо выбирается начальное
приближение уо(х). Для решения этой проблемы используется метод
продолжения по параметру совместно с методом погружения [9]. Эти методы
позволяют получать достаточно хорошее начальное приближение. Они
будут обсуждаться подробно во второй части книги, где решаются
нелинейные задачи механики деформируемого твердого тела методом
конечного элемента.
§18. Прогонка для обыкновенных
дифференциальных уравнений
Дифференциальная прогонка (119). Решение разностных уравнений методом
прогонки (123).
Как следует из предыдущего исследования, для того чтобы пoлy^
чить численно корректное решение жесткой краевой задачи,
необходимо построить такой алгоритм, который бы не содержал интегрирования
задачи Коши в направлении быстрого роста функции. Этого можно
добиться методом прогонки.
Дифференциальная прогонка. Рассмотрим метод прогонки на
примере уравнения второго порядка общего вида
ey"=p(x)y + F(x) е<1, р(х) > 0, х G [0,1] (18.1)|

120
Гл.З Методы решения краевых задач для систем ОДУ
при краевых условиях также общего вида
?/ = aoJ/+AU (18J
у' = aNy + (3N\X=1
Левое условие, продолженное в каждую точку х отрезка [0,1], выделяет
семейство интегральных кривых уравнения (18.1). Это семейство
удовлетворяет линейному уравнению первого порядка, которое можно записать
так:
у\х) = а{х)у + р(х) (18.3)
где а(х) и (3(х) неизвестные функции, подлежащие определению из
уравнения (18.1).
Соотношение (18.3) можно трактовать как перенос левого гранич^
ного условия в любую точку отрезка 0 ^ х ^ 1, в том числе и в
точку х = 1. Подставляя (18.3) в (18.1), получаем уравнения для
определения а{х) и (3(х)
у" = а'у + у'а + 0 = а'у + а2у + (3а + /3' = ^у + F^
е е\
Приравнивая нулю коэффициенты при у и свободные члены два
дифференциальных уравнения для определения а(х) и (3(х)
еф' + ар) = F(x) /3(0) = Д, (18.5)
Начальные условия (18.4)-(18.5) следуют из условия (18.2) при х = 0.
Алгоритм решения задачи (18.1)-(18.2) состоит из следующих этапов.
1) Решая уравнения Риккати (18.4) находим функцию а(х)\
2) Затем из уравнения (18.5) находим функцию (3{х).
3) Зная а{х) и /3{х) из (18.3), получаем дополнительное граничное
условие на правом конце х = 1
l/(l) = a(l)j/(l)+/?(l) (18.6)
которое вместе со вторым условием (18.2) при х = 1 дает хорошо
обусловленную систему 2-х уравнений для определения у(1) и у'{1)\

§18. Прогонка для ОДУ
121
4) При найденном значении у(1) на правом конце интегрируем урав^
нение (18.3) справа налево. Это так называемая обратная прогонка\
Прямая прогонка состоит в интегрировании уравнений для а{х)\
и (3(х) слева направо.
Таким образом, краевая задача (18.1)-(18.2) сводится к решению
трех задач Коши для 3-х уравнений первого порядка (18.4), (18.5)
и (18.3).
Заметим, что некорректно решать задачу Коши для уравнения
2го порядка (18.1) с двумя найденными условиями на правом конце
при х = 1, в обратном направлении, т.к. такая задача содержит
быстрорастущее решение.
Дадим обоснование корректности решения для всех трех задач
Коши, идея которого состоит в том, что для каждой из них интегрирование
ведется в направлении убывания решения так, что ошибка может только
затухать при изменении ж, [5].
Доказательство. Рассмотрим уравнение (18.4)
а\х) + а2{х) = —, а(0) = a0l х2 = р(х) > 0 (18.7)
при постоянном р(х) = к2 уравнение (18.7) допускает точное решение!
i
0
— ао
^а i
/1хо
зее"1/2
thz
1
X
(b)
Рис. 18.1.
(а)

122
Гл.З Методы решения краевых задач для систем ОДУ
II К ( \ К л ( Н ^ \
Ы>-щ, a(x) = ^cth^x + C1)
Ы<^, a(x) = -^th(^x + C2) (18.8)|
п l^ к + aps1/2 1 х + аде1/2 J
Ci = - In ш > О, С2 = -- In — < О
2 -х + аде ' 2 ус — а^е1'2
I. а) Пусть olq ^ 0, тогда в зависимости от знака неравенства |а0| < x/z1^2-,
поведение решения а{х) меняется, в первом случае оно возрастает, а во втором убывает,
но всегда а{х) ^ 0 для всех х (рис. 18.1а). Докажем, что ошибка ба(х) затухает.
Действительно, 5а(х) — погрешность решения уравнения (18.7) удовлетворяет
уравнению в вариациях:
5а{х)' + 2а{х) 6а(х) = О
и при а(х) ^ О при возрастании х величина |£а(ж)| убывает.
б) Решение второй задачи (18.5) при а(х) ^ 0 затухает с ростом х, поэтому и|
ошибка также затухает.
в) Наконец, третье уравнение (18.3) относительно у интегрируется в направле-|
нии убывания х и при а(х) ^ 0 у(х) и 5у(х) также убывают.
П. Пусть теперь а0 ^ 0, тогда, поскольку, |а0| «С н/е1/2 справедливо второе
выражение (18.8), показанное на рисунке 18.15.
а(х) = -^ th (^ + С2) , где - -^ < а0 ^ 0
Точка пересечения х0 кривой а(х) с осью х будет больше нуля и, следовательно J
а(х) на отрезке 0 ^ х ^ х0 будет отрицательной, ошибка 5а(х) на этом интервале
будет возрастать. Оценим это возрастание
-2 f a(x)dx
5а'{х) + 2а(х) 5а{х) = 0 5а = 5(0)е °
т.к. а[х) монотонная функция на рассматриваемом отрезке, то max 5a будет при
х = Xq. Вычислим интеграл в показателе экспоненты
а(х) dx = J —= cth I —=x + С2
, х 1 к + а0ф
- m cth | —=х — - In ■=
^е 2 х - apyje
-In
V~*

§18. Прогонка для ОДУ
123
т.е. 8а(х) мало отличается от ёа(0). Только, если ^ ~ £~1//2, тогда возможны
неприятности, но это противоречит исходному предположению, что все параметры
задачи 0(1), кроме е. На остальной части интервала х0 ^ х ^ 1 а(х) > О и ошибка 5а(ж),
как было показано выше, затухает. Очевидно, что такая же оценка будет
справедлива и для ошибок 5f3(x) и 5у(х). Таким образом корректность метода прогонки для
решения двухточечной краевой задачи (18.1)-(18.2) доказана.
Решение разностных уравнений методом прогонки. Перейдем
теперь к рассмотрению разностной аппроксимации краевой задачи для
произвольного уравнения второго порядка и перенесем рассмотренный
метод на решение разностных уравнений.
Разностная аппроксимация любого дифференциального уравнения
2-го порядка с гладкими коэффициентами
А{х)и" + В(х)и' + С{х)и = F(x)
при краевых условиях: х = 0 и(0) = (/?, х = 1 и(1) = ф будет выглядеть
следующим образом
anun-i + Ьпип + спип+г = fn (18.9)|
где номер внутренних точек п = 1... (7V — 1) х Е [0,1], а на концах
отрезка заданы краевые условия
х = О, щ = (р, х = 1, и^ = Ф (18.10)|
Уравнения (18.9)-(18.10) представляет собой систему линейных
уравнений с трехдиагональной матрицей, решение которой можно получить
методом прогонки в виде рекуррентных соотношений.
Первый интеграл уравнения (18.9) (аналог уравнения (18.3))
запишется в виде
ип = Ln+i2un+l + Кпц (18.11)
где Ln+i и Kn+i — коэффициенты, которые определяются из
уравнения (18.9). Начальные значения Ьп+\ и Кп+\ находим из краевых
условий (18.10)
П = О Щ = ^1/2^1 + Kl/2

124
Гл.З Методы решения краевых задач для систем ОДУ
из (18.10) следует, что Lx/2 = 0 if 1/2 = ^
1 ai{Ll/2ui + К i/2) + Ьгщ + ciu2 = /i|
wi
Ci^2
/l - ^ 1/2^1
^3/2 :
Cl
«1^1/2 + h
aiL1/2 + 61 aiL1/2 + bi
Si - Кфсц
aiLi/2 + 61
if.
3/2 :
(18.12)|
L.
n+1/2 '■
К
n+1/2 :
fn — ifn-l/2an
anLn-\/2 + bn UnLn-i/2 + bn
Формулы (18.12) для определения коэффициентов являются
разностными аналогами уравнений (18.4) и (18.5). Прямая прогонка заключается
в определении коэффициентов Ln+i/2 и ifn+i/2 по рекуррентным
формулам (18.12). Обратная прогонка состоит в определении ип по
найденному un+i с помощью формулы (18.11). При п = N — 1 из (18.11) находим
Ujsf-i — LN_ij2 Un + if JV-1/2 = LN_ii<2 ф + ifjv_i/2|
и далее для п = (N — 2)... 1 получаем последовательно u^-2, • • • щ\
Нетрудно видеть, что более общий вид краевых условий (18.2) при^
водит к решению системы 2-х уравнений при х = 1 относительно
неизвестных un и ujsr-i- После чего следует обратная прогонка.
Для счетной устойчивости из принципа максимума (см. §19) в
частном случае одной переменной следует, что должны выполняться условия
ап > 0, Ьп < 0,
> 0, К + сп| ^ \Ьп\
(18.13)
В этом случае ошибка округления при реализации метода не будет
возрастать.
Метод прогонки так же, как и метод Гаусса, является методом
факторизации. По-существу он представляет своеобразную реализацию
метода Гаусса в виде явных рекуррентных формул, получающихся
благодаря простому трехдиагональному виду матрицы А системы разностных
уравнений (18.9).

§18. Прогонка для ОДУ
125
Может показаться, что метод прогонки пригоден только для
решения краевой задачи для обыкновенного дифференциального уравнения
2го порядка. Однако это не так, решение этой задачи является составным
элементом решения широкого круга задач уравнений в частных
производных. Покажем это на примерах.
Прогонка для уравнения теплопроводности. Рассмотрим начально-
краевую задачу для уравнения теплопроводности
Ж = 0: «(О, *) = #*)
~i = ^^2 х = 1: u[l,t) = m (18.14j
t = 0 : и(х, 0) = щ(х)
Аппроксимируя функцию u(x,t) дискретными значениями ип{х) по
переменной t и полагая || = -—W~u w? получим N обыкновенных
дифференциальных уравнений второго порядка по переменной х\
^т(^гТ -ип+1(х) = -ип(х) п = 1...м\
с краевыми условиями
х = 1 ип+1(1, tn+l) = ф{1п+1) ип+\0, tn+l) = фп+1]
которая решается методом дифференциальной прогонки.
Если же для уравнения (18.12) взять простейшую неявную схему]
то получим
„.п+1 п „.п+1 _ 07/™+1 i п+1
Uk ~Uk _ „2Uk+l ZUk +%-!
— УС
2 / 2 \ 2 (18.15)
КТ п+1 (л О \ п+1 КТ п+1 _ п
~f^Tuk+i - 1 х + 2~^рГ J uk + ~f^Tuk-i - ~ик
Уравнение (18.15) и граничные условия (18.14) имеют при
фиксированном п вид уравнений (18.9)-(18.10), и решаются последовательно для
каждого временного слоя п + 1 при известном решении для слоя п
(рис 18.2) методом разностной прогонки (18.11)-(18.12).
Очевидно, что условия счетной устойчивости (18.13) для
коэффициентов уравнения (18.15) выполняются.

126
Гл.З Методы решения краевых задач для систем ОДУ
н
1
1
1
f~~V~^
г я
У\
.
А
•
J .
WWW
к~\ к к+\
А
W
А
ь ~_i_i
W \w ч i 1
А Т
Н
V
•
Т h
п
—►
о
1 х
Рис. 18.2.
§19. Решение краевых задач для эллиптических
уравнений
Прогонка для уравнения теплопроводности (125). Уравнение Пуассона (126).
Принцип максимума для разностных уравнений второго порядка (129).
Устойчивость уравнения Пуассона (130). Диагональное преобладание (130). Решение
уравнения Пуассона методом матричной прогонки (132). Метод разделения
переменных Фурье (136).
Уравнение Пуассона. Перейдем к разностным уравнениям краевых
задач для уравнений эллиптического типа. Простейшим примером таких
уравнений является уравнение Пуассона к решению которого сводится
большое число задач математической физики
д2и д2и
дх2 ду2
f(x,y)
Краевые условия должны быть следующего вида (xi £ сЮ)
либо условия Дирихле
либо условия Неймана
либо смешанного типа
и = Lp{x,
ди
дп
(f(Xi)
ч ди
a(xi)— + b(xi)u = <p(xi)
(19.1}
(19.2)|
(19-3)
(19.4)1

§19. Решение краевых задач для эллиптических уравнений
127
где п — нормаль к контуру <9П области П, а(а^), b{xi)) ф{хъ) ~ функции,
заданные на д£1.
Наконец, могут быть заданы на разных частях поверхности dVt
разные условия одного из трех перечисленных типов.
Рассмотрим разностную аппроксимацию задачи (19.1)-(19.4).
Введем для простоты квадратную сетку с шагом /i, которой покроем
область Vt. Введем сеточную функцию u(xk,ym) = Щ,т, определенную в
узлах сетки к = 0,1...К,т = 0,1... М.
тх\
о — внутренние
узлы
■• • - граничные
узлы
(а)
(т-\,п)
• (аи,и+1)
(т,п)
(т+1,п)\
(т,п-\)
(б)
Рис. 19.1.
Аппроксимация уравнения (19.1) запишется на простейшем пяти^
точечном шаблоне заменой вторых производных центральными
разностями второго порядка (рис. 19.1)
л ик+1,т ^клп ~т~ Uk—im
h2
. ^к,т+1 ^^к.т т ^клп—\ £ /1Г. г\
+ тъ = fk,m (19.5)|
Уравнение (19.5) определено только для узлов шаблонов, полностью
находящихся внутри области £1.
Аппроксимацию краевых условий проведем для простоты в случае
прямоугольной области. Например, на границах x = 0hi = 1 краевое
условие (19.4) можно аппроксимировать левой односторонней разностью

128
Гл.З Методы решения краевых задач для систем ОДУ
У\
f-[
А
и-
1 _
U
..--<
т—с
м
> «-
L1 -1
\ 1
L
■1-i-J
с
1и+1 А
Р А У
1 ▼ Л
г т
^ 1
» ^-
2 .
"I
1 ' 1
Т* 1
1
>
>
>
>
>
1
V 1
1 1
__J X
►1
Рис. 19.2.
на границе х = 0 и правой на границе х = 1.
V>l,m — и0,п
«0,7
^К,г
^,7
h
- uK-i,
h
Tl , 7
0,1... M
(19.6
Аналогично аппроксимируются краевые условия при у = 0 и у = /.
Порядок аппроксимации дифференциального уравнения (19.1)
О (/г2), а краевых условий только 0(h). Чтобы получить второй
порядок и для краевых условий нужно сдвинуть сетку по х и по у на h/ 2,
ввести законтурные узлы вдоль контура <9П, как показано на рис. 19.2 и
ввести аппроксимацию, например по ж, в виде
Щ,т + Щ,т
«1/2,:
^1,т — Щ,г
+ 6
1/2,?
^1/2,7
(19.7
h ' ~1/z'"i 2
где величины с индексом k = | точно соответствуют значениям на
контуре х = 0. Предполагается, что решение и(х,у) имеет четыре
непрерывные производные, это требует двух непрерывных производных
функции f(x,y) в Q и непрерывности функции ц) на контуре dQ\

§19. Решение краевых задач для эллиптических уравнений
129
Принцип максимума для разностных уравнений второго
порядка. Устойчивость разностной схемы (19.5)-(19.7) может быть доказана, исходя из
принципа максимума, справедливого для дифференциального эллиптического
уравнения (19.1). Из принципа максимума следует, что значение гармонической
функции в некоторой точке (ж, у) равно ее среднему значению, взятому по окружности
с центром в этой же точке. Принцип максимума справедлив также для разностной
аппроксимации (19.5).
Пусть имеем дифференциальное уравнение второго порядка эллиптического
типа
-Е^-(а<а)(г)^)+«(г)« = /(г). ^€0 (19.8
Разностная аппроксимация (19.8) на прямоугольной сетке с шагом ha по
координате ха (а = 1, 2) имеет вид
гл f \ ( (1) uk+l,l~Uk,l (l)Uk,l — Uk-l,l .
Qh(u) = ~ [al+1,1 ^ " - аК1^—^ +
. (2) Щ,1+1 — Щ,1 (2)Uk,l - Uk,l-l\ f /inrJ
+ak,i+i h~2 ak,i /^ ) + Qk,iXk,iuktl = fk,i (19.9)
ak+1j = a(xk+1 - 0.5/ii, y) q(xk, yi) = qkyl
au+1 = a(xk, yl+1 - 0.bh2) f(xk, уг) = fkj
Перепишем разностное уравнение (19.9) в условном виде
Qh(uktl) = Aukii + ^ B{xkth 0^(0 (19-10)|
где £ — номера узлов шаблона Ni равного полному шаблону N без центральной
точки (fc, /).
А и В коэффициенты, удовлетворяющие условиям
:n"i
A(x)>0, B{x)>0, C(x) = A-^2B(x,£)>0, хеП (19.11)
Для решения разностного уравнения (19.10) справедлив следующий принцип
максимума: если сеточная функция и^ удовлетворяет однородным граничным условиям
Uh{xi) = 0 Xi Е <9Г2, а правая часть f(x*i) ^ 0 (либо f{x*i) ^ 0), то u{xj) ^ 0 (u(x*i) ^ 0)
и соответственно max \uh(xi)\ или min \uh(xi)\ достигается на границе (xi G дО)\
С помощью принципа максимума доказывается устойчивость и сходимость]
разностных схем вида (19.9)-(19.10).

130
Гл.З Методы решения краевых задач для систем ОДУ
Устойчивость уравнения Пуассона. Докажем, что если (Аи)к^т > 0 во
всех внутренних узлах, то max 1г&,т достигается в граничном узле.
к,т
Доказательство. Из условия (Au)k,m > 0, находим
Uk,m < Т (Uk-l,m + Uk,m-1 + Щ-1,т+1 + Щт+х) (19.12)
Откуда видно, что значение во внутренней точке действительно меньше, чем среднее
значение в узлах на описанной окружности.
Таким же образом получаем из условия (Аи)к^т < 0, что тпшик^т находится)
на границе.
В теории уравнений в частных производных на основании принципа
максимума доказано, что для решения задачи Дирихле для уравнения Пуассона имеет место
следующая оценка
NKIMI + WII (19-13]
где ||и|| = max |и|, R — радиус круга, охватывающего область Г2, функции / и ip в
правых частях уравнения (19.1) и граничных условиях (19.2)-(19.4).
Аналогичную оценку можно получить для разностного уравнения (19.5), рав-|
номерную по шагу сетки h
№<imi + wii (i9.i4)i
где \\u\\h = max \uk,m\.
Оценка (19.14) означает непрерывную зависимость решения разностного
уравнения (19.5) при условии (19.2) от правой части задачи. Из нее следует устойчивость
и сходимость с порядком равным порядку аппроксимации
Km - Ukim\\h < Crf + C2R2h2 (19.15)1
где Uk^m — решение разностного уравнения, Uk,m ~ проекция точного решения на
сетку /с, га, С\ — оценка нормальной производной U(x,y) на границе, С2 — оценка
четвертых производных U(x,y).
Диагональное преобладание. Из разностного уравнения (19.5)
следует, что коэффициент, стоящий на диагонали матрицы системы (19.5)
по модулю не меньше суммы модулей всех остальных коэффициентов
недиагональных членов. Это свойство есть следствие принципа
максимума и называется диагональным преобладанием. Оно очень важно и
его стараются обеспечить при построении разностный схем и для более

§19. Решение краевых задач для эллиптических уравнений
131
сложных уравнений, например для уравнения со смешанной производной
ихх + иху + иуу = f(x, у) (19.16)1
Аппроксимация смешанной производной со вторым порядком точности
по четырем узлам на девятиточечном шаблоне
1
(и
ху)к,т '
Ah2
,га+1 — ик+1,т-1 ~ ик-1,т+1 + uk-l,m-lj
дает уравнение с коэффициентами, приведенными на рис. 19.3а.
Видно, что при такой аппроксимации диагональное преобладание
нарушается. На рис. 19.36 производная построена с первым порядком точности
74
Г
1 i
> •
1 1
-4 |
> 1
V4 l
(а)
*-[ 7
'l |-
+
1 <
+1
1 <
1—1
'-4 |
> 1
-v4
+ = +
1 •
t1 )~3 1
' • '
(б)
Рис. 19.3.
по точкам, соединенным стрелками. Диагональное преобладание имеет
место и шаблон состоит всего из 4х точек. Коэффициенты
вычисляются как сумма коэффициентов для уравнения Пуассона на пятиточечном
шаблоне, а для смешанной производной на шаблоне, показанном
пунктирными стрелками.

132
Гл.З Методы решения краевых задач для систем ОДУ
Следует предупредить, что отсутствие диагонального преобладав
ния еще не является достаточным условием для неустойчивости схемы.
Решение уравнения Пуассона методом матричной прогонки.
Метод прогонки можно обобщить и применить для решения
эллиптических задач. Покажем это на примере уравнения Пуассона.
Рассмотрим обобщенную задачу Неймана для прямоугольника, показанного на
рис. 19.2.
uxx + Uyy = f(x,y) х,у е П
ди or ч , ^ .о (19Л?)|
граничные условия —— = Ь{х,у) и + q{x,у) х^у Е oil
Введем прямоугольную сетку с координатами узлов (т, п) и примем
следующие обозначения:
т = -1,0,1,... М п = -1,0,1,... iV
Система координат (х, у) показана на рис. 19.2. Внутренние точки с
координатами хт = (т + |) /г, уп = (п + |) к отмечены светлыми кружками.
Узловые значения и(хт, уп) = итп. Законтурные точки расположены на
расстоянии /г/2 и к/2 от границ прямоугольника и отмечены темными
кружочками. Они отвечают номерам (—1, М) и (—1, TV), соответственно,
в направлениях х и у.
Разностные уравнения для задачи (19.17) во внутренних точках
будут
^т+1,п ^^т,п ~т~ Urn—1,п ^т,п+1 ^^т,п ~т~ Um,n—l г
К* + jfe2 " /m'n (19.18)|
m = 0,l...M-l, n = 0,l...iV-l
записываются только для внутренних узлов, их число равно M7V, в
крайние уравнения входят значения в фиктивных законтурных точках
т=— 1,М, п = — 1,7V. Записываем граничные условия (19.17) через
центральные разности в точках контура dtt. Производная берется в
направлении внешней нормали к контуру dft.

§19. Решение краевых задач для эллиптических уравнений
133
Краевые условия на линиях х = 0 и х = а будут
—'— - = £(0, уп)—^-^ + <?(0, уп)
= £(а, Уп) + q(a, yn)
Краевые условия на линиях у = 0 и у = Ъ такие
—^— - = S(xmi 0)—^ — + д(жт, 0)
Um,N — UmN-i UmN + UmN-l , / ,ч
T = S(xm, b) + q(xm, b)
n = 0,...N-l
(19.19)
m = 0....M-1
(19.20)
Число уравнений (19.19) равно 27V, а уравнений (19.20) равно 2М,
т.к. для угловых точек уравнения вообще не требуются. Число
неизвестных в задаче (5 = (М + 2)(7V + 2) — 4. Общее число уравнений
а = М • N + 2(М + AT") равно числу неизвестных а = (3.
Введем в рассмотрение TV-мерные векторы
um{umfi • • • um,N-i) га = -1, 0,... М|
Qo(qo--qN-i) где qn = q(0, yn)
Qm{Qo-"Qn-i) где Qn = q(cL,yn)
компонентами которых являются значения функций только во
внутренних узлах. Значения функций в узлах, отмеченных темными
кружочками в направлении у в них не входят. Тогда условия (19.19) на границах
х = 0 и х = а можно записать в векторном виде
h В 2 + Qm
где А и В — диагональные матрицы.
Неизвестные в законтурных точках в направлении у образуют
векторы u_i и Им- Таким образом число векторов равно М + 2. В
уравнение (19.18) входят три компонента вектора ит и по одной
компоненте йт+\ и йт-\. Представим уравнение (19.18) в матричном виде, для

134
Гл.З Методы решения краевых задач для систем ОДУ
этого величины ит^_\ и ит^, которые не входят в вектор ит,
должны быть исключены из уравнений (19.18) с помощью уравнений (19.20)
при у = 0 и у = Ь.
^m,-l — ит,0
S(xm, 0)—-— - + q(xm, 0)
= S(xmj b) h q{xm, b
Тогда
um-i
= l + |S(sTO,0)
1-|5(жго,0)
1 + |5(xm, b)
um$ + ■
<k
2o(:rm, (J)
nNk
\-\S{xm,b) \-\S[xm,b)
Эти выражения подставляем в уравнения (19.18), содержащие узлы в
направлении у, и получаем
V"m+l,n ^^т,п ~г Urn—l,n
Ап цш,п+1 ^У1 Рп ;^гтг,гг "Г ^п %,п-1 _ /. _ т^ fl Q 9~П
~т 7 2 Jm,n "ran ^J-У-^-И
где п = 0 ... 7V — 1
A„ -
1 О^п^ 7V-2
0 n = N-l
0 п = 0
1 п >0
Мп
1 + 1^
2(1-|^)
0 F
л \kaN
-l -г 2^т
2 (1 - |^)
« (1 2^т)
к (1 2^mj
П = 0
n = iV-l
Нерегулярный вид имеют только уравнения при n = 0nn = N— 1.
Запишем уравнение (19.21) в матричном виде. В (19.21) входят
векторы um-i йт, йт+1, по одной компоненте первого и третьего вектора

§19. Решение краевых задач для эллиптических уравнений
135
и три компоненты ит^п_\^ umj i£m,n+i второго вектора ит. В матричном
виде (19.21) запишется так
1
/i2
(йщ+1 - 2ВШ йт + um-i) + dn
О
где матрица 2ВШ трехдиагональная следующего вида:
|i + (i-^)f -й о о ... о
2B„
'к2
1 + (1-мГ)р
"А;2
О ••• О
о
о
о о
к2
l + (l-^-i)F
Вектор dm = -h2(fm- Fm).
Таким образом, имеем следующую систему уравнений:
u-i - щ
u-i + щ
А h <
/о
/г '" 2
ит+1 ~ 2l3mum + ^m-1 + ^m = 0
йм-йм-i ^йм + йм-i , л
: -Ь г (
т = 0...М-1
(19.22)|
^о
/i 2
Система (19.22) (М + 2)-х матричных уравнений решается мето-\
дом прогонки. Действуя по схеме, которая применялась в случае
скалярного уравнения, перегоняем систему граничных условий при х = 0 на
правый конец х = а с помощью промежуточного интеграла^
u-i = K_1/2U0 + Z-i/2
Из первого уравнения (19.22) находим
'Е_ АЛ
/i ~ 2
'Е А\-1/Е А\ -* /Е Ал
w-i
Е А\ - ^
К
1/2 :
Е А
7Г + ¥
/_
1/2
wm-l — ^т-1/2ит + 'т-1/2
f 1 - 2Втит + К + 1
т— 1/2 +
^ш = Кш+1/2^т+1 + ^т+1/2
К
■т+1/2
= (2ВП
К
■т-1/2)
-1
т+1/2 :
К
-1
,(*■
,+1/2 ^т-1/2
+ fiL

136 Гл.З Методы решения краевых задач для систем ОДУ
при т = М получаем систему 2N уравнений относительно компонент
векторов им и йм+i
йм-1 = ^М-1/2% + УМ
Доказательство корректности матричной прогонки следует из
условия, что ||ХШ|| < 1, [6].
Решение задачи (19.17) с использованием матричной прогонки
требует М + 1 обращения матриц размерами М х N вместо решения систе-
мы из а = М х N + 2(М + N) уравнений.
Оценим выигрыш, который получается благодаря прогонке. При
М = N = 50 число уравнений а = 2500 + 200 = 2700, решение которых
методом Гаусса пропорционально а3 операций: а3 = 54 • 503. При
прогонке требуется М + 1 = 51 обращений матрицы 50 х 50, что
составляет ~ 503 операций: А = 51 х 503, т.е. выигрыш более чем в 542 раз:
а3 ^ 542Ai. Ясно, что направление прогонки должно быть выбрано так,
чтобы было N < М, тогда приходится обращать матрицу меньшей
размерности. Приведенная оценка конечно грубая, но дает качественную
оценку эффективности метода.
Метод разделения переменных Фурье. Перейдем теперь к
методам эффективного решения полученной системы сеточных уравнений.
Системы линейных алгебраических уравнений, которые получают^
ся из разностной аппроксимации дифференциальных уравнений
эллиптического типа (19.9) обладают определенной спецификой. Число
уравнений этой системы как правило очень велико и равно числу
используемых узлов сетки. Её матрица — слабо заполненная с большим
числом нулей и имеет ленточную структуру, состоящую из диагональных и
нескольких над и поддиагональных элементов.
Для решения таких систем целесообразно использовать
специальные, а не стандартные методы линейной алгебры, учитывающие специ-

§19. Решение краевых задач для эллиптических уравнений
137
фику системы уравнений. В частности для решения простых уравнений
математической физики второго порядка вида (19.8) широко
используется весьма эффективный метод разделения переменных Фурье.
Аналогичный метод, применим и к разностной аппроксимации этих
уравнений (19.9). Рассмотрим его на примере решения задачи Дирихле для
разностного уравнения Пуасона (19.5) на прямоугольнике при нулевых
краевых условиях
и = 0 при х, у £ dft
О ^ х ^ /i, 0 ^ у ^ 1% при х, у £ Q
После разделения переменных решение сводится к задаче на
собственные значения для разностного уравнения второго порядка по
переменной х
г4(ж) - Xkv\x) = О, г;0 = 0, vN = 0 (19.23)|
где обозначено
v(x + h, у) — v(x, у) v(x, у) — v(x — h, y)\
Vx = '■, 5 ^х = -
h h
_ Vx ~Vx
^XX 7
h
соответственно vx — правая, vx — левая и, vxx — вторая разностная
производная по х.
Решение уравнения (19.23) ищется в виде v = Vo ехр(шкх).
Подставляя, получаем собственные значения и собственные функции
4 . 9 fiikhi\ /дл. ч [2 . (irkx\ ,«„„*]
Xk = rin{-w)' v{)^=4hsin{~v) (19-241
k = 1,2...N-1
где hi — шаг сетки по ж, a hi — шаг сетки по у.
Решение щ(х, у) уравнения (19.5) будем искать в виде разложения
по собственным функциям у^к\х)
N-1
uh(x,y) = ^k\x)w^(y) (19.25)1
к=1

138
Гл.З Методы решения краевых задач для систем ОДУ
где w^k\y) — решение неоднородного уравнения по переменной у при
нулевых краевых условиях, для определения которого получаем систему
разностных уравнений
1 («,«1>в - 2<)„ + w<£>_ 1>n) - А« <>„ = /<*> (ут,п) (19.26j
при условиях WqI = 0 и wN' п = О, где f^k\y) — коэффициенты Фурье
правой части уравнения Пуассона
N-1
f(kKy) = YlHy,x)vik\x)hi (1Э.27]|
к=1
Итак, метод Фурье заключается в нахождении собственных
значений и собственных функций разностной задачи по переменной х (по
ней коэффициенты постоянны), вычислению коэффициентов Фурье по
формуле (19.27), нахождению коэффициентов w^ в разложении (19.25)
с помощью решения краевой задачи (19.26) методом скалярной
прогонки. Вычисление коэффициентов Фурье (19.27) осуществляется
эффективным численным методом быстрого преобразования Фурье при
затратах пропорциональных числу операций q = (^(A^A^logiVi). См. [13, 20].
Для задач с постоянными коэффициентами можно использовать
преобразование Фурье по обеим переменным х и у, раскладывая по
собственным функциям двумерного разностного оператора Lh(x,y)\
Численное решение эллиптических уравнений представляет собой
хорошо разработанный раздел вычислительной математики, который
будет рассмотрен во второй части книги при использовании метода
конечных элементов (КЭ). Методу КЭ в последнее время отдается
предпочтение по сравнению с методом конечных разностей, особенно при решении
задач для трехмерных областей сложной геометрии.

§20. Жесткие краевые задачи
139
§20. Жесткие краевые задачи
Жесткие системы дифференциальных уравнений (139). Обобщенный метод
начальных параметров (141). Ортогональная прогонка (143).
Жесткие системы дифференциальных уравнений. В §18
рассмотрено решение жесткой краевой задачи для линейного
дифференциального уравнения второго порядка. Проанализировано в чем состоят
трудности такой задачи и изложена идея метода их преодоления.
Перейдем к обобщению этого метода решения на общий случай жестких
краевых задач для систем линейных дифференциаьных уравнений.
Пусть дана система линейных уравнений общего вида
^-А(х)и = а(х) (20.1
ах
Краевые условия заданы в следующем виде
х = 0 [1г,й(0))=щ г = 1,2... к к<п
x = b lli,u(b)\=ai г = к-\-1...п
где г — номер краевого условия; круглыми скобками обозначено
скалярное произведение векторов. Заданы к условий на левом и (п — к)
условий на правом конце интервала [0,6]. Векторы й, U и а
принадлежат n-мерному векторному пространству Rn. Матрицу А для простоты
будем считать постоянной, хотя все последующие построения остаются
справедливыми и для А = А (ж), если элементы матрицы А(х) медленно
меняющиеся функции.
Определение Систему (20.1) будем называть жесткой, если спектр
матрицы А можно разделить на три характерные части, показанные
на рис. 20.1, удовлетворяющие следующим условиям
ReX^ ^ —L г = 1, 2 ... к~ отрицательная жесткая\
Re\~l ^ L г = 1, 2 ... к+ положительная жесткая (20.3)
|Aj| < Z г = 1, 2 ... га и мягкая часть спектра\

140
Гл.З Методы решения краевых задач для систем ОДУ
х; -/ / а;
//////////} ( Г ) (чччччччччччН
-L L
Рис. 20.1.
где к+ + к~ + т = п.
Деление это достаточно условное, т.к. четких границ между ча^
стями спектра может не быть. Важно, что существуют такие области
спектра положительных и отрицательных значений, для которых будет
выполняться условие ЪЬ ^> 1, а Ы ~ 0(1). Здесь Ъ — отрезок
интегрирования. Соответственно общее решение уравнения (20.1) тоже можно
представить состоящим из трех частей
к~ к+ т
и(х) = £ СгСЗГ ех~х + £ сД+ ех-х + £ ctCJt ел'ж (20.4^
г=1 г=1 г=1
где (Лг+, (Лг~, UJi собственные векторы матрицы А, отвечающие
соответствующим собственным значениям А+, Лг~ и Л^ в трех разных частях
спектра.
Будем рассматривать такие краевые задачи (20.1)-(20.2), для
которых решение ограничено
\\й(х)\\<С(\\а(х)\\ + \\а\\) (20.5)
где справа стоят нормы правых частей уравнения (20.1) и граничных
условий (20.2).
Подчеркнем, что это вполне определенный класс краевых задач,
особенность которого состоит в том, что несмотря на то, что в общем
решении (20.4) присутствуют быстро растущие компоненты, для
которых exp(L6) ^> 1, тем не менее решение краевой задачи имеет порядок
и(х) ~ 0(1). Решение не всех краевые задачи будут удовлетворять
условию ограниченности (20.5), если спектр матрицы А отвечает
условиям (20.3), и соответственно не все задачи (20.1)-(20.2) могут быть
корректно решены вычислительными методами.

§20. Жесткие краевые задачи
141
Будет ли задача корректно разрешимой численно или нет для
заданной системы уравнений зависит от числа краевых условий (20.2) на
правом и левом концах интервала и от соотношения этого числа с числом
точек жесткого спектра принадлежащих его правой и левой частям.
Необходимыми условиями корректной разрешимости задачи (20.1)-
(20.2) являются следующие
1)к^к~, 2)п-к^к+ (20.6)|
Число краевых условий к на левом конце должно быть не
меньше, чем число быстроубывающих решений при возрастании ж, и число
краевых условий на правом конце п — к тоже должно быть не меньше,
чем число убывающих решений при уменьшении х. Действительно,
если это условие будет нарушено, то тогда найдется ненулевое решение,
удовлетворяющее к однородным условиям на левом конце и убывающее
с ростом х. Ему будет отвечать краевое условие на правом конце и
тогда его интегрирование будет происходить справа налево в направлении
быстрого роста, а это приведет к быстрому росту малого возмущения.
Напомним, что подобная ситуация возникала при интегрировании
уравнения 2-го порядка методом начальных параметров в §18 и приводила к
плохообусловленной системе алгебраических уравнений на другом конце
интервала, при решении которой происходила потеря точности из-за
сокращения значащих знаков при вычитании двух больших близких по
величине чисел. Хотя сама исходная краевая задача (15.1) или (18.1)-(18.2)
были корректно разрешимыми, удовлетворяющими условию (20.6).
Обобщенный метод начальных параметров. Обозначим
семейство решений (20.1), удовлетворяющих только левым к краевым
условиям (20.2) через R~{x). Это многообразие решений при фиксированном х
является линейным (п — А:)-мерным подпространством n-мерного
пространства Rn(x) всех решений (20.1). Аналогично семейство решений,
удовлетворяющее правым (п — к) условиям, можно считать
многообразием Я+(ж), которое будет подпространством размерности к\

142
Гл.З Методы решения краевых задач для систем ОДУ
Явное выражение, например, для R~(x) можно получить следую^
щим образом.
1) Находим частное решение неоднородного уравнения (20.1),
удовлетворяющее неоднородным левым условиям (20.2). Для этого находим
какое-либо решение v°(0) Е Rni удовлетворяющее к левым линейным
уравнениям (20.2). При таких начальных значениях ti°(0) решаем
задачу Коши для уравнений (20.1) и находим vP(x).
2) Находим (п — к) линейно-независимых решений гГ(х), г =\
(1,2.. .п — к) однородного уравнения (20.1), удовлетворяющих
однородным левым краевым условиях (20.2). Решаем (п — к) задач Коши с
начальными условиями на левом конце, удовлетворяющих системе к
уравнений.
(£,й*(0)) =0, i = l...к
и дополненных для га-ой функции 1^(0) так, что г^+1(0) ... и^+п(0)
компоненты равны нулю, кроме г4+ш(0) равной единице.
( -^ W
гГ(0)= и1(0),гг2(0),...^(0),0...0, 1,0, ■■■01
\ п-к /|
т = 1,... п — к
Явное представление векторного многообразия R~(x) будет иметь вид
суммы
п—к
В.-{х) = гРл{х) + ^Ст^{х) (20.7)
т—1
В такой же форме представим и R+(x)
п
R+(x) = v°u(x) + J2 С^) (20-81
г=п—к+1
где и^(х) удовлетворяет условиям на левом, a v?n(x) — на правом конце
интервала, d произвольные константы.
Пересечение этих многообразий дает решение исходной задачи
(20.1)-(20.2)
и(х) = R-(x)C\R+(x)

§20. Жесткие краевые задачи
143
Таким способом можно решить задачу (20.1)-(20.2), если спектр
матрицы А не жесткий. Это обобщенный способ начальных параметров,
который рассматривался выше в §15 при решении уравнения второго
порядка.
Очевидно, что такой алгоритм можно осуществить численно.
Приравнивая полученные многообразия R~(x*) и R+(x*) в некоторой точке
х = ж*, получаем систему п уравнений для определения п произвольных
постоянных d (г = 1... п). Однако, непосредственно применять его для
решения жестких задач нельзя, т.к. не устранены причины, которые
приводят к плохой обусловленности системы уравнений интегрирование
задач Коши для быстрорастущих функций в направлении их роста.
Рассмотрим способ, который позволяет обойти эту трудность.
Ортогональная прогонка. Многообразия вида (20.7)-(20.8) состоят
из отдельных решений задач Коши, которые включают в себя сильно
растущие как вправо, так и влево решения вида eLx и e~Lx. Как
отмечалось выше они приводят к плохой обусловленности алгоритма начальных
параметров. Между тем, сами многообразия R±(x) при выполнении
условий корректной разрешимости (20.6) являются устойчивыми. Например,
R~{x) при малой вариации в правой части щ + Soli уравнений (20.2)
изменится на R~(x) + 5R~(x), где SR~(x) — малая величина.
\\SR-\\ <С\\5аг\\ С = 0(1)
хотя для составляющих его элементов константа С может быть большой
С = exp(6L) > 1.
Поэтому корректный метод надо строить так, чтобы оперировать
не с решениями отдельных задач Коши, а с устойчивыми
многообразиями R~(x) и R+(x). Как это можно сделать практически при численном
решении задачи?
Чтобы построить многообразие R~(x), которое представляет собой
гиперплоскость размерности (п — к) = г в пространстве Rn(%) нужно
задать точку на плоскости И° и базис из г линейно-независимых векто-

144
Гл.З Методы решения краевых задач для систем ОДУ
ров ei, бГ2, • • • еГ1 лежащих в этой плоскости. Тогда R (х) будет
множеством точек и Е Rn вида
г
и(х) = И°{х) + ^ С* ei(x) (20.9)|
г=1
где Ci — произвольные константы. Способ вычисления гл° и ei(x) был
указан выше. Но все дело в том, что базис ei(x), который в начальной
точке х = 0 выбирается ортогональным и хорошо определяет
многообразие R~(0), при интегрировании по х портится (сплющивается), т.е. углы
между этими векторами быстро уменьшаются и базис вырождается, а
многообразие (20.9) определяется неустойчивым образом и малые
возмущения в начальных условиях быстро растут и приводят к большим
ошибкам 5и. С вырождением базиса можно бороться, проводя
периодически (через небольшие промежутки) ортогонализацию базиса.)
Как мы видели выше, появление неустойчивых элементов в
многообразии R~(x) связано не просто с большими значениями А^, а с
произведениями \Ъ (Ь — отрезок интегрирования). Базис начинает портится,
только если достаточно велики eXiAx, где Ах отрезок интегрирования
отсчитываемый от точки, в которой проведена преддущая ортогонали-
зация базиса. Поэтому, беря отрезок интегрирования Ах достаточно
малым, можно сохранить базис близким к ортогональному. Проводя после
каждого такого промежутка интегрирования ортогонализацию базиса,
можно на всем интервале [0, Ь] иметь дело с вполне хорошим базисом.
В этом и состоит идея ортогональной прогонки по С. К. Годунову\
Подчеркнем, что при таком подходе мы уже имеем дело не с отдельными
несвязанными задачами Коши на большом отрезке интегрирования, как
это имело место в методе начальных параметров, а одновременно с целым
многообразием R~(x), его базисом ё^(х) и решением vP(x).
Приведем весь алгоритм ортогональной прогонки
последовательно. Перепишем краевые уравнения (20.2) в новых обозначениях. Введем
полную систему ортонормированных векторов /г~(0), г = 1, 2 ... га,
дополняя ортогональную систему ^-мерных векторов Z^(0), входящих в левые

§20. Жесткие краевые задачи
145
условия (20.2) до полной ортогонализированной системы, как
указывалось выше. Записываем левые к условий (20.2) через /гг(0)
(/г(о), й°(0)) =аг г = 1...к (20.10)
Интегрируем (п — к + 1) задач Коши, решения которых обозначим
y°(x),u1(x),.. .iT"fe(x), где а°(х) — обозначено решение неоднородного
уравнения (20.1). Определим его начальными условиями
к
Й°(0) = J>ir(0)
г=1
тогда удовлетворяться условия (20.10). Остальные (п — к) векторов гр(х)
определяются данными Коши для однородной системы
«J'(0) = C+J-(0) j = l,2...n-k
Тогда (20.7) дает при произвольных d все решения (20.1),
удовлетворяющие левым краевым условиям и является представлением R~(x]
п—к
R-(1\x)=H°(x1) + Y,CJuf)(x^>
3=1
которое можно использовать на отрезке 0 ^ х ^ А, где А такое, что|
ЦАЦА = LA« 0.1-1
При х\ = А многообразие R~{x) представлено формулой ви^
да (20.7). Заменяем его новым представлением. Для этого проведем
ортнормирование системы векторов й1(х\), U2(xi)... un~k(xi) и получим
вектора j\ (xi), Й> (xi).. .rn_k(xi). Новое представление той же
гиперплоскости в той же точке х\, будет
п—к
R~(xl) = if(xl)+Y/CJif)(xl)

146
Гл.З Методы решения краевых задач для систем ОДУ
После этого можно вместо решения неоднородной задачи Щ{х\)\
которое является некоторой точкой на гиперплоскости R~{x), найти
другую точку u°(xi), которая была бы ближе расположена к искомому
решению, чем точка vP(xi), так чтобы это расстояние было бы 0(1).
п—к
и°(Ж1) = й°оп) - £ (A**) V) f (20Л1|
3
Так как щ{х\) успела за промежуток интегрирования Ai несколько
отклониться и выйти из действительного положения гиперплоскости
R~(x). Поэтому каждый шаг в процессе ортонормированной прогонки
должен начинаться с восстановления ортонормированного базиса и
смещения начальной точки в новое положение на гиперплоскости.
После этого цикл вычислений повторяется на новом интервале Д^
Одновременно может проводится точно такой же процесс
интегрирования справа налево. Это делается для того, чтобы уменьшить
интервал интегрирования как при интегрировании слева направо, так и
справа налево. В результате мы перегоняем R~{x) и R+(x) в некоторую
точку, например, среднюю точку интервала х* = 1/2, оба многообразия.
Приравнивая их находим единственную точку их пересечения. Это дает
нам п линейных уравнений для определения п произвольных
постоянных Cj
п—к п—к
и°_Ы + ^Cjv?_{x*) = u°.(x#) + ^C/<(s*) (20Л2)|
i=i 3=1
Заметим, что положение точки ж* можно изменять и, следовательно,
получать разные системы уравнений (20.12), решения которых при
корректном расчете должны мало отличаться друг от друга. Такое
сравнение может служить проверкой точности полученного решения.
После определения постоянных интегрирования по формулам
(20.7)-(20.11) вычисляется искомое решение на всем интервале
интегрирования. Таким образом краевая задача (20.1)-(20.2) полностью решена.

§21. Упражнения
147
Жесткие задачи возникают при решении задач механики, в
которых процесс зависит от нескольких масштабов, сильно отличающихся
друг от друга. Например, в задаче могут фигурировать макромасштаб L
и микромасштаб /, связанный со структурой материала, отношение
которых 1/L <С 1. Либо при рассмотрении тонкостенных конструкций, когда
толщина тела h существенно меньше двух других размеров длины L:
h/L < 1. В задачах, зависящих от времени, могут сильно различаться
характерные времена одновременно протекающих процессов, благодаря
чему также появляется в уравнениях малый параметр т/t <С 1. Класс
таких задач очень широк и в механике и в физике и их корректное решение
численными методами крайне важно.
§21. Упражнения
1. Построить схему четвертого порядка аппроксимации задачи
Дирихле для уравнения Пуассона на 9-титочечном шаблоне 3x3. Нужно
дать аппроксимацию не расширяя 9-титочечного шаблона]
Указание: четвертую производную по ж и у в остаточном члене
смешанной производной 4-ого порядка выразить, используя уравнение
Пуассона.
Для центральной разности имеем
ui+i — 2щ + щ-\ д2и h2 д4'
+ 1-^ + 0(К) (21.1
Ы 8x1 \2dxi
Используя исходные уравнения Аи = /, находим
д4и д4и d2f д4и д4и d2f\
дх4 дх2ду2 дх2' ду4 дх2ду2 ду2\

148
Гл.З Методы решения краевых задач для систем ОДУ
Подставляя в правую часть (21.2), находим
U 12 \дх*ду2) + 12 \дх2 + ду2
д (хх) _ 9 А (хх) I Д (жж)
*J /,2
где Л(жж) вторая центральная разность по переменной х\
Расписывая смешанную производную на шаблоне 3x3, получаем^
схему 4-огопорядка на 9-титочечном шаблоне
2. Построить разностную схему с диагональным преобладанием для
уравнения
^хх ~г tLXy -\- Uyy J yX, у)
и исследовать ее аппроксимацию (рис. 21.1).
ХУ = Q^2~ \ик+1,т+1 — ик,т+1 ~ ик+1,т + ик,т) =
2h2
(Щ,т — Щт-1 — Uk-l,m + Щ-1,т-Ц
• 1
tl
• 1
1 Т
-4 t
> 4
+
/2
>2
'2
-3
/2
'v2
Ч
Рис. 21.1. Сумма шаблонов основной схемы и второй смешанной производной
3. Построить разностную схему для бигармонического уравнения
ААи = / на 13-титочечном шаблоне, показанном на рис. 21.2\

§21. Упражнения
149
(т-2,п)
(т+2,п)
(т,п-2)
Рис. 21.2. 13-титочечный шаблон
4. Решить двухточечную краевую задачу
у" = 10у3 + Зу + х2 0 < х < 1
2/(0) = 0, у{1) = 1
используя метод пристрелки. Предложить разумный способ при^
ближенного определения начального значения у;(0). Составить про^
грамму решения задачи Коши. Найти у'(0) при котором условие на
правом конце будет удовлетворено с точностью до 1%.
5. Решить краевую задачу из упражнения 4 методом конечных
разностей. Полученную нелинейную систему алгебраических уравнений
решить итерационным методом. В качестве начального приближе^
ния принять линейную функцию между граничными значениями.
Получить последовательность решений для п = 3, 7,15, где п — чи<^
л о разбиений отрезка [0,1]. Построить графики последовательных^
приближений.
6. Решить двухточечную краевую задачу для уравнения
У" = - (1 + еУ), 0 < * < 1
2/(0) = 0,
3/(1) =
методом конечных разностей. Полученную нелинейную систему ал^
гебраических уравнений решить итерационным методом. Получить

150
Гл.З Методы решения краевых задач для систем ОДУ
последовательность решений для п = 7,15, где п — число разбиений
отрезка [0,1]. Построить графики последовательных приближений.
7. Кривая провисания каната описывается системой четырех
нелинейных ОДУ
у[ = cos(y3)
2/2 = sin(2/s)
cos(y3) -sin(2/3) |sin(2/3)|
Уъ
У4
y'4 = sin(y3) - cos(y3) | cos(y3)|
где 2/i (t) и г/2 (i) — горизонтальные и вертикальные координаты точ^
ки каната, yz{t) угол касательной каната с горизонтальной осью,
2/4 (£) напряжения в канате, 0 ^ t ^ 1 — длина дуги кривой.
Используя метод конечных разностей определить кривую
провисания каната при граничных условиях
ЙО)
м
о
о
№
0
о
V 1 /
8. Используя метод пристрелки определить кривую провисания каната
(см. упражнение 7) при граничных условиях
№>
м
о
о
Vi/
т
0.50
0
V 1 /
9. Прогиб горизонтального стержня, опертого по обеим концам и
нагруженного продольной и поперечной силами описывается
уравнением второго порядка с переменными коэффициентами
У
Л (-t2 - 1) у,
-1 < t < 1

§21. Упражнения
151
при граничных условиях
г/(-1) = о, i/(i) = о
Собственные значения и собственные функции этой задачи
определяют частоты и моды колебаний стержня. Используя конечно^
разностную дискретизацию получить алгебраическую задачу для
определения собственных значений и собственных векторов.)
10. Составить уравнения для задачи растяжения упругого стержня
с переменной площадью сечения S = Soe~x'1. Решить полученные
уравнения методом прогонки при граничных условиях и(0) = uq\
и и(1) = щ. Модуль упругости стержня Е = 2 • 105 МПа, длина
I = 20 см, 5о = 1 см2> щ = щ = 2 • 10~2 см.
11. Решить задачу из упражнения 10, когда стержень растягивается
силами Р = 1 МПа, приложенными на концах стержня. Показать, что
решение этой задачи не однозначно и определяется с точностью до
перемещения стержня как жесткого целого. Задать дополнительные
граничные условия для однозначного определения перемещения.
12. Определить методом стрельбы перемещение в составном упругом
стержне длиной 21. Площадь поперечного сечения стержня на ин^
тервале 0 ^ х ^ / переменна S = So e~xll'. Модуль упругости
постоянен Е = Eq. На интервале / < х ^ 21 площадь постоянна S = So е-1,
а модуль упругости изменятся по закону Е = Еое~х^1. Стержень
растягивается с постоянной скоростью v = Vq. Второй конец
стержня неподвижен, vo = 1 см/час, остальные параметры те же, что в
упражнении 10.
13. Уравнение равновесия Л яме для осесимметричной трубы,
находящейся под действием внутреннего давления имеет вид
d2ur 1 dur ur
dr2 r dr r2
Граничные условия в напряжениях будут
arr\r=b= -Р, (Тгг\г=а = 0

152
Гл.З Методы решения краевых задач для систем ОДУ
где Ъ — внутренний, а — внешний радиус. Материал трубы упругий:
Е = 2 • 105 МПа, v = 0.3, р = 20 МПа, а = 5 см, Ъ = 2 см.|
Аппроксимировать уравнения разностной схемой второго порядка
на равномерной сетке из пяти ячеек. Граничные условия записать
относительно перемещений и аппроксимировать производные одно^
сторонними разностями. Полученную систему разностных
уравнений решить методом прогонки.
14. Показать, что метод прогонки эквивалентен методу факторизации1
трехдиагональной матрицы уравнения Аи = f на верхне и нижне
треугольные матрицы В и С
А = ВС, Bv = f, Cu = v
где В и С — двух диагональные матрицы, решения этих систем
уравнений соответствуют прямой и обратной прогонкам.
1 Факторизацией называется представление оператора (матрицы) в виде более простых
операторов.

Литература
1. А.А.Абрамов, Вариант метода прогонки // ЖВМиМФ. т.1, No.2,
1964.
2. Н. Бахвалов, Н. Жидков, Г. Кобельков, Численные методы М-СПб:
Физматлит. 2002. 632с.
3. Ю.Б.Брызгалов, В.Н. Кукуджанов, Численное моделирование
разрушения и локализации деформаций при импульсном разрушениии\
стрежней // Математическое моделирование. 2001. т. 13, №6. с.99-
103.
4. В.В. Воеводин, Вычислительные основы линейной алгебры. Теория и
алгоритмы. М.: Наука. 1977.
5. И.М. Гельфанд, О.В. Лукациевский Метод "прогонки" для решение
разностных уравнений. Дополнение II в книге С.К.Годунов, B.C. Ря^
бенький. Введение в теорию разностных схем М.: Физматгиз. 1962.
339с.
6. С.К. Годунов, B.C. Рябенький, Разностные схемы (введение в
теорию). М.: Наука. 1973. 400с.
7. Э.И. Григолюк, В.И. Шалашилин Проблемы нелинейного
деформирования. Метод продолжения по параметру в нелинейных задачах ме\
ханики твердого деформируемого тела. М.: Наука. Гл. ред.
Физматлит, 1988. 231с.
8. Н.Н. Калиткин, Численные методы. М.: Наука. 1978.
153

154
Литература
9. В.Н. Кукуджанов, Численные методы решения нелинейных задач
механики деформируемого твердого тела. Учебн. пособие. М.: МФ^
ТИ. 1990. 95с.
10. В.Н. Кукуджанов, Разностные методы решения задач механики
деформируемых тел. Учебн. пособие. М.: МФТИ. 1992. 123с.
11. В.Н. Кукуджанов, Распространение волн в упуговязкопластических
материалах с диаграммой общего вида. // Изв. РАН. МТТ. 2001.
№ 5. с.96-111.
12. Г.И.Марчук, В.В. Шайдуров, Повышение точности решений
разностных схем. М.: Наука. 1979. 319с.
13. Г.И.Марчук, Методы расщепления. М.: Наука. 1988. 263cJ
14. Б.Е. Победря, Численные методы теории упругости и
пластичности. М.: Изд-во Мск. ун-та, 1981. 344с.
15. Ю.В. Ракитский, С.М.Устинов, И.Г. Черноруцкий, Численные
методы решения жестких систем. М.: Наука. 1979. 208с.
16. Р. Рихтмайер, К. Мортон Разностные методы решения краевых
задач. М.: Мир. 1972. 418с.
17. Б.Л. Рождественский, Н.Н. Яненко Системы квазилинейных
уравнений. Издание второе. М.: Наука. 1978. 688с.
18. А.А. Самарский, Теория разностных схем. М.: Наука. 1971. 552с.
19. А.А. Самарский, П.Н. Вабищевич Е.А. Самарская, Учебное пособие.
Задачи и упражнения по численным методам. М.: Эдиториал.
УРСС. 2000. 207с.
20. Р.П. Федоренко, Введение в вычислительную физику. Учебное
пособие для ВУЗов. М.: Изд-во Моск.Физ.-тех.Инст-та. 1994. 528cJ

Предметный указатель
аппроксимация
дифференциального оператора
разностным, 10
конечно-разностная, 10
погрешность, 14
порядок, 14
согласованная граничных условий,
схемы высокого порядка, 23
условная, 54
диагональное преобладание, 130
жесткая
система ОДУ, 139
задача
жесткая, 115
Коши для ОДУ, 21
краевая
жесткая, 111
задачи
краевые, 140
итерационный
процесс
нестационарный, 75
стационарный, 75
итерационный метод
внешний, 96
внутренний, 96
двуступенчатый, 95
двухслойный, 74
двухшаговый, 74
Зейделя, 76
нелинейный, 96
метод верхней релаксации, 76
Ньютона нестационарный, 98
Ньютона-Рафсона, 92
модифицированный, 94,
одношаговый, 74
простой итерации, 89
секущих, 95
трехслойный, 74
Якоби, 76
краевой эффект, 112
матрица
верхнетреугольная, 701
Гессе, 101, 102
нижнетреугольная, 70|
нормальная, 46
Якоби, 29, 92
матрица перехода, 45
метод
Гаусса, 68
Гаусса с выбором главного элемента,
71
квадратного корня, Т2\
квазилинеаризации, 11^
начальных параметров, ХЩ
обобщенный, 141
погружения, 119
пристрелки, 115, 117
продолжения по параметру, 119|
сопряженных градиентов, 10^
спуска, 100
наискорейшего, 102
покоординатного, 101
установления, 77
Фурье, 137
Холецкого, 72
155

156
невязка, 75, 79
область
достижимости, 89
устойчивости, 33
обусловленные матрицы
плохо, 65
хорошо, 65
оператор
интерполирования, 12
сдвига, 44
сжатия, 88
оптимальный выбор шага, 98
паразитное решение разностного
уравнения, 19
переобуславливание, 105
погранслой, 29, 112
полином Чебышева, 84
преобразование
Фурье, 45
принцип
максимума, 129
неподвижной точки, 88
уменьшения невязки, 98
принцип замороженных коэффициентов,
51
прогонка, 73, 135
дифференциальная, 119
матричная, 132
обратная, 121, 124
ортогональная, 143
по Годунову, 144
прямая, 121, 124
промежуточный слой, 40
решение
быстрорастущее, 114
сетка, 11
равномерная, 21
сеточная функция, 10
норма, 11
сжимающее отображение, 88
Предметный указатель
система
жесткая, 29
показатель жесткости, 30|
сингулярно-возмущенная, 33|
скорость сходимости
асимптотическая, 81
оптимальная, 81
средняя, 81
спектральный радиус, 46
схема
Адамса, 23, 25
Гира неявные, 31
Лакса, 38
устойчивость, 47
Лакса-Вендроффа, 40
устойчивость, 49
неявная, 41
Рунге-Кутты, 26, 27
Эйлера
неявная, 22
предиктор-корректор, 261
явная, 22
экономичная, 23
явная, 41
сходимость
квадратичная, 87, 93
линейная, 87
разностной схемы, 17
уравнение
волновое, 37
Пуассона, 126
условие
согласованности, 40
устойчивости
Куранта, 47
условия
граничные
Дирихле, 126
Неймана, 126
смешанные, 126
устойчивость
Л-устойчивость, 30, 31|

Предметный указатель
157
А (а) -устойчивость, 31
абсолютная, 51
краевых задач, 55
разностной схемы, 17
схема Гира, 32
условие Неймана, 46
условная, 51
факторизация, 70
формула
экстрополяционная Ричардсона, 15
функция
выпуклая, 101
число обусловленности, 64
шаблон
2-х слойный, 38
многослойный, 38
сетки, 15
шаг
выбор для волнового уравнения, 58
выбор по неявной схеме, 56
экстраполяция, 24