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

Численные расчёты с SciPy

Численные расчёты с SciPy

В этом уроке мы переходим от "ручных" реализаций численных методов к готовым алгоритмам библиотеки SciPy. Мы разберем структуру SciPy, численное интегрирование (включая quad и интегрирование по точкам), численное решение ОДУ через solve_ivp, а также базовую интерполяцию и оптимизацию. В конце будут задачи с решениями и чек-лист.

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

Содержание

1. Цели урока

  • Понимать, зачем нужен SciPy и чем он дополняет NumPy.
  • Знать основные подпакеты SciPy и где искать нужные методы.
  • Уметь считать интегралы через scipy.integrate.quad и читать оценку погрешности.
  • Уметь интегрировать табличные данные через trapezoid и simpson.
  • Уметь решать ОДУ (задачу Коши) через scipy.integrate.solve_ivp и читать результат.
  • Уметь делать простую интерполяцию через scipy.interpolate.interp1d.
  • Понимать основы оптимизации через scipy.optimize.minimize и решения уравнений через scipy.optimize.root.
Что особенно важно запомнить: большинство функций SciPy работает с NumPy-массивами. Если вы не уверены в shape и dtype, сначала проверьте их.
^ К оглавлению

2. Контекст: что такое SciPy и как он связан с NumPy

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

  • SciPy -- библиотека Python для научных вычислении: интегрирование, оптимизация, интерполяция, линеиная алгебра, обработка сигналов и т.д.
  • NumPy -- библиотека для массивов и базовых численных операции. SciPy обычно использует numpy.ndarray как основнои формат данных.
  • подпакет (subpackage) -- "раздел" библиотеки, в котором собраны методы по теме. Пример: scipy.integrate.
  • API -- набор функции и правил их использования: какие аргументы, что возвращается.
  • численныи метод -- алгоритм, который дает приближенное решение (с погрешностью).
  • оценка погрешности -- дополнительное число, которое показывает, насколько методу "кажется", что ответ может отличаться от истинного.

Версия SciPy и почему это важно

Функции SciPy развиваются: появляются новые, старые могут устаревать (deprecated) и удаляться. Поэтому полезно уметь проверять версию.

Определение: deprecated (устарело) означает: "пока работает, но использовать не рекомендуется, и позже может быть удалено".
import scipy

print(scipy.__version__)
# 1.11.3
Python
Важно: если у вас выводится не 1.11.3, это не ошибка. Значит, у вас другая версия, и некоторые детали (предупреждения, вывод, доступность функций) могут отличаться.

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

Запомните: SciPy не заменяет NumPy. NumPy -- это массивы и базовые операции. SciPy -- это "продвинутые алгоритмы", которые обычно принимают и возвращают массивы NumPy.
import numpy as np
import scipy

a = np.array([1, 2, 3], dtype=np.float64)

print(a)
# [1. 2. 3.]

print(a.dtype)
# float64
Python
^ К оглавлению

3. Структура SciPy: основные подпакеты и типовой импорт

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

Ниже самые часто используемые подпакеты SciPy (не единственные, но базовые):

  • scipy.special -- специальные математические функции.
  • scipy.integrate -- интегрирование и численное решение ОДУ.
  • scipy.optimize -- оптимизация и решение уравнении (root-finding).
  • scipy.interpolate -- интерполяция.
  • scipy.fft -- преобразования Фурье.
  • scipy.signal -- обработка сигналов.
  • scipy.linalg -- линеиная алгебра (часто шире, чем numpy.linalg).
  • scipy.sparse -- разреженные матрицы и алгоритмы.
  • scipy.sparse.csgraph -- алгоритмы на разреженных графах.
  • scipy.spatial -- пространственные структуры и алгоритмы.
  • scipy.stats -- математическая статистика и распределения.
  • scipy.ndimage -- обработка изображении (n-мерных массивов).
  • scipy.io -- ввод-вывод научных форматов данных.
Определение: разреженная матрица -- матрица, где большинство элементов равно 0, поэтому ее выгодно хранить в специальном формате, чтобы экономить память и ускорять операции.

Как правильно импортировать

Типовые варианты импорта:

  • импортировать конкретные функции: from scipy.integrate import quad
  • импортировать подпакет: from scipy import integrate
from scipy.integrate import quad
from scipy import integrate

print(quad is integrate.quad)
# True
Python

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

Ошибка 1: перепутать SciPy и NumPy (вызвать функцию не там)

Например, искать solve_ivp в numpy.

import numpy as np

try:
    print(np.solve_ivp)
