Введение в научный Python-1. Введение в научный Python
Скачать 6.2 Mb.
|
v u a , где a – квадратная матрица (массив) коэффициентов, а v – вектор (массив) правых частей системы. >>> a=np.array([[2, -3],[3,1]]) >>> v=np.array([1,7]) >>> u=la.solve(a,v); print(u) [ 2. 1.] То же решение можно получить с использованием обратной матрицы. >>> a1=la.inv(a) >>> print(np.matmul(a1,v)) [ 2. 1.] В модуле имееется функция solve_banded(...), которая решает уравнение v u a с матрицей коэффициентов a, имеющей ленточную структуру. У ленточных матриц все ненулевые элементы сосредоточены на главной диагонали и нескольких побочных. Для таких матриц нет необходимости хранить все коэффициенты системы, поскольку большинство из них равно нулю и, кроме того, имеются эффективные алгоритмы решения таких систем уравнений. На следующем рисунке слева показан пример ленточной матрицы, а справа – ее «неленточная» форма (прямоугольная матрица, данные которой используются в вычислениях). Неленточная форма строится по следующему правилу. Каждая диагональ исходной матрицы представляется строкой прямоугольной матрицы. Строки, представляющие верхние побочные диагонали, имеют ведущие нули. Количество нулей равно «номеру» побочной диагонали, отсчитываемому от главной диагонали. Строки, представляющие нижние диагонали дополняются нулями справа. Для решения уравнение v u a с ленточной матрицей коэффициентов a используется функция solve_banded((l,u),A,v[,...]). Здесь (l,u) – кортеж, в котором l является количеством ненулевых нижних побочных диагоналей, а u – количеством ненулевых верхних побочных диагоналей. A – матрица коэффициентов в «неленточной» форме (описано выше), v – вектор правой части системы. 99 Пример 1. Решим краевую задачу 0 1 , 0 0 , 1 y y x y методом конечных разностей (ее точное решение известно 2 1 x x y ). В соответствии с этим методом область непрерывного изменения аргумента заменяется дискретным множеством (сеткой), а дифференциальное уравнение заменяется разностным. В нашем примере сеткой будет разбиение отрезка [0,1] точками 1 ..., , 2 , 1 , 0 , n i h i x i , где 1 / 1 n h – шаг сетки. Вторая производная приближается выражением 2 1 1 2 h y y y x y i i i i , где i i x y y . Приближенная задача формулируется следующим образом: найти дискретную функцию (массив) 1 1 0 ..., , , n y y y y , которая удовлетворяет системе уравнений. 0 2 ,..., 2 , 1 , 1 2 0 1 2 1 1 0 n i i i y n i h y y y y Здесь первое и последнее уравнение выражают краевые (граничные) условия. В матричной записи эта система имеет вид: 0 0 1 0 0 0 0 0 0 0 0 0 1 2 1 0 0 0 0 0 0 0 0 1 2 1 0 0 0 0 0 0 0 0 0 0 0 0 1 2 1 0 0 0 0 0 0 0 0 1 2 1 0 0 0 0 0 0 0 0 0 1 2 2 2 2 1 2 3 2 1 0 h h h h y y y y y y n n n Матрица коэффициентов ленточная, и систему уравнений можно решать с помощью функции solve_banded(...) . Решение задачи можно условно разделить на три части. Вначале создаем матрицу коэффициентов в «неленточной» форме. >>>import numpy as np >>>import scipy.linalg as la >>>n=9 # количество узлов сетки # формирование матрицы >>>a0=np.ones(n); >>>a1=np.ones(n)*( -2) >>>a=np.vstack((a0,a1,a0)); >>>a[0][0]=0; a[0][1]=0 >>>a[1][0]=1; a[1][-1]=1 >>>a[2][-1]=0; a[2][-2]=0 >>>print(a) 100 [[ 0. 0. 1. 1. 1. 1. 1. 1. 1.] [ 1. -2. -2. -2. -2. -2. -2. -2. 1.] [ 1. 1. 1. 1. 1. 1. 1. 0. 0.]] Затем формируем вектор правой части системы. >>>h=1/(n-1) >>>v=np.ones(n)*(h**2) >>>v[0]=0 >>>v[n-1]=0 >>>print(v) [ 0. 0.015625 0.015625 0.015625 0.015625 0.015625 0.015625 0.015625 0. ] Теперь решаем систему уравнений. >>>y=la.solve_banded((1,1),a,v); >>>print(y) [ 0. -0.0546875 -0.09375 -0.1171875 -0.125 -0.1171875 -0.09375 -0.0546875 0. ] Для сравнения строим график дискретного решения i i y x , (множество найденных нами точек), и кривой 2 1 x x y , представляющей точное решение. >>>import matplotlib.pyplot as plt >>>xx=np.linspace(0,1,101) >>>yy=xx*(xx-1)/2 >>>plt.plot(xx, yy, '-r', linewidth=3) # кривая точного решения >>>x=np.linspace(0,1,n); >>>plt.plot(x, y, 'o', markersize=10) # точки дискретного решения >>>plt.grid(True) >>>plt.show() Пояснение графических функций вы найдете в следующей главе. Вкратце это выглядит так. Вначале импортируется графический модуль matplotlib.pyplot. Затем с помощью функции plot() рисуются точки и кривая. Функции grid(True) рисует координатную сетку. Функция show() показывает график в отдельном окне. ■ 101 Команда r=lstsq(a,v[,...]) решает систему линейных уравнений v x a методом наименьших квадратов. Решением считается вектор x, для которого евклидова норма вектора x a v минимальна. Функция возвращает кортеж, первым элементом которого является массив решений, вторым элементом – невязка, третьим идет ранг матрицы коэффициентов, а затем массив сингулярных чисел. Т.о. массив решений находится в r[0]. Пример 2. Проведем аппроксимирующую прямую через n точек. Линейная аппроксимация состоит в отыскании коэффициентов a и b уравнения b x a y таких, чтобы все экспериментальные точки лежали максимально близко к аппроксимирующей прямой. В качестве меры близости множества точек n i i i y x 1 , и прямой обычно принимают выражение n i i i b x a y 1 2 , которое и следует минимизировать. Эту задачу можно переформулировать так: решить методом наименьших квадратов переопределенную систему линейных уравнений n i y b x a i i ,..., 2 , 1 относительно неизвестных a и b. Систему можно записать в матричном виде. n n y y y b a x x x 2 1 2 1 1 1 1 Пусть координаты точек равны 1 , 3 , 5 1 , 2 , 5 0 , 1 , 1 , 0 . Следующий код создает массив A коэффициентов системы, и вектор y ее правой части. >>> x = np.array([0, 1, 2, 3]) >>> y = np.array([-1, 0.5, 1.5, 1]) >>> A = np.vstack([x, np.ones(len(x))]).T >>> A array([[ 0., 1.], [ 1., 1.], [ 2., 1.], [ 3., 1.]]) Теперь решаем систему уравнений методом наименьших квадратов. >>> a, b = la.lstsq(A, y)[0] >>> print(a,b) # коэффициенты a и b уравнения прямой 0.7 -0.55 Поясним решение графически. >>> import matplotlib.pyplot as plt >>> plt.plot(x, y, 'o', markersize=10) >>> plt.plot(x, a*x + b, 'r' , linewidth=3) >>> plt.xlim(-0.1,3.1);plt.ylim(-1.1,2.6);plt.grid(True) >>> plt.show() 102 Здесь вначале импортируется графический модуль matplotlib.pyplot. Затем с помощью функции plot() рисуются исходные точки и аппроксимирующая прямая. Функции xlim(), ylim()задают диапазон графика по X и Y координатам, а функция grid(True) рисует координатную сетку. Функция show() показывает график в отдельном окне. ■ Произведение Кронекера матриц (массивов) вычисляется с помощью функции kron(a, b). >>> u=np.array([[2,3]]) >>> v=np.array([[5,10,15]]).T >>> la.kron(u,v) array([[10, 15], [20, 30], [30, 45]]) Функция tril(m[,k=0]) возвращает матрицу с элементами, расположенными выше k–ой диагонали, равными нулю. Например, k=1 соответствует первой верхней побочной диагонали, а k=-1 обнуляет все элементы, которые расположены выше нижней побочной диагонали. >>> la.tril([[1,2,3,4],[4,5 ,6,7],[7,8,9,10],[10,11,12,13]], -1) array([[ 0, 0, 0, 0], [ 4, 0, 0, 0], [ 7, 8, 0, 0], [10, 11, 12, 0]]) Аналогично, функция triu(m[,k=0]) возвращает матрицу со всеми нулевыми элементами, которые расположены ниже k–ой диагонали. Функция w,v=eig(a[,...]) вычисляет собственные значения и собственные векторы квадратной матрицы a. Собственные значения возвращаются в векторе w, и повторяются в нем столько раз, какова их кратность. Собственные векторы возвращаются в матрице v, каждый столбец которой представляет собственные вектор, отвечающий соответствующему собственному значению в векторе w. >>> w, v = la.eig(np.diag((1, 2, 3))) >>> print(w); print(v) [ 1. 2. 3.] [[ 1. 0. 0.] [ 0. 1. 0.] [ 0. 0. 1.]] 103 >>> w, v = la.eig(np.array([[5,4],[8,9]])) >>> w,v (array([ 1.+0.j, 13.+0.j]), array([[-0.70710678, -0.4472136 ], [ 0.70710678, -0.89442719]])) Таким образом, в этом примере собственному числу w[0]=1 отвечает собственный вектор, элементы которого стоят в первом столбце матрицы v. >>> v[:,0] array([-0.70710678, 0.70710678]) Легко проверить, что собственные векторы нормированы. >>> print(v[:,1]) # 2 – й собственный вектор, отвечает с.ч. w[1]=13 [-0.4472136 -0.89442719] >>> la.norm(v[:,1]) 1.0 Для решения системы дифференциальных уравнений x A x , где x – вектор, t d d x x , A – матрица: n x x 1 x , n n n n a a a a A 1 1 1 1 надо найти собственные числа n i i ,..., 1 матрицы A. Если все собственные числа простые, то каждому i соответствует решение t i i i e C v , где i C – произвольная постоянная, i v – собственный вектор, соответствующий этому i Пример 3. Решим задачу Коши для системы дифференциальных уравнений y x y y x x 4 3 2 , 5 0 , 2 0 y x Составим матрицу коэффициентов и найдем ее собственные числа. >>> A=np.array([[2,1],[3,4]]) >>> w, v = la.eig(A); print(w); print(v) [ 1.+0.j 5.+0.j] [[-0.70710678 -0.31622777] [ 0.70710678 -0.9486833 ]] Отбросим, если нужно, нулевые мнимые части собственных чисел. >>> w=w.real; print(w) [ 1. 5.] Общее решение можно записать в виде ) ] 1 [ exp( ] 1 :, [ ) ] 0 [ exp( ] 0 :, [ 1 0 t w C t w C v v , где v[:,0] и v[:,1] являются собственными векторами. Если сюда подставить начальные значения t=0 и 0 0 0 , 0 y y x x , то мы приходим к задаче решения линейной алгебраической системы уравнений 0 0 1 0 y x C C v , 104 где v – матрица, столбцы которой состоят из собственных векторов, т.е. это матрица v, которую вернула функция eig( ). Решим эту систему уравнений. >>> r=np.array([3,5]) # вектор правых частей >>> c=la.solve(v,r); print(c) [-1.41421356 -6.32455532] Построим вектора i i i i t y y t x x , , содержащие значения функций x(t), y(t) решения в дискретные моменты времени (в точках некоторого вектора t). >>> t=np.linspace(0,1,51) # значения независимого аргумента >>> et=np.exp(np.outer(w,t)) >>> cet=c.reshape((2,1))*et # cet=(c*et.T).T >>> x,y=np.dot(v,cet) # значения функций x(t), y(t) в точках вектора t Строим графики решений. >>> import matplotlib.pyplot as plt >>> plt.plot(t, x, ' -b', t,y,'-r', linewidth=3) # построение графиков >>> plt.grid(True) # отображение сетки >>> plt.show() # открыть графическое окно ■ В модуле scipy.linalg есть много других функций, предназначенных, для решения задач линейной алгебры, когда матрицы имеют специальных вид. Кроме того, имеются функции предназначенные для представления матриц через другие матрицы, например, в произведение верхне– и нижне– треугольных. Имеется много функций, предназначенных для конструирования матриц, например, блочно – диагональных, матриц Ханкеля, Гильберта и других. В модуле имеются также матричные функции. Например, функция 0 ! expm n n n M M не совпадает с поэлементной функцией M exp >>> la.expm(M) array([[ 1.54308063, 1.17520119], [ 1.17520119, 1.54308063]]) Сравните >>> np.exp(M) array([[ 1. , 2.71828183], [ 2.71828183, 1. ]]) 105 Имеются матричные функции, такие как cosm(M) (матричный косинус), sinm(A) (матричный синус), sqrtm(A[, ...]) матричный корень квадратный и другие. Кроме модуля scipy.linalg, в пакете NumPy имеется аналогичный модуль numpy.linalg. Многие его функции идентичны функциям модуля scipy.linalg. Однако в numpy.linalg есть несколько функций, которых нет в scipy.linalg. Например, функция numpy.linalg.matrix_power(M, n) возводит матрицу M в степень n. >>> M=np.array([[1,2],[2, -1]]) >>> M2=np.linalg.matrix_power(M, 2) >>> print(M2) [[5 0] [0 5]] >>> M=np.array([[1,2,3],[4,5, -1],[2,5,-3]]) >>> np.linalg.matrix_power(M, 2) array([[15, 27, -8], [22, 28, 10], [16, 14, 10]]) Все перечисленные выше функции выполняли матричные операции с массивами. Если же вы хотите использовать для матричных операций привычные обозначения, тогда вы должны преобразовать массив в матричный тип. Для создания матриц в NumPy предназначен специальный класс matrix, объекты которого всегда являются двумерными. Вместо matrix можно использовать сокращение mat. В отличие от массивов (объектов класса ndarray), для объектов типа matrix операция „*‟ обозначает матричное произведение, а „**‟ – операцию матричного возведения в степень. Атрибут matrix.T возвращает транспонированную матрицу, а атрибут matrix.I – обратную матрицу. Для создания объектов matrix можно использовать специальный синтаксис, когда конструктору матрицы передается текстовая строка. В ней элементы каждой строки матрицы разделяются пробелами, а строки – символом „;‟ (точка с запятой). >>> A = np.matrix('1 1;1 1');print(A) [[1 1] [1 1]] >>> B=np.mat([[1, -2],[-1,1]]);print(B) [[ 1 -2] [-1 1]] >>> A+B matrix([[ 2, -1], [ 0, 2]]) >>> A*B matrix([[ 0, -1], [ 0, -1]]) 106 >>> B*A matrix([[-1, -1], [ 0, 0]]) >>> B**2 matrix([[ 3, -4], [-2, 3]]) >>> B.T matrix([[ 1, -1], [-2, 1]]) >>> B.I # вычисление обратной матрицы matrix([[-1., -2.], [-1., -1.]]) >>> C=B.getI(); print(C) # вычисление обратной матрицы [[-1. -2.] [-1. -1.]] >>> A = np.arange(12) >>> A.shape = (3,4) >>> print(A) [[ 0 1 2 3] [ 4 5 6 7] [ 8 9 10 11]] >>> M = np.mat(A.copy()) # конструктору матрицы можно передавать массив >>> M.T matrix([[ 0, 4, 8], [ 1, 5, 9], [ 2, 6, 10], [ 3, 7, 11]]) Многие функции, предназначенные для работы с матрицами, расположены в модуле numpy.matlib . >>> import numpy.matlib as ml >>> ml.zeros((3,4)) matrix([[ 0., 0., 0., 0.], [ 0., 0., 0., 0.], [ 0., 0., 0., 0.]]) Здесь мы привели только основные процедуры, предназначенные для решения задач линейной алгебры. Имеется еще много других функций и методов, не описанных нами. Вы можете познакомиться с ними самостоятельно по справочной системе. |