Чисельне диференціювання та інтегрування: від скінченних різниць до квадратури Гаусса

1. Похибка зрізання проти похибки округлення

Під час чисельного наближення похідної чи інтеграла конкурують два незалежні джерела похибки:

• Похибка зрізання: похибка від заміни справжньої математичної операції (границя при h→0) скінченним наближенням із кроком h. Зменшення h зменшує цю похибку — зазвичай як степінь h.
• Похибка округлення: похибка від використання арифметики з рухомою комою скінченної точності. Зменшення h підсилює округлення, бо воно вимагає віднімання двох майже рівних чисел (катастрофічного скорочення).

Оптимальний крок балансує між цими двома: надто великий — і домінує зрізання; надто малий — і домінує округлення. Для подвійної точності IEEE 754 (ε_mach ≈ 2.2 × 10⁻¹⁶) оптимальний h для центральної різниці першої похідної приблизно дорівнює:

За межею h_opt межу похибки задає округлення, і її не можна зменшити далі зменшенням h. Це фундаментальне обмеження — щоб отримати кращі похідні, потрібна або арифметика вищої точності, або математично інший підхід (автоматичне диференціювання, комплексний крок).

2. Формули скінченних різниць

Усі формули скінченних різниць виводяться з рядів Тейлора. Розклавши f(x + h) навколо x:

Формули першої похідної

Центральна різниця досягає точності другого порядку, бо провідний член похибки (від h²·f'') скорочується при відніманні f(x+h) − f(x−h). Саме тому центральні різниці є кращими щоразу, коли доступні обчислення функції з обох боків x.

Друга похідна

Формула другої похідної використовує три точки й має похибку O(h²). Зверніть увагу, що коефіцієнт похибки — 1/12 × f⁽⁴⁾, тож функції з великими четвертими похідними (як-от sin(100x)) страждають більше. Ця формула всюдисуща у фізиці: вона реалізує дискретний 1D-лапласіан і є основою розв'язувачів PDE на скінченних різницях.

Мішані частинні похідні

Цей 4-точковий хрестоподібний шаблон має похибку O(h²) і використовується для обчислення матриць Гессе для оптимізації та в розв'язувачах задач пружності.

3. Екстраполяція Річардсона

Екстраполяція Річардсона — потужний прийом для підвищення порядку точності будь-якого наближення скінченними різницями без виведення нової формули. Ідея: якщо відомо, що похибка наближення D(h) має вигляд

то, обчисливши при h та h/2, ми можемо усунути провідний член h²:

Це фактично підвищило центральну різницю O(h²) до O(h⁴) — так само, як 5-точковий шаблон — за допомогою лише двох дешевших обчислень. Процес можна повторювати рекурсивно (інтегрування Ромберга для інтегралів або ітеративний Річардсон для похідних), щоб досягти довільно високого порядку.

Інтегрування Ромберга

Інтегрування Ромберга застосовує екстраполяцію Річардсона до складеного правила трапецій:

На практиці інтегрування Ромберга збігається настільки швидко, що навіть R[4][4] (потребуючи лише 17 обчислень функції) досягає точності подвійної точності для більшості гладких функцій. Це основний метод для одновимірних інтегралів на замкнених інтервалах.

4. Диференціювання вищого порядку та з комплексним кроком

Шаблони вищого порядку через невизначені коефіцієнти

Будь-який шаблон скінченних різниць можна вивести систематично. Щоб знайти k-ту похідну з n точками, запишіть:

Коефіцієнти cᵢ знаходять, вимагаючи, щоб формула була точною для многочленів 1, x, x², ..., x^(n−1) (тобто розв'язуючи систему Вандермонда). Сайт findiff автоматизує це. Для робочого коду надійно реалізують це findiff.coefficients() зі scipy або numpy.gradient() з numpy.

Диференціювання з комплексним кроком

Прийом комплексного кроку — мабуть, найбільш недооцінений прийом чисельного диференціювання:

Вирішальна перевага: оскільки ми беремо уявну частину, немає віднімання майже рівних чисел. Похибка округлення суто від обчислення f у комплексній точці, а не від скорочення. Ми можемо використати h = 10⁻²⁰ й отримати майже точні перші похідні, обмежені лише власною точністю обчислення функції. Єдина вимога: функцію f має бути можливо реалізувати над комплексними числами (легко робиться в Python, Julia, C++ з std::complex).

5. Правила інтегрування Ньютона-Котеса

Правила Ньютона-Котеса інтегрують, підганяючи многочлен через рівновіддалені точки й інтегруючи многочлен точно.

Складені правила

Одне правило Сімпсона на [a, b] використовує 3 точки. Складене правило Сімпсона ділить [a, b] на n (парне число) підінтервалів ширини h = (b−a)/n і застосовує Сімпсона на кожну пару:

Чому не Ньютон-Котес вищого порядку?

Для n ≥ 8 ваги Ньютона-Котеса містять від'ємні значення — тобто формула віднімає додатні значення функції, що призводить до катастрофічного скорочення та поганих числових властивостей. Це феномен Рунге для інтегрування. Натомість для високої точності слід використовувати квадратуру Гаусса (яка хитро розташовує вузли в нерівномірних позиціях, щоб цього уникнути) або складений Сімпсон/трапеції з екстраполяцією Річардсона (Ромберг).

6. Квадратура Гаусса

Тоді як Ньютон-Котес фіксує позиції вузлів і знаходить ваги, квадратура Гаусса спільно оптимізує обидва. n-точкове правило Гаусса точне для многочленів степеня до 2n − 1 — удвічі точніше за n-точкове правило Ньютона-Котеса. Усі ваги додатні (уникає скорочення), а вузли — корені ортогональних многочленів.

Квадратура Гаусса-Лежандра

Для інтегралів на [−1, 1] правило Гаусса-Лежандра використовує вузли xᵢ = корені многочлена Лежандра Pₙ(x) та ваги:

Щоб інтегрувати на [a, b] замість [−1, 1], застосуйте заміну змінних x = (b−a)/2 · t + (a+b)/2:

Правила Гаусса-Кронрода (G7K15)

Багато адаптивних інтеграторів (QUADPACK, scipy.integrate.quad) використовують 15-точкове правило Гаусса-Кронрода. Воно вкладає 7-точкове правило Гаусса-Лежандра як підмножину: 7 наявних вузлів повторно використовуються, а 8 нових вузлів додаються. Це дає водночас високоточний результат (15 точок) та оцінку похибки (|G15 − G7|) без додаткових обчислень функції порівняно з обчисленням G15 з нуля. Оцінена похибка керує адаптивним поділом.

7. Адаптивне інтегрування та контроль похибки

Адаптивні алгоритми автоматично концентрують обчислення функції там, де підінтегральна функція змінюється найбільше, разюче підвищуючи ефективність для негладких чи гостропікових функцій.

Адаптивне правило Сімпсона

Найпростіша адаптивна схема оцінює похибку на [a, b], порівнюючи результат Сімпсона на всьому інтервалі із сумою двох результатів на півінтервалах:

Якщо оцінка похибки < допуск × (b−a) / (b−a_total), рекурсуйте на кожній половині; інакше прийміть S₁ + S₂ як локальний результат. Цей процес завершується за O(n log n) обчислень функції для ліпшиц-неперервних підінтегральних функцій.

Реалізація на JavaScript

// Адаптивне інтегрування за Сімпсоном
function adaptiveSimpson(f, a, b, tol = 1e-9) {
  function simpson(f, a, b) {
    const m = (a + b) / 2;
    return (b - a) / 6 * (f(a) + 4 * f(m) + f(b));
  }
  function recurse(a, b, tol, whole, fa, fm, fb, depth) {
    const m = (a + b) / 2;
    const lm = (a + m) / 2, rm = (m + b) / 2;
    const flm = f(lm), frm = f(rm);
    const left  = (m - a) / 6 * (fa + 4 * flm + fm);
    const right = (b - m) / 6 * (fm + 4 * frm + fb);
    const delta = left + right - whole;
    if (depth >= 50 || Math.abs(delta) <= 15 * tol)
      return left + right + delta / 15;
    return recurse(a, m, tol/2, left,  fa, flm, fm, depth+1)
         + recurse(m, b, tol/2, right, fm, frm, fb, depth+1);
  }
  const fa = f(a), fb = f(b), fm = f((a+b)/2);
  const whole = (b - a) / 6 * (fa + 4 * fm + fb);
  return recurse(a, b, tol, whole, fa, fm, fb, 0);
}

// Приклади використання:
console.log(adaptiveSimpson(Math.sin, 0, Math.PI));  // ≈ 2.0
console.log(adaptiveSimpson(x => 1/x, 1, 100));      // ≈ ln(100) ≈ 4.6052
    

Адаптивний Гаусс-Кронрод (у стилі scipy)

Робочі інтегратори, як-от scipy.integrate.quad, використовують правило G7K15 зі стратегією поділу на основі черги з пріоритетами (підінтервал із найбільшою похибкою ділиться першим). Вони також обробляють:

• Нескінченні межі через заміну (наприклад, x = 1/t для [1, ∞))
• Інтегровні особливості через правила Гаусса-Якобі чи Гаусса-Лагерра
• Осцилюючі підінтегральні функції через методи інтегралів Фур'є (QUADPACK QAWO/QAWS)

8. Застосування в симуляції

Обчислення градієнта у фізичних рушіях

Багатьом фізичним рушіям потрібні градієнти функцій енергії чи потенціалу. Для простих аналітичних потенціалів (пружини, гравітація) кращими є точні похідні. Для складних потенціалів чи потенціалів-«чорних скриньок» (силові поля нейронних мереж, табличні потенціали) використовують скінченні різниці або похідні з комплексним кроком.

У молекулярній динаміці міжатомна сила — це F = −∇U(r). Потенціал Леннарда-Джонса:

Чисельне інтегрування в розв'язувачах ОДР

Зв'язок між інтегруванням і розв'язанням ОДР: метод Ейлера — це просто правило прямокутників за лівим кінцем, застосоване до інтегральної форми ОДР:

RK4 — це точно правило Сімпсона, застосоване до підінтегральної функції f(y)! Ця еквівалентність пояснює, чому RK4 досягає точності четвертого порядку (так само, як Сімпсон) і чому йому потрібні 4 обчислення функції на крок.

Аналіз чутливості

У симуляції на основі оптимізації (наприклад, обернена кінематика, фізично-інформовані нейронні мережі) нам потрібні похідні виходів за параметрами. Для малої кількості параметрів скінченні різниці працюють. Для великої кількості параметрів (нейронні мережі з мільйонами ваг) потрібне зворотне поширення (автоматичне диференціювання в зворотному режимі) — скінченні різниці були б у O(n) разів дорожчими. Комплексний крок займає золоту середину: точний до машинної точності, без потреби в AD-інфраструктурі, працює для помірної кількості параметрів (до ~1000 параметрів).