except Exception as e:
    print("ERROR:", type(e).__name__)
# ERROR: AttributeError
Python

Ошибка 2: назвать свой файл scipy.py и "сломать" импорт

Это частая проблема новичков: Python импортирует ваш файл вместо библиотеки. Признак: странные ошибки импорта.

Практическии вывод: не называите свои файлы как numpy.py, scipy.py, random.py и т.д.

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

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

Интегрирование по функции: quad

Определения:
  • определенныи интеграл -- число \(\int_a^b f(x)\,dx\).
  • quad -- функция SciPy для численного интегрирования на отрезке \([a,b]\).
  • оценка абсолютной ошибки -- число, которое quad возвращает вторым элементом, как оценку \(|ошибка|\).

Пример: посчитаем интеграл

$$ \int_0^3 (x^2 - x)\,dx $$

Аналитически это просто, и можно проверить:

$$ \int_0^3 (x^2 - x)\,dx = \left[\frac{x^3}{3} - \frac{x^2}{2}\right]_0^3 = 4.5 $$

from scipy.integrate import quad

val, err = quad(lambda x: x**2 - x, 0, 3)

print(round(val, 4))
# 4.5

print(f"{err:.2e}")
# 5.37e-14
Python
Как читать вывод: первое число -- значение интеграла, второе -- оценка абсолютной ошибки. Если ошибка очень маленькая (например 1e-14), результат обычно очень близок к истинному.

Интегрирование по точкам: trapezoid, simpson

Иногда у вас нет формулы \(f(x)\), а есть только табличные данные: массивы \(x\) и \(y\), где \(y_i \approx f(x_i)\). Тогда используют методы "по точкам".

Определение: табличные данные -- набор точек \((x_i, y_i)\), обычно полученный из измерении или симуляции.

Пример: приблизим тот же интеграл для \(f(x)=x^2-x\) на \([0,3]\), но через точки.

import numpy as np
from scipy.integrate import trapezoid, simpson

x = np.linspace(0.0, 3.0, num=7)
y = x**2 - x

val_trap = trapezoid(y, x)
val_simp = simpson(y, x)

print(np.round(x, 3))
# [0.  0.5 1.  1.5 2.  2.5 3. ]

print(np.round(val_trap, 6))
# 4.625

print(np.round(val_simp, 6))
# 4.5
Python
Практическии вывод: на гладких функциях метод Симпсона часто дает очень точныи результат уже на грубой сетке.

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

Ошибка 1: забыть, что quad возвращает два значения

Если вы ожидаете только число, можно случайно получить кортеж и потом "сломать" код.

from scipy.integrate import quad

res = quad(lambda x: x, 0, 1)

print(res)
# (0.5, 5.551115123125783e-15)

print(res[0])
# 0.5
Python

Ошибка 2: передать в quad функцию, которая возвращает массив

quad ожидает скалярныи результат для каждого x. Если вернуть массив, можно получить TypeError.

import numpy as np
from scipy.integrate import quad

def bad_f(x):
    return np.array([x, x**2])

try:
    print(quad(bad_f, 0, 1))
except Exception as e:
    print("ERROR:", type(e).__name__)
# ERROR: TypeError
Python
^ К оглавлению

5. Численное решение ОДУ: solve_ivp

Что такое IVP и ОДУ

Определения:
  • ОДУ (обыкновенное дифференциальное уравнение) -- уравнение, где неизвестная функция зависит от одного аргумента, например \(y(t)\).
  • задача Коши (initial value problem, IVP) -- ОДУ плюс начальное условие \(y(t_0)=y_0\).
  • solve_ivp -- функция SciPy для решения IVP численно.

Общии вид IVP:

$$ \frac{dy}{dt} = f(t, y), \quad y(t_0)=y_0 $$

Пример: dN/dt = N, N(0)=10

Рассмотрим уравнение роста популяции:

$$ \frac{dN}{dt} = N, \quad N(0)=10 $$

Его аналитическое решение:

$$ N(t) = 10 e^t $$

import numpy as np
from scipy.integrate import solve_ivp

# Правая часть: dN/dt = N
def rhs(t, N):
    return N

t_eval = np.array([0.0, 0.5, 1.0, 1.5, 2.0])

sol = solve_ivp(rhs, t_span=(0.0, 2.0), y0=[10.0], t_eval=t_eval)

print(sol.success)
# True

print(sol.t)
# [0.  0.5 1.  1.5 2. ]

