© Ф.Пинежанинов, П.Пинежанинов

В начало
Зеркало: http://pinega.da.ru/

• Вступление
• Квадратичный функционал в однородных координатах
• Алгоритм учета связи
• Алгоритм конгруэнтного преобразования матриц
• Преобразование в симметричных матрицах
• Прагматика
• Заключение
• Литература

Вступление

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

Введя интерполяцию этих функций, например, используя метод конечных элементов, функционал можно отобразить на функцию от многих переменных. Пусть функция от многих переменных, организованных в вектор   - искомая стационарная точка. Используя формулу Тейлора для многих переменных функцию можно разложить в ряд:. В стационарной точке левая часть равняется 0 , а правая может служить условием для поиска . Первый член в разложении определяет  линейную функцию или гиперплоскость, и не имеет стационарных точек или любая точка в вырожденном случае является стационарной. Таким образом, минимальное число членов разложения, которое позволяет искать стационарную точку 2. Именно из этого простого факта следует исключительное значение, которое имеют квадратичные функционалы. Задачей определения конкретного вида  вектора   и матрицы , обычно занимаются более конкретные разделы науки: механика, физика, экономика и так далее.

Мы же будем рассматривать, пока общие вопросы. Обычно мы знаем из этих дисциплин как построить функционал в общем виде, но "дьявол прячется в мелочах" и обычно какая-то информация о некоторых переменных нам известна. То есть искать стационарную точку мы должны уже для другого функционала, который учитывает уже известную нам информацию. Аналогичная ситуация может возникнуть и в системах линейных (для соблюдения аналогии)  уравнений: мы знаем как получить систему для общего случая, и ее надо привязать к известной информации о части неизвестных. В этой статье рассматривается, как действовать в том случае если дополнительная информация может быть выражена в виде линейных равенств, в которые входят неизвестные более общей постановки и предлагается эффективный алгоритм, позволяющий модифицировать систему линейных уравнений с учетом имеющейся информации. Все это конечно теория, но как говорил кто-то их великих: "нет ничего практичнее хорошей теории". А теперь к делу...   

Квадратичный функционал в однородных координатах

Напомним теорему из анализа: Пусть A - положительно определенная матрица, если   имеет решение  , то это решение дает минимум функционалу  справедливо и обратное: если существует    ,из области определения матрицы  A, который дает функционалу J(u) минимальное значение, то  является решением уравнения  .

Отметим тонкость, эта теорема справедлива, если A -положительно определенная матрица, в противном случае, если всего лишь симметричная, речь должна идти о стационарном, а не о минимальном значении.

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

Для начала убедимся, что функционал можно записать в однородных координатах     

               => {{0}}

Видим, что разность двух записей функционала равна нулю и проверим, что производная от функционала, если приравнять ее нулю даст нам линейное уравнение, а, приравняв градиент функционала нулю, получим систему уравнений для определения стационарной точки.  При этом под A можем понимать квадратную матрицу, а под u и f векторы согласованных с ней размерностей, в нужных  по контексту местах это строки или столбцы.  
  =>      

Проверим это еще в покомпонентном виде:

=>
{{0}} 

Видим, что выражения для функционала совпадают, проверим и для производных:

=>{{0}}  {{0}}

Приравняв полученные выражения к нулю, получим  необходимую систему линейных уравнений. Из полученных выражений также видно, что требование симметричности матрицы принципиально для справедливости теоремы.

Покажем на примере, что положительная определенность матрицы принципиальна для существования минимума:   

Для матрицы, не обладающей диагональным преобладанием,  точнее не являющейся положительно определенной, хотя и симметричной, это стационарная точка, так как функционал не имеет минимума или максимума:     

Дополнительные условия на переменные

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

Пусть нам что-то известно об одной из переменных, например -  v.  В рамках линейных функций наиболее общая ситуация v=a u +b w +c,  приравнивая  a и(или) b  и(или) c  нулю, можем получать частные случаи. В матричном виде это можно записать:

, а еще удобнее в виде произведения двух матриц  . Проверим это:

 =>

Свойства матрицы S

Если обратить матрицу S, то можно обнаружить интересное свойство,  если оно подтвердится для любых связей, то можно легко вычислять обратную матрицу по правилу  :

 => => => True

=>=>True.

Этот результат важен в теоретическом аспекте, так как доказывает существование обратной матрицы для S. У Воеводина и Кузнецова приведен следующий результат: "12.61. Если A>0  (?,<,?) , то для любой невырожденной матрицы B имеем ”.

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

