Введение в предмет

Исторически, математика появилась в первую очередь как ответ на необходимость измерять те или иные величины окружающего мира и производить над ними вычисления. Однако оказалось, что некоторые, даже сравнительно простые вычислительные задачи, не так легко решить. Например, вавилонская глиняная табличка YBC 7289 (1800-1600 до н.э.) содержит в себе сравнительно точное приближение иррационального числа $\sqrt{2} \approx 1.414222$, которое является длиной диагонали квадрата с единичной стороной. Для нахождения этой аппроксимации вавилоняне использовали одну из вариаций численного метода решения алгебраического нелинейного уравнения $x^2 = 2$, который сейчас известен как метод простой итерации (он будет пройдет в курсе). Таким образом уже древние вавилоняне столкнулись с фундаментальной проблемой, которая породила всё направление вычислительной математики -- сформулировать математическую задачу еще не значит ее решить. Проходя сквозь века математическая наука совершенствовалась и научилась формулировать свои модели не только в виде алгебраических уравнений, но и в виде обыкновенных дифференциальных уравнений, уравнений в частных производных, интегральных уравнений и проч. Развивались и методы нахождения аналитического решения математических задач -- например, Эйлер изобрел метод решения однородных линейных дифференциальных уравнений, а Фурье обнаружил, что некоторые линейные уравнений в частных производных могут быть решены с помощью бесконечных рядов косинусов и синусов. За вычетом нескольких исключительных случаев (например, корни многочленов до 4 степени включительно или уравнения Бернулли), математикам удавалось находить точные аналитические решения только для линейных уравнений. Однако, как показала практика, многие интересные физические явления описываются нелинейными уравнениями. Так, например, динамика вязкой несжимаемой жидкости описывается уравнением Навье--Стокса, которое за счет члена $\bm{u} \cdot \nabla\bm{u}$, где $\bm{u}$ – поле скоростей потока, является нелинейным уравнением в частных производных. Для этого уравнения, дополненного уравнением непрерывности, не только не найдено аналитического решения в общем виде, но даже не доказано, что оно в принципе существует и является гладким для гладких начальных условий (это доказательство является одной из "задач тысячелетия", сформулированных институтом Клэя). В таком случае ученым приходится полагаться исключительно на приближенные решения, полученные либо аналитически с помощью упрощения оригинального уравнения (такой подход используется, например, в асимптотических методах или слабонелинейном анализе), либо с помощью численных методов. Очевидно, что упрощение оригинального уравнения сильно ограничивает валидность полученного решения, что и мотивирует активное использование численных методов при решении современных задач математической физики.

Величайшие ученые своего времени были озадачены вопросами численного решения математических задач – по мере прохождения курса мы встретим фамилии Ньютона, Лагранжа, Эйлера, Гаусса и многих других известных математиков. Нельзя не сказать о вкладе русских ученых в развитие численных методов – полиномы Чебышёва, подпространство Крылова, метод Галеркина, метод Годунова хорошо известны современным ученым по всему миру. Глобально вычислительные методы можно разбить на три больших группы (по крайней мере в рамках нашего курса):

1. методы приближения (аппроксимации) дискретных данных, принадлежащих некоторой сложной, часто недоступной, функции, к сравнительно простым аналитическим функциям;
2. численные методы линейной и нелинейной алгебры.
3. вероятностные численные методы;

Первая группа исторически развивалась раньше двух других и была нацелена на нахождение решений нелинейных уравнений с помощью аппроксимации решения набором простых функций, что сводило нелинейную задачу к более простой линейной (под "сложной" функцией подразумевается сильно нелинейная или зашумленная функция). Линейная задача чаще всего выражена с помощью СЛАУ, и подавляющая часть методов второй группы посвящена точным и приближенным методам решения СЛАУ. Приближенные методы решения СЛАУ появились вследствие того, что с каждым днем размерность вычислительных задач $N$ увеличивалась, и решение СЛАУ, например, методом Гаусса не всегда было реализуемо ввиду его большой сложности $O(N^3)$. Увеличение размерности задач так же мотивировало создание многих методов из третьей группы. Например, решение задачи многих тел в квантовой механике сводится к вычислению значения интеграла, кратность которого может достигать нескольких тысяч. В таких случаях вместо классических квадратур численного интегрирования используется метод Монте-Карло, где подынтегральная функция вычисляется в $M$ случайных точках, сгенерированных с помощью заданного закона распределения, после чего их значения суммируются. Это позволяет уменьшить зависимость сложности от размерности задачи.

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

