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

Методы многомерной оптимизации: градиентный спуск и метод Ньютона

Поиск наилучшего решения в мире множества вариантов

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

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

Градиентный спуск: движение по наискорейшему склону

Самым известным и концептуально простым алгоритмом является метод градиентного спуска (Gradient Descent), предложенный еще Огюстеном Луи Коши. Градиент — это вектор, указывающий направление наискорейшего возрастания функции. Следовательно, антиградиент (вектор, взятый со знаком минус) указывает направление наискорейшего убывания.

Алгоритм стартует из случайной начальной точки. На каждом шаге мы вычисляем вектор градиента в текущей точке и делаем небольшой шаг в противоположном направлении. Размер этого шага определяется параметром, который называется «скоростью обучения» (learning rate). Главная проблема классического градиентного спуска — сильная зависимость от этого параметра. Если шаг слишком велик, алгоритм начнет перепрыгивать через минимум и расходиться; если слишком мал — сходимость растянется на долгие часы. Кроме того, в так называемых «овражных» функциях градиентный спуск начинает совершать зигзагообразные колебания поперек оврага, крайне медленно продвигаясь к его дну.

Метод Ньютона и учет кривизны пространства

Чтобы ускорить сходимость и избежать зигзагов в оврагах, необходимо использовать информацию не только о наклоне функции (первой производной), но и о ее кривизне. Эту информацию предоставляет матрица вторых производных, называемая матрицей Гессе (или гессианом). Использование гессиана лежит в основе многомерного метода Ньютона.

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

Подробнее

Численная линейная алгебра: вычисление собственных значений и векторов

Скрытые резонансы: проблема спектрального анализа

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

Физический и инженерный смысл этой задачи колоссален. Собственные значения описывают частоты свободных (резонансных) колебаний механических систем: от гитарной струны и подвески автомобиля до мостов и небоскребов. Совпадение внешней нагрузки с этими частотами вызывает разрушительный резонанс. В квантовой механике собственные значения оператора Гамильтона определяют дискретные уровни энергии атомов (отсюда происходит термин «спектр»). В теории графов и анализе социальных сетей анализ спектра матриц смежности позволяет находить скрытые кластеры и самых влиятельных участников.

Степенной метод для поиска доминирующего значения

Аналитически собственные значения находятся путем решения так называемого характеристического уравнения (поиска корней полинома, степень которого равна размерности матрицы). Однако, как мы знаем из теоремы Абеля-Руффини, для матриц размером 5x5 и больше точных формул корней не существует. Поэтому в реальности собственные значения всегда ищут только итерационными численными методами.

Самым простым подходом является степенной метод (итерации фон Неймана). Он позволяет найти наибольшее по модулю (доминирующее) собственное значение. Идея поразительно проста: мы берем случайный начальный вектор и начинаем многократно умножать его на нашу матрицу. На каждом шаге вектор будет поворачиваться в сторону собственного вектора, соответствующего максимальному собственному значению. После достаточного количества итераций (и периодической нормализации вектора, чтобы числа не ушли в бесконечность) процесс сходится к искомому главному собственному вектору. Именно этот простой принцип лежит в основе первоначального алгоритма PageRank, созданного Ларри Пейджем и Сергеем Брином для ранжирования веб-страниц в поисковой системе Google, где матрицей выступает граф ссылок всего интернета.

Алгоритм QR-разложения для нахождения всего спектра

Степенной метод хорош, если нам нужно только одно максимальное значение. Но что делать, если инженеру нужен весь спектр частот (все собственные значения)? Для этого в 1961 году Джоном Фрэнсисом и Верой Кублановской был независимо разработан QR-алгоритм, который считается одним из 10 самых важных численных алгоритмов 20-го века.

Алгоритм основан на факторизации (разложении) матрицы. На каждом шаге исходная матрица A представляется в виде произведения ортогональной матрицы Q и верхнетреугольной матрицы R (A = QR). Затем эти матрицы умножаются в обратном порядке: строится новая матрица A1 = RQ. Удивительный математический факт состоит в том, что новая матрица A1 имеет точно такие же собственные значения, что и исходная матрица A (они подобны). Если повторять этот процесс итеративно (находя QR-разложение и перемножая в обратном порядке), то матрица постепенно стремится к диагональному (или блочно-диагональному) виду, где все ее собственные значения просто стоят на главной диагонали. С применением предварительного приведения к форме Хессенберга и стратегий сдвигов, QR-алгоритм работает феноменально быстро и стабильно, являясь сегодня ядром библиотек линейной алгебры, таких как LAPACK, встроенных в Python, MATLAB и R.

