Введение в вычислительную физику

ФЕДОРЕНКО Р.П.

ВВЕДЕНИЕ 

В ВЫЧИСЛИТЕЛЬНУЮ 

ФИЗИКУ

2008

ФИЗТЕХОВСКИЙ УЧЕБНИК

Второе, исправленное и дополненное издание

Под редакцией и с дополнениями А.И. Лобанова

УДК 519.63(075.8)
ББК 22.31
Ф33

Федоренко Р. П.
Ф33
Введение в вычислительнуюфизику: Учебное пособие
для вузов / Р. П. Федоренко / Под ред. А. И. Лобанова. —
2-е изд., испр. и доп. — Долгопрудный: Издательский
Дом <Интеллект>, 2008. — 504 с.
ISBN 978-5-91559-011-2

Книга посвящена описанию методов приближенного решения задач
математической физики, возникающих в различных областях. Изложение основных понятий и средств численного анализа доводится
до описания специальных алгоритмов решения важных прикладных задач, разработка которых продолжается в настоящее время.
Приближенные решения сложных задач получаются как общими
средствами вычислительной математики, так и специфическими для
данного узкого класса задач приемами, которые позволяют обходить
существенные трудности в современной вычислительной работе и
делают расчеты посильными для ЭВМ.
Для студентов и аспирантов факультетов прикладной математики и физико-технических специальностей вузов с достаточно
высоким уровнем преподавания математики, а также для научных
работников, специализирующихся в области применения численных методов в научных исследованиях.
ББК 22.31
УДК 519.63(075.8)

ISBN 978-5-91559-011-2
© 2008, Федоренко Р. П.
© 2008, ООО Издательский Дом
<Интеллект>, оригинал-макет,
оформление

ОГЛАВЛЕНИЕ

Предисловие ко второму изданию
. . . . . . . . . . . . . . . . . . . . . . . .
4
Предисловие автора
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5

Ч А С Т Ь
П Е Р В А Я

ОСНОВЫ ВЫЧИСЛИТЕЛЬНОЙ МАТЕМАТИКИ
. . . . . . . . . . . . . . .
9

§
1. Решение систем нелинейных уравнений . . . . . . . . . . . . . . . . . .
9
§
2. Численное дифференцирование
. . . . . . . . . . . . . . . . . . . . . .
23
§
3. Интерполяция функций
. . . . . . . . . . . . . . . . . . . . . . . . .
26
§
4. Вычисление определенных интегралов
. . . . . . . . . . . . . . . . . .
45
§
5. Численное интегрирование задачи Коши для систем обыкновенных дифференциальных уравнений
. . . . . . . . . . . . . . . . . . . . . . . . . . .
54
§
6. Абстрактная форма приближенного метода
. . . . . . . . . . . . . . . .
60
§
7. Исследование сходимости методов Рунге—Кутты
. . . . . . . . . . . . . .
65
§
8. Приближенное решение краевых задач для систем обыкновенных дифференциальных уравнений
. . . . . . . . . . . . . . . . . . . . . . . . . . .
78
§
9. Метод дифференциальной прогонки
. . . . . . . . . . . . . . . . . . . .
86
§ 10. Прогонка в разностной задаче Штурма—Лиувилля
. . . . . . . . . . . .
91
§ 11. Численное интегрирование задачи Коши для уравнений с частными производными
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
97
§ 12. Спектральный признак устойчивости
. . . . . . . . . . . . . . . . . . .
111
§ 13. Метод переменных направлений
. . . . . . . . . . . . . . . . . . . . .
130
§ 14. Решение эллиптических задач методом сеток
. . . . . . . . . . . . . . .
137

Ч А С Т Ь
В Т О Р А Я

ПРИБЛИЖЕННЫЕ МЕТОДЫ ВЫЧИСЛИТЕЛЬНОЙ ФИЗИКИ
. . . . . . . .
175

