====== Вычислительные задачи ====== Для решение более содержательных математических задач существует модуль ''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) {{:calc:filt.png?400|}}
===== Сглаживание апериодического сигнала ===== ==== Подготовка данных ==== 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')
==== Результат сглаживания ==== {{ :calc:figure_1-9.png |}}
===== Интерполяция ===== Будем рассматривать случаи получения функции в неявном виде. ==== В случае одной переменной ==== Метод ''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() {{:calc:kde.png?400|}}
===== Анализ спектра ===== Будем разлагать значения сигнала на сетке в спектр с помощью быстрого преобразования Фурье. 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.