Различия

Показаны различия между двумя версиями страницы.

--- calc:math [2021/01/18 11:43]
127.0.0.1 внешнее изменение
+++ calc:math [2023/03/14 15:39] (текущий)
root
@@ Строка 1: / Строка 1: @@
 <div slide>
 ====== Вычислительные задачи ======
-===== Функция на сетке =====
-Важная особенность ''np.array'' — возможность обходиться при вычислении значения функции без циклов.
-<sxh python>
-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.])
-</sxh>
-При этом в функции в качестве ''x'' и ''y'' передаются массивы (никаких векторных операций в стиле языка R в Python нет). Об этом надо помнить при написании функции. Также необходимо чтобы массивы имели одинаковый ''shape''.
-</div><div slide>
-===== Некоторые операции на сетках =====
-^ Метод ^ Назначение ^
-| np.diff | "Дифференциал" от массива по выбранной оси |
-| np.gradient | "Градиент" от N-мерного массива |
-| np.trapz | "Интеграл" от массива по выбранной оси |
-| np.interp | Одномерная кусочно-линейная интерполяция |
-</div><div slide>
-==== Примеры ====
-=== diff ===
-<sxh python>
-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])
-</sxh>
-Полученных значений меньше на 1 чем исходных. Можно получать "диференциал" ''n'' раз и по определенной оси многомерного массива (аргумент ''axis'').
-</div><div slide>
-=== gradient ===
-<sxh python>
-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])
-</sxh>
-В одномерном случае аналогичен ''diff'', но возвращает столько же значений, сколько исходных данных.
-</div><div slide>
-=== trapz ===
-<sxh python>
-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
-</sxh>
-Интегрирует вдоль оси (аргумент ''axis'') с шагом ''dx'' или принимает два массива: ''y'' и ''x'' (отсчеты по оси, по умолчанию 1). Интегрирования выполняется методом трапеций.
-</div><div slide>
-=== interp ===
-<sxh python>
-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])
-</sxh>
-В качестве аргументов принимает: координату точек в которых надо рассчитать значения, координату и значение в известных точках.
-</div><div slide>
-====== Модуль scipy  ======
 Для решение более содержательных математических задач существует модуль  ''scipy''<sup>pip</sup>.
@@ Строка 133: / Строка 11: @@
 import scipy.optimize
-def foo(x): return 4*x+6
+def foo(x):
+    y = 4*x+6
+    return y
+# Найдем решение foo(x) = 0
+# с начальным приближением x = 0
 optobj = scipy.optimize.root(foo, 0)
@@ Строка 154: / Строка 36: @@
     return (f1,f2)
+# Найдем решение foo(x,y) = 0
+# с начальным приближением x = 5, y = 5
 optobj = scipy.optimize.root(func, (5,5))
@@ Строка 172: / Строка 56: @@
     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])
@@ Строка 185: / Строка 72: @@
 <sxh python>
 import scipy.linalg
+# Решаем систему Ax=B
 # Переопределенная система
@@ Строка 195: / Строка 84: @@
 B = array([ 0. ,  1. ,  2. ,  0.1])
-p, res, rnk, s = scipy.linalg.lstsq(a,b)
+x, res, rnk, s = scipy.linalg.lstsq(A,B)
-#p == array([ 0.07299771, -0.23112712,  0.16181287])
+#x == array([ 0.07299771, -0.23112712,  0.16181287])
 </sxh>
@@ Строка 220: / Строка 109: @@
 Поиск происходит методом наименьших квадратов. Важный аргумент ''bounds'' принимает кортеж из двух списков. Содержит минимальные и максимальные ограничения для значений аргументов.
+</div><div slide>
+===== Фильтрация сигнала =====
+<sxh python>
+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)
+</sxh>
+{{:calc:filt.png?400|}}
 </div><div slide>
@@ Строка 296: / Строка 212: @@
 <sxh python>
 from scipy.interpolate import RegularGridInterpolator
+import numpy as np
+import matplotlib.pyplot as plt
 def func(x, y):
@@ Строка 305: / Строка 223: @@
 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))
+#Нерегулярная сетка
+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)
-func(-0.5,1)
+# точки
-# 1.0
+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')
-foo([-0.5,1])
+# результат интерполяции
-# array([ 1.08])
+plt.imshow(foo((grid_x, grid_y)), extent=[-5,5,-2,2], origin='lower')
 </sxh>
@@ Строка 349: / Строка 286: @@
 z, ss = UK.execute('grid', gridx, gridy)
 </sxh>
+</div><div slide>
+==== Построение функции плотности распределения ====
+<sxh python>
+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()
+</sxh>
+{{:calc:kde.png?400|}}
 </div><div slide>