• в методе конечного элемента (МКЭ), используемого для численного решения уравнений в частных производных, функции формы аппроксимируют решение уравнения в элементе с помощью полиномов Лагранжа или Эрмита (классический МКЭ) или полиномов Чебышёва (метод спектральных элементов);
• при обучении нейронных сетей, локальная оптимизация сильно нелинейной целевой функции ошибки происходит с помощью вариации метода градиентного спуска (стохастический градиентный спуск).

Для успешного прохождения курса представляется важным освежить воспоминания из прошлых курсов:

• математический анализ:
• понятия предела, непрерывности, производной и интеграла Римана;
• теорема о среднем значении;
• теорема о промежуточном значении;
• теорема Ньютона-Лейбница;
• ряды Тейлора;
• линейная алгебра:
• понятия линейного нормированного пространства, матрицы, определителя, обратной матрицы;
• собственные числа и собственные вектора матрицы;
• свойства положительно определенных и симметричных матриц;
• простейшие методы решения систем линейных алгебраических уравнений (СЛАУ);

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

1. Численные методы. Н.Н. Калиткин. Главная редакция физико-математической литературы изд-ва <<Наука>>, М., 1978.
• комментарий: оптимальный учебник для первого знакомства с численными методами (многие примеры и описания в курсе взяты из него).
2. Численные методы / Н.С. Бахвалов, Н.П. Жидков, Г.М. Кобельков. -- 7-е изд. -- М. : БИНОМ. Лаборатория знаний, 2017. -- 636 с. : ил. -- (Классический университетский учебник).
• комментарий: полезен для тех случаев, когда необходимо разобраться со строгим теоретическим обоснованием численных методов.
3. Элементы теории функций и функционального анализа. А.Н. Колмогоров, С.В. Фомин. Главная редакция физико-математической литературы изд-ва <<Наука>>, М., 1976.
• комментарий: удобен для подробного теоретического знакомства с функциональными пространствами, интегралом Лебега и тригонометрическими рядами.

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

1. Numerical Analysis, Ninth Edition, Richard L. Burden and J. Douglas Faires.
• комментарий: полный, подробный и вместе с тем доступный курс численного анализа.
2. Analysis of Numerical Methods, Eugene Isaacson and Herbert Bishop Keller (1994).
• комментарий: курс численного анализа, ориентированный на строгое доказательное обоснование свойств численных методов.
3. Introductory Functional Analysis with Applications, Erwin Kreyszig.
• комментарий: курс функционального анализа с подробным разбором примеров использования функциональных пространств в различных прикладных задачах.
4. The Elements of Real Analysis, Robert G Bartle.
• комментарий: углубленный курс математического анализа функций вещественного переменного.
5. Essential Topology, Martin D Crossley.
• комментарий: книга доступно излагает базовые понятия топологии, которая, будучи наукой о непрерывных преобразованиях, представляется полезной для освоения в контексте курса численных методов.
6. Private lecture notes by Prof Mark Kelmanson and Dr Kevin Gouder.
• комментарий: эти лекции отсутствуют в открытом доступе, однако соответствующие конспекты представляются ценными и могут быть предоставлены по запросу.

Моделирование объектов исследования

