Вычислительные задачи

Для решение более содержательных математических задач существует модуль scipy^pip.

Решение уравнений

В случае одной переменной

import scipy.optimize

def foo(x):
    y = 4*x+6
    return y

# Найдем решение foo(x) = 0
# с начальным приближением x = 0
optobj = scipy.optimize.root(foo, 0)

if optobj.success:
    print("x = ", optobj.x)

Совместное решение

import scipy.optimize

def func(d):
    x,y = d
    f1 = 2*x + 3*y
    f2 = -x**2 + 3
    return (f1,f2)

# Найдем решение foo(x,y) = 0
# с начальным приближением x = 5, y = 5
optobj = scipy.optimize.root(func, (5,5))

if optobj.success:
    print("x = ", optobj.x)
# [ 1.73205081 -1.15470054]

Минимизация функции

import scipy.optimize

def f(c):
    x, y = c
    return 0.5*(1 - x)**2 + (y - x**2)**2


# Выполним минимизацию f(x,y) -> min
# с начальным приближением x = 2, y = -1
optobj = scipy.optimize.minimize(f, [2, -1])

if optobj.success:
    print("x = ", optobj.x)
# [ 0.99999991,  0.99999979]

Приближенное решение системы линейных уравнений

import scipy.linalg

# Решаем систему Ax=B

# Переопределенная система
m = [[  1,   3,   4],
     [  2,  -3,   1],
     [ -3,  -4,   8],
     [1.1, 2.9, 4.1]]
A = np.array(m)

B = array([ 0. ,  1. ,  2. ,  0.1])

x, res, rnk, s = scipy.linalg.lstsq(A,B)

#x == array([ 0.07299771, -0.23112712,  0.16181287])

Функция возвращает кортеж, важнейшие поля: значение, сумма ошибок и эффективный ранг матрицы,

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

import scipy.optimize

def foo(x, a, b, c):
    return a*x**2+b*x+c

xdata = [0,2,4,6]
ydata = [6.1,9.9,29,65]

popt, pcov = scipy.optimize.curve_fit(foo, xdata, ydata)
# popt == array([ 2.0125, -2.285 ,  6.18  ])

Поиск происходит методом наименьших квадратов. Важный аргумент bounds принимает кортеж из двух списков. Содержит минимальные и максимальные ограничения для значений аргументов.

Фильтрация сигнала

import numpy as np
from scipy import signal

# Исходный сигнал
data = ...

#Фильтр Баттерворта

fs = 5 # Частота сигнала герцах
fc = 0.05 # Частота отсечения
# Частота в единицах отсчетов
# (одна для НЧ и ВЧ или кортеж из двух для полосового и режекторного фильтра)
w = fc / (fs / 2)
# Синтез НЧ фильтра 5-го порядка 
b, a = signal.butter(5, w, 'low')
filtered = signal.filtfilt(b, a, data)

Сглаживание апериодического сигнала

Подготовка данных

import numpy as np

# Исходный зашумленный сигнал
rys = np.array(...)

# Окно сглаживания
window_len = 45

