22.4 Пример использования разреженных матриц в реальных задачах

Распространенное применение разреженных матриц — решение моделей конечных элементов. Модели конечных элементов позволяют численно решать дифференциальные уравнения в частных производных, у которых нет замкнутых решений, обычно из-за сложной формы области.

Для мотивации этого применения рассмотрим краевую задачу для уравнения Лапласа. Эта система может моделировать скалярные потенциальные поля, такие как тепло или электрический потенциал. Дана среда Ω с границей ∂Ω. В каждой точке на ∂Ω известны граничные условия, и нам нужно вычислить потенциал в Ω. Граничные условия могут задавать потенциал (условие Дирихле), его нормальную производную через границу (условие Неймана) или взвешенную сумму потенциала и его производной (условие Коши).

В тепловой модели мы хотим рассчитать температуру в Ω и знаем граничную температуру (условие Дирихле) или тепловой поток (из которого мы можем вычислить условие Неймана, разделив на теплопроводность на границе). Аналогично, в электрической модели мы хотим рассчитать напряжение в Ω и знаем граничное напряжение (Дирихле) или ток (условие Неймана после деления на электрическую проводимость). В электрической модели часто большая часть границы является электрически изолированной; это условие Неймана с током, равным нулю.

Простейшие модели конечных элементов разделят Ω на простые симплексы (треугольники в 2D, пирамиды в 3D). В качестве примера 3D рассмотрим цилиндрический резервуар с жидкостью с маленьким непроводящим шаром из проекта EIDORS¹¹. Эта модель разработана для отражения применения электроимпедансной томографии, где к такому резервуару прикладываются токи для отображения внутреннего распределения проводимости. Для описания геометрии КЭ имеем матрицу вершин nodes и симплексов elems.

Следующий пример создает простую прямоугольную 2D электрически проводящую среду с наложенными 10 В и 20 В на противоположных сторонах (граничные условия Дирихле). Все остальные ребра электрически изолированы.

node_y = [1;1.2;1.5;1.8;2]*ones(1,11);
   node_x = ones(5,1)*[1,1.05,1.1,1.2, ...
             1.3,1.5,1.7,1.8,1.9,1.95,2];
   nodes = [node_x(:), node_y(:)];

   [h,w] = size (node_x);
   elems = [];
   for idx = 1:w-1
     widx = (idx-1)*h;
     elems = [elems; ...
       widx+[(1:h-1);(2:h);h+(1:h-1)]'; ...
       widx+[(2:h);h+(2:h);h+(1:h-1)]' ];
   endfor

   E = size (elems,1); # No. of simplices
   N = size (nodes,1); # No. of vertices
   D = size (elems,2); # dimensions+1

Это создает матрицу N×2 nodes и матрицу E×3 elems со значениями, которые определяют треугольники конечных элементов:

nodes(1:7,:)'
    1.00 1.00 1.00 1.00 1.00 1.05 1.05 …
    1.00 1.20 1.50 1.80 2.00 1.00 1.20 …

  elems(1:7,:)'
    1    2    3    4    2    3    4 …
    2    3    4    5    7    8    9 …
    6    7    8    9    6    7    8 …

Используя конечно-элементную модель первого порядка, мы аппроксимируем распределение электрической проводимости в Ω как постоянное на каждом симплексе (представленном вектором conductivity). Основываясь на геометрии конечных элементов, мы сначала вычисляем матрицу системы (или жесткости) для каждого симплекса (представленная как 3×3 элементы на диагонали элементной матрицы системы SE). На основе SE и матрицы связей N×DE C, представляющей связи между симплексами и вершинами, вычисляется глобальная матрица связей S.

## Element conductivity
  conductivity = [1*ones(1,16), ...
         2*ones(1,48), 1*ones(1,16)];

  ## Connectivity matrix
  C = sparse ((1:D*E), reshape (elems', ...
         D*E, 1), 1, D*E, N);

  ## Calculate system matrix
  Siidx = floor ([0:D*E-1]'/D) * D * ...
         ones(1,D) + ones(D*E,1)*(1:D) ;
  Sjidx = [1:D*E]'*ones (1,D);
  Sdata = zeros (D*E,D);
  dfact = factorial (D-1);
  for j = 1:E
     a = inv ([ones(D,1), ...
         nodes(elems(j,:), :)]);
     const = conductivity(j) * 2 / ...
         dfact / abs (det (a));
     Sdata(D*(j-1)+(1:D),:) = const * ...
         a(2:D,:)' * a(2:D,:);
  endfor
  ## Element-wise system matrix
  SE = sparse(Siidx,Sjidx,Sdata);
  ## Global system matrix
  S = C'* SE *C;

Матрица системы действует как проводимость S в законе Ома S * V = I. Основываясь на граничных условиях Дирихле и Неймана, мы можем найти напряжения в каждой вершине V.

## Dirichlet boundary conditions
  D_nodes = [1:5, 51:55];
  D_value = [10*ones(1,5), 20*ones(1,5)];

  V = zeros (N,1);
  V(D_nodes) = D_value;
  idx = 1:N; # vertices without Dirichlet
             # boundary condns
  idx(D_nodes) = [];

  ## Neumann boundary conditions.  Note that
  ## N_value must be normalized by the
  ## boundary length and element conductivity
  N_nodes = [];
  N_value = [];

  Q = zeros (N,1);
  Q(N_nodes) = N_value;

  V(idx) = S(idx,idx) \ ( Q(idx) - ...
            S(idx,D_nodes) * V(D_nodes));

Наконец, для отображения решения мы показываем каждое вычисленное значение напряжения на оси z для каждой вершины симплекса. См. Рисунок 22.6.

elemx = elems(:,[1,2,3,1])';
  xelems = reshape (nodes(elemx, 1), 4, E);
  yelems = reshape (nodes(elemx, 2), 4, E);
  velems = reshape (V(elemx), 4, E);
  plot3 (xelems,yelems,velems,"k");
  print "grid.eps";

Рисунок 22.6: Пример модели конечных элементов, показывающий треугольные элементы. Высота каждой вершины соответствует значению решения.