Решение задач вычислительной математики с использованием NumPy

Основы NumPy для численных вычислений

NumPy — базовая библиотека экосистемы Python для работы с массивами и выполнения быстрых численных вычислений. В вычислительной математике мы почти всегда работаем с векторами, матрицами, сетками, выборками измерений и табличными данными. NumPy даёт единый инструмент для всех этих случаев: многомерный массив ndarray и набор операций над ним.

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

> Официальная документация: NumPy документация

Что такое ndarray и почему он важен

numpy.ndarray — это:

Многомерный массив фиксированного типа данных (все элементы одного dtype).

Эффективное хранение (как правило, в непрерывном блоке памяти).

Быстрые операции благодаря реализации на C и векторизации.

Вычислительная математика постоянно использует объекты, которые естественно представляются массивами:

вектор значений функции на сетке

матрица коэффициентов системы линейных уравнений

набор наблюдений (строки) и признаков (столбцы)

Справка по объекту массива: numpy.ndarray

Установка и импорт

Чаще всего NumPy импортируют как np.

Создание массивов

Из Python-списков

Важно: при создании NumPy выберет общий тип данных (dtype), который сможет вместить все элементы.

Специальные конструкторы

Равномерные сетки: arange и linspace

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

Практический совет:

np.linspace(a, b, n) удобен, когда важно число узлов.

np.arange(a, b, h) удобен, когда важен шаг.

Форма массива: shape, размерность: ndim, количество элементов: size

Ключевые характеристики массива:

shape — кортеж размеров по осям.

ndim — число осей (измерений).

size — общее число элементов.

!Визуальная интуиция формы массива и осей

Типы данных: dtype

NumPy хранит элементы в одном типе данных, например int64, float64.

Почему это важно для численных задач:

float64 обычно безопаснее для вычислений (меньше ошибок округления, чем у float32).

целочисленные типы могут приводить к нежелательному целочисленному делению или переполнению, если не контролировать dtype.

Индексация и срезы

Одномерные массивы

Двумерные массивы

Важный момент: срезы — это представления (views)

Срез часто не копирует данные, а даёт вид на исходный массив.

Если нужна независимая копия:

Векторизация: вычисления без явных циклов

Вычислительная математика выигрывает от того, что операции применяются сразу ко всем элементам.

Пример: посчитать на сетке.

Почему это лучше циклов for:

обычно быстрее

код короче и ближе к математической записи

меньше шансов ошибиться с индексами

Broadcasting: как NumPy “растягивает” размеры

Broadcasting — правило, по которому NumPy выполняет поэлементные операции над массивами разных форм, если их размеры совместимы.

Классический пример: прибавить вектор к каждой строке матрицы.

Интуиция: v как бы “повторяется” для каждой строки.

Официальное руководство: NumPy Broadcasting

!Наглядное объяснение, как работает broadcasting

Практический приём: “сделать столбец” через добавление оси.

Универсальные функции (ufunc) и агрегирования

Поэлементные функции

NumPy предоставляет быстрые поэлементные функции:

Агрегирования

Агрегирования сводят массив к меньшей размерности: сумма, среднее, максимум.

axis=0 означает “вдоль строк”, то есть результат по столбцам.

axis=1 означает “вдоль столбцов”, то есть результат по строкам.

Матричные операции и линейная алгебра

В вычислительной математике очень часто встречаются выражения вида .

— матрица коэффициентов (двумерный массив)

— вектор неизвестных (одномерный массив)

— правая часть (одномерный массив)

Умножение: поэлементное и матричное — это разные операции

Оператор @ соответствует матричному произведению.

Скалярное произведение и норма

Скалярное произведение двух векторов и часто записывают как . В NumPy:

Евклидова норма (длина) вектора :

Пояснение символов:

— длина (норма) вектора в евклидовом смысле

— число элементов вектора

— -й элемент вектора

— суммирование по индексам от 1 до

— квадратный корень

В NumPy это обычно делают так:

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

Для квадратной системы :

Справочник по линейной алгебре: numpy.linalg

Практический совет:

не используйте np.linalg.inv(A) @ b для решения системы, если можно np.linalg.solve(A, b) — решение через solve обычно устойчивее и быстрее.

Генерация данных: numpy.random

Для моделирования, Монте-Карло и тестирования алгоритмов нужны случайные данные.

default_rng предпочтительнее старого интерфейса, потому что даёт явный генератор и более предсказуемое управление случайностью.

Чтение и запись массивов

Текстовые файлы

Двоичный формат NumPy (быстро и без потери точности)

Ошибки, на которые часто наталкиваются

Путаница между * и @:

- * — поэлементное умножение - @ — матричное умножение

Неожиданные изменения исходного массива из-за view (срез без .copy()).

Несовместимые формы при broadcasting.

Смешивание целых и вещественных типов без контроля dtype.

Минимальный чек-лист перед численным экспериментом

Определите, что вы моделируете: вектор, матрицу, сетку.

Проверьте shape, dtype и оси для агрегирования.

Используйте векторизацию и broadcasting вместо циклов.

Для линейной алгебры применяйте np.linalg.solve, np.linalg.norm, оператор @.

Для воспроизводимости случайных экспериментов используйте default_rng(seed).

Что дальше в курсе

Дальнейшие темы вычислительной математики (аппроксимация, численное дифференцирование и интегрирование, решение систем и оптимизация) будут опираться на навыки из этой статьи: уверенное владение формой массивов, векторизацией, broadcasting и базовой линейной алгеброй в NumPy.

Погрешности, устойчивость и численная точность

В прошлой статье вы освоили базовый инструмент вычислительной математики в Python — массив numpy.ndarray: формы (shape), типы (dtype), векторизацию, broadcasting и базовую линейную алгебру (@, np.linalg.solve). Следующий шаг — понять, почему даже при правильной формуле и правильном коде численный ответ может отличаться от “математически точного”, и как сделать вычисления надёжными.

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

погрешность — насколько вычисленное значение отличается от истинного;

устойчивость алгоритма — насколько алгоритм усиливает или подавляет неизбежные ошибки округления;

численная точность — ограничения представления чисел (в первую очередь чисел с плавающей точкой) и итоговая точность результата.

Откуда берутся ошибки в численных вычислениях

Обычно источники ошибок такие:

Погрешность модели: исходная математическая модель упрощает реальность.

Погрешность метода (аппроксимации): вы заменяете непрерывную задачу дискретной (например, сеткой linspace), интеграл — суммой, производную — разностью.

Ошибки округления (floating-point rounding): компьютер хранит числа приближённо, и операции выполняются с округлением.

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

Абсолютная и относительная погрешность

Пусть истинное значение — , а вычисленное (приближённое) — .

Абсолютная погрешность:

Где:

— истинное значение (то, к чему стремимся);

— вычисленное значение;

— модуль (в одномерном случае);

