Численные методы

Численные (вычислительные) методы — методы решения математических задач в численном виде^[1].

Представление как исходных данных в задаче, так и её решения — в виде числа или набора чисел.

Многие численные методы являются частью библиотек математических программ^[2]. В системе подготовки инженеров технических специальностей являются важной составляющей.

Основами для вычислительных методов являются:

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

Вавилонские математики (II тысячелетие до н. э.) разработали для извлечения квадратного корня особый численный метод^[3].

Вавилонская глиняная табличка YBC 7289 из вавилонской коллекции Йельского университета была создана между 1800 и 1600 годами до н. э. и демонстрирует √2 и √2/2 соответственно в шестидесятиричной системе счисления: 1;24,51,10 и 0;42,25,35 на квадрате, пересечённом двумя диагоналями^[4]. (1;24,51,10) по основанию 60 соответствует 1,41421296, что является правильным значением с точностью до 5 десятичных знаков: ${\displaystyle 1+24/60+51/60^{2}+10/60^{3}=1{,}41421296.}$

Сложные проценты вычислялись с помощью метода линейной интерполяции.

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

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

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

В древней Греции Евдокс (метод исчерпывания Евдокса), Архимед (метод верхних и нижних интегральных сумм ), и Герон Александрийский (правило извлечения кубического корня)применяли численные методы к решению некоторых задач. ^[5]

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

В Новое время численные методы продолжали развиваться: Николай Кузанский нашел формулу приближенного спрямления окружност, Рене Декарт решая задачу по преломлению света нашел овалы Декарта, Бернулли предложил метод изоклин. Исаак Ньютон предложил метод решения алгебраических уравнений (метод Ньютона).

Все задачи вычислительной математики решаются в следующей последовательности^[6]:

1. Исходная математическая задача заменяется другой задачей — вычислительным алгоритмом. Основными требованиями к вычислительному алгоритму являются: высокая точность, устойчивость и экономичность. При переходе к дискретной модели появляется погрешность аппроксимации, а при реализации вычислений — погрешность округления, поэтому для реальных вычислительных алгоритмов проводится анализ погрешностей и устойчивости вычислительного алгоритма^[2]. В современной науке для решения задач прикладной математики формулируется математическая модель в терминах интегральных и дифференциальных уравнений функций непрерывного аргумента. Переход от континуальной к дискретной математической модели осуществляется заменой функций непрерывного аргумента функциями дискретного аргумента. В получившихся конечно-разностных уравнениях интеграл и производная представлены конечной суммой и разностным отношением, соответственно^[2]. Получившаяся модель представляет собой систему алгебраических уравнений, для решения которой с определённой точностью составляется вычислительный алгоритм, который реализуется на вычислительных машинах^[2][7]. При решении больших систем необходимо вычислять собственные значения и вектора матриц, сводить нелинейные системы уравнений к линейным. Для некоторых задач (нейронная физика, физика плазмы, экономика) модель строится непосредственно на статистической выборке или на крупных объектах. Кроме того, строятся нерегулярные системы, для которых численные методы сочетаются с теорией графов. Отдельный класс представляют некорректно поставленные задачи^[2].
2. Вычислительный алгоритм содержит параметр ${\displaystyle N}$, которого нет в исходной задаче;
3. Выбором этого параметра ${\displaystyle N}$ можно добиться любой близости решения второй задачи к решению первой. Для многих важных классов задач разработаны разнообразные численные методы решения. По способу дискретизации численные методы делятся на проекционные и конечно-разностные, по способу решения — на прямые и итерационные. В методах конечных разностей ставится задача определить значения функции на дискретном множестве точек, в то время как в проекционных методах функция представлена линейной комбинацией элементов. При этом дискретная функция также может рассматриваться как линейная комбинация полиномов. Прямые методы решения обладают слабой устойчивостью, в то время как итерационные методы более устойчивы и обеспечивают быструю сходимость^[2].
4. Неточная реализация алгоритма, вызванная округлениями при вычислениях, не меняет существенно его свойств. Необходимо помнить, что вычислительная машина выполняет только четыре основных арифметических операции^[8]. Точность решения при этом должна быть несколько выше ожидаемой точности физического эксперимента^[9]. При определении критериев и условий роста погрешности долгое время не принималась во внимание погрешность округления. Необходимость гарантированных оценок точности реальных вычислений привела к возникновению интервального анализа. Оптимальным алгоритмом считается алгоритм с минимальной погрешностью или с минимальным числом операций при заданной погрешности. При этом разрабатывается теория параллельных вычислительных алгоритмов^[2].

