Перейти к содержимому
3.6

Численные методы исчисления в Python.

Численные методы исчисления в Python.

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

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

Содержание

1. Цели урока

  • Понимать разницу между символьным и численным подходом к решению задач.
  • Знать 3 основных источника погрешностеи: неустранимая, методическая, машинная.
  • Уметь приближенно находить производную в точке правои и центральнои разностными формулами.
  • Уметь приближенно вычислять определенныи интеграл методами прямоугольников, трапеций и Симпсона.
  • Понимать, что такое интерполяция, экстраполяция и аппроксимация, и чем они отличаются.
  • Уметь сделать линеиную интерполяцию и интерполяцию Лагранжа на простом примере.
  • Понимать, как выбирать шаг (h) и число разбиении (n), и как проверять результат.
Что особенно важно запомнить: уменьшить шаг или увеличить число разбиении недостаточно. Слишком маленькии шаг может усилить машинную погрешность, и результат станет хуже.
^ К оглавлению

2. Символьные и численные расчеты. Источники погрешностеи

Новые термины и определения

  • символьные вычисления -- подход, где формулы хранятся и преобразуются как символы. Результат обычно "точныи": функция, выражение, матрица без округления (например \(\sqrt{2}\)).
  • численные расчеты -- подход, где мы получаем приближенное число (например 1.4142...), часто быстрее и проще, но с погрешностью.
  • погрешность -- разница между истинным значением и приближенным. Часто используют абсолютную погрешность \(|x_{true} - x_{approx}|\).
  • шаг (обычно \(h\)) -- маленькое положительное число, которое задает "насколько мы сдвигаемся" по аргументу в разностных формулах.
  • порядок точности -- характеристика метода: как быстро уменьшается ошибка при уменьшении шага (или увеличении числа разбиении).

Три вида погрешностеи

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

Покажем машинную погрешность на классическом примере: в десятичном виде кажется, что \(0.1 + 0.2 = 0.3\). Но в двоичном представлении float числа 0.1 и 0.2 не представимы точно, поэтому возникает маленькая разница.

a = 0.1 + 0.2
b = 0.3

print(a)
# 0.30000000000000004

print(a == b)
# False

Частая путаница

Запомните: "численныи ответ отличается от аналитического" не означает, что код неправильныи. Это может быть нормальная методическая или машинная погрешность. Важно уметь оценить, насколько отличие допустимо.
^ К оглавлению

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

Правая разностная производная

Определения:
  • производная в точке \(x_0\) -- скорость изменения функции около точки.
  • разностная производная -- приближение производнои через значения функции в соседних точках.

Правая разностная производная (приближение):

$$ f'(x_0) \approx \frac{f(x_0 + h) - f(x_0)}{h} $$

Определение: первыи порядок точности означает, что при уменьшении шага \(h\) в 10 раз ошибка обычно уменьшается примерно в 10 раз (линеино).
f = lambda x: x**3

x0 = 2
h = 0.1

df = (f(x0 + h) - f(x0)) / h
print(f"{df:.2f}")
# 12.61

Для функции \(f(x)=x^3\) аналитическая производная:

$$ f'(x) = 3x^2 $$