§ 15. Спектральная задача Штурма—Лиувилля
. . . . . . . . . . . . . . . . .
175
§ 16. Главная спектральная задача для краевых задач математической физики
. .
184
§ 17. Жесткие системы обыкновенных дифференциальных уравнений
. . . . . .
201
§ 18. Жесткие линейные краевые задачи
. . . . . . . . . . . . . . . . . . . .
234
§ 19. Осреднение быстрых вращений
. . . . . . . . . . . . . . . . . . . . . .
252
§ 20. Одномерные уравнения газовой динамики и их численное интегрирование
272
§ 21. Нелинейное уравнение теплопроводности
. . . . . . . . . . . . . . . . .
299
§ 22. Реализация разностной схемы для уравнений газовой динамики с теплопроводностью
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
311
§ 23. Приближенное решение двумерных задач газовой динамики
. . . . . . . .
330
§ 24. Приближенное интегрирование уравнения Власова
. . . . . . . . . . . .
362
§ 25. Некорректные задачи и их приближенное решение
. . . . . . . . . . . .
376
§ 26. Поиск минимума
. . . . . . . . . . . . . . . . . . . . . . . . . . . .
392
§ 27. Дифференцирование функционалов
. . . . . . . . . . . . . . . . . . . .
415
§ 28. Задачи оптимального управления
. . . . . . . . . . . . . . . . . . . . .
434
§ 29. Вариационные задачи механики с недифференцируемыми функционалами
449
§ 30. Псевдодифференциальные уравнения
. . . . . . . . . . . . . . . . . . .
466
§ 31. Метод конечных суперэлементов
. . . . . . . . . . . . . . . . . . . . .
478

Список литературы
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
494
Библиографический комментарий
. . . . . . . . . . . . . . . . . . . . . . . .
499

ПРЕДИСЛОВИЯ

ПРЕДИСЛОВИЕ КО ВТОРОМУ ИЗДАНИЮ

Вашему вниманию предлагается второе издание книги Радия
Петровича Федоренко <Введение в вычислительную физику>. Со времени появления первого издания прошел довольно значительный срок.
Книга приобрела популярность и среди студентов, и среди вычислителей-практиков. Она сразу обратила не себя внимание и своеобразным
авторским подходом, и стилем изложения. Дело здесь не только в том,
что в основу положены лекции, которые читались студентам факультета
проблем физики и энергетики, а затем — общей и прикладной физики
МФТИ. Когда создавался и <обкатывался> предлагаемый курс, учебным
планом был предусмотрен вводный (семестровый) курс вычислительной математики (поэтому первая — элементарная — часть книги для
современного читателя может показаться неожиданно короткой). Затем
следовало продолжение — годовой курс вычислительной физики, являющийся мостиком от элементарного, модельного изложения к реальным
физическим задачам. Сейчас они могут показаться довольно простыми,
но они замечательно иллюстрируют все трудности, которые и сейчас с
той или иной степенью успеха преодолеваются при решении сложных
нелинейных задач. А мост этот перекинут с гладкого и неоднократно
пройденного берега линейных задач на берег нелинейный, высокий, обрывистый и не всегда гладкий.
Материал книги (особенно второй части) во многом опирается на оригинальные работы автора. Отметим, что Радий Петрович, кроме огромного
числа решенных задач, предложил три замечательных подхода, развиваемых
сейчас во всем мире. Во-первых, это концепция гибридных разностных
схем, в дальнейшем развившаяся в TVD-схемы, другие методы монотонизации разностных схем высокого порядка аппроксимации. Во-вторых,
многосеточный итерационный метод решения уравнений эллиптического
типа, за разработку которого Радий Петрович вместе с Н. С. Бахваловым,
Г. П. Астраханцевым и В. В. Шайдуровым в 2004 году был удостоен Государственной премии Российской Федерации. В-третьих, это развиваемый в
настоящее время метод граничных суперэлементов. Вместе с тем, в книге делается основной упор именно на идейную сторону методов, на связь
математики (иногда весьма сложной) с физической постановкой задач.
При подготовке настоящего издания сохранены содержание и структура изложения, принятые автором. Исправлены некоторые опечатки и