Свойства матрицы G

Для дальнейшего важны три свойства матрицы G:

1.    При умножении любой матрицы на G слева (G M) происходит обнуление строк матрицы M, соответствующих нулевым значениям на диагонали G;

2.   При умножении любой матрицы на G справа ( M G) происходит обнуление столбцов матрицы M, соответствующих нулевым значениям на диагонали G;

3. Для любой  матрицы  справедлива формула G.M==G.(E - G + G.M).

Покажем справедливость этих утверждений на модельных матрицах.    

=>=> True

Это было небольшое полезное отвлечение. Вернемся к функционалу и подставим в него матрицу  связей  (или условий - это синонимы).

Трансформация функционала с учетом матрицы связей

=>

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

Выручает в этой ситуации третье свойство матрицы G:     G.M==G.(E - G + G.M), которое позволяет записать матрицу A в следующем виде:

=>   

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

 =>=>0

Алгоритм учета связи

Результаты, полученные в предыдущем пункте можно получить и чисто алгебраическим путем. Систему уравнений AU=B запишем в однородных координатах в виде однородного векторного уравнения .Умножим его слева на некоторую матрицу  и вектор неизвестных заменим по правилу , тогда получим новое уравнение  из решения которого легко найти U. Если М будет удовлетворять условиям для матрицы связей, то есть в ней можно выделить две группы неизвестных - независимых и зависящих от  независимых, то матрицу M можно разложить на две M=S G. Если матрица А симметричная, то результаты полностью совпадут с полученными ранее. 

Другим путем такие же соотношения можно получить, если функционал рассматривать как функцию от неизвестных, а зависимые неизвестные рассматривать как функции от независимых переменных. Применяя правило дифференцирования сложной функции, в случае линейной зависимости получим такие же соотношения, но в случае нелинейных связей уже можем выделить нелинейные члены.

 Для практического использования технологии учета связей важно изучить свойства матрицы    для поиска эффективного алгоритма, так как эти матрицы в МКЭ, сильно разрежены и большие и прямое использование матричных произведений приведет к затратам времени соизмеримым с разложением матрицы, не говоря уж о дополнительных затратах памяти. Наиболее интересен случай симметричных матриц, поэтому будем изучать его.

Тривиально, но убедимся в тождественном отображении:

Все в порядке и посмотрим, что происходит, если условие - равенство двух переменных, а G=E:

Отметим, что пока у нас IdentityMatrix[Length[G]]-G тождественный нуль, а  G единичная матрица и для того, чтобы суть дела была более наглядна, опустим их:

=>=> 

Что получилось? 3 строка у нас зависит от 1, то есть зависимый столбец прибавили к независимому столбцу и получили матрицу Am1. Теперь подействуем слева транспонированной  матрицей S на полученную матрицу:

=>

=>=>True

Ну а тут что произошло? "Элементарно Ватсон" сказал бы Шерлок Холмс -" теперь уже  строку, соответствующую  зависимой переменной, прибавили к независимой и получили то же самое, что и при прямом вычислении .

 Интересно, почему Шерлок Холмс, метод, основанный на изучении деталей и применении аналитического продолженния для получения общей картины, называет дедуктивным? Наверное шутка Конан- Дойля... .

 Отметим, что при втором подходе, в операциях участвуют только строки и столбцы, входящие в связь. Это очень важное свойство для вычислений!

Можно проделать и наоборот, сначала строки, а потом столбцы:

=>

=>

=>True

Об этом можно было бы и догадаться, вспомнив об ассоциативном законе для произведения матриц, но лучше один раз увидеть, чем семь раз услышать.  Далее посмотрим уже для реальной матрицы G:

=>True

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

=>
=>

Видим, что при  все совпадает с предыдущим. Посмотрим, как это скажется на алгоритме:

=>=>

Зависимый столбец умножился  на множитель и как обычно добавился к независимому столбцу. И далее

=>

=>

Теперь уже зависимую строку умножили на множитель и добавили к независимой строке. Как обычно сравним результаты с прямым произведением матриц и убедимся, что структура действий и  основные свойства  алгоритма сохраняются:

=>True

Алгоритм конгруэнтного преобразования матриц

Замена переменных, которую  по существу мы осуществляли с помощью выражения SX, привела к матрице   . В линейной  алгебре есть для этого случая определение: матрица называется конгруэнтной  матрице A. В геометрии, например, фигуры конгруэнтны, если существуют движения, при которых фигуры можно совместить без деформаций фигур.  О интересных свойствах конгруэнтных преобразований можно почитать, например, у Хорна и Джонсона. Выше мы внимательно понаблюдали, что происходит с конгруэнтной матрицей, при специальном выборе матрицы S, и обнаружили, что происходит умножение какой- либо строки из A на коэффициент из S и сложение с другими строками или  строкой, аналогично и со столбцами.