— величина ошибки в тех же единицах, что и .

Относительная погрешность:

Где:

— ошибка “в долях” от масштаба самого значения;

знаменатель показывает, что одна и та же абсолютная ошибка может быть существенной для малых величин и несущественной для больших.

Практический смысл:

абсолютная погрешность важна, когда есть естественная шкала (например, “ошибка не больше ”);

относительная погрешность важна, когда важен процент/доля (“ошибка не больше от значения”).

Числа с плавающей точкой: почему float64 не “вещественные числа”

Большинство вычислений в NumPy по умолчанию выполняются в float64. Это формат IEEE 754 двойной точности: число хранится приближённо, потому что имеет ограниченное число бит.

Часть десятичных дробей (например, ) не представима точно в двоичной системе.

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

Справочный материал:

IEEE 754

numpy.finfo

Машинный эпсилон и пределы типа

Для типа float64 полезно знать две характеристики:

машинный эпсилон eps — характерный масштаб относительной ошибки округления около 1;

максимум/минимум — границы диапазона.

Пример в NumPy:

Практический вывод:

если вы ожидаете точность “на уровне ” в обычных расчётах float64, это нереалистично;

если вы работаете с очень большими/очень малыми числами, возможны переполнение (inf) и потеря значащих разрядов.

!Наглядно показывает, почему дроби вроде 0.1 не хранятся точно

Накопление ошибок: почему сумма зависит от порядка

Ошибки округления особенно заметны, когда:

вы складываете числа очень разного масштаба;

выполняете длинные суммы (например, интегрирование суммированием);

многократно применяете шаговый метод (итерации, рекурсии).

Пример: “большое + маленькое”

Если число очень большое, добавление очень маленького может вообще не изменить результат (маленькое “не помещается” в мантиссу при данном масштабе).

Это не ошибка NumPy, а ограничение представления float64.

Суммирование и dtype

Если вы суммируете массив float32, результат может быть заметно хуже, чем при суммировании в float64.

Практический вывод:

контролируйте dtype на входе и при агрегированиях;

если точность важна, часто разумно аккумулировать в float64 даже при хранении данных в float32.

См. документацию: numpy.sum

Катастрофическая потеря значимости (cancellation)

Одна из самых неприятных ситуаций — вычитание близких чисел. В результате остаётся малая разность, но значащие цифры “съедаются” округлением.

Пример: для малых

При малых значения очень близки к 1, и выражение вычисляется как разность почти равных чисел.

Часто более устойчивой формой является:

Пояснение символов:

— аргумент;

и — тригонометрические функции;

означает квадрат значения синуса;

равенство показывает два математически эквивалентных способа вычисления, но численно они могут вести себя по-разному.

Численный пример:

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

Специальные функции для “малых величин”

NumPy предоставляет функции, специально сделанные для повышения точности в критичных зонах:

np.expm1(y) вычисляет точнее, чем np.exp(y) - 1 при малых ;

np.log1p(y) вычисляет точнее, чем np.log(1 + y) при малых ;

np.hypot(a, b) устойчиво вычисляет без лишнего переполнения/потери точности.

Ссылки:

numpy.expm1

numpy.log1p

numpy.hypot

Обусловленность задачи и устойчивость алгоритма

Важно разделять:

обусловленность задачи: насколько сама математическая задача чувствительна к малым изменениям входных данных;

устойчивость алгоритма: насколько конкретный алгоритм усиливает неизбежные ошибки округления.

Интуиция на примере линейной системы

Мы часто решаем (в прошлой статье вы делали это через np.linalg.solve).

Если матрица плохо обусловлена, то даже маленькая ошибка в или может сильно изменить решение . Это свойство задачи, а не “плохой код”.

Для квадратной матрицы часто используют число обусловленности:

Пояснение символов:

— матрица;

— обратная матрица (если существует);

— норма (способ измерить “размер” матрицы);

— число обусловленности: чем оно больше, тем потенциально сильнее усиливаются относительные ошибки.

В NumPy число обусловленности оценивают так:

Документация: numpy.linalg.cond

Устойчивость: почему solve лучше, чем inv

Даже если задача решаема, выбор реализации влияет на точность.

Плохая практика для решения :

сначала найти обратную: np.linalg.inv(A);

затем умножить на b.

Хорошая практика:

использовать np.linalg.solve(A, b).

Причины:

solve использует численно устойчивые разложения (в типичных случаях LU с частичным выбором главного элемента);

вычисление обратной матрицы обычно дороже и часто усиливает ошибки.

Ссылки:

numpy.linalg.solve

numpy.linalg.inv

!Различие между обусловленностью задачи и устойчивостью алгоритма

Практические приёмы точных вычислений в NumPy

Ниже набор приёмов, которые особенно полезны в вычислительной математике.

Следите за dtype

- используйте float64 как безопасный стандарт для расчётов; - явно приводите типы при чтении данных и перед вычислениями.

Масштабируйте данные

- избегайте ситуаций, когда в одной формуле встречаются числа порядка и без необходимости; - нормализация часто улучшает численное поведение.

Переписывайте формулы, чтобы избегать вычитания близких чисел

- используйте тождества (как с ); - применяйте log1p, expm1, hypot.

В линейной алгебре выбирайте устойчивые операции

- для используйте np.linalg.solve, а не np.linalg.inv(A) @ b; - проверяйте np.linalg.cond(A) для диагностики “опасных” матриц.

Сравнивайте числа с допусками, а не через ==

- используйте np.isclose и np.allclose.

Документация:

numpy.allclose

numpy.isclose

Что означает допуск в allclose

По умолчанию близость проверяется примерно по правилу:

Пояснение символов:

и — сравниваемые числа (или массивы поэлементно);

— абсолютная разница;

atol — абсолютный допуск (важен около нуля);

rtol — относительный допуск (важен при больших масштабах);

— масштаб второго аргумента, относительно которого измеряется относительная часть.

Практический совет:

если ожидаемые значения близки к нулю, подбирайте atol осознанно;

если величины большие, критичен rtol.

Частые ловушки в вычислительных экспериментах

Незаметный переход к float32

- например, данные пришли как float32, и дальше весь расчёт идёт с меньшей точностью.

Смешивание целых и вещественных типов

- особенно важно следить за типом при чтении/создании массивов.

Неверные ожидания от “точности печати”

- то, что NumPy печатает мало знаков, не означает, что значение хранится именно так; - для диагностики можно управлять выводом:

Документация: numpy.set_printoptions

Итог

Числа с плавающей точкой — приближённые, и ошибки округления неизбежны.

Ошибки могут накапливаться и усиливаться, особенно при вычитании близких чисел и при плохой обусловленности.

Устойчивость алгоритма и обусловленность задачи — разные вещи, но обе критичны.