Подробнее

Метод Монте-Карло: стохастический подход к вычислениям и интегрирование

Случайность на службе строгой математики

Численные методы, которые мы обсуждали до сих пор, являются детерминированными: они следуют строгим формулам, и при одинаковых входных данных всегда дают абсолютно идентичный результат с точностью до бита. Однако в середине 20-го века, во время работы над Манхэттенским проектом, физики Станислав Улам, Джон фон Нейман и Николас Метрополис предложили радикально иной подход, названный кодовым словом «Монте-Карло» в честь знаменитого казино. Этот метод использует генерацию случайных чисел для решения сугубо детерминированных математических задач.

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

Проклятие размерности в задачах интегрирования

Где же метод Монте-Карло проявляет себя лучше всего? Его главная ниша — вычисление кратных интегралов высокой размерности. В финансовой математике (ценообразование опционов), статистической физике или квантовой механике часто возникают интегралы, зависящие от десятков, сотен и даже тысяч переменных.

Если мы попытаемся использовать классические детерминированные методы (например, метод Симпсона) для вычисления D-мерного интеграла, мы столкнемся с «проклятием размерности». Если на каждую ось мы поместим по 10 узлов сетки, то для 3-мерного объема потребуется 10^3 = 1000 узлов. Но для 100-мерного интеграла потребуется 10^100 вычислений функции — число, превышающее количество атомов во Вселенной. Ни один суперкомпьютер с этим не справится. Ошибка классических методов зависит от размерности пространства. Метод Монте-Карло же обладает уникальным свойством: скорость его сходимости (убывания ошибки) равна O(1/sqrt(N)), где N — количество испытаний, и эта скорость абсолютно не зависит от размерности пространства! Это делает его единственным рабочим инструментом для задач высокой размерности.

Алгоритмы уменьшения дисперсии и цепи Маркова (MCMC)

Главный недостаток базового метода Монте-Карло — его медленная сходимость. Из-за корня в формуле сходимости, чтобы уменьшить ошибку в 10 раз, нужно увеличить количество испытаний в 100 раз. Чтобы бороться с этим, математики разработали изощренные техники снижения дисперсии (variance reduction techniques).

Одним из мощнейших инструментов является метод выборки по значимости (importance sampling). Вместо того чтобы разбрасывать случайные точки равномерно по всему объему, алгоритм концентрирует выборку в тех областях, где подынтегральная функция вносит наибольший вклад в итоговый результат. Другим прорывом стали алгоритмы Монте-Карло по схеме марковских цепей (Markov Chain Monte Carlo, MCMC), в частности алгоритм Метрополиса-Гастингса. В этом случае случайные блуждания не являются независимыми: следующая точка выбирается на основе положения предыдущей, формируя интеллектуальный поиск в многомерном пространстве параметров. Алгоритмы MCMC произвели революцию в байесовском машинном обучении, искусственном интеллекте и расшифровке генома, доказав, что контролируемая случайность — один из самых мощных инструментов познания.

Подробнее

Решение жестких систем обыкновенных дифференциальных уравнений (ОДУ)

Феномен жесткости: когда классические методы дают сбой

При численном моделировании многих химических, биологических и электронных систем исследователи сталкиваются с загадочным и крайне неприятным вычислительным явлением, получившим название «жесткость» (stiffness). Представьте себе химический реактор, в котором протекают две параллельные реакции: одна завершается за миллисекунды (например, горение или взрыв), а вторая длится часами (медленное окисление или диффузия). Эти процессы описываются системой обыкновенных дифференциальных уравнений (ОДУ).

