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

56
mport sympy as sm sm.var(’x у’)
p = sm.Plot(visible=False) p[l] = x**2 + y**2
p[1] .style = ’solid’ p [2] = - x**2 - y**2 p [2].style = ’wireframe’ p.showO
Рис. 3.53 получен как копия экрана (горячая клавиша F8) при установленном пакете Python Imaging Library (PIL)⁷⁴.
Разреженные матрицы (SciPy.sparse)
Разреженные матрицы, которые характеризуются большим числом нулевых элементов, возникают при численном решении краевых задач для уравне­ний с частными производными. В этом случае при использовании конечно­разностных или конечно-элементных аппроксимаций имеется зависимость ре­шения дискретной задачи только от небольшого числа соседей. Вычисли­тельные алгоритмы линейной алгебры должны учитывать отмеченную спе­цифику разреженных матриц. Здесь мы отметим возможности пакета SciPy (модуль sparse).
В качестве иллюстративного примера рассмотрим разностную задачу Ди­рихле для уравнения Пуассона в единичном квадрате. Функция -н(х), х = (xi,x₂), удовлетворяет в П = {х | 0 < х_а < 1, а = 1,2 уравнению
д
П

²и д²и дх\ дх²
и однородным граничным условиям
'и(х) = 0, х G дП.
* Для приближенного решения введем равномерную разностную сетку с шагом h (Nh = 1) по каждому направлению. Разностное решение i = 0,1,..., TV, j = 0,1,..., TV определяется как решение сеточной задачи
Vi—ij Уг+ij Ui,j—i yij+1 fi,j> ^ 1,2, ...,7V 1, j 1,2, ...,7V 1,
u_id = 0, i = 1,7V, j = 0,7V,
где f_id = h² f(ih_y jh).
Сформулируем матричную задачу для определения разностного решения во внутренних узлах сетки. Нумерацию узлов проведем но правилу
j = 1,2,...,7V - 1 : i = 1,2, ...,7V - 1 : v_k = y_ij} к = i + (N - 1 )j.
Общее число неизвестных равно (N — I)². В представлении
Av — b
разреженная матрица имеет диагональную структуру (ненулевые элемет содержат только пять диагоналей). Вместо 0(N⁴) нужно хранить толь 0(N²) ненулевых элементов.
Для работы с разреженными матрицами используются специальные форм ты хранения матриц. Простейший из них COOrdinate format (С00) связан использованием трех массивов. Первый из них содержит значения ненулевь элементов (aij ф 0), второй — номер строки (г), а третий — номер столб) (г), т.с. храниться тройка (а^-,г, j).
Для выполнения основных операций над разреженными матрицами (слож ние, умножение, транспонирование, решение систем линейных уравнений) б лее удобен разреженный строчный формат — Compressed Sparse Row forma (CSR). В этом случае значения ненулевых элементов и номера столбцов хр нятся по строкам в двух массивах: ац (массив АЛ) и j (массив JA). Трет! массив (обозначим его IA) отвечает за номер строки: его г-й элемент указ] вает, с какой позиции в массивах элементов и столбцов (АА и JA) начинает г-я строка. В силу этого IA(i + 1) - IA(г) есть число элементов в г-й строк
Имеются также варианты разреженного формата для хранения по столбце — Compressed Sparse Column format (CSC), по блокам, для симметричных ленточных матриц. Основные форматы поддерживаются в модуле sparse п кета SciPy.
Первая проблема при работе с разреженными матрицами связана с задание матрицы в нужном формате. Мы можем начать с простейшего способа, к торый связан с заданием матрицы как полной и последующей конвертации нужный формат. Аналогично стандартному массиву для инициализации ра реженной матрицы используется lil\_matrix(), после этого можно задава^г отдельные ненулевые элементы. Для больших матриц нужно ориентирован¹ ся на явное задание матрицы в разреженном формате.
Для работы с Compressed Sparse Row форматом в модуле sparse предназначе! класс csr\_matrix(), для Compressed Sparse Column — csc\_matrix(), ДJ DlAgonal format — dia\_matrix(), для COOrdinate format — coo\_matrix().
from scipy import sparse import matplotlib.pyplot as pit import time def poisson_lil(n):
A = sparse.lil_matrix((n*n, n*n)) for j in range(n):
for i in range(n): k = i + n*j A [k, k] = 4
i



Рис. 3.54 Портрет разностного оператора Лапласа

f i > 0:
A[k,k-1] = -1 if i < n-1:
A[k,k+1] = -1 if j > 0:
A [k,k-n] = -1 if j < n-1:
A[k,k+n] = -1
return A
def poisson_coo(n):