Предисловие автора
5

неточности первого издания. Наибольшим изменениям подверглась первая часть книги. Вторая часть, основанная на личных результатах Радия
Петровича, осталась практически без изменений, за исключением дополнений по жестким системам ОДУ.
Желаем Вам приятного чтения и освоения идей курса и некоторых
методов вычислений.

Заведующий кафедрой вычислительной математики МФТИ
член-корреспондент РАН А. С. Холодов

ПРЕДИСЛОВИЕ АВТОРА

Предлагаемая вниманию читателя книга написана на основе двух курсов лекций, в течение ряда лет читавшихся студентам
Московского физико-технического института. Им соответствуют две части книги. Первая часть содержит основы вычислительной математики
(такой курс слушают студенты всех факультетов). Вторая часть соответствует годовому курсу вычислительной физики (на факультете общей и
прикладной физики).
Почему книга называется <Введение в вычислительную физику>, а
не <Методы вычислительной математики>, например? Это объясняется
характером будущей работы слушателей. Для них вычислительная математика в первую очередь будет инструментом научных исследований, а
не их предметом. Методы приближенныхвычисленийизлагаютсяв книге не
как самостоятельнаянаучная дисциплина,а как набор средств,позволяющих
продвинуться в исследовании тех или иных прикладных проблем физики,
химии, аэромеханики и т. п. Это соответствует характеру образования,
получаемого в Московском физико-техническом институте, и научному
стилю Института прикладной математики им. М. В. Келдыша. Работа автора
в этом институте определила его понимание науки, называемой <вычислительная математика>, и нашла отражение как в содержании книги,
так и в характере изложения.
Не следует думать, что физик-вычислитель обречен лишь на пассивное использование средств, развиваемых математиками. История развития методов приближенных вычислений ясно показывает большую,
можно сказать, определяющую, роль решения именно частных прикладных задач, для которых известные методы оказываются неэффективными
в силу каких-то специфических особенностей. Наличие важных приложений оправдывает выделение такого (неестественно вырожденного с
общематематической точки зрения) класса задач в самостоятельный объект, заслуживающий отдельного углубленного изучения, а привлечение

Предисловие автора

содержательной интуиции и неформализованных знаний той прикладной области, в которой возникла задача, помогает понять ее специфику
и разработать эффективный метод решения. Эти же знания существенно используются для контроля приближенных решений задачи. В этой
связи с практикой — сила физика, позволяющая ему часто решающим
образом влиять на развитие вычислительной математики. Однако в этом
же и его слабость: нередко такой специалист воспринимает свою задачу
слишком обособленной, понятной только ему и не имеющей никакого
отношения к общематематической теории.
В современных методах приближенных вычислений можно выделить
методы, имеющие широкое применение и уже ставшие достоянием математической теории, и методы, развитые для специальных, но важных
в приложениях классов задач. Этому делению и соответствуют две части
книги. Первая часть по содержанию близка к традиционным курсам численных методов, однако отбор материала, внимание, уделяемое тем или
иным вопросам, и характер изложения определяются в первую очередь
местом, которое эти вопросы занимали в практике автора и его коллег по Институту прикладной математики. В частности, относительно
небольшое место отведено таким сильно развитым разделам, как теория интерполяции и квадратурные формулы, а вычислительные методы
общей линейной алгебры совсем не отражены в книге. Это объясняется обилием стандартных программ и руководств, отражающих развитие
теории соответствующих разделов вычислительной математики.
В ходе изложения мы не стремимся к максимальной общности и
безупречной строгости формулировок. Современный стиль изложения
математических результатов требует четкой и полной формулировки всех
используемых в доказательстве предположений о свойствах встречающихся функций. Среди них можно выделить две характерные группы.
Первая группа условий строго оговаривает свойства общего, типичного
характера, отсутствие которых с прикладной точки зрения является исключением, редко встречающимся вырождением. Вторая группа условий
выделяет специальный, частный случай, рассмотрение которого оправдано наличием важных и интересных приложений. Именно такие условия
мы считаем необходимым выделять, обсуждать и комментировать. Условия первой группы обычно используются <неявно>. К ним, в частности,
относятся предположения о гладкости функций. Вместо строгого оформления таких условий в тексте часто используется термин <гладкая функция>,
означающий функцию, ограниченную вместе со своими производными того порядка, который используется в выкладках. При этом предполагается,
что функция является гладкой в той части пространства, с которой мы имеем дело при решении задачи (т. е. <там, где нам это нужно>).
Такое отступление от педантичного стиля современной математической
литературы
представляется
соответствующим
духу
<вычислительной