Символически задача поиска неизвестной величины записывается в виде ${\displaystyle y=A(x)}$. Для отыскания ${\displaystyle y}$ в вычислительной математике используют одну или несколько замен пространств, в которых определены величины ${\displaystyle x}$, ${\displaystyle y}$, или функции ${\displaystyle A}$, чтобы сделать вычисления более удобными. Получившаяся новая задача ${\displaystyle {\bar {y}}={\bar {A}}({\bar {x}})}$ должна иметь решение, близкое к решению исходной задачи. Например, при вычислении интеграла ${\displaystyle \int _{a}^{b}f(x)dx}$, непрерывную функцию на отрезке ${\displaystyle [a,b]}$ можно всегда заменить полиномом ${\displaystyle P(x)}$, для которого интеграл легко определяется; или же заменить интеграл конечной суммой ${\displaystyle \sum _{i=1}^{n}{f(x_{i})\Delta x_{i}}}$ и решать получившуюся задачу. Для того чтобы осуществить подобную замену, необходимо отыскать конечное множество элементов, хорошо аппроксимирующих основное пространство. Последнее условие накладывает ограничения на метрическое пространство. Основным ограничением является наличие ${\displaystyle \epsilon }$-сети, из которого вытекает компактность пространства в себе и сепарабельность. Вместе с тем, это ограничение не является обязательным. Современные методы функционального анализа позволяют выбрать метрические пространства, наиболее подходящие условиям задачи^[10].

При использовании численных методов возникает несколько видов погрешностей. При приближении одного числа другим возникает погрешность округления, погрешность связанная с неточными начальными данными называется неустранимой, кроме того, в связи с заменой исходной задачи на приближённую существует погрешность метода. Полная погрешность при этом складывается из погрешности метода и погрешности вычислений, иными словами, вместо уравнения ${\displaystyle y=A(x)}$ решается уравнения ${\displaystyle {\bar {\bar {y}}}={\bar {\bar {A}}}({\bar {\bar {x}}})}$, точность решения которого определяется по формуле^[11]

${\displaystyle y-{\bar {\bar {y}}}=(y-{\bar {y}})+({\bar {y}}-{\bar {\bar {y}}})}$

Для определения величины погрешности пользуются понятиями абсолютной и относительной погрешности, а также предельной абсолютной и относительной погрешности, при этом теория погрешностей определяет изменение величин погрешностей при различных арифметических действиях^[12]. Наряду с методами точной оценки погрешностей, в результате которых определяются предельные величины погрешностей, используют статистические методы, позволяющие определить возможность достижения отдельных погрешностей^[13], а также учитывают математические характеристики случайных ошибок, связанных с отклонением от заданных условий опыта, когда по нескольким результатам измерения физической величины определяется её приближённое значение^[14].

Для получения значения функции ${\displaystyle f(x)}$, заданной таблицей значений, на промежуточных значениях аргумента строят приближённую функцию ${\displaystyle \varphi (x)}$, которая в заданных точках ${\displaystyle x_{0},\dots ,x_{m}}$, которые называются узлами интерполирования, принимает значения ${\displaystyle f(x_{0}),\dots ,f(x_{m})}$, а в остальных точках принадлежат области определения функции. Чаще всего приближённая функция строится в виде алгебраического многочлена, включающего первые ${\displaystyle n+1}$ элементов линейно независимой системы. На практике в качестве элементов линейно независимой системы используют последовательность степеней ${\displaystyle x}$: ${\displaystyle 1,x,x^{2},\dots }$, тригонометрических функций: ${\displaystyle 1,\sin x,\cos x,\sin 2x,\dots }$, показательных функций: ${\displaystyle 1,e^{\alpha _{1}x},e^{\alpha _{2}x},\dots }$^[15].

Для построения интерполирующей функции в таком случае необходимо решить систему из ${\displaystyle m+1}$ уравнений с ${\displaystyle n+1}$ неизвестными. На получившуюся матрицу системы накладываются определённые условия: ранг матрицы должен быть равен ${\displaystyle m+1}$, а ${\displaystyle n\geq m}$ — чтобы гарантировать условие линейной независимости, ${\displaystyle n=m}$ — чтобы решение задачи было однозначным, определитель матрицы ${\displaystyle \Delta \neq 0}$ — чтобы существовало решение и притом единственное^[16]. Построение интерполяционного многочлена Лагранжа ${\displaystyle L_{n}(x)}$ является базовым методом решения подобного рода задач, очень ресурсоёмким и трудно расширяемым^[17].

