numpy¶
Пакет numpy предоставляет $n$-мерные однородные массивы (все элементы одного типа); в них нельзя вставить или удалить элемент в произвольном месте. В numpy реализовано много операций над массивами в целом. Если задачу можно решить, произведя некоторую последовательность операций над массивами, то это будет столь же эффективно, как в C или matlab - львиная доля времени тратится в библиотечных функциях, написанных на C.
Одномерные массивы¶
from numpy import (array,zeros,ones,arange,linspace,logspace,
float64,int64,sin,cos,pi,exp,log,sqrt,abs,
nan,inf,any,all,sort,hstack,vstack,hsplit,
delete,insert,append,eye,fromfunction,
trace,diag,average,std,outer,meshgrid)
Можно преобразовать список в массив.
a=array([0,2,1])
a,type(a)
print печатает массивы в удобной форме.
print(a)
Класс ndarray имеет много методов.
set(dir(a))-set(dir(object))
Наш массив одномерный.
a.ndim
В $n$-мерном случае возвращается кортеж размеров по каждой координате.
a.shape
size - это полное число элементов в массиве; len - размер по первой координате (в 1-мерном случае это то же самое).
len(a),a.size
numpy предоставляет несколько типов для целых (int16, int32, int64) и чисел с плавающей точкой (float32, float64).
a.dtype,a.dtype.name,a.itemsize
Индексировать массив можно обычным образом.
a[1]
Массивы - изменяемые объекты.
a[1]=3
print(a)
Массивы, разумеется, можно использовать в for циклах. Но при этом теряется главное преимущество numpy - быстродействие. Всегда, когда это возможно, лучше использовать операции над массивами как едиными целыми.
for i in a:
print(i)
Массив чисел с плавающей точкой.
b=array([0.,2,1])
b.dtype
Точно такой же массив.
c=array([0,2,1],dtype=float64)
print(c)
Массив, значения которого вычисляются функцией. Функции передаётся массив. Так что в ней можно использовать только такие операции, которые применимы к массивам.
def f(i):
print(i)
return i**2
a=fromfunction(f,(5,),dtype=int64)
print(a)
a=fromfunction(f,(5,),dtype=float64)
print(a)
Массивы, заполненные нулями или единицами. Часто лучше сначала создать такой массив, а потом присваивать значения его элементам.
a=zeros(3)
print(a)
b=ones(3,dtype=int64)
print(b)
Функция arange подобна range. Аргументы могут быть с плавающей точкой. Следует избегать ситуаций, когда $(конец-начало)/шаг$ - целое число, потому что в этом случае включение последнего элемента зависит от ошибок округления. Лучше, чтобы конец диапазона был где-то посредине шага.
a=arange(0,9,2)
print(a)
b=arange(0.,9,2)
print(b)
Последовательности чисел с постоянным шагом можно также создавать функцией linspace. Начало и конец диапазона включаются; последний аргумент - число точек.
a=linspace(0,8,5)
print(a)
Последовательность чисел с постоянным шагом по логарифмической шкале от $10^0$ до $10^1$.
b=logspace(0,1,5)
print(b)
Массив случайных чисел.
from numpy.random import random,normal
print(random(5))
Случайные числа с нормальным (гауссовым) распределением (среднее 0, среднеквадратичное отклонение 1).
print(normal(size=5))
Операции над одномерными массивами¶
Арифметические операции проводятся поэлементно.
print(a+b)
print(a-b)
print(a*b)
Скалярное произведение
a@b
print(a/b)
print(a**2)
Когда операнды разных типов, они пиводятся к большему типу.
i=ones(5,dtype=int64)
print(a+i)
numpy содержит элементарные функции, которые тоже применяются к массивам поэлементно. Они называются универсальными функциями (ufunc).
sin,type(sin)
print(sin(a))
Один из операндов может быть скаляром, а не массивом.
print(a+1)
print(2*a)
Сравнения дают булевы массивы.
print(a>b)
print(a==b)
c=a>5
print(c)
Кванторы "существует" и "для всех".
any(c),all(c)
Модификация на месте.
a+=1
print(a)
b*=2
print(b)
b/=a
print(b)
Так делать можно.
a+=i
А так нельзя.
i+=a
При выполнении операций над массивами деление на 0 не возбуждает исключения, а даёт значения np.nan или np.inf.
print(array([0.0,0.0,1.0,-1.0])/array([1.0,0.0,0.0,0.0]))
nan+1,inf+1,inf*0,1./inf
Сумма и произведение всех элементов массива; максимальный и минимальный элемент; среднее и среднеквадратичное отклонение.
b.sum(),b.prod(),b.max(),b.min(),b.mean(),b.std()
x=normal(size=1000)
x.mean(),x.std()
Функция sort возвращает отсортированную копию, метод sort сортирует на месте.
print(sort(b))
print(b)
b.sort()
print(b)
Объединение массивов.
a=hstack((a,b))
print(a)
Расщепление массива в позициях 3 и 6.
hsplit(a,[3,6])
Функции delete, insert и append не меняют массив на месте, а возвращают новый массив, в котором удалены, вставлены в середину или добавлены в конец какие-то элементы.
a=delete(a,[5,7])
print(a)
a=insert(a,2,[0,0])
print(a)
a=append(a,[1,2,3])
print(a)
Есть несколько способов индексации массива. Вот обычный индекс.
a=linspace(0,1,11)
print(a)
b=a[2]
print(b)
Диапазон индексов. Создаётся новый заголовок массива, указывающий на те же данные. Изменения, сделанные через такой массив, видны и в исходном массиве.
b=a[2:6]
print(b)
b[0]=-0.2
print(b)
print(a)
Диапазон с шагом 2.
b=a[1:10:2]
print(b)
b[0]=-0.1
print(a)
Массив в обратном порядке.
b=a[len(a):0:-1]
print(b)
Подмассиву можно присвоить значение - массив правильного размера или скаляр.
a[1:10:3]=0
print(a)
Тут опять создаётся только новый заголовок, указывающий на те же данные.
b=a[:]
b[1]=0.1
print(a)
Чтобы скопировать и данные массива, нужно использовать метод copy.
b=a.copy()
b[2]=0
print(b)
print(a)
Можно задать список индексов.
print(a[[2,3,5]])
Можно задать булев массив той же величины.
b=a>0
print(b)
print(a[b])
2-мерные массивы¶
a=array([[0.0,1.0],[-1.0,0.0]])
print(a)
a.ndim
a.shape
len(a),a.size
a[1,0]
Атрибуту shape можно присвоить новое значение - кортеж размеров по всем координатам. Получится новый заголовок массива; его данные не изменятся.
b=linspace(0,3,4)
print(b)
b.shape
b.shape=2,2
print(b)
Поэлементное и матричное умножение.
print(a*b)
print(a@b)
print(b@a)
Умножение матрицы на вектор.
v=array([1,-1],dtype=float64)
print(b@v)
print(v@b)
Внешнее произведение $a_{ij}=u_i v_j$
u=linspace(1,2,2)
v=linspace(2,4,3)
print(u)
print(v)
a=outer(u,v)
print(a)
Двумерные массивы, зависящие только от одного индекса: $x_{ij}=u_j$, $y_{ij}=v_i$
x,y=meshgrid(u,v)
print(x)
print(y)
Единичная матрица.
I=eye(4)
print(I)
Метод reshape делает то же самое, что присваивание атрибуту shape.
print(I.reshape(16))
print(I.reshape(2,8))
Строка.
print(I[1])
Цикл по строкам.
for row in I:
print(row)
Столбец.
print(I[:,2])
Подматрица.
print(I[0:2,1:3])
Можно построить двумерный массив из функции.
def f(i,j):
print(i)
print(j)
return 10*i+j
print(fromfunction(f,(4,4),dtype=int64))
Транспонированная матрица.
print(b.T)
Соединение матриц по горизонтали и по вертикали.
a=array([[0,1],[2,3]])
b=array([[4,5,6],[7,8,9]])
c=array([[4,5],[6,7],[8,9]])
print(a)
print(b)
print(c)
print(hstack((a,b)))
print(vstack((a,c)))
Сумма всех элементов; суммы столбцов; суммы строк.
print(b.sum())
print(b.sum(0))
print(b.sum(1))
Аналогично работают prod, max, min и т.д.
print(b.max(0))
print(b.min(1))
След - сумма диагональных элементов.
trace(a)
Линейная алгебра¶
from numpy.linalg import det,inv,solve,eig
det(a)
Обратная матрица.
a1=inv(a)
print(a1)
print(a@a1)
print(a1@a)
Решение линейной системы $au=v$.
v=array([0,1],dtype=float64)
print(a1@v)
u=solve(a,v)
print(u)
Проверим.
print(a@u-v)
Собственные значения и собственные векторы: $a u_i = \lambda_i u_i$. l - одномерный массив собственных значений $\lambda_i$, столбцы матрицы $u$ - собственные векторы $u_i$.
l,u=eig(a)
print(l)
print(u)
Проверим.
for i in range(2):
print(a@u[:,i]-l[i]*u[:,i])
Функция diag от одномерного массива строит диагональную матрицу; от квадратной матрицы - возвращает одномерный массив её диагональных элементов.
L=diag(l)
print(L)
print(diag(L))
Все уравнения $a u_i = \lambda_i u_i$ можно собрать в одно матричное уравнение $a u = u \Lambda$, где $\Lambda$ - диагональная матрица с собственными значениями $\lambda_i$ по диагонали.
print(a@u-u@L)
Поэтому $u^{-1} a u = \Lambda$.
print(inv(u)@a@u)
Найдём теперь левые собственные векторы $v_i a = \lambda_i v_i$ (собственные значения $\lambda_i$ те же самые).
l,v=eig(a.T)
print(l)
print(v)
Собственные векторы нормированы на 1.
print(u.T@u)
print(v.T@v)
Левые и правые собственные векторы, соответствующие разным собственным значениям, ортогональны, потому что $v_i a u_j = \lambda_i v_i u_j = \lambda_j v_i u_j$.
print(v.T@u)
Преобразование Фурье¶
a=linspace(0,1,11)
print(a)
from numpy.fft import fft,ifft
b=fft(a)
print(b)
Обратное преобразование Фурье.
print(ifft(b))
Интегрирование¶
from scipy.integrate import quad,odeint
from scipy.special import erf
def f(x):
return exp(-x**2)
Адаптивное численное интегрирование (может быть до бесконечности). err - оценка ошибки.
res,err=quad(f,0,inf)
print(sqrt(pi)/2,res,err)
res,err=quad(f,0,1)
print(sqrt(pi)/2*erf(1),res,err)
Дифференциальные уравнения¶
Уравнение осциллятора с затуханием $\ddot{x} + 2 a \dot{x} + x = 0$. Перепишем его в виде системы уравнений первого порядка для $x$, $v=\dot{x}$:
$\frac{d}{dt} \begin{pmatrix}x\\v\end{pmatrix} = \begin{pmatrix}v\\-2av-x\end{pmatrix}$
Решим эту систему численно при $a=0.2$ с начальным условием $\begin{pmatrix}x\\v\end{pmatrix}=\begin{pmatrix}1\\0\end{pmatrix}$
a=0.2
def f(x,t):
global a
return [x[1],-x[0]-2*a*x[1]]
t=linspace(0,10,1000)
x=odeint(f,[1,0],t)
Графики координаты и скорости.
from matplotlib.pyplot import plot
%matplotlib inline
plot(t,x)
Точное решение для координаты.
b=sqrt(1-a**2)
x0=exp(-a*t)*(cos(b*t)+a/b*sin(b*t))
Максимальное отличие численного решения от точного.
abs(x[:,0]-x0).max()