Конечной целью прикладной математики, а следовательно и вычислительной математики как ее подраздела, является моделирование некоторых объектов и процессов окружающего мира (далее – объектов исследования). В математике под моделированием подразумевается описание самых важных и существенных свойств объекта исследования в математических понятиях. Допустим, например, что объектом исследования является футбольный мяч. Если из всех его свойств нам больше всего интересна геометрическая форма, то его моделью будет сфера как математический объект. Однако, если нас интересует его деформация при сжатии руками (в этом случае объектом исследования формально является процесс сжатия мяча, а не сам мяч), то модель будет сформулирована как система дифференциальных и алгебраических уравнений статической задачи линейной теории упругости. В обоих случаях полученная модель называется математической моделью объекта исследования. Математические модели зачастую настолько сложны, что для них не существует способа нахождения аналитического (точного) решения ("сложными" могут быть как разрешающие уравнения, так и граничные условия). В таком случае может ли сама математическая модель быть объектом исследования? Безусловно, так как она в свою очередь тоже может быть смоделирована. К примеру, система уравнений упругости, приведенная на рисунке выше, может быть дискретизирована с помощью метода конечных элементов. Полученная модель будет представлять собой алгоритм по сборке СЛАУ и саму СЛАУ. Мы будем называть подобную модель дискретной моделью. Принципиально, точное решение такой модели может быть найдено, однако для этого потребуется неоправданно большое количество времени и людей. В современном мире для нахождения решений дискретных моделей используются компьютеры, что вынуждает трансформировать дискретную модель в компьютерную программу, называемую также компьютерной моделью. Результирующая диаграмма подобного моделирования показана на рисунке \ref{fig:modelling_diagram}. Заметим, что решения всех трех моделей в общем случае не равны и могут лишь сходиться друг к другу при определенных обстоятельствах.

Приближенный анализ

Очевидно, что упрощения, осуществляемые на каждом шаге моделирования, вносят определенную неточность в решение относительно изначального объекта исследования, которую называют погрешностью. Принято выделять четыре вида погрешностей: неустранимая погрешность, погрешность математической модели, погрешность метода и вычислительная погрешность. Дадим им определения. Неустранимой погрешностью называют неточность при задании исходных данных. Допустим, вы бросаете шар с Пизанской башни. В независимости от того, как вы моделируете законы природы, для того, чтобы рассчитать скорость шара в момент касания поверхности Земли, вам необходимо знать высоту, с которой вы отпускаете шар в свободное падение. В случае, если вы измеряете высоту рулеткой с миллиметровой шкалой, погрешность измерения составляет $O(10^{-3})$ метра. Тогда, даже если вы способны идеально точно предсказывать скорость шара, погрешность в $O(10^{-3})$ метра устранить не удастся. Погрешностью математической модели называют неточность при описании реального объекта математическими понятиями. Математическая модель, описывающая движение тела в свободном падении, выражается простым уравнением Ньютона

\begin{equation} \label{eq:newton_law} m \frac{d^2 x}{dt^2} = -m g \end{equation}

где $m$ – масса тела, $g$ – ускорение свободного падения, $x(t)$ – высота тела в момент времени $t$. Применительно к земной атмосфере, погрешность этой математической модели возникает в результате пренебрежения сопротивлением воздуха. Более точной моделью была бы:

\begin{equation} \label{eq:newton_law_better} m \frac{d^2 x}{dt^2} = -m g - A \text{sign} \left(\dfrac{dx}{dt}\right) \dfrac{dx}{dt}^2 \end{equation}

Очевидно, что погрешность растет пропорционально квадрату скорости, что делает уравнение \eqref{eq:newton_law} неприменимым, скажем, для моделирования движения спускаемого аппарата сквозь атмосферу. Погрешностью метода называют неточность при замене математической модели приближенной. Уравнения \ref{eq:newton_law} и \ref{eq:newton_law_better} имеют аналитические решения, однако в случае более сложной правой части нам необходимо будет использовать один из численных методов решения обыкновенных дифференциальных уравнений (ОДУ), например, метод Рунге-Кутта, который мы будем проходить в курсе. В таком случае погрешностью метода будет разница между точным решением и решением уравнения, дискретизированного методом Рунге-Кутта. Вычислительной погрешностью называют погрешность математических операций, производимых компьютером. С математической точки зрения вычислительная погрешность появляется в результате того, что при программировании численного метода, алгебраические структуры которого определены для поля вещественных чисел $\mathbb{R}$, мы неявно заменяем их на алгебраические структуры, определенные для поля рациональных чисел $\mathbb{Q}$. Это связано с формой представления (формально, приближения) вещественных чисел в памяти компьютера. Стандарт IEEE 754-2008 гласит, что 64-битное представление вещественного числа состоит из последовательно расположенных бита знака $s$, 11 бит экспоненты (порядка) $c$ и 52 бит мантиссы $f$, а само значение числа вычисляется как где $f = \displaystyle\sum_{i = 1}^{52} 2^{-b_i}$ и $b_i \in \{0, 1\}$ – $i$-ый бит мантиссы. Возникающую погрешность при вычислительных операциях часто называют погрешностью округления.