В NumPy есть практические инструменты для контроля точности: np.finfo, устойчивые функции (log1p, expm1, hypot), диагностика (np.linalg.cond) и правильные методы (np.linalg.solve).

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

Решение систем линейных уравнений и матричные разложения

Системы линейных уравнений вида — один из центральных объектов вычислительной математики. Практически любая задача аппроксимации, численного решения дифференциальных уравнений, оценки параметров или оптимизации в итоге сводится к линейной алгебре.

В первой статье курса вы научились работать с массивами NumPy и вызывать базовые операции линейной алгебры. Во второй статье — увидели, что точность результата зависит не только от формулы, но и от устойчивости и обусловленности. Здесь мы объединим это в практический инструментарий:

как решать правильно и диагностировать качество решения

как решать переопределённые системы (наименьшие квадраты)

зачем нужны матричные разложения (QR, Холецкого, SVD, собственные значения) и когда какое применять

Основная документация для этого модуля: numpy.linalg

Постановка задачи

Запись означает:

— матрица коэффициентов размера

— вектор неизвестных длины

— вектор правой части длины

Если , система называется квадратной. Если дополнительно существует единственное решение, NumPy чаще всего решает её через np.linalg.solve.

Пояснение символов:

— таблица чисел (матрица), где каждая строка — одно уравнение

— неизвестные, которые мы ищем

— заданные значения в правой части

знак равенства означает, что после умножения матрицы на вектор мы должны получить

!Наглядная геометрия объектов в уравнении Ax=b

Как решать квадратную систему в NumPy

Базовый способ: np.linalg.solve

Для квадратной матрицы и вектора используйте:

Ссылка: numpy.linalg.solve

Ключевые свойства:

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

внутри используются разложения (обычно LU с выбором главного элемента), но вам не нужно реализовывать их вручную

Почему не надо решать через обратную матрицу

Математически можно написать , но вычислительно это почти всегда хуже.

Плохая практика:

Хорошая практика:

Причины:

вычисление inv обычно дороже, чем решение системы

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

Ссылка: numpy.linalg.inv

Проверка качества решения: невязка и норма

Даже если solve вернул ответ, полезно проверить насколько хорошо он удовлетворяет исходной системе.

Определим вектор невязки:

Пояснение символов:

— найденное численное решение

— левая часть системы при подстановке найденного

— исходная правая часть

— насколько «не сошлось» по каждому уравнению

Чтобы получить одно число-оценку, берут норму невязки, например евклидову:

Пояснение символов:

— евклидова норма: «длина вектора»

— размер ошибки в единицах

В NumPy:

Ссылка: numpy.linalg.norm

Практическая интерпретация:

маленькая означает, что вы хорошо решили именно эту задачу

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

Обусловленность и чувствительность решения

Если матрица плохо обусловлена, малые ошибки в или могут привести к большим изменениям в . Стандартная диагностика — число обусловленности.

В NumPy:

Ссылка: numpy.linalg.cond

Интуитивно:

если cond(A) близко к 1, задача обычно «здоровая»

если cond(A) очень велико, решение может быть чувствительным к округлению и шуму в данных

Практический совет, связанный с темой точности:

проверяйте dtype и старайтесь решать в float64

думайте о масштабах: иногда нормализация столбцов/строк существенно улучшает поведение

Переопределённые системы и метод наименьших квадратов

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

Идея наименьших квадратов: подобрать , чтобы сделать минимальной.

— вектор ошибок по всем уравнениям

минимизация нормы означает «в среднем» приблизить все уравнения как можно лучше

В NumPy это делается через:

Ссылка: numpy.linalg.lstsq

Что возвращает lstsq:

x — найденное решение (в смысле наименьших квадратов)

residuals — сумма квадратов невязок (может быть пустой в некоторых случаях)

rank — оценка ранга матрицы

s — сингулярные числа (важны для диагностики)

!Геометрический смысл МНК как проекции

Практические замечания:

lstsq полезен для регрессии, калибровок, восстановления параметров по шумным данным

ранг (rank) и сингулярные числа (s) помогают понять, не является ли задача вырожденной или почти вырожденной

Зачем нужны матричные разложения

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

Ниже — самые важные разложения, доступные в numpy.linalg, и их типичные применения.

QR-разложение

QR-разложение представляет матрицу как:

Пояснение символов:

— исходная матрица

— матрица с ортонормированными столбцами (геометрически: поворот/отражение без изменения длины)

— верхнетреугольная матрица (удобна для решения систем обратной подстановкой)

В NumPy:

Ссылка: numpy.linalg.qr

Где QR полезен:

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

построение ортонормированных базисов

Разложение Холецкого

Если матрица симметрична и положительно определена, то есть:

(симметрия)

для любого ненулевого вектора верно (положительная определённость)

то существует разложение Холецкого:

Пояснение символов:

— нижнетреугольная матрица

— транспонированная

такое разложение обычно быстрее и численно приятнее общего случая

В NumPy:

Ссылка: numpy.linalg.cholesky

Где Холецкий полезен:

решения систем с матрицами ковариаций

методы оптимизации второго порядка, где возникает симметричная положительно определённая матрица

многие задачи МНК через нормальные уравнения (но сами нормальные уравнения могут ухудшать обусловленность, поэтому часто предпочитают QR или SVD)

SVD: сингулярное разложение

SVD — одно из самых мощных и устойчивых разложений. Оно представляется так:

Пояснение символов:

— ортонормированная матрица (поворот/отражение в пространстве выходов)

— диагональная матрица сингулярных чисел (неотрицательные числа, отражающие «силу» направлений)

— транспонированная ортонормированная матрица (поворот/отражение в пространстве входов)

В NumPy:

Ссылка: numpy.linalg.svd

Где SVD полезен:

диагностика ранга и почти линейной зависимости столбцов

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

построение псевдообратной матрицы (когда это действительно нужно)

Практическое правило:

QR обычно быстрее

SVD обычно устойчивее и информативнее для диагностики

Собственные значения и собственные векторы

Для квадратной матрицы встречается задача найти такие и , что:

Пояснение символов:

— квадратная матрица

— собственный вектор (ненулевой)

— собственное значение

равенство означает, что действие на меняет только масштаб (на ), но не направление

В NumPy:

Ссылка: numpy.linalg.eig

Где это применяется:

устойчивость динамических систем

спектральные методы

анализ квадратичных форм

Для симметричных матриц обычно лучше использовать специализированную функцию:

Ссылка: numpy.linalg.eigh

eigh использует свойства симметрии и обычно работает точнее и быстрее.

Как выбрать разложение или метод

