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

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

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

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

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

В следующем примере создаётся простая прямоугольная двумерная электропроводящая среда, на противоположные стороны которой подаются напряжения 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 ...

Используя МКЭ первого порядка, мы аппроксимируем распределение электропроводности в Omega постоянной величиной на каждом симплексе (представленной вектором 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: Пример модели конечных элементов с треугольными элементами. Высота каждой вершины соответствует значению решения.