Так как процесс моделирования реального объекта является последовательным (сначала составляется математическая модель, затем дискретная модель, а после компьютерная), логично заключить, что погрешность, вносимая на очередном этапе моделирования не должна быть меньше, чем погрешность предыдущих этапов. Вернемся к примеру с бросанием тела с Пизанской башни. Допустим, мы используем такой численный метод для интегрирования уравнений движения тела, что погрешность вычисляемой высоты тела получилась $O(10^{-5})$ метра. Очевидно, что при погрешности линейки, с помощью которой мы измеряли начальную высоту, равной $O(10^{-3})$ метра эти усилия окажутся напрасными. Определим, что мы формально называем погрешностями. Абсолютной погрешностью приближенного значения $a^{*}$ называют величину $\Delta(a^{*})$, которая определена как

\begin{equation} \Delta(a^{*}) = |a - a^{*}|, \end{equation}

где $a$ – точное значение. Число $a$ записывают с учетом абсолютной погрешности в следующей форме: Относительной погрешностью приближенного значения $a^{*}$ называют величину $\delta(a^{*})$, которая определена как

\begin{equation} \delta(a^{*}) = \left|\frac{a - a^{*}}{a}\right|, \end{equation}

где $a$ – точное значение. Число $a$ записывают с учетом относительной погрешности в следующей форме:

Некоторые понятия функционального анализа}

Понятия абсолютной и относительной погрешности определены для некоторых приближенных величин, которые неявно предполагаются принадлежащими одному из числовых множеств (чаще всего $\mathbb{R}$ или $\mathbb{Q}$). Однако понятие близости распространяется на куда более широкий класс множеств и, в частности, на множества функций, что особенно важно при доказательстве сходимости численных методов. Соответствующей обобщенной мерой "близости" называется метрика. Множество $X$ называется метрическим пространством, если на нем определена функция $\rho: X \times X \longrightarrow \mathbb{R}$, называемая метрикой или расстоянием, для которой выполняются следующие аксиомы:

1. $\rho(x_1, x_2) \geq 0$,
2. $\rho(x_1, x_2) = 0 \Longleftrightarrow x_1 = x_2$,
3. $\rho(x_1, x_2) = \rho(x_2, x_1)$,
4. $\rho(x_1, x_3) \leq \rho(x_1, x_2) + \rho(x_2, x_3) \text{~(аксиома треугольника)}$,

Метрическим пространством же чаще называют пару $(X, \rho)$. Последовательность $\{x_n\}_{n=1}^{\infty}$ метрического пространства $(X, \rho)$ называется сходящейся к элементу $x \in X$, если для нее верно

$$ n \to \infty \Longrightarrow \rho(x_n, x) \to 0. $$

Последовательность $\{x_n\}_{n=1}^{\infty}$ метрического пространства $(X, \rho)$ называется фундаментальной, если $\forall n > 1, \epsilon > 0\ \exists k(\epsilon): \rho(x_n, x_m) < \epsilon, m > k$ Фундаментальную последовательность также называют сходящейся в себе последовательностью или последовательностью Коши. Метрическое пространство называется полным, если любая последовательность его элементов сходится к элементу того же пространства. Сложность при анализе вычислительной погрешности в частности возникает из-за того, что множество рациональных чисел $\mathbb{Q}$ является неполным. К примеру, последовательность сходится к $e \in \mathbb{R}$.