| Задача | Рекомендуемый инструмент | Требования к матрице | Типичный смысл | |---|---|---|---| | Квадратная система | np.linalg.solve | квадратная, невырожденная | Стандартное решение | | Переопределённая система (МНК) | np.linalg.lstsq | Любая форма | Приближение с минимальной невязкой | | Быстрое решение для SPD матриц | np.linalg.cholesky + подстановка | Симметричная положительно определённая | Эффективность и устойчивость | | Ортогонализация, МНК через разложение | np.linalg.qr | Любая форма | Стабильное преобразование | | Диагностика ранга, плохо обусловленные задачи | np.linalg.svd | Любая форма | Самая информативная диагностика | | Спектральный анализ | np.linalg.eig / np.linalg.eigh | Квадратная (для eigh — симметричная) | Направления и масштабы |

Практический чек-лист перед запуском вычислений

Приведите данные к float64, если точность важна.

Решаете квадратную систему — начинайте с np.linalg.solve.

Всегда проверяйте невязку np.linalg.norm(A @ x - b).

Оцените np.linalg.cond(A), если решение кажется «подозрительным».

Если уравнений больше, чем неизвестных, используйте np.linalg.lstsq.

Если матрица симметричная и положительно определённая, подумайте о np.linalg.cholesky.

Для диагностики вырождения и ранга используйте np.linalg.svd.

Итог

Решение в NumPy следует начинать с np.linalg.solve, а не с вычисления обратной матрицы.

Качество решения удобно контролировать через невязку и её норму.

Плохая обусловленность может сделать задачу чувствительной к округлению и шуму — это ключевая идея из темы численной точности.

Матричные разложения (QR, Холецкого, SVD, спектральные методы) — это не «теория ради теории», а практические инструменты устойчивых и быстрых вычислений.

Собственные значения и задачи линейной алгебры

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

после решения систем и контроля невязки вы умеете проверять качество численного результата;

после обсуждения погрешностей и обусловленности вы понимаете, что спектральные задачи тоже могут быть чувствительными к округлению;

после обзора матричных разложений вы уже видели eig и eigh, а теперь разберём когда и зачем их применять.

Основные ссылки на документацию NumPy:

numpy.linalg.eig

numpy.linalg.eigh

numpy.linalg.svd

numpy.linalg.norm

Что такое собственные значения и собственные векторы

Пусть дана квадратная матрица размера . Пара называется собственным значением и собственным вектором матрицы , если выполняется:

Пояснение символов:

— квадратная матрица (линейное преобразование)

— ненулевой вектор (направление в пространстве)

— число (масштаб)

равенство означает: при действии вектор не меняет направление, а только масштабируется в раз

Геометрическая интерпретация:

вектор — это особое направление, которое матрица не «поворачивает» (в идеальном случае), а растягивает/сжимает

значение говорит, насколько сильно происходит растяжение (если ) или сжатие (если ) вдоль этого направления

!Геометрический смысл уравнения Av = λv

Спектр матрицы и почему он важен в вычислениях

Набор всех собственных значений матрицы называют спектром. В вычислительной математике спектр важен потому что:

он определяет поведение итерационных процессов (сходимость, скорость затухания/роста)

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

он связан с квадратичными формами и минимумами/максимумами функций

он лежит в основе методов понижения размерности (например, PCA через собственные значения ковариационной матрицы)

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

Как искать собственные значения в NumPy

Общий случай: np.linalg.eig

Для произвольной квадратной матрицы применяют:

Интерпретация результата:

w — массив собственных значений

V — матрица, где V[:, i] — собственный вектор , соответствующий w[i]

Практическое замечание:

собственные векторы могут быть масштабированы произвольно: если — собственный вектор, то тоже собственный вектор для любого ненулевого

Симметричные матрицы: np.linalg.eigh

Если матрица симметрична, то есть , используйте:

Почему это важно:

для симметричных матриц собственные значения гарантированно вещественные

алгоритмы для eigh учитывают симметрию, обычно дают более точный результат и работают быстрее

Пояснение символов для условия симметрии:

— матрица

— транспонированная матрица

равенство означает, что элементы зеркальны относительно главной диагонали

Проверка качества найденных собственных пар

Как и в задаче решения , полезно проверять результат численно.

Для каждой пары определим невязку (ошибку уравнения собственного вектора):

Пояснение символов:

— результат действия матрицы на вектор

— тот же вектор, умноженный на число

— вектор, который должен быть близок к нулю, если пара вычислена хорошо

На практике проверяют норму невязки:

Пояснение символов:

— евклидова норма (длина вектора)

малое значение означает, что равенство выполняется хорошо

Пример проверки:

Практический смысл:

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

Диагонализация и спектральное разложение

Диагонализация (общая идея)

Если у матрицы есть полный набор линейно независимых собственных векторов, её можно представить как:

Пояснение символов:

— исходная матрица

— матрица, составленная из собственных векторов по столбцам

— диагональная матрица (на диагонали стоят собственные значения)

— обратная матрица к

Почему это полезно:

многие операции с матрицами упрощаются (например, степень матрицы)

поведение системы часто можно понять через значения на диагонали

Важное ограничение:

не всякая матрица диагонализируема (особенно при кратных собственных значениях)

даже если диагонализация существует, вычислительно она может быть чувствительной к ошибкам

Симметричный случай: ортогональное разложение

Для симметричной матрицы есть особенно удобная форма:

Пояснение символов:

— матрица собственных векторов, причём столбцы ортонормированы

— транспонированная

ортонормированность означает: столбцы взаимно перпендикулярны и имеют длину 1

Практический плюс:

ортонормированная матрица не усиливает ошибки так, как может усиливать плохая в общей формуле

Собственные значения в типовых задачах вычислительной математики

Устойчивость дискретной динамики

Рассмотрим линейную модель:

Пояснение символов:

— состояние системы на шаге (вектор)

— матрица перехода

— состояние на следующем шаге

Ключевая идея:

если по модулю все собственные значения меньше 1, траектории обычно затухают к нулю

если есть собственное значение с модулем больше 1, появляются направления роста

Часто используют спектральный радиус:

Пояснение символов:

— собственные значения

— модуль (важно, если значения комплексные)

максимум берётся по всем собственным значениям

В NumPy:

Документация: numpy.linalg.eigvals

!Спектральный радиус и устойчивость через расположение собственных значений

Квадратичные формы, минимум и максимум

Во многих задачах оптимизации и оценки качества встречается выражение вида:

Пояснение символов:

— вектор

— транспонированный вектор (строка)

— симметричная матрица (типично так в постановках)

результат — число

Если симметрична, то собственные значения описывают диапазон значений этой формы на единичной сфере :

минимальное значение связано с минимальным собственным значением

максимальное — с максимальным собственным значением

Это полезно для понимания:

насколько «круто» или «плоско» ведёт себя функция возле минимума

почему некоторые направления оптимизации сходятся медленно

PCA и ковариационные матрицы

Классическая постановка: у вас есть данные, вы строите ковариационную матрицу , которая симметрична и неотрицательно определена. Собственные значения показывают:

сколько дисперсии объясняет каждая главная компонента