Следующим шагом является введение понятия разделённой разности ${\displaystyle k}$-го порядка на базе отношений разности значения функции в соседних узлах к расстоянию между узлами, которая в силу своего определения обладает рядом полезных свойств, в частности разделённые разности порядка ${\displaystyle k}$ от многочлена степени ${\displaystyle n}$ имеют степень ${\displaystyle n-k}$, то есть разности порядка ${\displaystyle n}$ постоянны, а разности более высокого порядка равны ${\displaystyle 0}$^[18]. Разделённые разности позволяют переписать интерполяционный многочлен Лагранжа в виде, более удобном для вычислений. Новая формула носит название интерполяционного многочлена Ньютона^[19], в случае равных промежутков формула значительно упрощается^[20]. С использованием разделённых разностей строятся интерполяционные формулы Гаусса, Стирлинга, Бесселя, Эверетта^[21]. В общем случае разделённые разности сначала убывают с повышением порядка, а затем начинают снова расти, иными словами, нет смысла использовать разности высоких порядков в вычислениях^[22]. При этом возникает вопрос сходимости интерполяционного процесса, для решения которого привлекаются различные методы математического анализа^[23].

При решении практических задач необходимо многократно вычислять значения заданной функции, что в общем случае является ресурсоёмкой операцией. Возникает необходимость нахождения функции наилучшего равномерного приближения^[24]. Для приближения функции в линейном нормированном пространстве образуют подпространство размерности ${\displaystyle n+1}$ всевозможных линейных комбинаций, для которых опеределена норма и существует её точная нижняя грань. Элемент, в котором эта грань достигается называют элементом наилучшего приближения, или проекцией^[25]. Можно доказать что в подпространстве всегда существует элемент наилучшего приближения^[26], а при условии строгой нормированности пространства такой элемент является единственным^[27]. В пространстве непрерывных функций с нормой

${\displaystyle \lVert f\rVert =\sup _{x\in {[a,b)}}|f(x)|}$

также существует элемент наилучшего приближения^[28], но условием его единственности является наличие не более ${\displaystyle n}$ различных нулей обобщённого многочлена на отрезке (Многочлены Чебышёва)^[29].

Теория функций применима к системе степенных функций, так как она является системой Чебышёва на любом отрезке^[30]. Согласно теореме Вейерштрасса, при увеличении размерности подпространства (${\displaystyle n\to \infty }$) разность между проекцией и заданной функцией стремится к нулю^[31]. Порядок этого приближения зависит от структурных особенностей функции, его можно определить с помощью многочленов Бернштейна^[32]. Система тригонометрических функций также обладает свойствами системы Чебышёва на отрезке ${\displaystyle [0,2\pi )}$, для неё также разность между проекцией и заданной функцией стремится к нулю^[33].

Несмотря на показанное существование многочлена наилучшего приближения, способов его точного построения не существует. Вместо этого используют несколько способов приближённого построения многочленов наилучшего равномерного приближения^[34].

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

${\displaystyle \delta ={\sqrt {\int _{a}^{b}p(x)[f(x)-\phi (x)]^{2}dx}},}$

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

Уравнение вида ${\displaystyle y=A(x)}$, определённое на функциональном пространстве, может содержать операторы дифференцирования и интегрирования, для которых невозможно найти точное решение. Методы численного дифференцирования и интегрирования основаны на интерполяции^[36].

Производную основной функции считают приближённо равной производной интерполирующей функции, при этом производная остаточного члена интерполяционной формулы может быть велика, особенно для производных высших порядков^[37]. Формулы численного дифференцирования во многом основаны на непосредственном дифференцировании интерполяционных формул Ньютона^[38], Гаусса, Стирлинга и Бесселя^[39], построенных на распределённых разностях, но есть и безразностные формулы. В частности, когда для численного дифференциала используется непосредственно формула Лагранжа для равных промежутков^[40], метод неопределённых коэффициентов и другие^[41].