Предисловие автора
7

физики>. В свое время, прослушав аккуратный университетский курс
обыкновенных дифференциальных уравнений, в котором теорема существования была изложена со всеми необходимыми предположениями,
автор вынес впечатление, что решения этих уравнений если и существуют,
то только в виде редкого счастливого исключения и на очень малом интервале времени, который иногда можно и продолжить. Таков эффект
стиля, при котором все ограничения перечисляются <на равных правах> и не комментируются с некой не очень-то понятной и однозначной
<прикладной> точки зрения.
Достаточно большое внимание уделяется прикладному комментарию
к некоторым теоремам. Этим формируется своеобразное <прикладное
мировоззрение> читателя, его будущие взаимоотношения с теоретическими исследованиями. Дело в том, что алгоритм приближенного решения
сложной задачи математической физики практически никогда не бывает
строго обоснованным в том смысле, который придает этому слову математика. Полезные теоретические результаты, как правило, относятся
к выделенным из него идеализированным фрагментам. Использование
строгих результатов в практической вычислительной работе — своего
рода искусство, в котором результат оправдывает средства. Это характеризует <вычислительную физику> как науку в значительной мере
экспериментальную. Ее взаимоотношения с чистой теорией достаточно
сложны и неоднозначны.
Вторая часть книги точнее соответствует содержанию термина <вычислительная физика>. В ней собраны описания методов приближенного решения частных задач, имеющих, однако, важные области приложения
в современной науке и технике. Каждый параграф посвящен одной
из таких задач. Принят следующий способ изложения. Вначале дается
замкнутая математическая формулировка задачи, указывается ее <прикладное происхождение>. Физическая терминология используется для
<оживления> изложения, но никакого физического обоснования постановки задачи не проводится — это дело физика, а не вычислителя.
При этом указываются те особые обстоятельства, которые делают задачу нестандартной, требующей разработки специальных вычислительных
методов.
Затем описывается метод приближенного решения, оказавшийся достаточно эффективным. Основное внимание уделяется именно тем деталям метода, которые учитывают специфическую нестандартность данной
задачи и которым метод обязан своей эффективностью. Попутно обсуждаются те трудности, с которыми сталкиваются при стандартном
подходе к задаче (формально не только возможном, но иногда даже
строго обоснованном). В некоторых случаях приводятся и обсуждаются
характерные численные результаты. Стандартные детали вычислительной методики описываются бегло, а иногда и совсем опускаются.

Предисловие автора