print(np.round(sol.y[0], 3))
# [10.    16.487 27.183 44.817 73.891]
Python
Как читать результат:
  • sol.success -- успешно ли решено.
  • sol.t -- массив значении t.
  • sol.y -- массив значении решения. Если у вас 1 функция, то ее значения лежат в sol.y[0].

Сравним с аналитическим решением \(10 e^t\) на тех же точках:

import numpy as np

t = np.array([0.0, 0.5, 1.0, 1.5, 2.0])
true = 10.0 * np.exp(t)

print(np.round(true, 3))
# [10.    16.487 27.183 44.817 73.891]
Python

Устойчивость и шаг

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

Практическии вывод: если вы видите странные колебания, "взрыв" значении или очень медленную работу, возможно:

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

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

Ошибка 1: перепутать порядок аргументов в rhs

solve_ivp вызывает функцию как rhs(t, y). Если вы написали rhs(y, t), ошибка может быть тихой: код отработает, но решит другое уравнение.

import numpy as np
from scipy.integrate import solve_ivp

# Неправильно: аргументы перепутаны по смыслу
def rhs_wrong(N, t):
    return N  # на самом деле вернется t, потому что первый аргумент будет t

t_eval = np.array([0.0, 2.0])

sol = solve_ivp(rhs_wrong, t_span=(0.0, 2.0), y0=[10.0], t_eval=t_eval)

print(sol.success)
# True

print(np.round(sol.y[0], 3))
# [10. 12.]
Python

12 на конце означает, что решено примерно \(dN/dt = t\), то есть \(N(t)=10 + t^2/2\), а не экспонента.

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

Чаще всего нужно писать y0=[10.0], а не y0=10.0. Покажем тип ошибки безопасно.

from scipy.integrate import solve_ivp

def rhs(t, y):
    return y

try:
    sol = solve_ivp(rhs, t_span=(0.0, 1.0), y0=10.0)
    print(sol.success)
except Exception as e:
    print("ERROR:", type(e).__name__)
# ERROR: TypeError
Python
^ К оглавлению

6. Интерполяция, оптимизация, численная производная: где это в SciPy

Интерполяция: interp1d

Определения:
  • интерполяция -- восстановление значения функции между известными точками.
  • interp1d -- функция SciPy для одномернои интерполяции по точкам.
  • вне диапазона -- x меньше минимального или больше максимального из известных x.

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

import numpy as np
from scipy.interpolate import interp1d

x = np.array([1.0, 2.0])
y = np.array([1.0, 4.0])

f = interp1d(x, y, kind="linear")

print(float(f(1.5)))
# 2.5
Python

Оптимизация: minimize и root

Определения:
  • оптимизация -- поиск минимума или максимума функции.
  • целевая функция -- функция, которую мы минимизируем или максимизируем.
  • minimize -- поиск минимума (обычно численно).
  • root -- поиск корня уравнения \(g(x)=0\).
  • начальное приближение -- стартовая точка поиска (x0).

Минимизация простои функции:

$$ f(x) = (x-2)^2 $$

import numpy as np
from scipy.optimize import minimize

def f(x):
    return (x[0] - 2.0)**2

res = minimize(f, x0=np.array([0.0]))

print(np.round(res.x, 6))
# [2.]

print(np.round(res.fun, 6))
# 0.0
Python

Поиск корня уравнения:

$$ g(x) = x^2 - 2 = 0 $$

import numpy as np
from scipy.optimize import root

def g(x):
    return x**2 - 2.0

sol = root(g, x0=1.0)

print(sol.success)
# True

print(float(np.round(sol.x[0], 6)))
# 1.414214
Python

Численная производная: что изменилось в SciPy

Определения:
  • численная производная -- приближение производнои через значения функции (разностные формулы).
  • scipy.misc.derivative -- старая функция SciPy для численного дифференцирования (в центральнои разностнои форме), которая была помечена как устаревшая.

Важно знать про изменения: scipy.misc.derivative объявлена устаревшеи (deprecated) с SciPy 1.10.0 и полностью удалена в SciPy 1.12.0. Поэтому в учебных проектах лучше не строить код вокруг нее. Если вам нужна численная производная, используите специализированные инструменты или актуальные модули SciPy (в зависимости от версии) и проверяйте документацию.

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

Ошибка 1: интерполяция вне диапазона без явного разрешения

По умолчанию interp1d может падать, если x вне диапазона. Поиграем безопасно через try/except.

import numpy as np
from scipy.interpolate import interp1d

x = np.array([0.0, 1.0])
y = np.array([0.0, 1.0])

f = interp1d(x, y, kind="linear")