Математически система называется жесткой, если ее решения содержат компоненты с кардинально различающимися скоростями затухания. Если мы попытаемся решить такую систему классическим явным методом (например, популярным методом Эйлера или методом Рунге-Кутты 4-го порядка), мы столкнемся с катастрофой. Для того чтобы обеспечить вычислительную устойчивость и не дать быстрым компонентам «взорвать» решение, шаг интегрирования должен быть меньше характерного времени самого быстрого процесса. В результате, чтобы промоделировать медленный процесс длительностью в часы, компьютеру придется сделать триллионы крошечных миллисекундных шагов, что займет месяцы машинного времени и приведет к огромному накоплению ошибок округления.

Неявные методы и А-устойчивость

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

На первый взгляд это кажется усложнением: теперь на каждом шаге нужно запускать метод Ньютона и решать системы линейных уравнений. Однако неявные методы обладают потрясающим свойством — абсолютной устойчивостью (или А-устойчивостью). Для таких методов размер шага интегрирования ограничивается только требованиями к точности описания физического процесса, а не условиями математической устойчивости. Медленный процесс можно интегрировать огромными шагами, полностью игнорируя сверхбыстрые затухающие колебания, не боясь при этом «взрыва» алгоритма. Неявный метод Эйлера является простейшим представителем этого класса.

Формулы дифференцирования назад (BDF) и методы Гира

Простейшие неявные методы обладают низким порядком точности (первым или вторым). Для точного моделирования жестких систем требуются методы более высоких порядков. Американский математик Чарльз Гир в начале 1970-х годов разработал семейство жестко-устойчивых алгоритмов, известных как формулы дифференцирования назад (Backward Differentiation Formulas, BDF).

Вместо того чтобы аппроксимировать подынтегральную функцию правой части, методы BDF строят интерполяционный полином по нескольким последним вычисленным точкам самого решения, а затем дифференцируют этот полином и приравнивают его к значению правой части в новой, еще не вычисленной точке. Это позволяет создавать высокоточные многошаговые неявные схемы (вплоть до 6-го порядка). Современные программные библиотеки (например, LSODE или решатели ode15s в MATLAB) используют адаптивные алгоритмы на основе методов BDF. Они автоматически распознают жесткость системы, на лету меняя как размер шага, так и порядок метода, обеспечивая робастное и эффективное решение сложнейших многомасштабных задач химической кинетики и схемотехники.

Подробнее

Метод конечных элементов (МКЭ): фундаментальные концепции и приложения

От сложных геометрических форм к простым элементам

Метод конечных разностей, при всех его достоинствах, сталкивается с колоссальными трудностями при попытке смоделировать физические процессы в областях со сложной криволинейной геометрией. Построить прямоугольную сетку внутри фюзеляжа самолета или человеческого черепа так, чтобы она точно повторяла границы, практически невозможно. Для решения этой фундаментальной проблемы в середине 20-го века инженерами-прочнистами был разработан метод конечных элементов (МКЭ), который сегодня является абсолютным стандартом в промышленном проектировании систем (CAD/CAE).

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

Базисные функции и принцип локальности

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

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

Сборка глобальной матрицы жесткости и решение СЛАУ

Процесс вычислений в МКЭ строго формализован. Сначала вычисляются так называемые локальные матрицы жесткости для каждого отдельного элемента. Они описывают, как элемент сопротивляется деформации или пропускает тепло. Затем эти локальные матрицы «сшиваются» (агрегируются) в огромную глобальную матрицу жесткости всей конструкции. Этот процесс топологической сборки легко алгоритмизируется и выполняется компьютером автоматически на основе таблицы связности элементов.

После учета граничных условий (закрепления деталей, приложения внешних сил) задача сводится к решению гигантской системы линейных алгебраических уравнений (СЛАУ). Современные коммерческие пакеты МКЭ (такие как ANSYS, Abaqus, COMSOL) решают СЛАУ размерностью в десятки и сотни миллионов неизвестных, используя специализированные итерационные решатели (например, метод сопряженных градиентов с предобуславливанием) на суперкомпьютерных кластерах. МКЭ позволил инженерам отказаться от дорогостоящих натурных испытаний и перейти к виртуальному прототипированию, кардинально изменив облик современной промышленности.

Подробнее

Спектральные методы решения дифференциальных уравнений: от Фурье до Чебышева