В случае интегрирования, само определение интеграла говорит о возможности его замены интегральной суммой, но этот приём обладает медленной сходимостью и мало пригоден. Интеграл от основной функции считают приближённо равным интегралу от интерполирующей функции и в дальнейшем используют интерполяционные формулы с кратными узлами^[42]. Использование в качестве подынтегрального выражения интерполяционного многочлена Лагранжа для равных промежутков приводит к формулам Ньютона — Котеса^[43] и её частным случаям, формуле трапеций, когда кривая подынтегрального выражения заменяется хордой и интеграл равен площади трапеции, и формуле Симпсона, когда кривая подынтегрального выражения заменяется параболой, проходящей через три точки^[44]. Отказавшись от требования равных промежутков с помощью интерполяционного многочлена Лагранжа можно получить более точные формулы численного интегрирования, в частности формулы Гаусса^[45], формулы Эрмита^[46], формулы Маркова^[47], формулы Чебышёва^[48]. Квадратурные процессы, построенные на интерполяционных формулах Гаусса, всегда сходятся, в то время как формулы Ньютона — Котеса этим свойствам в общем случае не обладают^[49].

Существуют и другие способы численного интегрирования, основным из которых является использование формул Эйлера, в которых замена переменных и последующее интегрирование по частям приводят к формуле численного интегрирования трапецией и поправочного члена, к которому повторно применяется замена переменных и интегрирование по частям. В общем случае формула Эйлера использует в качестве коэффициентов числа и многочлены Бернулли^[50]. Вопрос применения того или иного метода численного интегрирования зависит от таких факторов, как вычислительные средства, требуемая точность, способ задания подынтегральной функции. Для ручных вычислений рекомендуется использовать формулы, содержащие разности, в то время как при автоматических вычислениях — безразностные формулы, в особенности формулы Гаусса^[51].

Для приближённого вычисления кратных интегралов повторно применяют формулы численного интегрирования однократных интегралов, при этом в зависимости от особенностей функции для разных интегралов можно использовать разные формулы. При использовании данного метода необходимо вычислять подынтегральную функцию в большом числе точек, поэтому целесообразно использовать формулы Гаусса и Чебышёва, которые являются более точными^[52]. Другим способом является замена подынтегральной функции интерполяционным многочленом от двух или несколько переменных^[53]. Люстерник и Диткин предложили использовать формулы Маклорена для приближённого вычисления кратного интеграла^[54]. Вместе с тем, при увеличении кратности интеграла резко растёт число точек, для которых необходимо знать значения подынтегральной функции, чтобы пользоваться методами, основанными на интерполяции. Для вычисления кратных интегралов чаще пользуются вероятностными методами Монте-Карло, при этом необходимость получения равновозможных последовательностей создаёт дополнительные погрешности, которые трудно оценить^[55].

Существует две группы методов решения систем линейных алгебраических уравнений: точные методы позволяют с помощью конечного числа операций получить точные значения неизвестных и включают преобразование системы к простому виду и решение упрощённой системы; методы последовательных приближений на основе начальных приближений позволяют получить «улучшенные» приближённые значения, для которых следует последовательно повторить операцию «улучшения»; методы Монте-Карло позволяют на основании математического ожидания случайных величин получить решение системы^[56].

Известный из школьного курса алгебры метод исключения позволяет свести матрицу системы к диагональному или треугольному виду^[57]. Схема исключения Гаусса с выбором главного элемента, который необходим чтобы уменьшить вычислительную погрешность, включает прямой ход (собственно процесс исключения) и обратный ход (решение системы с треугольной матрицей)^[58]. Её компактный вариант используется для определения обратной матрицы, что может быть полезно если в системе линейных уравнений меняется только правая часть^[59] и для вычисления определителей^[60]. Схема Жордана позволяет облегчить обратный ход^[61], а в схеме без обратного хода, которая основана на преобразовании клеточной матрицы ${\displaystyle {\begin{pmatrix}A&b\\-I&0\end{pmatrix}}}$, последний и не требуется^[62]. Условие симметричности матрицы позволяет сделать ряд упрощений и воспользоваться методом квадратного корня, в котором матрица системы представляется как произведение нижней треугольной матрицы на транспонированную по отношению к ней матрицу, в котором элементы треугольных матриц определяются по формулам через произведения элементы первоначальной матрицы (при отсутствии условия положительно определённой матрицы некоторые формулы могут содержать мнимые элементы), а система затем решается в два этапа через решение вспомогательных систем построенных на треугольных матрицах^[63]. Существуют также метод ортогонализации, основанный на свойствах скалярного произведения^[64], метод сопряжённых градиентов, при котором строится вспомогательная функция, которая образует семейство эллипсоидов с общим центром и для которой необходимо найти вектор, при котором она принимает минимальное значение^[65]. Для матриц высокого порядка применяют метод разбиения на клетки, когда задачу сводят к решению задач для матриц низших порядков^[66].

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