try:
    print(float(f(2.0)))
except Exception as e:
    print("ERROR:", type(e).__name__)
# ERROR: ValueError
Python

Ошибка 2: перепутать "minimize" и "root"

minimize ищет минимум, root ищет корень уравнения. Это разные задачи. Если вы решаете \(g(x)=0\), используите root.

import numpy as np
from scipy.optimize import minimize

def g(x):
    return x[0]**2 - 2.0

# minimize найдет минимум g(x), а не корень g(x)=0
res = minimize(g, x0=np.array([1.0]))

print(np.round(res.x, 6))
# [0.]

print(np.round(res.fun, 6))
# -2.0
Python

Видно: минимум получился в x=0, но это не корень уравнения \(x^2-2=0\).

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

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

Блок 1: интегрирование и ОДУ

Задача 1: quad с чтением погрешности

Посчитайте \(\int_0^3 (x^2-x)\,dx\) через quad. Выведите значение и оценку ошибки.

from scipy.integrate import quad

val, err = quad(lambda x: x**2 - x, 0, 3)

print(round(val, 4))
# 4.5

print(f"{err:.2e}")
# 5.37e-14
Python

Задача 2: trapezoid и simpson по точкам

Возьмите 7 точек на [0,3] и посчитайте интеграл \(x^2-x\) по точкам.

import numpy as np
from scipy.integrate import trapezoid, simpson

x = np.linspace(0.0, 3.0, num=7)
y = x**2 - x

print(np.round(trapezoid(y, x), 6))
# 4.625

print(np.round(simpson(y, x), 6))
# 4.5
Python

Задача 3: solve_ivp для dN/dt=N

Решите IVP \(\frac{dN}{dt}=N\), \(N(0)=10\) на [0,2] и выведите значения в точках 0, 0.5, 1, 1.5, 2.

import numpy as np
from scipy.integrate import solve_ivp

def rhs(t, N):
    return N

t_eval = np.array([0.0, 0.5, 1.0, 1.5, 2.0])

sol = solve_ivp(rhs, t_span=(0.0, 2.0), y0=[10.0], t_eval=t_eval)

print(sol.success)
# True

print(np.round(sol.y[0], 3))
# [10.    16.487 27.183 44.817 73.891]
Python

Блок 2: интерполяция и оптимизация

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

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

import numpy as np
from scipy.interpolate import interp1d

x = np.array([1.0, 2.0])
y = np.array([1.0, 4.0])

f = interp1d(x, y, kind="linear")

print(float(f(1.5)))
# 2.5
Python

Задача 5: minimize для (x-2)^2

Найдите минимум функции \(f(x)=(x-2)^2\). Выведите найденное x и значение функции.

import numpy as np
from scipy.optimize import minimize

def f(x):
    return (x[0] - 2.0)**2

res = minimize(f, x0=np.array([0.0]))

print(np.round(res.x, 6))
# [2.]

print(np.round(res.fun, 6))
# 0.0
Python

Задача 6: root для x^2-2=0

Найдите корень уравнения \(x^2-2=0\) рядом с x0=1.

import numpy as np
from scipy.optimize import root

def g(x):
    return x**2 - 2.0

sol = root(g, x0=1.0)

print(sol.success)
# True

print(float(np.round(sol.x[0], 6)))
# 1.414214
Python
^ К оглавлению

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

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

+/-НавыкПроверка
Понимаю роль SciPy Могу объяснить, чем SciPy дополняет NumPy
Ориентируюсь в подпакетах Знаю, где искать integrate, optimize, interpolate, linalg, stats
Проверяю версию SciPy Могу вывести scipy.__version__ и понять, почему это важно
Использую quad Могу посчитать интеграл и распаковать (val, err)
Читаю оценку ошибки quad Понимаю, что второе число это оценка абсолютной ошибки
Интегрирую по точкам Могу применить trapezoid(y,x) и simpson(y,x)
Решаю ОДУ через solve_ivp Могу задать rhs(t,y), t_span и y0, и получить sol.t и sol.y
Понимаю смысл t_eval Могу получить решение в заданных точках, а не только во внутренних шагах solver
Делаю интерполяцию interp1d Могу построить f(x) по точкам и вычислить f(1.5)
Отличаю minimize от root Понимаю, что minimize ищет минимум, а root ищет решение g(x)=0
Слежу за deprecated функциями Знаю, что часть функций может устаревать, и умею проверять документацию
^ К оглавлению
Конспект:
Четверг, 04 июня 2026
Численные расчёты с SciPy