какие направления в данных наиболее информативны

Для такой матрицы почти всегда применяют np.linalg.eigh.

Связь собственных значений с SVD

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

Пусть SVD матрицы имеет вид:

Пояснение символов:

— ортонормированная матрица левых сингулярных векторов

— диагональная матрица сингулярных чисел

— транспонированная матрица правых сингулярных векторов

Тогда для матрицы справедливо:

Пояснение символов:

— симметричная матрица (даже если несимметрична)

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

собственные значения равны , где — сингулярные числа

Практический вывод:

если цель — оценить ранк, обусловленность или найти устойчивое приближение, часто лучше использовать np.linalg.svd

если матрица симметрична и вас интересуют именно её собственные значения, используйте np.linalg.eigh

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

Связь со второй статьёй курса (про погрешности и устойчивость) здесь особенно важна.

Почему собственные значения могут быть чувствительными

Собственные значения и особенно собственные векторы могут сильно меняться при малых возмущениях матрицы, если:

у матрицы есть близкие или кратные собственные значения

матрица несимметрична и имеет «плохие» спектральные свойства

Практические советы:

предпочитайте float64 для спектральных задач

если матрица симметрична, используйте eigh, а не eig

проверяйте невязку для диагностики

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

Сортировка и интерпретация результата

eig и eigh могут возвращать значения в порядке, который не соответствует вашей задаче.

Типичные стратегии:

для PCA сортируют по убыванию собственных значений

для устойчивости дискретной системы смотрят на максимум

Пример сортировки:

Если значения комплексные (часто при eig):

сортировка по np.abs(w) может быть более осмысленной, если вас интересует рост/затухание

Минимальный практический чек-лист

Проверьте, симметрична ли матрица.

Если симметрична, используйте np.linalg.eigh.

Если несимметрична, используйте np.linalg.eig или np.linalg.eigvals (если векторы не нужны).

Контролируйте качество через невязку np.linalg.norm(A @ v - lam * v).

Для задач устойчивости используйте спектральный радиус max(abs(eigvals)).

Если задача чувствительна или плохо обусловлена, рассмотрите переход к SVD (np.linalg.svd) как более устойчивому инструменту диагностики.

Итог

Собственные значения и векторы описывают направления, которые матрица масштабирует без изменения направления.

Для симметричных матриц используйте np.linalg.eigh: это и точнее, и быстрее, и гарантирует вещественные собственные значения.

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

Спектр напрямую связан с устойчивостью итераций, квадратичными формами и методами анализа данных.

Связь с SVD помогает выбрать устойчивый инструмент, если спектральная задача численно «капризна».

Интерполяция, аппроксимация и МНК

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

Этот модуль опирается на предыдущие статьи курса:

из основ NumPy вам нужны массивы, broadcasting и векторизация

из темы численной точности — понимание, что алгоритм может быть чувствителен к масштабам и обусловленности

из темы систем линейных уравнений — навыки работы с np.linalg.solve, np.linalg.lstsq, нормами и невязкой

Главные инструменты NumPy для этого модуля:

np.interp для простой 1D линейной интерполяции

np.polyfit и np.polyval для полиномиальной аппроксимации и полиномиальной МНК

np.linalg.lstsq для общего МНК в произвольном базисе

Документация:

numpy.interp

numpy.polyfit

numpy.polyval

numpy.linalg.lstsq

numpy.linalg.norm

Интерполяция и аппроксимация: в чём разница

Пусть у нас есть данные .

Интерполяция строит функцию так, чтобы она точно проходила через заданные точки:

Пояснение символов:

— значения аргумента (узлы), где известны данные

— значения функции в узлах

— построенная интерполирующая функция

равенство означает, что мы требуем совпадения в каждом узле

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

Практическое правило выбора:

если данные точные, табличные и “должны совпасть” в узлах, чаще берут интерполяцию

если данные измеренные и содержат шум, чаще берут аппроксимацию, обычно через МНК

!Сравнение поведения интерполяции и аппроксимации на шумных данных

Линейная интерполяция в NumPy

Самый простой и часто достаточный вариант интерполяции в 1D — кусочно-линейная.

Идея: между соседними узлами и значение определяется по формуле прямой.

np.interp

np.interp делает именно это: линейную интерполяцию по отсортированным узлам.

Практические замечания:

x должен быть возрастающим (если нет — отсортируйте данные)

за пределами диапазона x функция по умолчанию возвращает крайние значения y (это не экстраполяция)

Если вам нужна экстраполяция или гладкие сплайны, обычно используют SciPy, но в рамках курса на NumPy полезно уверенно владеть линейной интерполяцией как базовой.

Полиномиальная интерполяция и матрица Вандермонда

Один из классических подходов интерполяции: найти полином степени , проходящий через точек.

Полином степени имеет вид:

Пояснение символов:

— значение полинома в точке

— неизвестные коэффициенты

— степень аргумента, задаёт базисные функции

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

Как построить Вандермондову матрицу в NumPy

После этого полином можно вычислять на сетке, используя векторизацию:

Численные проблемы полиномиальной интерполяции

Полиномиальная интерполяция на равномерной сетке при высокой степени часто ведёт к сильным колебаниям на концах интервала. Это известная проблема, связанная не только с математикой, но и с численной устойчивостью.

!Демонстрация эффекта Рунге: высокие степени на равномерных узлах дают колебания

Практические выводы:

интерполяционный полином высокой степени — не универсальный инструмент

даже если система решается, матрица Вандермонда может быть плохо обусловлена

для устойчивости часто переходят к кусочным моделям (например, кусочно-линейной интерполяции) или к аппроксимации

Аппроксимация и метод наименьших квадратов

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

Пусть есть модель, линейная по параметрам:

Пояснение символов:

— аппроксимирующая функция

— выбранные базисные функции, которые вы задаёте (например, )

— коэффициенты, которые нужно подобрать

сумма означает, что мы комбинируем базисные функции с весами

Для данных строят матрицу дизайна :

строка соответствует точке

столбец соответствует базисной функции

элемент

Тогда вектор предсказаний равен , а ошибка по всем точкам — .

Задача МНК

Метод наименьших квадратов (МНК) выбирает , минимизирующий сумму квадратов ошибок:

Пояснение символов:

— вектор неизвестных коэффициентов

— вектор ошибок (невязка) по всем точкам

— евклидова норма, то есть “длина” вектора ошибок

минимум означает: мы ищем такие коэффициенты, при которых общая ошибка минимальна

Решение МНК через np.linalg.lstsq

np.linalg.lstsq решает задачу МНК устойчивыми методами на основе матричных разложений.

Как интерпретировать возвращаемые значения:

c — найденные коэффициенты

residuals — сумма квадратов невязок (может быть пустой в некоторых случаях)

rank — оценка ранга матрицы A (важно для диагностики вырождения)

