Практикум j практическое примщенше численных методов

О 5 10 15 20 25 30 35


Рис. 3.29 Модули коэффициентов Фурье


it.xlabel(^,$k$^>)
pit.ylabel(’ $ Ix_kI $^y)
pit.figure(2)
yl = fftpack.ifft(x)
pit.plot(t, y, label=⁾exact’)
pit.plot(t, yl, label=’direct-invers’)
pit.xlabel(^,$t$⁾)
pit.legend(loc=l)
pit.show()
Для двумерного преобразования Фурье используются функции fft2() и ifft2(), а для Фурье-преобразований многомерных массивов имеются функ­ции fftnO и ifftn().
В модуле fftpack предусмотрена также возможность работы с сеточными функциями, которые принимают только вещественные значения — функции rfft() и irfftO для прямого и обратного преобразований Фурье соответ­ственно.
На основе дискретного преобразования Фурье реализованы дополнительные операции дифференцирования и интегрирования периодических сеточных функций (функция diff ()), некоторые исевдодифференциальные операции-
Р
0.5


7Г
\


exact
direct-invers


0.4 -


0.3


0.2 -


0.1 - /


о.о


-о.


l.o


0.2


0.4 0.6
t


0.8


1.0

ис. 3.30 Сеточные функции
Поддерживается также возможность быстрого вычисления свертки двух се­точных функций (функция convolve ()). Имеется также возможность исполь­зования прямого и обратного преобразований Гильберта (функции hilbertO и ihilbert()), которые тесно связаны с преобразованиями Фурье.
Линейная алгебра
Пакет SciPy предоставляет богатые возможности по решению задач линей­ной алгебры. Как мы уже отмечали выше, эти задачи представлены в паке­те N urn Ру, здесь же мы коротко опишем инструменты, предоставляемые мо­дулем linalg. Этот модуль содержит функции для вычисления норм матриц и векторов, обращения и факторизации матриц, решения систем линейных уравнений, позволяет находить собственные значения и собственные вектора, вычислять функции матриц. Для решения задач линейной алгебры с разре­женными матрицами, которые особенно важны при научных вычислениях, предназначен модуль sparse.linalg.
В SciPy, как и пакете NumPy, мы можем работать как с массивами, так и с матрицами. При задании матриц поддерживается синтаксис MATLAB (ини­циализация по строке). Напомним базовые операции над матрицами.


Рис. 3.31 Портрет матрицы А



Для вычисления определителя матрицы имеется функция det(). Для нахо­ждения нормы матрицы или вектора предназначена функция normQ. Второй аргумент этой функции принимает значения 1, 2 и inf (2 по умолчанию) и соответствует выбору нормы вектора в 1_р при р = 1,2, оо соответственно и подчиненной нормы матрицы. Здесь inf означает бесконечность (определено в пакете N urn Ру).
Для примера рассмотрим матрицу, которая соответствует разностной про­изводной второго порядка с обратным знаком па интервале 0 < х < 1 для функций, обращающихся в ноль на концах. Используется равномерная сетка с шагом h внутренние узлы ащ i = 0,1,..., т — 1 и h — l/(m+l). Для модель­ной матрицы А размерности гп х гп нарисуем портрет, сосчитаем определи­тель и различные нормы матрицы, проведем разложение Холсцкого, найдем обратную матрицу и получим решение уравнения Ay = Ь при b = 1.
import numpy as np
f



Рис. 3.32 Матрица U (А = UU*)


\

rom scipy import linalg import matplotlib.pyplot as pit m = 25
b = 1. / (m+1)
^A = np.mat(np.zeros((m, m), ’float¹)) for i in range(1, m):
A[i-l,i] = - 1. / h**2 A[i,i] = 2. / h**2 A [i, i-1] = - 1. / h**2 ^A[0,0] = 2. / h**2 ^A[m-l,_m-l] = 2. / h**2 Pit.matshow(A)
Ptint ’det:’, linalg.det(A)
P^rint ’norm-1:’, linalg.norm(A, 1)
p



Рис. 3.33 Портрет матрицы D = А ¹