Поставим задачу следующим образом: разработать алгоритм, который позволял бы вычислять  путем умножения строк из A на элементы из S и сложения с некоторыми  строками  и аналогичных действий со столбцами. При этом желательно, чтобы в процессе работы осуществлялось отображение A в себя, то есть, дополнительная память для    не выделялась. То есть, проведем геометрическую аналогию: фигура вращается на месте до тех пор, пока не сориентируется определенным образом.

Понятно, что для любых S это сделать не удастся, так как мы будем менять строки и столбцы, а, следовательно, и коэффициенты матрицы A.   Выше мы наблюдали, что, по крайней мере, для матриц со свойством   возможно удастся, а к ним мы сможем сводить условия для линейных связей. Нелишним будет проверить и другие "хорошие" матрицы.

 Автор немного поэкспериментировал с программированием и написал программу, которая позволяет это делать.

Алгоритм получился такой: последовательно проходим два раза матрицу S сначала по строкам, затем по столбцам, при каждом проходе действия аналогичные, поэтому опишем только строковые действия. Строку в A c номером i   умножаем на  и прибавляем к строке A с номером j и так по всем элементам S, а в случае  замещаем строку.

Еще предусмотрим анализ матрицы G для обнуления строк и столбцов и размещения 1 на диагонали для сохранения ранга матрицы. Для того, чтобы подчеркнуть характер взаимодействия с A доступ к ней, осуществляется через  функции  взятия строки (столбца) , записи строки (столбца) и функции обнуления. Это сделано для удобства переноса модельной функции на реальные задачи, где вспомогательные функции будут совсем другими даже для не разреженных матриц. Для сравнения используется функция, реализующая прямое произведение матриц. Описание описанием, но без живого кода все равно не обойтись. Вот он:

True

Этот результат говорит о том, что если матрица связей нижняя треугольная, то алгоритм работает аналогично прямому перемножению матриц. Из алгебры известно, что, используя линейные комбинации  и перестановки строк  систему линейных уравнений всегда можно представить в виде нижней треугольной, а учитывая, что наша матрица связана с условиями и обычно очень разрежена, то эта операция не займет много времени.   Найдем еще другие типы матриц, для которых алгоритм работает аналогично прямому произведению матриц:  

True

Из этих экспериментов можно сделать вывод, что в алгоритме допустимо использовать ненулевую  строку или столбец  с единицей на диагонали для исключения произвольной степени свободы. Это свойство может быть полезно при исключении внутренних мод или внутренних узлов при вычислении матрицы элемента, если брать в качестве строки соответствующую строку матрицы элемента,  деленную на коэффициент для нормировки. Последовательное применение этого алгоритма позволит исключить внутренние моды.

Рассмотрим еще один вариант матрицы S, для которого алгоритм работает  правильно. Пусть множество переменных разбивается на два - множество независимых переменных и множество зависимых, которые зависят только от независимых переменных. Заодно проверим справедливость ранее установленного свойства таких матриц связей  :

=>True

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

Для тех, кого не устраивает инженерный уровень строгости распространения на   , приведем абстрактно-математические рассуждения.  Квадратные матрицы порядка n образуют кольцо относительно операций сложения и умножения. Кольцо над   изоморфно кольцу над , другими словами эти кольца неразличимы по своим алгебраическим свойствам, и, следовательно, результаты, полученные на одном кольце в силу изоморфизма справедливы и для другого кольца. В силу произвольности n алгоритм будет работать на любых конечных матрицах.

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

Преобразование в симметричных матрицах

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

True

Проверим это на  матрицах из классов рассмотренных для несимметричного алгоритма и, учитывая важность этой технологии, проверим циклически по всем возможным вариантам. Сначала для матриц с одной неединичной  строкой или столбцом и единицами  на диагонали:

True

Теперь рассмотрим наиболее важный для метода конечных элементов вариант, когда зависимые неизвестные зависят только от независимых, а не друг от друга. Этого в непротиворечивых постановках задач всегда можно сделать. Заодно  еще раз проверим справедливость ранее установленного свойства таких матриц связей  :

True

Распространение на  осуществляется аналогично несимметричному случаю. Автор и люди в организациях, куда он передавал свои программы, пользуются похожим алгоритмом около 20 лет для решения трехмерных задач теории упругости, и проблем не возникало.