s — сингулярные числа (дают понимание обусловленности)

Связь с предыдущими статьями:

величина np.linalg.norm(A @ c - y) — это та же идея контроля невязки, что и для

если rank меньше числа столбцов, задача плохо определена: базисные функции почти линейно зависимы на ваших данных

Полиномиальная аппроксимация через np.polyfit

Для полиномиальной аппроксимации в NumPy есть удобный высокоуровневый интерфейс.

Практические замечания:

np.polyfit решает задачу МНК для полиномов

при больших степенях результат может стать численно нестабильным из-за масштабов степеней

часто помогает масштабирование аргумента, например перевод к диапазону около

Документация:

numpy.polyfit

numpy.polyval

Диагностика качества аппроксимации

Оценивать модель стоит как минимум по двум направлениям.

Невязка и её нормы

Вектор невязки:

Пояснение символов:

— предсказанные значения в точках

— исходные наблюдения

— ошибки по точкам

Часто используют:

как общий масштаб ошибки

RMSE, то есть корень из среднего квадрата ошибки

RMSE в терминах компонентов выглядит так:

Пояснение символов:

— число точек данных

— ошибка на точке

— суммирование по всем точкам

деление на — усреднение

квадрат и корень возвращают ошибку в тех же единицах, что и

Переобучение

Если вы увеличиваете сложность модели, ошибка на обучающих данных обычно уменьшается, но качество предсказания между точками или на новых данных может ухудшиться.

Признаки переобучения в численных экспериментах:

аппроксимирующая кривая начинает “ловить шум” и становится рваной

коэффициенты становятся очень большими по модулю

небольшое изменение данных сильно меняет коэффициенты

Это напрямую связано с темой обусловленности: “почти линейно зависимый” базис делает задачу МНК чувствительной.

Численная устойчивость: что важно помнить

Не решайте МНК через нормальные уравнения без необходимости

Иногда задачу МНК записывают как решение системы:

Пояснение символов:

— транспонированная матрица

— квадратная матрица размера “число параметров”

правая часть агрегирует вклад всех точек

В вычислениях эта форма может ухудшать обусловленность. Поэтому на практике лучше использовать np.linalg.lstsq, который опирается на устойчивые разложения.

Масштабирование аргумента

Если вы строите базис вроде , то при больших степени растут очень быстро, и матрица A становится плохо обусловленной.

Практические приёмы:

перенос и масштабирование:

выбор базиса, который лучше ведёт себя численно для вашего интервала

Проверяйте ранг и сингулярные числа

Если np.linalg.lstsq возвращает низкий rank или очень широкий разброс s, это сигнал, что:

параметры плохо идентифицируются по вашим данным

базисные функции почти линейно зависимы

решение будет нестабильным к шуму

Как выбирать метод в практике

| Ситуация | Цель | Рекомендуемый подход в NumPy | |---|---|---| | Табличная функция без шума, нужно значение между узлами | Совпасть в узлах и интерполировать между ними | np.interp | | Данные шумные, нужен общий тренд | Сгладить шум и получить параметры | np.linalg.lstsq или np.polyfit | | Нужна гибкость, но без высоких степеней | Избежать колебаний и плохой обусловленности | кусочная модель, простые базисы, контроль ранга |

Итог

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

Аппроксимация через МНК позволяет устойчиво подбирать параметры модели по шумным данным.

В NumPy базовые инструменты: np.interp для линейной интерполяции, np.polyfit для полиномов, np.linalg.lstsq для общего МНК.

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

Численное дифференцирование и интегрирование

Численное дифференцирование и интегрирование отвечают на две типовые ситуации вычислительной математики:

у вас есть формула , но вам нужно получить производную или интеграл на сетке;

у вас есть только табличные данные , и вы хотите оценить производную или площадь под кривой.

Эта тема напрямую опирается на предыдущие модули курса:

из основ NumPy — создание сеток np.linspace, векторизация и работа с осями axis;

из темы погрешностей — понимание компромисса между ошибкой аппроксимации и ошибкой округления;

из МНК и аппроксимации — идея локальной аппроксимации данных для устойчивых производных;

из линейной алгебры — нормы, невязки и оценка качества результата.

Важные ссылки на документацию NumPy:

numpy.diff

numpy.gradient

numpy.trapezoid

numpy.trapz

Сетка и постановка задачи

Обычно мы считаем, что функция задана на узлах сетки:

узлы

значения

Самый частый случай — равномерная сетка:

Пояснение символов:

— -й узел сетки

— первый узел

— индекс узла, целое число

— шаг сетки (расстояние между соседними узлами)

В NumPy равномерную сетку обычно строят так:

Численное дифференцирование

Идея конечных разностей

Производная — это предел отношения приращений. Численно мы заменяем предел на вычисление по соседним узлам.

#### Прямая разность (вперёд)

Пояснение символов:

— производная в узле

и — значения функции в соседних узлах

— шаг сетки

знак означает приближение, а не точное равенство

Точность этой формулы обычно первого порядка по шагу: ошибка убывает примерно как при уменьшении .

#### Обратная разность (назад)

Смысл тот же, но используются узлы слева.

#### Центральная разность

Пояснение символов:

и — значения справа и слева от узла

— расстояние между узлами и

Центральная разность обычно даёт второй порядок точности: ошибка убывает примерно как .

!Наглядное сравнение прямой, обратной и центральной разностей

Практика в NumPy: np.diff

np.diff вычисляет разности соседних элементов. Для прямой разности на равномерной сетке:

Замечания:

np.diff(y) имеет длину len(y) - 1, поэтому результат “привязан” к промежуткам.

если сетка неравномерная, деление на один h уже неверно.

Практика в NumPy: np.gradient

np.gradient удобен тем, что:

возвращает массив той же длины, что и y;

использует центральные разности внутри интервала и односторонние на границах;

умеет работать с неравномерными сетками, если передать массив x.

Пример:

Вторая производная

Часто нужна вторая производная (например, в задачах физики и оптимизации). Стандартная центральная схема:

Пояснение символов:

— вторая производная в точке

— центральный вес

— квадрат шага сетки

Реализация на равномерной сетке:

Ошибка шага: почему “меньше h” не всегда лучше

У численной производной есть два конкурирующих источника ошибки:

ошибка аппроксимации (мы заменили точную производную разностью)

ошибка округления (вычитание близких чисел и ограниченная точность float64)

Особенно опасно выражение вида при очень малом : происходит потеря значащих разрядов.

Практический подход:

начинайте с умеренного шага (например, – для “гладких” функций около масштаба 1)

проверяйте устойчивость результата при изменении шага

используйте float64

Пример эксперимента “ошибка от шага”:

Обычно вы увидите картину:

на “больших” доминирует аппроксимация (ошибка падает при уменьшении )

на слишком малых начинает доминировать округление (ошибка снова растёт)

