22.4 Пример из реальной жизни, использующий разреженные матрицы

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

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

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

Простейшие модели конечных элементов разделят Омега на простые симплексы (треугольники в 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";