Различия

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

--- calc:math [2021/01/29 11:47]
root
+++ 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.])
-S, P = np.meshgrid(s,p,indexing='ij')
-v = f (S,P)
-# array([[ 33.,  45.,  57.],
-#       [ 54.,  66.,  78.],
-#       [ 93., 105., 117.]])
-</sxh>
-</div><div slide>
-В случае простой подстановки ''f (s,p)'' будет массивы будут подставлены поэлементно. В случае подстановки с использованием ''meshgrid'' из массивов будет построена сетка, в узлах которой будут пары элементов и функция будет рассчитана для каждого элемента сетки.
-===== Дифференциальные операторы на сетке =====
-Для расчетов производных на сетке будем использовать модуль ''findiff''<sup>pip</sup> ([[https://findiff.readthedocs.io|документация]]).
-<sxh python>
-import numpy as np
-from findiff import FinDiff
-x = np.linspace(0, 10, 100)
-dx = x[1] - x[0]
-f = np.sin(x)
-d_dx = FinDiff(0, dx, 1)
-df_dx = d_dx(f)
-</sxh>
-</div><div slide>
-Объект ''FinDiff'' реализует дифференциальный оператор, аргументы: ''0'' — по первой оси, ''dx'' — размер шага, ''1'' — первая производная.
-''df_dx'' — массив ''np.array'' со значениями производной ''f'' в узлах сетки.
-</div><div slide>
-==== Дифференциальные операторы модуля findiff  ====
-^ Название ^ Значение ^ Запись ^ Вход ^ Выход ^
-| FinDiff | дифференциал | d/dx | скалярное поле | скалярное поле |
-| Gradient | градиент | ∇f | скалярное поле | скалярное поле |
-| Laplacian | оператор Лапласа | ∇<sup>2</sup> f | скалярное поле | скалярное поле |
-| Divergence | дивергенция | ∇· f | векторное поле | скалярное поле |
-| Curl | ротор | ∇×f | векторное поле | векторное поле |
-</div><div slide>
-<sxh python>
-import numpy as np
-from findiff import Gradient, Divergence, Laplacian, Curl
-x = np.linspace(0, 10, 100)
-dx = x[1]-x[0]
-y = np.linspace(0, 10, 100)
-dy = y[1]-y[0]
-z = np.linspace(0, 10, 100)
-dz = z[1]-z[0]
-X, Y, Z = np.meshgrid(x, y, z, indexing='ij')
-f = np.sin(X) * np.cos(Y) * np.sin(Z)
-grad = Gradient(h=[dx, dy, dz])
-grad_f = grad(f)
-laplace = Laplacian(h=[dx, dy, dz])
-laplace_f = laplace(f)
-g = np.array([f, 2*f, 3*f])
-div = Divergence(h=[dx, dy, dz])
-div_g = div(g)
-curl = Curl(h=[dx, dy, dz])
-curl_g = curl(g)
-</sxh>
-</div><div slide>
-====== Модуль scipy  ======
 Для решение более содержательных математических задач существует модуль  ''scipy''<sup>pip</sup>.
@@ Строка 106: / Строка 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)
@@ Строка 127: / Строка 36: @@
     return (f1,f2)
+# Найдем решение foo(x,y) = 0
+# с начальным приближением x = 5, y = 5
 optobj = scipy.optimize.root(func, (5,5))
@@ Строка 145: / Строка 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])
@@ Строка 158: / Строка 72: @@
 <sxh python>
 import scipy.linalg
+# Решаем систему Ax=B
 # Переопределенная система
@@ Строка 168: / Строка 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>
@@ Строка 193: / Строка 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>
@@ Строка 269: / Строка 212: @@
 <sxh python>
 from scipy.interpolate import RegularGridInterpolator
+import numpy as np
+import matplotlib.pyplot as plt
 def func(x, y):
@@ Строка 278: / Строка 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>
@@ Строка 322: / Строка 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>