!Типичная зависимость ошибки производной от шага сетки

Производная по данным: шум и устойчивость

Если содержат шум, то дифференцирование его усиливает. Простой признак проблемы:

интеграл “сглаживает” шум

производная “разгоняет” шум

Практические стратегии:

если данные шумные, сначала выполните аппроксимацию (например, МНК), затем дифференцируйте аппроксимирующую модель;

используйте локальную аппроксимацию полиномом (идея из МНК): на окне из нескольких точек подберите коэффициенты полинома и возьмите производную этого полинома.

Связь с предыдущим модулем:

в локальной аппроксимации вы снова решаете задачу МНК через np.linalg.lstsq и контролируете невязку.

Численное интегрирование

Интеграл как площадь и как сумма

Определённый интеграл можно понимать как площадь под графиком. Численно мы заменяем его суммой по промежуткам сетки.

Метод трапеций

Если на каждом промежутке заменить функцию прямой, то площадь — это трапеция:

Пояснение символов:

— соседние узлы

— значения функции в узлах

— ширина трапеции

дробь — средняя “высота”

Суммируя по всем промежуткам, получаем приближение интеграла на .

Интегрирование в NumPy: np.trapezoid

np.trapezoid(y, x) реализует метод трапеций и умеет работать с неравномерными x.

Полезные детали:

если x не задан, NumPy считает шаг равным 1 по индексу (это удобно для дискретных сумм, но опасно, если вы забыли передать реальную сетку);

интегрировать можно по оси axis, что важно для матриц и многомерных сеток.

np.trapz

Функция np.trapz исторически широко использовалась для метода трапеций. В современных версиях NumPy предпочтительнее np.trapezoid, но trapz может встречаться в коде и материалах.

Кумулятивный интеграл (интеграл “до точки”)

Иногда нужен не один интеграл на всём интервале, а функция накопленной площади:

Пояснение символов:

— накопленный интеграл до узла

— первый узел сетки

В NumPy можно собрать это вручную через трапеции и np.cumsum:

Теперь F[k] приближает интеграл от x[0] до x[k].

Интегрирование табличных данных

Если у вас есть измерения :

метод трапеций — базовый и часто разумный выбор;

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

Связь с модулем про интерполяцию:

если сетка редкая, иногда сначала строят более плотную сетку и линейно интерполируют y через np.interp, а затем интегрируют на плотной сетке.

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

Оценка ошибки на практике

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

сходимость при сгущении сетки: сравнить результат на n и на 2n узлах

согласованность операторов: например, если вы посчитали как интеграл от , то производная должна быть близка к (с учётом ошибок)

нормы ошибок: например, через np.linalg.norm

Когда сравниваете числа с плавающей точкой, используйте допуски:

np.isclose

np.allclose

Документация: numpy.allclose

Интегрирование и дифференцирование по осям массивов

NumPy-вычисления часто многомерные: например, у вас значения на сетке по времени и пространству.

Типовые приёмы:

производная по времени: np.gradient(F, t, axis=0)

производная по пространству: np.gradient(F, x, axis=1)

интеграл по пространству: np.trapezoid(F, x, axis=1)

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

Частые ошибки и чек-лист

Типовые проблемы:

забыли передать x в np.gradient или np.trapezoid и получили производную/интеграл “по индексам”

слишком маленький шаг для производной и потеря точности из-за вычитания близких чисел

шумные данные и “взрыв” производной

неравномерная сетка, а формулы применены как для равномерной

Мини-чек-лист перед вычислением:

Проверьте dtype и используйте float64, если точность важна.

Убедитесь, что сетка x отсортирована и соответствует y.

Для производной по данным начните с np.gradient(y, x).

Для интеграла по данным начните с np.trapezoid(y, x).

Проведите проверку сходимости: сравните результат при сгущении сетки.

Для шумных данных сначала аппроксимируйте (МНК), затем дифференцируйте модель.

Итог

Численные производные строятся конечными разностями; центральные разности обычно точнее, но выбор шага ограничен округлением.

np.gradient — практичный инструмент для производных на равномерных и неравномерных сетках.

Численное интегрирование по данным чаще всего начинают с метода трапеций; в NumPy это np.trapezoid.

Дифференцирование усиливает шум, интегрирование сглаживает — поэтому для табличных шумных данных полезны методы аппроксимации и МНК из предыдущего модуля.

Нелинейные уравнения и основы оптимизации

Во многих задачах вычислительной математики линейная алгебра — это лишь внутренний механизм вычислений. Реальная постановка часто выглядит иначе:

найти корень нелинейного уравнения ;

решить систему нелинейных уравнений ;

минимизировать функцию ошибки (оптимизация) .

Эта статья связывает предыдущие модули курса в единую картину:

из темы про погрешности и устойчивость нам нужны корректные критерии остановки и понимание, почему «хорошая формула» может быть численно плохой;

из темы про линейные системы и разложения нам нужен инструмент для шага Ньютона: решение систем с матрицей Якоби;

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

В рамках курса мы делаем акцент на решениях, которые можно реализовать на чистом NumPy (без SciPy), и на диагностике качества результата.

Нелинейное уравнение: что значит «найти корень»

Задача: найти такое , что .

— заданная функция (возвращает число);

— неизвестное число, которое мы подбираем;

— значение, к которому хотим привести функцию.

На практике корень ищут итерациями: строят последовательность , которая должна приближаться к решению.

Что важно определить заранее

Есть ли гарантии существования корня на интервале?

Нужна ли гарантированная сходимость или важнее скорость?

Как измерять «достаточно близко»?

Последний пункт особенно важен из-за чисел с плавающей точкой: сравнение с нулём через == обычно плохая идея. В NumPy для сравнения с допусками используют np.isclose и np.allclose.

Документация:

numpy.isclose

numpy.allclose

Критерии остановки: когда итерации можно прекратить

Обычно применяют комбинацию критериев.

Критерий по невязке: мало.

Критерий по шагу: мал.

Ограничение по числу итераций.

Типичная форма критерия по шагу:

Пояснение символов:

— текущее приближение;

— следующее приближение;

— модуль (абсолютное значение);

— абсолютный допуск (важен около нуля);

— относительный допуск (важен для больших масштабов).

Методы с гарантией: бисекция (деление пополам)

Если непрерывна на и знаки на концах разные, то корень внутри существует:

Пояснение символов:

— границы интервала;

— знак числа (минус/плюс);

условие означает, что функция «пересекает ноль» внутри.

Метод бисекции:

берём середину ;

выбираем половину интервала, где знак меняется;

повторяем.

Плюсы:

очень устойчивый;

гарантированно сходится при выполнении условия смены знака.

Минусы:

сходимость относительно медленная.

Пример реализации на NumPy:

Быстрые методы: Ньютон и секущие

Метод Ньютона для одного уравнения