# Формируем полуокна для предотвращения краевых эффектов
lead_0 = rys[0]
lead = rys[window_len//2-1:0:-1]
tail_0 = rys[-1]
tail = rys[-2:-window_len//2-2:-1]

Сглаживание

# Дополняем ряд данных полуокнами как четными функциями
ds = np.r_[lead, rys , tail]

# Дополняем ряд данных полуокнами как нечетными функциями
ds = np.r_[2*lead_0-lead, rys , 2*tail_0-tail]

# Выбираем оконную функцию
w = np.hanning(window_len)

# Выполняем свертку
rzlt = np.convolve(w/w.sum(), ds, mode='valid')

Результат сглаживания

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

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

В случае одной переменной

Метод interp1d возвращает функцию от одной переменной.

import scipy.interpolate
x = np.arange(0, 10)
y = 2*x**3-16*x**2
foo = scipy.interpolate.interp1d(x, y)

foo(0.2)
# array(-2.8000000000000003)

В случае нескольких переменных

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

from scipy.interpolate import RegularGridInterpolator
import numpy as np
import matplotlib.pyplot as plt

def func(x, y):
    return 8*x**2-y**3

xi = np.arange(-5,5,0.2)
yi = np.arange(-2,2,0.1)

grid_x, grid_y = np.meshgrid(xi, yi, indexing='ij')

#Прямой расчет на сетке
data = func(grid_x, grid_y)

# Гридинг нерегулярных данных

#Нерегулярная сетка
xs = np.random.uniform(-5, 5, 50)
ys = np.random.uniform(-2, 2, 50)

#Значение в точках
values = func(xs, ys)

# Линейная интерполяция на нерегулярной сетке
baz = scipy.interpolate.LinearNDInterpolator(np.rollaxis(np.array([xs,ys]), 1), values)
plt.imshow(baz((grid_x, grid_y)), extent=[-5,5,-2,2], origin='lower')

# Значения на регулярной сетке
data = scipy.interpolate.griddata(np.rollaxis(np.array([xs,ys]), 1), values, (grid_x, grid_y))

# Интерполяция
foo = RegularGridInterpolator((xi,yi),data)

# точки
plt.scatter(xs, ys)

# исходное поле
plt.imshow(func(grid_x, grid_y), extent=[-5,5,-2,2], origin='lower')

# полученное поле на регулярной сетке
plt.imshow(data, extent=[-5,5,-2,2], origin='lower')

# результат интерполяции
plt.imshow(foo((grid_x, grid_y)), extent=[-5,5,-2,2], origin='lower')

Кригинг

Отличие от других методов интерполяции — наилучшее линейное несмещенное предсказание промежуточных значений. В scipy криггинга нет, но он есть в пакете pykrige^pip.

from pykrige.uk import UniversalKriging
import numpy as np

# Значения в точках : x, y, value
data = np.array([[0.3, 1.2, 0.47],
                 [1.9, 0.6, 0.56],
                 [1.1, 3.2, 0.74],
                 [3.3, 4.4, 1.47],
                 [4.7, 3.8, 1.74]])

gridx = np.arange(0.0, 5.5, 0.5)
gridy = np.arange(0.0, 5.5, 0.5)

# Подготовка объекта кригинга
UK = UniversalKriging(data[:, 0], data[:, 1], data[:, 2],
                      variogram_model='linear',
                      drift_terms=['regional_linear'])

# Кригинг
z, ss = UK.execute('grid', gridx, gridy)

Построение функции плотности распределения

from scipy import stats

# data - массив измерений
# В данном примере одномерный, но многомерные тоже поддерживаются

kde = stats.gaussian_kde(data)

x_grid = np.linspace(min(data),max(data), 1000)

kde_vals = kde.evaluate(x_grid)

plt.hist(data, bins=50, density=True)
plt.plot(x_grid,kde_vals)
plt.show()

Анализ спектра

Будем разлагать значения сигнала на сетке в спектр с помощью быстрого преобразования Фурье.

from scipy.fftpack import fft

x = np.linspace(0.0, 2, 600)
y = np.sin(50.0 * 2.0*np.pi*x) + 0.5*np.sin(80.0 * 2.0*np.pi*x)
yf = fft(y)

# Наблюдаем пик частоты
print(np.abs(yf[98:103]))
# [  22.16734211   41.08774858  286.55468064   57.27069359   25.9449835 ]

Символьные вычисления

Помимо численных расчетов для Python существует достаточно простая система символьной математики sympy^pip.

from sympy import Symbol

x = sympy.Symbol('x')
y = sympy.Symbol('y')

# Важно! eq не значение выражения справа,
# а специальный объект типа sympy.equation
eq = y**2+x

eq.subs(x,2).subs(y,3).evalf()
# 11.0000000000000

eq2 = 1/( (x + 2)*(x + 1) )

sympy.apart(eq2, x)
# -1/(x + 2) + 1/(x + 1)

Примеры задач матанализа

import sympy
from sympy import oo, Symbol

x = Symbol("x")

# Предел
sympy.limit(sympy.sin(x)/x, x, 0)
# 1

# Производная
sympy.diff(sympy.tan(x), x)
# tan(x)**2 + 1

# Сумма ряда
sympy.summation(1/2**i, (i, 0, oo))
# 2

# Первообразная
sympy.integrate(2*x + sympy.sinh(x), x)
# x**2 + cosh(x)

Важно: при использовании символьных вычислений нужно использовать функции библиотеки sympy, а не стандартной библиотеки Python или NumPy/SciPy.