Различия
Показаны различия между двумя версиями страницы.
Предыдущая версия справа и слева Предыдущая версия Следующая версия | Предыдущая версия | ||
calc:math [2021/01/29 11:10] root |
calc:math [2023/03/14 15:39] (текущий) root |
||
---|---|---|---|
Строка 1: | Строка 1: | ||
<div slide> | <div slide> | ||
====== Вычислительные задачи ====== | ====== Вычислительные задачи ====== | ||
- | |||
- | ===== Функция на сетке ===== | ||
- | |||
- | Важная особенность '' | ||
- | |||
- | <sxh python> | ||
- | import numpy as np | ||
- | |||
- | def f(x,y): | ||
- | return x**2+2*x+6*y | ||
- | |||
- | s = np.linspace(1, | ||
- | p = np.linspace(5, | ||
- | |||
- | v = f (s,p) | ||
- | # v == array([ | ||
- | |||
- | S, P = np.meshgrid(s, | ||
- | v = f (s,p) | ||
- | # array([[ 33., 45., 57.], | ||
- | # [ 54., 66., 78.], | ||
- | # [ 93., 105., 117.]]) | ||
- | </ | ||
- | |||
- | </ | ||
- | |||
- | ===== Некоторые операции на сетках ===== | ||
- | |||
- | ^ Метод ^ Назначение ^ | ||
- | | np.diff | " | ||
- | | np.gradient | " | ||
- | | np.trapz | " | ||
- | | np.interp | Одномерная кусочно-линейная интерполяция | | ||
- | |||
- | </ | ||
- | |||
- | ==== Примеры ==== | ||
- | |||
- | === diff === | ||
- | |||
- | <sxh python> | ||
- | import numpy as np | ||
- | |||
- | x = np.linspace(0, | ||
- | # x == array([ 0. , 0.25, 0.5 , 0.75, 1. ]) | ||
- | |||
- | y = x**2 | ||
- | # y == array([ 0. , 0.0625, | ||
- | |||
- | dy = np.diff(y) | ||
- | # dy == array([ 0.0625, | ||
- | |||
- | 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 чем исходных. Можно получать " | ||
- | |||
- | </ | ||
- | |||
- | === gradient === | ||
- | |||
- | <sxh python> | ||
- | import numpy as np | ||
- | |||
- | x = np.linspace(0, | ||
- | # x == array([ 0. , 0.25, 0.5 , 0.75, 1. ]) | ||
- | |||
- | y = x**2 | ||
- | # y == array([ 0. , 0.0625, | ||
- | |||
- | dy = np.gradient(y) | ||
- | # dy == array([ 0.0625, | ||
- | |||
- | 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]) | ||
- | </ | ||
- | |||
- | В одномерном случае аналогичен '' | ||
- | |||
- | </ | ||
- | |||
- | === trapz === | ||
- | |||
- | <sxh python> | ||
- | import numpy as np | ||
- | |||
- | x = np.linspace(0, | ||
- | # x == array([ 0. , 0.25, 0.5 , 0.75, 1. ]) | ||
- | |||
- | s = np.trapz(x, | ||
- | # s == 0.5 | ||
- | </ | ||
- | |||
- | Интегрирует вдоль оси (аргумент '' | ||
- | |||
- | </ | ||
- | |||
- | === interp === | ||
- | |||
- | <sxh python> | ||
- | import numpy as np | ||
- | |||
- | x = np.linspace(0, | ||
- | # x == array([ 0. , 0.25, 0.5 , 0.75, 1. ]) | ||
- | |||
- | y = -x**2 | ||
- | # y = array([-0. | ||
- | |||
- | pts = [0.1, 0.9] | ||
- | |||
- | v = np.interp(pts, | ||
- | v == array([-0.025, | ||
- | </ | ||
- | |||
- | В качестве аргументов принимает: | ||
- | |||
- | </ | ||
- | |||
- | ====== Модуль scipy ====== | ||
Для решение более содержательных математических задач существует модуль | Для решение более содержательных математических задач существует модуль | ||
Строка 137: | Строка 11: | ||
import scipy.optimize | import scipy.optimize | ||
- | def foo(x): | + | def foo(x): |
+ | y = 4*x+6 | ||
+ | return y | ||
+ | # Найдем решение foo(x) = 0 | ||
+ | # с начальным приближением x = 0 | ||
optobj = scipy.optimize.root(foo, | optobj = scipy.optimize.root(foo, | ||
Строка 158: | Строка 36: | ||
return (f1,f2) | return (f1,f2) | ||
+ | # Найдем решение foo(x,y) = 0 | ||
+ | # с начальным приближением x = 5, y = 5 | ||
optobj = scipy.optimize.root(func, | optobj = scipy.optimize.root(func, | ||
Строка 176: | Строка 56: | ||
return 0.5*(1 - x)**2 + (y - x**2)**2 | return 0.5*(1 - x)**2 + (y - x**2)**2 | ||
+ | |||
+ | # Выполним минимизацию f(x,y) -> min | ||
+ | # с начальным приближением x = 2, y = -1 | ||
optobj = scipy.optimize.minimize(f, | optobj = scipy.optimize.minimize(f, | ||
Строка 189: | Строка 72: | ||
<sxh python> | <sxh python> | ||
import scipy.linalg | import scipy.linalg | ||
+ | |||
+ | # Решаем систему Ax=B | ||
# Переопределенная система | # Переопределенная система | ||
Строка 199: | Строка 84: | ||
B = array([ 0. , 1. , 2. , 0.1]) | 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, | + | #x == array([ 0.07299771, -0.23112712, |
</ | </ | ||
Строка 224: | Строка 109: | ||
Поиск происходит методом наименьших квадратов. Важный аргумент '' | Поиск происходит методом наименьших квадратов. Важный аргумент '' | ||
+ | |||
+ | </ | ||
+ | |||
+ | ===== Фильтрация сигнала ===== | ||
+ | |||
+ | <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, | ||
+ | filtered = signal.filtfilt(b, | ||
+ | |||
+ | </ | ||
+ | |||
+ | {{: | ||
</ | </ | ||
Строка 300: | Строка 212: | ||
<sxh python> | <sxh python> | ||
from scipy.interpolate import RegularGridInterpolator | from scipy.interpolate import RegularGridInterpolator | ||
+ | import numpy as np | ||
+ | import matplotlib.pyplot as plt | ||
def func(x, y): | def func(x, y): | ||
Строка 309: | Строка 223: | ||
grid_x, grid_y = np.meshgrid(xi, | grid_x, grid_y = np.meshgrid(xi, | ||
- | # | + | # |
data = func(grid_x, | data = func(grid_x, | ||
- | # Гридинг, где points — координаты точек, | + | # Гридинг |
- | # а values — значения в точках. | + | |
- | data = scipy.interpolate.griddata(points, | + | |
+ | # | ||
+ | xs = np.random.uniform(-5, | ||
+ | ys = np.random.uniform(-2, | ||
+ | |||
+ | # | ||
+ | values = func(xs, ys) | ||
+ | |||
+ | # Линейная интерполяция на нерегулярной сетке | ||
+ | baz = scipy.interpolate.LinearNDInterpolator(np.rollaxis(np.array([xs, | ||
+ | plt.imshow(baz((grid_x, | ||
+ | |||
+ | # Значения на регулярной сетке | ||
+ | data = scipy.interpolate.griddata(np.rollaxis(np.array([xs, | ||
+ | |||
+ | # Интерполяция | ||
foo = RegularGridInterpolator((xi, | foo = RegularGridInterpolator((xi, | ||
- | func(-0.5,1) | + | # точки |
- | # 1.0 | + | plt.scatter(xs, |
+ | |||
+ | # исходное поле | ||
+ | plt.imshow(func(grid_x, grid_y), extent=[-5,5,-2,2], origin=' | ||
+ | |||
+ | # полученное поле на регулярной сетке | ||
+ | plt.imshow(data, | ||
- | foo([-0.5,1]) | + | # результат интерполяции |
- | # array([ 1.08]) | + | plt.imshow(foo((grid_x, grid_y)), extent=[-5,5,-2,2], origin=' |
</ | </ | ||
Строка 353: | Строка 286: | ||
z, ss = UK.execute(' | z, ss = UK.execute(' | ||
</ | </ | ||
+ | |||
+ | |||
+ | </ | ||
+ | |||
+ | ==== Построение функции плотности распределения ==== | ||
+ | |||
+ | <sxh python> | ||
+ | from scipy import stats | ||
+ | |||
+ | # data - массив измерений | ||
+ | # В данном примере одномерный, | ||
+ | |||
+ | kde = stats.gaussian_kde(data) | ||
+ | |||
+ | x_grid = np.linspace(min(data), | ||
+ | |||
+ | kde_vals = kde.evaluate(x_grid) | ||
+ | |||
+ | plt.hist(data, | ||
+ | plt.plot(x_grid, | ||
+ | plt.show() | ||
+ | </ | ||
+ | |||
+ | {{: | ||
</ | </ |