Это старая версия документа!
Вычислительные задачи
Функция на сетке
Важная особенность np.array — возможность обходиться при вычислении значения функции без циклов.
import numpy as np
def f(x,y):
return x**2+2*x+6*y
s = np.linspace(1,7,3)
p = np.linspace(5,9,3)
v = f (s,p)
# v == array([ 33., 66., 117.])
При этом в функции в качестве x и y передаются массивы (никаких векторных операций в стиле языка R в Python нет). Об этом надо помнить при написании функции. Также необходимо чтобы массивы имели одинаковый shape.
Некоторые операции на сетках
| Метод | Назначение |
|---|---|
| np.diff | «Дифференциал» от массива по выбранной оси |
| np.gradient | «Градиент» от N-мерного массива |
| np.trapz | «Интеграл» от массива по выбранной оси |
| np.interp | Одномерная кусочно-линейная интерполяция |
Примеры
diff
import numpy as np
x = np.linspace(0,1,5)
# x == array([ 0. , 0.25, 0.5 , 0.75, 1. ])
y = x**2
# y == array([ 0. , 0.0625, 0.25 , 0.5625, 1. ])
dy = np.diff(y)
# dy == array([ 0.0625, 0.1875, 0.3125, 0.4375])
dx = np.diff(x)
# dx = array([ 0.25, 0.25, 0.25, 0.25])
y_dev = dy/dx
# y_dev == array([ 0.25, 0.75, 1.25, 1.75])
Полученных значений меньше на 1 чем исходных. Можно получать «диференциал» n раз и по определенной оси многомерного массива (аргумент axis).
gradient
import numpy as np
x = np.linspace(0,1,5)
# x == array([ 0. , 0.25, 0.5 , 0.75, 1. ])
y = x**2
# y == array([ 0. , 0.0625, 0.25 , 0.5625, 1. ])
dy = np.gradient(y)
# dy == array([ 0.0625, 0.125 , 0.25 , 0.375 , 0.4375])
dx = np.gradient(x)
# dx = array([ 0.25, 0.25, 0.25, 0.25, 0.25])
y_dev = dy/dx
# y_dev == array([ 0.25, 0.5 , 1. , 1.5 , 1.75])
В одномерном случае аналогичен diff, но возвращает столько же значений, сколько исходных данных.
trapz
import numpy as np
x = np.linspace(0,1,5)
# x == array([ 0. , 0.25, 0.5 , 0.75, 1. ])
s = np.trapz(x,dx=0.25)
# s == 0.5
Интегрирует вдоль оси (аргумент axis) с шагом dx или принимает два массива: y и x (отсчеты по оси, по умолчанию 1). Интегрирования выполняется методом трапеций.
interp
import numpy as np
x = np.linspace(0,1,5)
# x == array([ 0. , 0.25, 0.5 , 0.75, 1. ])
y = -x**2
# y = array([-0. , -0.0625, -0.25 , -0.5625, -1. ])
pts = [0.1, 0.9]
v = np.interp(pts, x, y)
v == array([-0.025, -0.825])
В качестве аргументов принимает: координату точек в которых надо рассчитать значения, координату и значение в известных точках.
Модуль scipy
Для решение более содержательных математических задач существует модуль scipypip.
Решение уравнений
В случае одной переменной
import scipy.optimize
def foo(x): return 4*x+6
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)
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
optobj = scipy.optimize.minimize(f, [2, -1])
if optobj.success:
print("x = ", optobj.x)
# [ 0.99999991, 0.99999979]
Приближенное решение системы линейных уравнений
import scipy.linalg
# Переопределенная система
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])
p, res, rnk, s = scipy.linalg.lstsq(a,b)
#p == 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
# Исходный зашумленный сигнал
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
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)
# Гридинг, где points — координаты точек,
# а values — значения в точках.
data = scipy.interpolate.griddata(points, values, (grid_x, grid_y))
foo = RegularGridInterpolator((xi,yi),data)
func(-0.5,1)
# 1.0
foo([-0.5,1])
# array([ 1.08])
Кригинг
Отличие от других методов интерполяции — наилучшее линейное несмещенное предсказание промежуточных значений. В scipy криггинга нет, но он есть в пакете pykrigepip.
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.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 существует достаточно простая система символьной математики sympypip.
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.