Альтернатива локальным методам сеток

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

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

Ряды Фурье для периодических задач

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

Для перехода между физическим и спектральным пространствами используется алгоритм быстрого преобразования Фурье (БПФ). Открытие БПФ в середине 20-го века Кули и Тьюки произвело настоящую революцию в вычислительной математике, снизив алгоритмическую сложность преобразования с квадратичной до логарифмической. Это позволило спектральным методам успешно конкурировать с разностными схемами на мощных суперкомпьютерах.

Полиномы Чебышева и Лежандра: решение непериодических задач

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

Узлы интерполяции в этих методах располагаются неравномерно: они сгущаются к краям расчетной области и разрежаются в центре. Такое распределение (узлы Гаусса-Чебышева или Гаусса-Лобатто) позволяет избежать феномена Рунге и обеспечивает так называемую спектральную (или экспоненциальную) точность. Это означает, что при увеличении количества узлов N погрешность убывает быстрее, чем любая конечная степень 1/N. Если решение является гладкой аналитической функцией, спектральные методы могут достичь машинной точности при использовании всего нескольких десятков узлов, тогда как разностным методам для этого потребовались бы миллионы точек. Это делает спектральные методы незаменимым инструментом в высокоточной вычислительной гидродинамике и квантовой химии.

Подробнее

Численное решение уравнений в частных производных: метод сеток

Повелители сложных систем: Уравнения математической физики

Если обыкновенные дифференциальные уравнения (ОДУ) описывают системы, зависящие только от одного параметра (например, от времени), то уравнения в частных производных (УЧП) описывают многомерные процессы, которые разворачиваются как во времени, так и в пространстве. Вся современная физика базируется на УЧП. Уравнение теплопроводности Фурье описывает нагрев и охлаждение деталей двигателя. Волновое уравнение описывает распространение звука, радиоволн и света. Уравнения Навье-Стокса управляют течением жидкостей и газов, позволяя проектировать аэродинамику самолетов и прогнозировать погоду. Уравнение Шредингера лежит в основе квантовой механики.

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

Дискретизация: превращение континуума в матрицу

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

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

Явные и неявные схемы: проблема устойчивости

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

Неявные схемы работают иначе: они связывают неизвестные значения на новом временном слое в единую систему уравнений. Чтобы сделать шаг по времени, компьютеру приходится решать гигантскую СЛАУ. Это требует гораздо больше вычислительных ресурсов на один шаг, но зато неявные схемы обладают абсолютной устойчивостью. Инженер может задать огромный шаг по времени (ограниченный лишь требуемой точностью отображения физики процесса), и вычисления не развалятся. Балансировка между этими методами, создание адаптивных сеток и распараллеливание на кластерах — суть современной вычислительной гидродинамики (CFD).

Подробнее

Метод наименьших квадратов (МНК) в обработке экспериментальных данных

Борьба с шумом: почему интерполяция здесь бессильна?

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

С этой задачей блестяще справляется Метод Наименьших Квадратов (МНК), впервые опубликованный Лежандром и строго обоснованный Карлом Фридрихом Гауссом. Этот алгоритм стал фундаментом математической статистики, эконометрики и современного машинного обучения (в частности, алгоритмов линейной и полиномиальной регрессии).

Математическая суть МНК

Идея метода заложена в его названии. Мы заранее выбираем вид функции (модель), которая, по нашему мнению, должна описывать процесс. Это может быть прямая линия (y = ax + b), парабола, экспонента или синусоида. Эта функция содержит неизвестные параметры (коэффициенты a и b), которые нам предстоит найти.

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

Сведение к системе нормальных уравнений

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

Эта система называется системой нормальных уравнений. Матрица этой системы обладает рядом полезных свойств (она симметрична и положительно определена), что делает ее решение устойчивым. Вычислив коэффициенты, мы получаем ту самую «наилучшую» кривую, которая проходит сквозь облако экспериментальных точек. МНК обладает высочайшей статистической эффективностью: согласно теореме Гаусса-Маркова, при определенных условиях (гомоскедастичность и отсутствие автокорреляции ошибок), оценки параметров, полученные методом наименьших квадратов, являются лучшими, линейными и несмещенными (BLUE).