rint ’norm-2:’, linalg.norm(A, 2) print ’norm-inflinalg.norm(A, np.inf) L = linalg.cholesky(A) pit.matshow(L)
В = linalg.inv(A)
pit.mat show(В, cmap=’gray’)
b = np.mat(np.ones((m), ’float’))
у = B*b.T
pit.figure(4)
x = np.linspace(h, m*h, m)
pit.plot(x, y, x, x*(l-x)/2, ’:’)
plt.xlabel(’$x$’)
pit.ylabel(’$y$’)
pit. showO

0.14


Рис. 3.34 Решение уравнения



:|||i||||||l|i^;l||||;|:||||||||||;;
На рис. 3.31 приведен портрет матрицы А (использовалась функция matshowO из пакета Matplotlib). Матрица А является симметричной и трехдиагональ­ной.
В пашем примере решение системы уравнений проведено с явным исполь­зованием обратной матрицы. Вторая возможность связана с применением функции solve () (вместо у = В*Ь.Т положим у = linalg. solve (А, Ь.Т)). Функция solve_banded() предназначена для решения систем уравнений с ленточной матрицей.
С учетом того, что матрица А симметрична и положительно определена мы представили ее в факторизованном виде А = LL* = U*U (ленгочность мно­жителей матрицы наблюдается на рис. 3.32). Обратная матрица (рис. 3.33), конечно, является полной. Сравнение с точным решением показано на рис. 3.34 — разностный оператор второй производной является точным для полиномов второй степени.

Для более общих матриц в модуле linalg имеются функции LU-разложения. Более того, для работы с любыми невырожденными матрицами реализова­но более общее разложение А = PLU — функция 1и(). Здесь L — нижняя треугольная матрица с диагональными элементами, равными единице, U верхняя треугольная матрица, а Р — матрица перестановок. Другие возмож­ности //{/-разложения реализованы в функциях lu_factor(), lu_solve() (см. документацию по модулю linalg).
Имеются также функции для других разложений матриц:

• 
  разложение Холецкого — функции choleskyO, cholesky_banded(), cho_factor(), cho_solve();
  
• 
  Q/2-разложепие — функция qr();
  
• 
  сингулярное (SVD) разложение — функции svd(), svdvalsO, diagsvdO, diagsvdO;
  
• 
  разложение Шура — функции schurO, rsf2csf () .
  

Использование метода наименьших квадратов проиллюстрируем па задаче подбора коэффициентов сглаживающего полинома. Пусть в точках Xi = ih, i = 0,1, ...,т известны значения функции у_и i = 0,1, ...,т. Эти данные будем приближать с помощью полинома z(a\x) = a₀ + aiX-|-...+a_na;ⁿ, причем п < тп. Коэффициенты а^к = 0,1, ...,п находятся из условия минимума
тп
^J(^a) = YK^z(^a’^x*)-y*)²-
i=1
В матричной записи имеем J(a) = ||Аа — Ь\\².
Реализация метода наименьших квадратов с использованием функции IstsqO дается следующей программой.
import numpy as np from scipy import linalg import matplotlib.pyplot as pit m = 20
n = 3
h = l./m
A = np.mat(np.zeros((m+1, n+1), ’float’))
x = np.zeros((m+1), ’float’) у = np.zeros((m+1), ’float’) yl = np.zeros((m+1), ’float’) for i in range(0, m+1): xx = i*h x[i] = xx
y[i] = np.exp(-xx**2) for k in range(0, n):
A[i,k] = xx**k
a, resid, rank, sigma = linalg.lstsq(A, y) for i in range(0, m+1):


Рис. 3.35 Аппроксимация полиномом