Материал второй части книги несет двойную нагрузку. Во-первых,
описываемые задачи достаточно интересны в приложениях и опыт их
успешного решения представляет прямой интерес в связи с задачами
именно этого типа. Во-вторых, разработка эффективного алгоритма частной задачи обычно связана с использованием приемов, имеющих более
широкое, выходящее за рамки данной задачи значение. Автор предпочитает знакомить читателя с такими приемами на примерах конкретных
задач, в которых они были использованы с большим эффектом. Есть
и другой путь — выделить эти приемы как отдельные самостоятельные
сюжеты, дать абстрактное описание ситуаций, в которых их применение
целесообразно. Подобный способ изложения представляется нам чуждым
духу вычислительной физики.
Отметим некоторые технические детали изложения. Текст книги разбит на параграфы, каждый из которых имеет свою нумерацию формул.
При ссылке на формулы другого параграфа используется двойной номер
(параграфа и формулы). Впрочем, автор стремился свести к минимуму
подобные ссылки. В тексте опускаются и библиографические ссылки.
Этот недостаток компенсируется библиографическим комментарием, тем
более необходимым, что во многих местах излагаются результаты, еще
не вошедшие прочно в учебную литературу и часто освещенные лишь в
журнальных, а то и препринтных публикациях. Курсивом в тексте выделены общеупотребительные термины вычислительной математики.
Оба курса, на основе которых написана эта книга, читались по
предложению академика О. М. Белоцерковского, много сделавшего для
внедрения компьютерных наук в <систему физтеха>. Пользуюсь случаем
высказать Олегу Михайловичу свою искреннюю благодарность.
Автор должен отметить и неоценимое влияние, оказанное на него
коллегами по Институту прикладной математики им. М. В. Келдыша.
Нет возможности упоминать их здесь, автор постарался должным образом отметить их вклад в развитие предмета книги в библиографическом
комментарии к списку литературы. Это будет полезно и для будущих
историков науки, которым рано или поздно предстоит изучать историю
становления отечественной вычислительной математики.

Ч А С Т Ь
ПЕ Р ВА Я
ОСНОВЫ ВЫЧИСЛИТЕЛЬНОЙ
МАТЕМАТИКИ

§ 1.
Решение систем нелинейных уравнений

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

f(x)=0,
(1)

имеяввиду,что x естьточка n-мерного пространства, т. е. x={x1, x2, . . . , xn},
а f — n-мерная вектор-функция, т. е. f ={f1, f2, . . . , fn}. Таким образом, (1) есть система n уравнений с n неизвестными:

f1(x1, x2, . . . , xn)=0,
. . . . . . . . .
fn(x1, x2, . . . , xn)=0.
(2)

Конечно, когда речь идет о нелинейных уравнениях в общем случае,
мы не имеем ни теорем о существовании, ни теорем о единственности
решения. Тем не менее, имея дело с системой (1), мы предполагаем,
что искомое решение существует. Оно, быть может, не единственно, и
методы, которые будут рассмотрены ниже, не имеют целью найти все
решения; обычно достаточно будет какого-то одного. Более того, мы
предполагаем, что из каких-то содержательных соображений нам известно примерное расположение этого решения, некоторая не очень большая
область, в которой оно находится. Таким образом, лучше говорить не
о решении системы нелинейных уравнений, а об уточнении имеющегося весьма грубого приближения к некоторому решению. Ниже это будет
должным образом конкретизировано.

Метод Ньютона. Основная вычислительная конструкция, применяемая для решения системы (1), по традиции приписывается Ньютону,
хотя теоретические исследования этого алгоритма были выполнены лет
на сто позже (Фурье, Коши). Основу метода составляет фундаментальная для вычислительной математики конструкция — метод итераций
(последовательных приближений) и линеаризация уравнений.
В методе Ньютона, начиная с некоторого начального приближения
x0, последовательно находятся точки x1, x2, . . . , xk, . . . таким образом,
что lim
k→∞ f(xk)=0, a lim xk=х∗, где х∗ — решение системы (1).

Часть первая. Основы вычислительной математики