Подробнее

Интерполяция и аппроксимация функций: многочлены Лагранжа и сплайны

Задача восстановления функции по точкам

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

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

Глобальная интерполяция: многочлен Лагранжа

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

Самым известным способом записи такого полинома является интерполяционный многочлен Лагранжа. Его формула строится элегантно и интуитивно понятно: это линейная комбинация базисных полиномов. Каждый базисный полином устроен так, что в одном конкретном узле сетки он равен единице, а во всех остальных узлах обращается в ноль. Многочлен Лагранжа прекрасен в теории, но имеет серьезный изъян на практике при большом количестве точек. Этот изъян известен как «феномен Рунге». Если мы попытаемся интерполировать функцию по большому количеству узлов, мы получим полином очень высокой степени (например, 20-й степени). Между узлами на краях интервала этот полином начнет демонстрировать дикие, высокоамплитудные осцилляции (раскачку), которые не имеют ничего общего с реальным поведением исходной физической величины.

Кубические сплайны: гибкость и гладкость

Чтобы избежать феномена Рунге и не использовать полиномы гигантских степеней, вычислительная математика перешла к кусочно-полиномиальной интерполяции, жемчужиной которой являются сплайны. Само слово «сплайн» пришло из инженерного черчения — так называли гибкую металлическую линейку, которую прижимали грузиками к узловым точкам чертежа, чтобы провести плавную кривую.

В математике сплайн — это функция, которая на каждом отдельном интервале между двумя соседними точками задается своим собственным полиномом (обычно третьей степени — кубический сплайн). Однако, в отличие от простого набора кривых, эти полиномы сшиваются в узлах особым образом. Накладываются жесткие условия не только на непрерывность самой функции (кривая не должна разрываться), но и на непрерывность ее первой производной (отсутствие изломов, гладкость) и второй производной (непрерывность кривизны). В результате получается идеально гладкая, визуально эстетичная и математически устойчивая кривая. Сплайны не подвержены осцилляциям Рунге и сегодня лежат в основе всей компьютерной графики (векторные шрифты, CAD-системы, 3D-моделирование).

Подробнее

Численное решение обыкновенных дифференциальных уравнений (ОДУ)

Задача Коши и потребность в численных методах

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

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

Метод Эйлера: простота и геометрический смысл

Исторически первым и самым простым численным методом решения ОДУ является метод Эйлера, предложенный великим Леонардом Эйлером еще в 18 веке. Этот метод основан на разложении функции в ряд Тейлора и удержании только первых двух членов разложения.

Геометрический смысл метода предельно нагляден. Пусть мы находимся в начальной точке графика. Дифференциальное уравнение дает нам значение производной (то есть тангенс угла наклона касательной) в этой точке. Метод Эйлера предлагает сделать небольшой шаг вдоль этой касательной. Мы смещаемся на расстояние шага интегрирования h по оси абсцисс, вычисляем новую ординату и принимаем полученную точку за новое начальное условие. Процесс повторяется. Главный недостаток метода Эйлера — его низкая точность. Погрешность метода на одном шаге пропорциональна квадрату шага, а глобальная погрешность на всем интервале — первой степени шага (метод первого порядка точности). Чтобы получить приемлемый результат, приходится брать микроскопически малый шаг, что ведет к накоплению катастрофических ошибок округления.

Семейство методов Рунге-Кутты

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

Идея методов Рунге-Кутты заключается в том, чтобы на каждом шаге интегрирования вычислить значения функции правой части (наклоны касательных) не только в начальной точке отрезка, но и в нескольких промежуточных точках (пробных точках) внутри текущего интервала. В классическом методе RK4 вычисляются четыре таких коэффициента. Итоговое приращение функции вычисляется как взвешенная сумма этих четырех коэффициентов. Метод RK4 обладает глобальной погрешностью порядка O(h^4). Это значит, что уменьшение шага всего в два раза приводит к уменьшению ошибки в 16 раз. Благодаря идеальному балансу между вычислительной сложностью и высокой точностью, метод RK4 стал стандартом де-факто («рабочей лошадкой») во многих программных пакетах инженерного анализа.

Подробнее