Прагматика

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

Предложенная технология, учета связей сокращает число неизвестных и упрощает ситуацию в отличие от альтернативного подхода, связанного с введением множителей Лагранжа. Фундаментальные свойства модифицированной матрицы не ухудшаются  по сравнению с исходной, по крайней мере положительно определенные матрицы таковыми и остаются.

Прямое перемножение матриц, для вычисления конгруэнтной матрицы, требует 0() умножений, ну в лучшем случае, при использовании алгоритма Штрассена  0(). А при использовании предлагаемой технологии  0(n* <среднее число элементов в связи > *<число связей>). При обычных постановках задач можно считать 0(n) и это при заполненных матрицах, а в случае применения технологий  для разреженных матриц ситуация еще улучшается.

Приведем некоторые примеры из полезных задач для расчета машин и механизмов. Будем использовать следующие обозначения:  - неизвестный вектор перемещений в i узле, S - фиксированная квадратная матрица,  - фиксированный вектор, связанный с i узлом.

Условия типа  позволяют моделировать ситуации, когда в некоторых узлах известны величины перемещений.

На условиях  , основан сам принцип  сшивки разбиений на подобласти в единую систему, обычно это делается автоматически при формировании  матрицы жесткости или эквивалентной ей, но можно и смоделировать отдельные детали и потребовать их соединения.  Если детали не соединены, то матрица будет распавшейся. Вот что при этом происходит:

В первом случае 3 неизвестная равняется 2, а во втором случае 4 неизвестная  равнялась 1. Понятно, что элементы матриц мы можем рассматривать как клетки и операции понимать как клеточные сложения. Прямое сложение элементов справедливо для любого количества  клеток, но матрица обязательно должна быть клеточно-диагональной или полученная из нее путем согласованной перестановки строк и столбцов. Или другими словами, граф матрицы должен быть распавшимся. Уж извините за “выкрутас”, но в разреженных матрицах без теории графов не разобраться, как это убедительно показали Джордж и Лю , а так же Писсанецки в своих прекрасных книгах.

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

Интересные возможности открывают и комбинации  модель скольжения, или   и вы легко поймете , что происходит при запресовке деталей с натягом d или оцените влияние сварки, если знаете усадку -d, о подобных подходах к сварке почитайте у Талыпова. О возможностях подобных комбинаций можно писать много, но я думаю читатель и сам знает или догадается о приложениях подобного подхода.

Условия с матрицей симметрии  , если i на одной поверхности разреза, а j сходственная точка на другой поверхности разреза циклически симметричного объекта, и S такова, что координаты точек одной поверхности разреза переводит в координаты другой поверхности,  при условии циклически симметричной нагрузки, позволяет легко изучать подобные объекты. А это практически все многопоточные  машины и механизмы: осевые турбомашины, электродвигатели, планетарные механизмы, муфты, шлицевые соединения и много других.

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

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

Заключение

Как это все использовать? Для метода конечных элементов желательна технология разреженных матриц. Как это должно выглядеть, можно почитать у Герберта Шилдта и в классике от Дональда Кнута.  Задача значительно упростится, если использовать STL - библиотеку стандартных шаблонов, которой мы обязаны Александру Степанову и менеджменту Дмитрия Ленкова, низкий им поклон за это.

Если Вы совсем ленивы или не любите заниматься изящным программированием ,то есть не любите "копаться в болтах и гайках и выковыривать орешки из шоколада" ,  то можете воспользоваться библиотеками для C++ от MATLAB или  MATCOM. Если соберетесь программировать, то удачи, а будут проблемы - пишите.

Литература

1. Воеводин В.В. Кузнецов Ю.А. Матрицы и вычисления: М. Наука.1984.
2. Галлагер Р. Метод конечных элементов.Основы:М. Мир.1984.
3. Джордж А.,Лю Дж. Численное решение больших разреженных систем уравнений. М.: Мир.1984.
4. Кнут Д. Искусство программирования для ЭВМ. Основные алгоритмы. М.: Мир.1976.
5. Писсанецки С.Технология разреженных матриц. М.: Мир,1988.
6. Талыпов Г.Б. Сварочные деформации и напряжения. Л.: Машиностроение.1973.
7. Хорн Р. Джонсон Ч. Матричный анализ. М.: Мир.1989.
8. Шилдт Г.Теория и практика C++. BHV-Санкт-Петербург.1996.  

Санкт-Петербург, май-июнь 2000

В начало