Итерация Ньютона строится так:

Пояснение символов:

— текущее приближение;

— следующее приближение;

— значение функции в текущей точке;

— производная функции в текущей точке;

дробь показывает, как касательная «подсказывает» поправку к .

Интуиция: мы заменяем график касательной в и берём точку пересечения касательной с осью .

!Геометрический смысл шага метода Ньютона

Плюсы:

очень быстрая сходимость около корня (часто говорят о квадратичной скорости).

Минусы:

требуется производная;

при плохом начальном приближении может расходиться;

если близка к нулю, шаг становится огромным и метод теряет устойчивость.

Пример реализации:

Если производной нет: метод секущих

Метод секущих приближает производную разностью:

Пояснение символов:

— два последних приближения;

числитель — разность значений функции;

знаменатель — разность аргументов;

отношение — оценка наклона секущей.

Плюсы:

не требует аналитической производной;

часто почти так же эффективен, как Ньютон.

Минусы:

нет полной гарантии сходимости без дополнительных условий;

может быть нестабилен при плохом старте.

Система нелинейных уравнений: многомерный Ньютон

Теперь — это вектор, а — вектор функций:

Пояснение символов:

— вектор неизвестных (NumPy-массив формы (n,));

— вектор невязок по уравнениям (массив формы (n,));

— нулевой вектор.

Якобиан и шаг Ньютона

В многомерном случае появляется матрица Якоби :

Пояснение символов:

— элемент матрицы Якоби в строке и столбце ;

— -я компонента вектор-функции ;

— -я компонента вектора ;

— частная производная.

Шаг Ньютона записывают через решение линейной системы:

Пояснение символов:

— текущее приближение;

— вектор невязки;

— Якобиан;

— поправка (вектор шага);

первое равенство — линейная система, которую нужно решить;

второе равенство — обновление приближения.

Ключевая связь с прошлым модулем: каждая итерация Ньютона — это решение линейной системы, поэтому используем np.linalg.solve, а не явную инверсию.

Документация:

numpy.linalg.solve

numpy.linalg.norm

Численный Якобиан (конечные разности)

Если аналитического Якобиана нет, его часто приближают конечными разностями. Для столбца можно использовать:

Пояснение символов:

— базисный вектор (все нули, кроме 1 на позиции );

— малый шаг для разностей;

числитель — изменение при малом изменении ;

деление на — оценка производной.

Практически нельзя делать «как можно меньше»: слишком маленький шаг усиливает ошибки округления (связь с темой про численную точность). Обычно разумно начинать с – для масштаба около 1.

Пример функции для численного Якобиана:

Ньютон для системы на NumPy

Комментарии к критериям:

np.linalg.norm(fx) — это «размер» невязки по всем уравнениям;

np.linalg.norm(dx) — насколько сильно меняется решение на текущем шаге;

множитель 1.0 + norm(x_new) добавляет масштабную устойчивость критерия.

Почему методы могут работать плохо: масштабирование и обусловленность

Даже «правильный» метод может вести себя плохо, если:

переменные имеют очень разные масштабы (например, одна компонента порядка , другая порядка );

Якобиан плохо обусловлен (аналогично плохой обусловленности матрицы в );

начальное приближение далеко от решения.

Практические приёмы:

приводите переменные к сопоставимым масштабам (нормализация);

проверяйте нормы и невязки, а не только факт «что-то сошлось»;

если Ньютон делает слишком большие шаги, часто помогает демпфирование шага: с .

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

Задача оптимизации в базовой форме:

Пояснение символов:

— неизвестный параметр (число или вектор);

— целевая функция (часто это ошибка, энергия, штраф);

означает «найти , при котором значение минимально».

Связь с нелинейными уравнениями:

если гладкая, то в точке минимума часто выполняется условие ;

это превращает оптимизацию в задачу решения системы нелинейных уравнений.

Градиент: что это и зачем

Градиент — это вектор производных по компонентам:

Пояснение символов:

— компоненты вектора ;

— насколько быстро меняется при изменении ;

весь вектор указывает направление наиболее быстрого роста .

Значит, для уменьшения функции логично идти в направлении против градиента.

!Интуиция градиентного спуска по линиям уровня

Градиентный спуск: базовый алгоритм

Итерация градиентного спуска:

Пояснение символов:

— текущее приближение;

— градиент в текущей точке;

— шаг (скорость обучения);

минус означает движение в сторону уменьшения .

Выбор — критичный момент:

слишком большой шаг может привести к расходимости;

слишком маленький — к очень медленной сходимости.

Простейшая реализация (градиент задаём функцией):

Численный градиент конечными разностями

Если аналитического градиента нет, можно использовать разности:

Пояснение символов:

— текущая точка;

— малый шаг;

— базисный вектор, меняющий только компоненту .

Пример:

Метод Ньютона для оптимизации: когда нужен второй порядок

Если доступна матрица вторых производных (Гессиан) , можно использовать шаг Ньютона для минимизации:

Пояснение символов:

— градиент;

— Гессиан (матрица вторых производных);

— шаг;

линейная система снова решается через np.linalg.solve.

Практический смысл:

градиентный спуск использует только направление;

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

Ограничения:

вычислять и хранить может быть дорого;

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

Важная связка с МНК

Из модуля про МНК вы уже знаете задачу:

Пояснение символов:

— матрица признаков;

— параметры;

— наблюдения;

— евклидова норма.

Это оптимизация, но квадратичная и выпуклая; для неё есть надёжные прямые методы (np.linalg.lstsq).

Нелинейная аппроксимация часто выглядит так:

Пояснение символов:

— параметры модели;

— невязка по -му измерению;

сумма квадратов — типичная функция ошибки.

Дальнейшее развитие этой темы (в более продвинутых курсах) приводит к методам Гаусса–Ньютона и Левенберга–Марквардта, но даже базовые инструменты из этой статьи уже позволяют решать небольшой класс задач на NumPy.

Мини-чек-лист для практики

Для 1D корня с гарантией используйте бисекцию, если можно выделить интервал со сменой знака.

Для быстрых корней используйте Ньютон, но контролируйте производную и делайте осмысленный выбор начальной точки.

Для систем: формулируйте как NumPy-вектор и решайте шаг Ньютона через np.linalg.solve.

Обязательно проверяйте невязку: в 1D это , в системе это .

В оптимизации начинайте с градиентного спуска и разумного масштаба переменных.

Если производных нет, используйте конечные разности, но подбирайте шаг осознанно.

Итог

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

Метод Ньютона связывает нелинейные уравнения с линейной алгеброй: на каждом шаге решается система с Якобианом.

Оптимизация тесно связана с решением уравнений через условие .

NumPy даёт все базовые инструменты для реализации: векторизацию, нормы (np.linalg.norm), решение систем (np.linalg.solve), сравнение с допусками (np.isclose).