Нужно только иметь в виду, что вычислительную математику интересует не только факт сходимости, но и скорость сходимости. Метод
Ньютона особенно ценен тем, что обеспечивает очень высокую, как говорят, <квадратичную> скорость сходимости (точный смысл этого термина
выяснится позже, после доказательства соответствующей теоремы).
Рассмотрим стандартный шаг итерационного процесса метода Ньютона.
Пусть имеется некоторое уже найденное приближение xk; следующее
приближение xk+1 ищем в виде xk+1=xk+δх, где δх — малая поправка, уточняющая xk. Для ее определения выпишем уравнение f(xk+δx)=0.
Само по себе оно не проще исходного уравнения (1), но, используя
предположение о малости δх, его можно линеаризовать, т. е. использовать разложение f по δх с точностью только до членов первого порядка:

f(xk+δx)=f(xk)+fх(xk)δх+O(∥δх∥2).

Пренебрегая членами O(∥δх∥2), получаем линеаризованное уравнение
для δх:

f(xk)+fx(xk)δx=0,
(3)

которое уже решается, и можно выписать его явное решение:

δx=−f−1
x (xk)f(xk).

Итак, алгоритм метода Ньютона (MH) имеет следующую форму:
1) имеется некоторое уже найденное приближение х;
2) вычисляются вектор f(x) и матрица fx(x);
3) решается система (линейных) уравнений (3);
4) пересчитывается приближенное решение х:=х+δх.
Далее процесс повторяется циклически до получения достаточно малой
величины ∥f(x)∥.
Прежде чем перейти к теоретическому исследованию, рассмотрим
некоторые связанные с методом вопросы.
1. Что такое fx? Это есть производная вектор-функции по векторному
аргументу. Точный смысл fх определяется первым членом ряда Тейлора

f(х+δx)≈f(x)+fx(x)δx.

Это — векторная форма записи отрезков ряда Тейлора для всех n
функций:

f1(x1, . . . , xn)+ f 1

x1
δx1+ f1

x2
δx2+. . .+ f1

xn
δxn,
. . . . . . . . . . . . . . . . . . . . . .

fn(x1, . . . , xn)+ f n

x1
δx1+ f n

x2
δx2+. . .+ fn

xn
δxn.

§ 1. Решение систем нелинейных уравнений
11

Отсюда следует расшифровка:

fx=

0
44444444444444444444444442

f 1
x1
. . .
f 1
xn
.
.
.
.
.
.
f n
x1
. . .
f n
xn

1
55555555555555555555555553

2. В теоретических оценках будем использовать явную форму

δx=−f−1
x f,

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

Рис. 1

задачи n сравнительно невелика. В дальнейшем мы встретимся с ситуациями, когда n настолько велико, что применение
стандартных методов решения линейных
систем невозможно или по крайней мере
нерационально. Обычно в таких ситуациях матрица fx имеет специальную структуру, и часто удается построить специальные
методы решения, существенно более эффективные, чем общие.
При n=1 метод Ньютона известен как
<метод касательных>. Этот термин поясняет рис. 1. Линеаризация состоит в том, что
кривая у=f(x) заменяется касательной, проведенной в уже найденной
точке xk, а в качестве следующего приближения xk+1 берется пересечение касательной с осью абсцисс. На рис. 1 показано несколько таких
итераций и, естественно, возникает предположение о том, что процесс
должен быть достаточно эффективным. Это и подтверждается более точным исследованием.

Сходимость метода Ньютона. Докажем теорему о квадратичной сходимости метода.
Т е о р е м а 1. Пусть х∗ — решение системы (1). Предположим, что
в некоторой окрестности х∗:
а) f(x) является гладкой функцией в том смысле, что существуют ее
производные до второго порядка и имеет место оценка ∥f xx(x)∥≤C2;
б) отображение х→f(x) равномерно невырождено в том смысле, что
f−1
x (x) существует и ограничена: ∥f−1
x (x)∥≤C1.
Тогда, если начальное приближение x0 достаточно близко к х∗, метод
Ньютона сходится и имеет квадратичную скорость сходимости.