В точке \(x_0=2\) получаем \(f'(2)=12\). При \(h=0.1\) у нас 12.61, то есть близко, но не точно.

f = lambda x: x**3

x0 = 2
h = 0.01

df = (f(x0 + h) - f(x0)) / h
print(f"{df:.2f}")
# 12.06

Центральная разностная производная

Центральная разностная производная использует две точки: \(x_0-h\) и \(x_0+h\).

$$ f'(x_0) \approx \frac{f(x_0 + h) - f(x_0 - h)}{2h} $$

Определение: второи порядок точности означает, что при уменьшении шага \(h\) в 10 раз ошибка обычно уменьшается примерно в 100 раз (квадратично).
f = lambda x: x**3

def dfdx_central(x0, h):
    return (f(x0 + h) - f(x0 - h)) / (2*h)

print(f"{dfdx_central(x0=2, h=0.1):.3f}")
# 12.010

print(f"{dfdx_central(x0=2, h=0.01):.3f}")
# 12.000

Типичные ошибки

Ошибка 1: взять слишком большои шаг

Если \(h\) слишком большои, вы измеряете "среднии наклон" на большом отрезке, а не локальную производную.

f = lambda x: x**3
x0 = 2

for h in [1.0, 0.5, 0.1]:
    df = (f(x0 + h) - f(x0)) / h
    print(h, round(df, 3))
# 1.0 19.0
# 0.5 15.25
# 0.1 12.61

Ошибка 2: взять слишком маленькии шаг и ухудшить результат из-за float

При очень маленьком \(h\) выражение \(f(x_0+h)-f(x_0)\) может терять значимые цифры (вычитание близких чисел). Это проявление машиннои погрешности.

f = lambda x: x**3
x0 = 2

for h in [1e-2, 1e-6, 1e-12]:
    df = (f(x0 + h) - f(x0)) / h
    print(h, df)
# 0.01 12.06009999999994
# 1e-06 12.000006000533403
# 1e-12 11.99928961026391
^ К оглавлению

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

Определения:
  • определенныи интеграл \(\int_a^b f(x)\,dx\) -- "площадь под кривои" (в простом геометрическом смысле, если \(f(x)\ge 0\)).
  • аппроксимация площади -- замена реальнои кривои простыми фигурами (прямоугольники, трапеции, параболы).
  • разбиение -- деление отрезка \([a,b]\) на \(n\) частей.

Метод прямоугольников

numerical_methods-1_ecd85.png

Идея: разбиваем \([a,b]\) на \(n\) равных частей, ширина шага:

$$ h = \frac{b-a}{n} $$

Один из вариантов (правые прямоугольники):

$$ \int_a^b f(x)\,dx \approx \sum_{i=1}^{n} f(x_i)\,h,\quad x_i = a + i h $$

Определение: первыи порядок точности для интегрирования здесь означает: ошибка обычно уменьшается примерно пропорционально \(1/n\) (при увеличении числа прямоугольников).

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

numerical_methods-2_40b98.png

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

$$ \int_a^b f(x)\,dx \approx \sum_{i=1}^{n} \frac{f(x_{i-1}) + f(x_i)}{2}\,h $$

Пример, где легко сравнить с аналитическим ответом:

$$ \int_{1}^{3} \frac{1}{x}\,dx = \ln 3 $$

Определение: натуральныи логарифм \(\ln(x)\) -- логарифм по основанию \(e\), где \(e\) -- основание натурального логарифма.
import math

f = lambda x: 1/x

def integrate_trapezoid(a, b, n):
    res = 0.0
    h = (b - a) / n
    for i in range(1, n + 1):
        x_left = a + (i - 1)*h
        x_right = a + i*h
        res += (f(x_left) + f(x_right)) / 2.0 * h
    return res

for n in [4, 8, 16, 32]:
    print(n, f"{integrate_trapezoid(1, 3, n):.4f}")

print("true", f"{math.log(3):.4f}")
# 4 1.1167
# 8 1.1032
# 16 1.0998
# 32 1.0989
# true 1.0986

Метод Симпсона

numerical_methods-3_4e258.png

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

Определение: четное n означает, что \(n\) делится на 2 без остатка. Пример: 10 четное, 11 нечетное.

Формула (один из стандартных вариантов):

$$ \int_a^b f(x)\,dx \approx \frac{h}{3}\left[f(x_0) + f(x_n) + 4\sum_{\text{нечетн } i} f(x_i) + 2\sum_{\text{четн } i,\ i\ne n} f(x_i)\right] $$

import math

f = lambda x: 1/x

def integrate_simpson(a, b, n):
    if n % 2 != 0:
        raise ValueError("n must be even")
    h = (b - a) / n
    s = f(a) + f(b)
    for i in range(1, n):
        x = a + i*h
        s += (4.0 if i % 2 == 1 else 2.0) * f(x)
    return s * h / 3.0

for n in [4, 8, 16, 32]:
    print(n, f"{integrate_simpson(1, 3, n):.6f}")

print("true", f"{math.log(3):.6f}")
# 4 1.100400
# 8 1.098758
# 16 1.098622
# 32 1.098613
# true 1.098612

Типичные ошибки

Ошибка 1: использовать метод Симпсона с нечетным n

Покажем ошибку безопасно через try/except.

def integrate_simpson(a, b, n):
    if n % 2 != 0:
        raise ValueError("n must be even")
    return 0.0

try:
    print(integrate_simpson(1, 3, 3))
except Exception as e:
    print("ERROR:", type(e).__name__)
# ERROR: ValueError

Ошибка 2: интегрировать функцию с разрывом на интервале без контроля

Например, \(f(x)=1/x\) на интервале [-1, 1] имеет разрыв в 0. Простые формулы могут выдавать мусор или падать.

f = lambda x: 1/x

def integrate_trapezoid(a, b, n):
    res = 0.0
    h = (b - a) / n
    for i in range(1, n + 1):
        x_left = a + (i - 1)*h
        x_right = a + i*h
        res += (f(x_left) + f(x_right)) / 2.0 * h
    return res

try:
    print(integrate_trapezoid(-1, 1, 4))
except Exception as e:
    print("ERROR:", type(e).__name__)
# ERROR: ZeroDivisionError
^ К оглавлению

5. Интерполяция, экстраполяция, аппроксимация

Определения и различия

Определения:
  • дискретные данные -- значения, известные только в конечном наборе точек (например измерения датчика).
  • интерполяция -- построение функции, которая проходит через заданные точки, и вычисление значений внутри диапазона известных \(x\).
  • экстраполяция -- вычисление значений вне диапазона известных \(x\). Обычно рискованнее.
  • аппроксимация -- построение функции, которая "хорошо приближает" данные, но не обязана проходить через все точки.
Запомните: интерполяция и аппроксимация это разные задачи. Интерполяция гарантирует попадание в точки, а аппроксимация часто специально "сглаживает" шум и может не попадать в каждую точку.

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

Пусть известны две соседние точки \((x_1, y_1)\) и \((x_2, y_2)\), и нужно оценить \(f(x)\) для \(x\) между ними. Линеиная интерполяция:

$$ f(x) \approx y_1 + \frac{y_2 - y_1}{x_2 - x_1}\,(x - x_1) $$

def linear_interp(x1, y1, x2, y2, x):
    return y1 + (y2 - y1) / (x2 - x1) * (x - x1)

# Пример: f(x)=x^2, но мы "делаем вид", что знаем только две точки:
# (1, 1) и (2, 4). Оценим f(1.5).
print(linear_interp(1, 1, 2, 4, 1.5))
# 2.5

Проверка с истиннои функциеи \(x^2\): \(1.5^2 = 2.25\). Линеиная интерполяция дала 2.5, потому что заменяет параболу прямои.

print(1.5**2)
# 2.25

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

Определения:
  • полином -- выражение вида \(a_0 + a_1 x + a_2 x^2 + \dots + a_n x^n\).
  • полином Лагранжа -- способ построить интерполяционныи полином по \(n\) точкам так, чтобы он проходил через все точки.

Формула полинома Лагранжа (идея):

$$ P(x)=\sum_{i=1}^{n} f(x_i)\,L_i(x) $$

где базисные полиномы:

$$ L_i(x)=\prod_{\substack{j=1 \\ j\ne i}}^{n}\frac{x-x_j}{x_i-x_j} $$

Ниже не строим "красивую формулу" полинома, а просто вычисляем значение \(P(x)\) в заданнои точке \(x\). Это практичныи подход.

def lagrange_value(xs, ys, x):
    # xs и ys -- списки одинаковои длины
    n = len(xs)
    total = 0.0
    for i in range(n):
        li = 1.0
        for j in range(n):
            if j == i:
                continue
            li *= (x - xs[j]) / (xs[i] - xs[j])
        total += ys[i] * li
    return total

# Возьмем три точки от f(x)=x^2: (0,0), (1,1), (2,4).
# По ним интерполяционныи полином Лагранжа будет ровно x^2.
xs = [0.0, 1.0, 2.0]
ys = [0.0, 1.0, 4.0]

print(lagrange_value(xs, ys, 1.5))
# 2.25

Типичные ошибки

Ошибка 1: пытаться "интерполировать" вне диапазона (получается экстраполяция)

Экстраполяция может сильно ошибаться. В примере ниже мы считаем значение при x=10, хотя точки только 0..2.

def lagrange_value(xs, ys, x):
    n = len(xs)
    total = 0.0
    for i in range(n):
        li = 1.0
        for j in range(n):
            if j == i:
                continue
            li *= (x - xs[j]) / (xs[i] - xs[j])
        total += ys[i] * li
    return total

xs = [0.0, 1.0, 2.0]
ys = [0.0, 1.0, 4.0]

print(lagrange_value(xs, ys, 10.0))
# 100.0

В этом примере получилось "правдоподобно", потому что исходная функция была \(x^2\). Но для реальных шумных данных так повезет не всегда.

Ошибка 2: дубли x-точек (деление на ноль)

В формуле Лагранжа есть деление на \(x_i-x_j\). Если \(x_i=x_j\), будет деление на ноль.

def lagrange_value(xs, ys, x):
    n = len(xs)
    total = 0.0
    for i in range(n):
        li = 1.0
        for j in range(n):
            if j == i:
                continue
            li *= (x - xs[j]) / (xs[i] - xs[j])
        total += ys[i] * li
    return total

xs = [0.0, 1.0, 1.0]  # дублируется 1.0
ys = [0.0, 1.0, 1.0]

try:
    print(lagrange_value(xs, ys, 0.5))
except Exception as e:
    print("ERROR:", type(e).__name__)
# ERROR: ZeroDivisionError
^ К оглавлению

6. Практические выводы: как выбирать шаг и как проверять результат

Шаг, число разбиении, порядок точности

Определения:
  • параметр точности -- то, чем вы управляете в численном методе: шаг \(h\) или число разбиении \(n\).
  • сходимость -- ситуация, когда при "улучшении сетки" (уменьшении \(h\) или увеличении \(n\)) результат приближается к истинному значению.

Главная идея: при уменьшении \(h\) и увеличении \(n\) методическая погрешность уменьшается, но машинная погрешность может начать доминировать. Поэтому часто есть "золотая середина".

Проверки здравого смысла

  • Проверка размерности и типа: что вы вообще считаете (число, вектор, матрица).
  • Проверка граничных случаев: что будет при очень малом \(h\) или очень большом \(n\).
  • Сравнение с аналитическим ответом на простом примере (если возможно).
  • Повтор вычисления другим методом (например трапеции и Симпсон) и сравнение.

Что будет дальше (SciPy)

Многие численные методы имеют "промышленную" реализацию в SciPy: оптимизация, интегрирование, интерполяция, решение ОДУ и т.д. После того как вы поняли смысл базовых формул на простых примерах, правильныи шаг -- использовать готовые и проверенные алгоритмы из SciPy.

Типичные ошибки

Ошибка 1: сравнивать результаты без одинакового формата вывода

Если вы выводите одно число с 2 знаками, а другое с 10, визуально можно сделать неправильныи вывод. Форматируите вывод одинаково.

x = 1.0986122886681098

print(f"{x:.2f}")
# 1.10

print(f"{x:.6f}")
# 1.098612

Ошибка 2: использовать очень сложныи метод без понимания ограничении

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

^ К оглавлению

7. Практика: задачи (с решениями)

Блок 1: производная и интеграл

Задача 1: правая разностная производная

Для \(f(x)=x^3\) оцените \(f'(2)\) при \(h=0.1\) и \(h=0.01\). Выведите оба результата.

f = lambda x: x**3
x0 = 2

for h in [0.1, 0.01]:
    df = (f(x0 + h) - f(x0)) / h
    print(h, f"{df:.2f}")
# 0.1 12.61
# 0.01 12.06

Задача 2: центральная разностная производная

Для \(f(x)=x^3\) оцените \(f'(2)\) при \(h=0.1\) и \(h=0.01\) центральнои формулои.

f = lambda x: x**3

def dfdx_central(x0, h):
    return (f(x0 + h) - f(x0 - h)) / (2*h)

print(f"{dfdx_central(2, 0.1):.3f}")
print(f"{dfdx_central(2, 0.01):.3f}")
# 12.010
# 12.000

Задача 3: метод трапеций для \(\int_1^3 \frac{1}{x}\,dx\)

Посчитайте приближение при n=4, 8, 16, 32 и сравните с \(\ln 3\).

import math

f = lambda x: 1/x

def integrate_trapezoid(a, b, n):
    res = 0.0
    h = (b - a) / n
    for i in range(1, n + 1):
        x_left = a + (i - 1)*h
        x_right = a + i*h
        res += (f(x_left) + f(x_right)) / 2.0 * h
    return res

for n in [4, 8, 16, 32]:
    print(n, f"{integrate_trapezoid(1, 3, n):.4f}")

print("true", f"{math.log(3):.4f}")
# 4 1.1167
# 8 1.1032
# 16 1.0998
# 32 1.0989
# true 1.0986

Блок 2: интерполяция и погрешности

Задача 4: линеиная интерполяция

По точкам (1,1) и (2,4) оцените значение при x=1.5.

def linear_interp(x1, y1, x2, y2, x):
    return y1 + (y2 - y1) / (x2 - x1) * (x - x1)

print(linear_interp(1, 1, 2, 4, 1.5))
# 2.5

Задача 5: интерполяция Лагранжа по 3 точкам

Возьмите точки (0,0), (1,1), (2,4) и вычислите значение при x=1.5. Должно получиться 2.25.

def lagrange_value(xs, ys, x):
    n = len(xs)
    total = 0.0
    for i in range(n):
        li = 1.0
        for j in range(n):
            if j == i:
                continue
            li *= (x - xs[j]) / (xs[i] - xs[j])
        total += ys[i] * li
    return total

xs = [0.0, 1.0, 2.0]
ys = [0.0, 1.0, 4.0]

print(lagrange_value(xs, ys, 1.5))
# 2.25

Задача 6: показать машинную погрешность на 0.1+0.2

Выведите сумму и сравнение с 0.3.

a = 0.1 + 0.2
b = 0.3

print(a)
# 0.30000000000000004

print(a == b)
# False
^ К оглавлению

8. Чек-лист самопроверки

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

+/-НавыкПроверка
Отличаю символьныи и численныи подход Могу объяснить, когда ответ получается "точным", а когда "приближенным"
Понимаю виды погрешностеи Могу привести примеры неустранимои, методическои и машиннои погрешности
Считаю правую разностную производную Могу применить \( (f(x_0+h)-f(x_0))/h \)
Считаю центральную разностную производную Могу применить \( (f(x_0+h)-f(x_0-h))/(2h) \)
Понимаю, что шаг нельзя уменьшать бесконечно Могу объяснить, почему при слишком маленьком h результат может ухудшиться
Считаю интеграл методом трапеций Могу реализовать сумму трапеций и увидеть улучшение при росте n
Знаю ограничение метода Симпсона Помню, что часто требуется четное n, и могу обработать ошибку
Отличаю интерполяцию, экстраполяцию, аппроксимацию Могу объяснить разницу и риски экстраполяции
Делаю линеиную интерполяцию Могу вычислить значение по формуле прямои между двумя точками
Делаю интерполяцию Лагранжа Могу вычислить значение полинома по нескольким точкам и избежать дублирующихся x
Понимаю машинную погрешность float Могу показать пример 0.1+0.2 и объяснить причину
^ К оглавлению

Ссылки по теме

Суббота, 25 апреля 2026
Численные методы исчисления в Python.