На протяжении всего курса мы практически всегда будем иметь дело с элементами линейных (векторных) пространств, т.е. множеств с определенными для них операциями сложения и умножения на число. Для линейных пространств в качестве метрики часто выбирают норму, что делает его линейным нормированным пространством. Линейным нормированным пространством называется пара $(X, ||\cdot||)$, где $X$ -- линейное пространство, а $||\cdot||: X \longrightarrow \mathbb{R}$ -- норма, удовлятворяющая следующим аксиомам:

1. $||x|| \geq 0$,
2. $||x|| = 0 \Longleftrightarrow x = \bm{0}$,
3. $||\lambda x|| = |\lambda| \cdot ||x||$,
4. $||x_1 + x_2|| \leq ||x_1|| + ||x_2||$,

где $x \in X$. Легко показать, что линейное нормированное пространство является метрическим пространством с метрикой $\rho(x_1, x_2) = ||x_1 - x_2||$. Банаховым пространством называется полное линейное нормированное пространство. Самым очевидным примером банахова пространства является $(\mathbb{R}, |\cdot|)$, где нормой является модуль числа. Так как банаховы пространства и соответствующие им нормы играют важную роль в численных методах, мы рассмотрим внимательно несколько самых важных банаховых пространств. Конечномерным нормированным пространством $l_p^{(n)}$ называется пара $(X, ||\cdot||_p)$, где $X$ – множество векторов $\bm{x} = (x_1, ..., x_n)$ в $n$-мерном линейном пространстве и норма определена функцией

\begin{equation} ||\bm{x}||_p = \left(\frac{1}{n}\sum_{i=1}^n |x_i|^p \right)^{1/p}. \end{equation}

Для простоты будем предполагать, что $n$-мерное линейное пространство определенно над полем вещественных или комплексных чисел. Эвклидово пространство $\mathbb{R}^n$ является частным случаем конечномерного нормированного пространства: $l_2^{(n)}$. Если индекс нормы опущен, то предполагается, что используется классическая эвклидова норма, т.е. $||\cdot|| = ||\cdot||_2$.

Важно отметить, что сходимость в одной из норм $||\cdot||_p$ гарантирует сходимость во всех остальных нормах этого типа. Если последовательность векторов $\{\bm{x_m}\}_{n=1}^{\infty}$ не сходится, но сходится последовательность то говорят о сходимости по направлению.

Логичным расширением случая конечномерных нормированных пространств является аналогичное пространство с бесконечной размерностью. Заметим, что в таком случае множества, образуемые векторами, должны оставаться счетными. Бесконечномерным нормированным пространством $l_p$ называется пара $(X, ||\cdot||_p)$, где $X$ -- множество векторов $\bm{x} = (x_1, x_2, ...)$, каждый из которых является в свою очередь счетным множеством, и норма определена функцией

\begin{equation} \label{eq:inf_norm_space} ||\bm{x}||_p = \lim_{n to \infty}\left(\frac{1}{n}\sum_{i=1}^n |x_i|^p \right)^{1/p}. \end{equation}

Следующий шаг состоит в переходе от счетных векторов к несчетным. Такой шаг не является излишним теоретизированием, так как функции, с которыми мы чаще всего имеет дело, являются бесконечномерными векторами с несчетным числом элементов. Это утверждение звучит контринтуитивно из-за терминологической путаницы, которая вносится понятием мерности пространства. Предположим, что у нас есть функция $d(x, y, z)$, вычисляющая расстояние от начала координат до точки в трехмерном пространстве: С точки зрения теории множеств, эта функция задана на трехмерном пространстве и отображает его на положительную вещественную ось: $f: \mathbb{R}^3 \longrightarrow \mathbb{R}_{+}$. Однако если мы рассматриваем саму функцию как элемент некоторого множества, то выясняется, что она является бесконечномерным вектором с несчетным числом элементов. Действительно, рассмотрим следующий $(n+1)$-мерный вектор: где $f_j = \sin \frac{2\pi j}{n}$ и $j = 0, \dots, n$. При $n \to \infty$ вектор $\bm{f}$ стремится к функции $\sin x$ на отрезке $[0; 2\pi]$, что демонстрируется на рисунке ниже. Учитывая, что отрезок $[0; 2\pi]$ является несчетным множеством, соответствующий бесконечномерный вектор, восстанавливающий функцию $\sin x$, так же будет несчетным множеством.