${\displaystyle {\overline {x}}_{k+1}=F_{k}({\overline {x}}_{0},\dots ,{\overline {x}}_{k}),}$

где ${\displaystyle F_{k}}$ — функция, которая зависит от матрицы системы, правой части, номера приближения и предыдущих приближений ${\displaystyle {\overline {x}}_{0},\dots ,{\overline {x}}_{k}}$, где ${\displaystyle {\overline {x}}_{0}}$ — начальный вектор. При этом считается, что метод имеет первый порядок, если функция зависит только от последнего из предыдущих приближений. В этом случае формула ${\displaystyle {\overline {x}}_{k+1}=B_{k}{\overline {x}}_{k}+{\overline {c}}_{k}}$ может быть записана в виде ${\displaystyle D_{k}{\overline {x}}_{k+1}+E_{k}{\overline {x}}_{k}={\overline {b}}}$, где ${\displaystyle D_{k}+E_{k}=A}$. Для удобства вычислений желательно использовать диагональную или треугольную матрицу ${\displaystyle D_{k}}$, которую будет удобно обратить. В зависимости от выбора этой матрицы методы называют полношаговыми и одношаговыми, соответственно^[67]. К линейным полношаговым методам относят простую итерацию^[68], метод Ричардсона^[69]; к линейным одношаговым методам — метод Зейделя^[70], релаксационный метод^[71]; к нелинейным методам — метод скорейшего спуска^[72].

Решения алгебраического уравнения ${\displaystyle f(z)=0}$, где в левой части находится функция действительного или комплексного аргумента, лежит в комплексной плоскости^[73]. Для его определения в первую очередь необходимо заключить каждый корень в достаточно малую область, то есть отделить его, для чего часто используют графические методы^[74]. Для действительных корней используют также обобщённое правило Декарта, теорему Штурма^[75], метод Фурье^[76]. Широкое применение нашёл метод квадратного корня, или метод Лобачевского^[77]. В его основной формулировке он применим к действительным корням^[78], далеко отстоящим друг от друга, но существуют обобщения как на комплексные^[79], так и на действительные равные или близкие корни^[80].

Итерационные методы решения алгебраических уравнений делятся на стационарные, когда функции ставится в соответствие другая функция с теми же корнями, не зависящая от номера итерации^[81], и нестационарные, когда функция может зависеть от номера итерации. К простейшим стационарным итерационным методам относят метод секущих (или метод линейного интерполирования) и метод касательных (или метод Ньютона), которые являются методами первого и второго порядка, соответственно. Комбинация этих методов, при которой последовательные приближения лежат по разные стороны от корня, позволяет достичь более быстрой сходимости^[82]. Метод Чебышева, основанный на разложении обратной функции по формуле Тейлора, позволяет построить методы более высоких порядков, обладающие очень быстрой сходимостью^[83]. Существуют также метод, основанный на теореме Кёнига^[84], и метод Эйткена^[85]. Для доказательства сходимости итерационных методов используется принцип сжатых отображений^[86].

• Метод конечных элементов

• Амосов А. А., Дубинский Ю. А., Копченова Н. В. Вычислительные методы для инженеров. — 1994.
• История математики, в трёх томах / Под редакцией А. П. Юшкевича. — М.: Наука, 1970—1972.
• Березин И. С., Жидков Н. П. Методы вычислений. — М.: Наука, 1962. — Т. 1.
• Березин И. С., Жидков Н. П. Методы вычислений. — М.: Наука, 1959. — Т. 2.
• Калиткин Н. Н. Численные методы. — М.: Наука, 1978.
• Численные методы : теория и практика : учебное пособие для бакалавров, для студентов высших учебных заведений, обучающихся по направлению подготовки «Математика. Прикладная математика» / У. Г. Пирумов, Гидаспов В.Ю., Иванов И.Э., Ревизников Д. Л., Стрельцов В.Ю., Формалев В.Ф.; Московский авиационный ин-т-нац. исслед. ун-т. — 5-е изд., перераб. и доп. — Москва : Юрайт, 2012. — 421 с. : ил., табл.; 22 см. — (Бакалавр. Базовый курс).; ISBN 978-5-9916-1867-0
• Рыжиков Ю. Вычислительные методы. — СПб.: BHV, 2007. — 400 с. — ISBN 978-5-9775-0137-8

• Научный журнал «Вычислительные методы и программирование. Новые вычислительные технологии»
• Материалы по численным методам
• Численные методы
• Вычислительные методы в прикладной математике, Международный журнал, ISSN 1609-4840