хх = i*h sum = 0.
for k in range(0, n): sum += a[k]*xx**k yl[i] = sum
pit.plot(x, y, 1аЬе1=’exact’)
pit.plot(x, yl, ’ + label=^,polinom^>) pit.xlabel( pit.legend() pit. showO
Рис. 3.34 демонстрирует точность аппроксимации функции е~^х2 полиномом третьей степени при 0 < х < 1. Вместо функции IstsqO можно использо­вать pinv() или pinv2(), которые позволяют найти псевдорешение системы линейных уравнений.
Приведем также пример нахождения собственных значений матриц. Будем искать собственные значения оператора центральной производной для пери­
одических сеточных функций:
У
Ау =
г+1 - Vi-1
2 h
при y_i+m = y_h h = l/m. Собственными значения кососимметричиой матрицы А являются чисто мнимые:
s
v£T
h

Рис. 3.36 Портрет матрицы А


in(27r/c/i),
Ниже представлена программа вычисления собственных значений этой мат­рицы (на рис. 3.36 показан портрет матрицы) с использованием функции eigvalsO модуля linalg. На рис. 3.37 дано сравнение точных и вычисленных собственных значений.
import numpy as np

20 25 30 35

Рис. З.ЗТ Точные и вычисленные собственные значения
from scipy import linalg import matplotlib.pyplot as pit m = 32 h = 1. / m
A = np.mat(np.zeros((m, m), ’float*))
.for i in range(1, m):
A[i-1,i] = 1. / (2*h)
A[i,i-1] = -1. / (2*h)
A[0,m-1] = -1. / (2*h)
A[m-1,0] = 1. / (2*h)
Pit.matshow(A)
Pit.figure(2)
l^a = linalg.eigvals(A)
1^1 = np.sort(la.imag)
■‘■l = np.linspace(0, m-1, m)
Pit.plot(ii, lal, label=*numerical *)
l^E = np.sort(l./h * np.sin(2*np.pi*ii*h)) Pit.plot(ii, laE, *label=*exact’)
Pit.Iegend(loc=0)
Pit.show()

Для нахождения как собственных значений так и собственных функций ис- пользуется функция eig(). Реализованы также алгоритмы решения спек­тральных задач для специальных классов матриц: функции eigvalshQ _и eigh() для комплексных эрмитовых и вещественных симметричных матриц, функции eigvals_banded() и eig_banded() для ленточных матриц.
В модуле linalg имеется поддержка работы с функциями матриц. Для вычис­ления матричной экспоненты





используются три функции: ехрш() — аппроксимация Паде, ехрт2() — спек­тральное разложение, ехртЗО — ряд Тейлора. Для вычисления логарифма имеется функция logmQ. Реализованы также матричный синус, косинус и тангенс (функции sinmQ, cosmO и tanmO), гиперболический синус, косинус и тангенс (функции sinhmO, coshmO и tanhmO). Для вычисления матричная sgn (sign) функции используется signmO, квадратного корня — sqrtmO. От­метим также возможность вычисления функций матриц (с помощью funmO), которые связываются с любой функцией Python.
Рассмотрим пример использования матричной экспоненты при решении за­дачи Коши:
^ + Ау = 0, 0 < t < Т,
2/(0) = и.
При постоянной (не зависящей от t) матрице А решение представляется в виде
y(t) = e~^Atu.
В программе матрица А соответствует второй разностной производной се­точных функций, обращающихся в ноль на концах интервала [0,1]. Точное решение дискретной и непрерывной задачи представляется в виде разложе­ния по собственным функциям у>_к = v^2sin(7rfca;), к = 1,2,....
import numpy as np from scipy import linalg import matplotlib.pyplot as pit n = 3 m = 19 t = 0.1
h = 1. / (m+1)
A = np.mat(np.zeros((m, m), ’float’)) for i in ranged, m):
A [i-1,i] = - 1. / h**2 A[i,i] = 2. / h**2
A



Рис. 3.38 Решение задачи Коши при t = 0.1

[i,i-1] = - 1. / h**2 A[0,0] = 2. / h**2 A[m-l,m-l] = 2. / h**2 x = np.zeros((m), ’float’)

• 
  = np.zeros((m), ’float’) for i in range(0, m):
  

xx = (i+l)*h x[i] = xx
y[i] = np.sin(np.pi*n*xx)
У0 = np.exp(-(np.pi*n)**2 * t) * у

6. 
   = np.exp(-4/h**2*(np.sin(np.pi*n*h/2))**2 * t) * у В = np.mat(linalg.expm(-A*t))
   

У2 = B*np.mat(у).T
Pit.plot(x, yO, label=’exact (continue)’)
Pit.plot(x, yl, label=’exact (dicret)’)
Pit.plot(x, y2, label=’numerical’)
Pit.xlabel(’$x$’)
Pit.ylabel(’$y$’)
Pit.Iegend(loc=0)
Pit.grid(True)
Pit.show()
\