В таком случае естественное обобщение предела суммы в норме \eqref{eq:inf_norm_space} до интеграла приводит к лебеговым пространствам $L_p$. Лебеговым пространством $L_p$ называется пара $(F, ||\cdot||_p)$, где $F$ -- множество функций $x(t)$, $p$-я степень которых интегрируема на отрезке $[a, b]$, и норма определена функционалом

\begin{equation} ||x(t)||_p = \left(\int_a^b |x(t)|^p dt \right)^{1/p}. \end{equation}

Пространство $L_2$ называют гильбертовым, а его норму $||\cdot||_2$ среднеквадратичной. Норму $||x(t)||_\infty = \max_{t \in [a; b]} |x(t)|$ называют равномерной или чебышевской. Несложно доказать, что для норм верны следующие соотношения: что означает, что из равномерной сходимости (т.е. сходимости в норме $||x(t)||_\infty$) следует сходимость в среднем (т.е. сходимость в норме $||x(t)||_p, p < \infty$), однако обратное не является в общем случае верным. На рисунке демонстрируется сходимость последовательности функций \newline $\{f_1(x), f_2(x), f_3(x), \dots\}$ к некоторой функции $f(x)$. На левом рисунке последовательность функций сходится равномерно и, следовательно, в среднем, в то время как на правом рисунке последовательность функций сходится только в среднем: "пик", формируемый $f_i(x)$ при $i \to \infty$, не окажет влияния на значение интеграла в норме $||\cdot||_2$, но при этом приведет к ненулевому значению $||f_i(x)||_\infty = \max_{x \in [a; b]} |f_i(x)|$. Вопрос равномерной и средней сходимости имеет прикладное значение, так как если некоторый численный метод сходится в среднем, но не сходится равномерно, это означает, что численное решение может включать в себя паразитное решение, например, в форме паразитных осцилляций. Сходимость в среднем гарантирует, что паразитные осцилляции имеют меру нуль (т.е. определены в одной точке и нигде больше) при бесконечно малом размере сетки. Однако, так как размер сетки всегда конечен, они появятся в численном решении задачи и могут сделать его неудовлетворительным. В курсе мы встретимся с примерами паразитных осцилляций.

Так как класс непрерывных функций играет важную роль в уравнениях математической физики, рассмотрим их отдельно. Пространством непрерывных функций $C[a, b]$ называется пара $(F, ||\cdot||_C)$, где $F$ – множество функций, непрерывных на отрезке $[a, b]$, а $||x||_C = \max_{t \in [a, b]}|x(t)|$ – равномерная норма. Пространство непрерывных функций также обозначают $C^0[a, b]$. В свою очередь пространство $C^p[a, b]$ называют пространство функций, $p$-я производная которых непрерывна. Гладкой называют функцию, имеющую непрерывные производные (порядок последней непрерывной производной часто опускают, подразумевая, что функция достаточно гладкая для решения данной задачи). Бесконечно гладкие функции называют аналитическими – такие функции можно представить в виде бесконечной суммы ряда Тейлора. Функция $x(t)$ называется равномерно-непрерывной на заданном отрезке, если $\forall \epsilon > 0\ \exists \delta = \delta(\epsilon): |x(t_1) - x(t_2)| \leq \epsilon, |t_1 - t_2| \leq \delta$. Функция $x(t)$ называется липшиц-непрерывной, если $\epsilon \leq K \delta$, что эквивалентно $|x(t_1) - x(t_2)| \leq K|t_1 - t_2|$, где $K$ называется константой Липшица. Если функция имеет ограниченную производную, т.е. $|x'(t)| \leq K$, то она липщиц-непрерывна, причем точная верхняя грань модуля производной равна константе Липшица $K$. Липшиц-непрерывные функции играют важную роль в теории дифференциальных уравнений, так как если функция $f$ в уравнении $x'(t) = f(t, x)$ является липшиц-непрерывной по переменной $x$ (т.е. имеет ограниченную частную производную по $x$), то из этого следует существование и единственность решения уравнения.