ДОКАЗАТЕЛЬСТВО СУЩЕСТВОВАНИЯ ПРЕДЕЛА ТОЧНОСТИ МКЭ ПРИ ПРИМЕНЕНИИ ОДИНАКОВЫХ ОДНОМЕРНЫХ ЛИНЕЙНЫХ КОНЕЧНЫХ ЭЛЕМЕНТОВ Текст научной статьи по специальности «Математика»

Вычислительные технологии, 2022, том 27, № 1, с. 39-51. © ФИЦ ИВТ, 2022 Computational Technologies, 2022, vol. 27, no. 1, pp. 39-51. © FRC ICT, 2022

ISSN 1560-7534 elSSN 2313-691X

ВЫЧИСЛИТЕЛЬНЫЕ ТЕХНОЛОГИИ

DOI: 10.25743/ICT.2022.27.1.004

Доказательство существования предела точности МКЭ при применении одинаковых одномерных линейных конечных элементов

П.М. Винник1'*, Т. В. Вин ник2, А. А. X Акимов3

балтийский государственный технический университет "ВОЕНМЕХ" им. Д. Ф. Устинова, 190005, Санкт-Петербург, Россия

2

университет), 190013, Санкт-Петербург, Россия

3

141701, Долгопрудный, Россия

*

Поступила 16 августа 2021 г., доработана 9 ноября 2021 г., принята в печать 10 декабря 2021 г.

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

Ключевые слова: метод конечных элементов, собственные числа матрицы жесткости, число обусловленности.

Цитирование: Винник П.М., Винник Т.В., Хакимов A.A. Доказательство существования предела точности МКЭ при применении одинаковых одномерных линейных конечных элементов. Вычислительные технологии. 2022; 27(1):39—51. DOI: 10.25743/ICT.2022.27.1.004.

Введение

Одним из основных численных методов решения задач механики сплошных сред является метод конечных элементов (МКЭ). В рамках этого метода задача сводится к решению системы линейных алгебраических уравнений (СЛАУ). Для построения матрицы СЛАУ по каждому конечному элементу (КЭ) е вычисляется матрица [К(е)], называемая локальной матрицей жесткости элемента е. Затем матрица СЛАУ [KG], называемая

Рис. 1. Одномерные конечные элементы Fig. 1. One-dimensional finite elements

глобальной матрицей жесткости, получается суммированием матриц одного порядка, каждая из которых строится вложением матрицы [К(е)] в нулевую матрицу большего порядка [1].

Матрица [KG] является симметричной [1], т.е. частным случаем нормальной [2], а также неотрицательно-определенной, т. е. ее собственные числа (СЧ) удовлетворяют неравенству

где N — размер матрицы, а число г зависит от размерности решаемой задачи (для одномерной задачи г =1).

Устойчивость решения СЛАУ относительно погрешностей исходных данных определяется числом обусловленности матрицы СЛАУ [3]. Спектральное число обусловленности С (А) для невырожденных нормальных матриц А равно отношению модуля наибольшего собственного числа к модулю наименьшего [2].

В настоящей статье построена аналитическая оценка спектрального числа обусловленности глобальной матрицы жесткости [КС] для одномерной задачи, решаемой МКЭ с одинаковыми одномерными линейными КЭ (рис. 1).

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

1. Глобальная матрица жесткости для цепочки одномерных КЭ и ее регуляризация

Для одномерного конечного элемента матрица [К(е)] имеет вид [6]

где величина Q = const зависит от длины конечного элемента и физических констант задачи и одинакова для всех КЭ. Для цепочки из (п — 1) одинаковых последовательных одномерных КЭ, из которых каждые два соседних соединены друг с другом в общем узле, глобальная матрица жесткости имеет вид

Ai ^ А2 ^ ... ^ А

•N-г > А N-г+1

Xn-i — А N — 0,

( 1

-1

[KG] = Q

0 0

V 0

-1 2

1

0 0 0

00 -1 0 21

0 0 0

0 0 0

0

0 0

0 -12 -10 0 0 -1 2 -1 0 0 0 -1 1

(1)

Для устранения сингулярности матрицы (1) нужно задать (в рассматриваемом одномерном случае) одно перемещение. Как показано в [1, 6] (вопросы технической организации самих вычислений МКЭ опущены), это соответствует вычеркиванию некоторой строки и столбца матрицы (1) с одним и тем же номером ] (считаем, что перемещение задано в узле с номером ]).

Так как в рассматриваемом одномерном случае можно задать перемещение в любом узле (т.е. вычеркнуть строку и столбец с любым одним и тем же номером), имеется всего п возможностей такого вычеркивания. На рис, 2 показаны пример вычеркивания и получившаяся после него матрица. Видно, что при вычеркивании любого диагонального элемента, кроме первого или последнего {] = 1 и ] = п), полученная матрица включает два блока одинаковой структуры, которые (с точностью до перенумерации строк и столбцов) имеют вид

1

1

[Ко] = Q

0

0 0

-1 2

1

0 0 0

00 -1 0 21

0 0 0

0 0 0

0

0 0

0 -12 -10 0 0 -12 -1 0 0 0 -1 2 /

(2)

т. е, отличаются от матрицы [КС] заменой последнего элемента диагонали на 2, Отметим, что ожидаемая невырожденность матрицы (2) очевидна, так как последователь-

1 -1 0 0 0 0 0

-1 2 -1 0 0 0 0

0 -1 2 -1 0 0 0

0 0 -1 2 -1 0 0

0 0 0 -1 2 -1 0

0 0 0 0 -1 2 -1

0 0 0 0 0 -1 1

1 -1 0 0 0 0

-1 2 0 0 0 0

0 0 2 -1 0 0

0 0 -1 2 -1 0

0 0 0 -1 2 -1

0 0 0 0 -1 1

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

Fig. 2. Deleting the diagonal element of the matrix (a) and the matrix obtained after such deletion (6)

a

ным прибавлением первой строки ко второй, получившейся второй — к третьей и т, д, имеем ёе^Хо] = 1.

В случае вычеркивания диагонального элемента матрицы [КС] с номер ом ] = 1 или ] = п полученная матрица состоит из одного блока вида (2),

Таким образом, если обозначить характеристический многочлен матрицы (2) размера к х к через Н^ (полагая Н0 = 1), то характеристический многочлен матрицы, получающейся при вычеркивании ^'-го диагонального элемента матрицы [КС] размера п х п, будет равен

— Hj—lHn—j. (3)

Так как различные многочлены Н^, как можно убедиться непосредственными вычислениями, имеют, вообще говоря, различные корни, регуляризация матрицы [КС] является неканонической (ее результат зависит от выбора задаваемого перемещения). Матрица (2) является частным случаем квазитеплицевых трехдиагональных матриц, Оценкам спектра и числа обусловленности квазитеплицевых трехдиагональных матриц посвящены работы [7-11], причем в работе [11] установлено асимптотическое поведение числа обусловленности ряда матриц указанного типа без аналитического вычисления их собственных чисел, и, в частности, для матрицы (2) установлено стремление числа обусловленности к бесконечности,

В [12] с целью устранения этой неканоничности регуляризации предложено определить спектральное число вырожденной матрицы [КС] указанного типа (нормальной и имеющей вещественные неотрицательные собственные числа) как

с ([КС]) = (4)

Лы-г

и рекомендовано использовать число С ([КС]) для оценки чисел обусловленности регу-ляризаций матрицы жесткости.

Таким образом, необходимо вычислить все собственные числа матриц вида (2) и вида (1),

2. Вывод рекуррентной формулы для характеристических многочленов матриц (1) и (2)

Рассмотрим матрицу

К

( -1 -1 -1 0 0 -1

0

0 0

0 0 0

0 -1 0

0 0

1

0 0 0

0 0 0

0

0 0

0 -1 0 -1 0 0 0 -1 0 -1 0 0 0 -1 * )

(5)

Если обозначить через I единичную матрицу соответствующего размера, то обе матрицы [KG]/Q - 21 и [K0]/Q - 21 будут иметь вид (5), только правый нижний элемент "*" будет равен (-1) для [KG]/Q - 21 и нулю для [K0]/Q - 21.

Вычислим характеристический многочлен матрицы К:

сЫ(х1 — К) = Д

1 + ж 1 0

0 0 0

1

х 1

0 0 0

0 1

х

0 0 1

0 0 0

0 0 0

1 0 0

х 1 0

0 0 0

1

х 1

0 0 0

0 1

х + *

Индекс п указывает здесь на размер определителя.

Считаем п > 3. Заметим, что Дп можно представить в виде

Д*

^о(ж) Со (ж) 0 0 0

1 X 10 0 0 0

0 1 X 1 0 0 0

0 0 ... 0 1 X 1 0

0 0 . . . 0 0 1 X 1

0 0 . . . 0 0 0 1 X + 1

0

0

где Ро(х) = ж+1, С0(ж) = 1. Вычитая из первой строки Дп вторую строку, домноженную па Ро(х), получаем

Д.

0 Со(х) — хРо(х) -^о(ж) 0

1 0

0 0 0

X 1

0 0 0

1

X

0 1

0 0 0

0 0 0

1 0

0 1 X 1

0 0 1 X

0 0 0 1

ад 0 0 0

X 1 0 0

1 X 1 0

0 0 0

0 1

X + 1

Раскладываем по первому столбцу и меняем знаки в первой строке

0 0 0

0 0 0

0 0. . . 0 1 X 1 0

0 0. . . 0 0 1 X 1

0 0. . . 0 0 0 1 X + 1

п— 1

Обозначаем Р1(х) = хРо(х) — Со(х), С1(х) = Ро(х). Продолжая точно так же, получаем

(х + *)Еп-2(х) — Сп-2(х),

Д,г

Рп-2(х) Сп-2(х) 1 х + *

п

71

П

где Сп(х) = Рп-1(х) и Рп(х) = хРп-1(х) - Сп-Ъ откуда

Ап = (х + *)Рп-2(х) - Рп-з(х).

Или

Лп =(х + *)Рп-2(х) - Рп-з(х) =

= (х + *) (хРп-з(х) - Рп-4) - (хРп-4 (х) - Рп-5) =

=х ((ж + *)Рп-3(х) - Рп-4(х)) - ((х + *)Рп-4(х) - Гп-5(х)).

Заменяя содержимое последних двух скобок в (7) по (6), устанавливаем

Лп = хЛп_ 1 - Лп_2.

(6)

(7)

Для матрицы [KG]/Q - 21 имеем Л2([КС]/Я - 21)

1 + х 1

1 1 + X

х2 + 2х = (х + 2)х,

(8)

(9)

Л3([КС]^ - 21)

1 + ж 1 0 1 ж 1 0 11 + ж

= ж3 + 2х2 - х - 2 = (х + 2)(х2 - 1). (10)

Из (8)—(10) следует, что при любом п ^ 2 многочлен Лn([KG]/Q - 21) делится па скобку (х + 2). В многочленах Лn([KG]/Q - 21) обозначим множитель при (х + 2) через Рп(х):

Лп([КС]/Я - 21) = (х + 2)Рп(х), (11)

т.е. Р2(х) = ж, Р3(х) = х2 - 1. Тогда из равенства (8) получаем

(х + 2)Рп = х(х + 2)Рп-1 - (х + 2)Рп-2, Рп хРп—1 Рп-2.

Для матрицы [K0]/Q - 21, обозначая Нп = Лn([K0]/Q - 21), имеем

(12)

Н2

1 + х 1

1 X

х + х — 1, Н3

1 + х 10 1 х 1 0 1 ж

= х3 + X2 - 2х- 1.

3. Конечные формулы для характеристических многочленов

Утверждение 1.

1, Если п = 2М ^ 2 — четное, то

N-1

Рп(х) = у(-1г-1-к с2^х2к+1.

к=0

(13)

2. Если п = 2Ы- 1 > 3

м-1

Рп(х) = ^(-1)М-1-к с2к

N+к-1

X

2к

(14)

к=0

нечетное, то

Доказательство выполним по индукции. Для п = 2 и п = 3 многочлены Рп(х) выписаны выше. Для четного п = 2N + 2 вычислим многочлен Рп(х) по формуле (12)

с учетом индукционного предположения:

N

Р2М+2(х) = х £ (_1^-кС

^кгЯк ^2к

Кк=0 ^ 1

х

2N+1

+ -кс

N+к

_2к+1

\ N-1

И _£<

/ к=0

(_1)N -1-к£ 2к+1х2к+1

^ 1

N+к

х

N-к ^2к+1 х2к+1 =

к=0

к=0

^ 1

„2 N +1

+ £(-1)

N-к„.2к+1{ + Ц!

х

2N+1

к=0 N-1

+ Е (_1)

(

( N + к)!

х~'" ( 1Г { рк)!^ _ к)! + (2к + 1)!(N-к _ 1)!

)

N-к„ 2к+1 ^ + к)!

х

к=0

(2 к)!( N _к_ 1)!\N _ к 2к + 1

11

+

Ш+1 + 1)N-кх2к+1 (N + к+ 1)! = ^ ( 1)N-кг2к+1 х2к+1

х + 2^( 1) х (2& + 1)!( N — &)! -¿—-А 1) ^N+1+кX

к=0

(2 к + 1)!(N _к)!

к=0

Аналогично для нечетного п = 2N + 1:

(N-1 \ N-1

£(_1^-1-кс^+к1х2к+11 _ £ к=0 ) к=0

N-2

„2 N , 1 — к

х

(_ -1-кС2к х ( 1) ^N-1+^

N-1

+ ^(_1^-1-кс2к+кх2к+2 _ £(_1Г-1-кс2к-1+кх2к _ (_1)

к=1

к=0

^ 1

Делая в первой сумме сдвиг суммирования т = к + 1, получаем

N-1 N-1

„2N I ( -тгу2т-1 „,2т , / 1 NN-к^2к

СN+т-1X ^ 2^(_1) С т=1 к=1

^ , 2к / (N + /С _ 1)! , (N + /С _ 1)!

х

х

х

+ (-1)

+ £(_!) "-кх* (

к=1 ^ N-1

N + Е(_1)

+

N-к х2к

(2 & _ 1)!(N _ к)! (2^)!^ _ 1) ( N + к _ 1)!

^-1+к

+ к _ N _1

11

х2к _ (_1^-1 =

у) _ (_1Г-1

к=1 ^ 1

(2 к _ 1)!(N _к _ 1)!\N _к 2к

N

)! & + 2&)

+ ^ _ (_ 1)

N-1

х

2 N

+ £(_ 1 Г

к=1

-кх2к (N + к)\ _ -1 х (2^)!^ _ к)! ( )

£(_1)ЛТ-кС

^кгч2к ^2к

N+к

х

к=0

□

Замечание. С учетом того, что для четных п = 2N выполнено N _ 1, а для нечетных п = 2N _ 1 выполнено

п 1 "2 N _ 1

2 2

п 1 "2 N _ 2"

2 2

N _ 1, сделав

замену индекса суммирования т = N _ 1 _ к ж учитывая равенство С(м = С^-е (для произвольных целых М ^ 0, 0 ^ I ^ М), для многочленов Рп(х) из утверждения 1 формулы (13), (14) можно переписать в едином виде:

' Ж

Рп (х) = £ (_1)тСт-т-1хП-2т-1. (15)

т=0

Утверждение 2.

Для п ^ 2 п

Нп = ^(-1)п-кс:-кк (2 + х)к. (16)

1п

к=0

Доказательство выполним по индукции. Для п = 2 и п = 3 многочлены Нп(х) выписаны выше. Заменим у = х + 2.

п п—1

Нп+Г=] хнп - Нп-1 = (у - 2) ^(-1)п-к^ук - ^(-1)п-1-кСп-_\+ккук =

к=0 к=0 п— 1

Т,(-1)П-кУк+1 - 2^(-1)П—кук - Т,(-1)П— 1—кСп—Иук к=0 к=0 к=0 п+1 п—1 п—1

= Т.(-1)^+1С^1у£-2 Т.(-1)П—кУк-2(-1)П-2уп-Т,(-1)П—1—кУк- (-1)п—1 = 1=1 к=1 к=1

п—1 п—1 п—1

= уп+1 - (2п- 1)уп+^(-^С^—у1 - 2 ^(-1)п— ук - ^(-^С1—— Ук+ 1=1 к=1 к=1 п— 1

+-1)<+ -2уп = уп+1 - (2п+ 1)уп + ^(-1)п+1—кук (С^ + 2С— -С—*) + (-1)п*1.

к=1

Отдельно преобразуем скобку в сумме:

гп—к+1 + 2Гп—к_ Гп—1 — к = (п + к - 1)! + 2(п + к)! (п + к - 1)! =

^п+к—1 + ^п+к ^П—1+к (п - к + 1у(2к - 2)! +(п - к)\(2к)\ (п - к - 1)\(2к)\ = = 2к)\ -)!к)\ {(2к - 1)2к + 2(п + к)(п +1 - к) - (п - к)(п - к + =

-—т-—(4к2-2к + 2п2 + 2кп+2п+2к - 2пк-2к2 - п2 + 2пк - к2 - п+к]

(2к)\(п + 1 - к)\\ /

(п + к 1)! I л12 , 0_2 , , о„ , т„2_ п2 _

(к2 + п2 + 2кп + п + к) = 7777^+^-^^ттт ((к + п)2 + п + к] V ) (2к)\(п +1 - у )

т+1—к

(2к)\(п +1 - к)\ ~п+к+1'

(2к)\(п +1 - к)\\ ) (2к)\(п +1 - к)

(п + к + = +1—к

1 \ . / ч . ^^^и^Л

Таким образом,

п+1

Нп+1 = ^(-1)п+1—к СппХИ (2 + х)к. □

к=0

4. Собственные числа глобальной матрицы жесткости [КО] и ее регуляризаций

Вычислим СЧ для матрицы [КС]. Из (11) и (15)

/М \

V1 {1)тПт

/ у V V Куп—т—1"

Лп([КС]/Я - 2 • I) = (х + 2)

п—2т—1

'п—т— 1 X

т=0

(17)

Для синуса кратного угла есть формула (см, формулу (11) в [13] ' ' [ ^ ] ' ' '

sin(n a) = sin а ^^ (—1)mC%- m-1 (2 cosa)

\n-2m-1

(18)

m=0

Из сравнения формул (17) и (18) следует, что все корни многочлена Дп — это

nk

х = -2, Xk = 2cos—, k = 1,

n

,n — 1.

Тогда собственными числами матрицы [KG] (1) являются

n k \

X = 0, Xk = Q( 2 + 2 cos — j

2

k = 1 , . . . , n 1 .

(19)

Для регуляризации матрицы [КС], полученной вычеркиванием j-тo диагонального элемента, рассмотрим сумму

2 - i 2N+1 \ (a sin(( N +1)a) sin( Na) V 2 V ^2

Í2N + 1 \ (a \ . (2N +1 \ -a cos — sin -a

\ 2 J \2J _ \ 2 J

sin a

sin a

sin a

С другой стороны, по формуле (18) получаем

2

sin ( 2)

(20)

sin(( N + 1)a) + sin( Na)

sin a

sin a

[ ^^ ]

Ё (— 1)mCmN+1-m-1 (2cos(^

a\ \ 2N+1-2m-1 2.

M D)

m=0

a 2 N-2 m N

Y,(—1)mC2N-m (2 cos (= £(—1)mC2mN-m (2 + 2 cos (a))

N m

m=0

2

k = N — m m = N — k

N

m=0

J2(—1)N-kcN-k (2 + 2 cos (a))k

(21) (22) (23)

k=0

Сравнивая (16), (20) и (21)—(23), устанавливаем, что корнями многочлена Нп являются

2n k

X = 2 cos

По (3) и (24) корни многочлена h

2 N + 1

k = 1,

N.

nj

X = 2 + 2 cos

2nk

/ ч , k = 1,... ,j — 1,

2(j — 1) + 1, , ,J , 2ш

х = 2 + 2 cos —-„ , г = 1,

(24)

(25)

2( n— ) + 1

, n —

5. Числа обусловленности регуляризаций глобальной матрицы жесткости

Если обозначить mj = тах(п _ ],] _ 1), то минимальным и максимальным корнями хттхтах^ много члена по (25) будут

X

min j

2 + 2 cos ■

2nm,-

X

2 + 2 cos ■

2n

2mj + 1' 2mj + 1'

поэтому спектральное число обусловленности матрицы, получающейся при вычеркивании j-тo диагонального элемента матрицы [ КС], будет равно

Cj =

х

1 + cos

2ж

шах,у

2mj + 1

X

mm, j

1 + cos j

2mj + 1

Из определения числа т^ и монотонного возраетания числа С^ при т^ € [1; следует, что для фиксированного размера матрицы п х п

2к

1 + cos

с,

regular

min Cj

2N + 1

16N 2

1 + cos

2nN

N n

2

где N

n

L2J

2М + 1

которое для нечетных п = 2М + 1 будет достигаться при ] = N + 1, а для четных п = 2М — при ] = N и ] = N +1,

В свою очередь, число обусловленности глобальной матрицы жесткости, определенное по формуле (4), по (19) будет равно

С ([KG])

2 + 2 cos — п

2 + 2 cos

ж(п — 1)

п

1 + cos —

_п_

•к

4п2 2'

п

1 — cos - п

п

Сравнивая С ([КС]), Сгеди1аг и произ вольные Сустанавливаем, что

С ([KG]) < С regular < min

1 ^ j ^ n,

j = N, j = N + 1

Gregul ar < С ([KG]) < min Cj

1 ^ j ^ n, j = N + 1

Cj

при п

2N,

при n = 2N + 1.

Причем асимптотически (как при п = 2М, так и при п = 2М + 1) они эквивалентны

16М 2

С ([КС])

Сп

regular „т 2

N^+ж N^+ж Л2

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

Из стремления спектральных чисел обусловленности к бесконечности и факта разбалтывания решения СЛАУ [3] при большом числе обусловленности следует, что существует предел точности решения задачи при помощи одинаковых одномерных конечных элементов.

Заметим, что матрицы (1) и (2) от матриц вида

а Ь 0 0 0 ... 0 0

(26)

b а b 0 0 0 0

0 b a b 0 0 0

0 0 . . . 0 b а b 0

0 0 . . . 0 0 b a b

0 0 . . . 0 0 0 b a

имеющих важное значение в методе сеток [14], существенно отличаются тем, что крайние элементы главной диагонали матрицы (26) совпадают с остальными элементами диагонали. Полный перечень собственных чисел и векторов для матриц (26) известен.

Выводы

1, Аналитически вычислены все собственные числа глобальной матрицы жесткости для одномерной задачи, решаемой МКЭ с одинаковыми одномерными линейными конечными элементами,

2, Вычислено наименьшее возможное спектральное число обусловленности всех ре-гуляризаций вырожденной глобальной матрицы жесткости,

3, Из стремления к бесконечности числа обусловленности (с ростом количества КЭ) следует существование предела точности решения для задач указанного класса.

Список литературы

[1] Зенкевич О. Метод конечных элементов в технике. М.: Мир; 1975: 542.

[2] Хорн Р., Джонсон Ч. Матричный анализ. М.: Мир; 1989: 656.

[3] Фаддеев Д.К., Фаддеева В.Н. Вычислительные методы линейной алгебры. СПб.: Лань; 2009: 736.

[4] Самарский A.A., Вабищевич П.Н. Численные методы решения задач конвекции-диффузии. М.: Книжный дом "ЛИБРОКОМ"; 2015: 248.

[5] Olshanskii М.А., Reusken A. Convergence analysis of a multigrid method for a convection-dominated model problem. SIAM Journal on Numerical Analysis. 2004; (42):1261—1291. DOI: 10.1137/S0036142902418679.

[6] Сегерлинд Л. Применение метода конечных элементов. М.: Мир; 1979: 392.

[7] Veerman J.J.P., Hammondy D.K., Baldiviesoz P.E. Spectra of tridiagonal matrices. Available at: https://arxiv.org/pdf/1801.04977.pdf (accessed November 02, 2021).

[8] Noschese S., Pasquini L., Reichel L. Tridiagonal Toeplitz matrices: Properties and novel applications. Available at: https://www.math.kent.edu/~reichel/publications/toep3.pdf

(accessed November 02, 2021).

[9] Alvarez-Nodarse R., Petronilho J., Quinteroc N.R. Spectral properties of certain tridiagonal matrices. Linear Algebra and Its Applications. 2012; (436):682-698.

[10] Годунов С.К., Малышев А.Н. О специальном базисе из приближенных собственных векторов с локализованными носителями для изолированного узкого кластера собственных значений симметричной трехдиагональной матрицы. Журнал вычислительной математики и математической физики. 2008; 48(7) :1156—1166.

[11] Батагаева Т.А., Казаков А.Л. Аналитическое и численное исследование числа обусловленности квазитеплицевой трехдиагональной матрицы с неограниченной размерностью. Вестник ИрГТУ. 2015; 4(99):122—130. Адрес доступа: https://www.elibrary.ru/item. asp?id=23411978.

[12] Иванов K.M., Винник П.М., Иванов В.Н. Численное моделирование разделительных процессов обработки давлением. Вестник Самарского государственного аэрокосмического университета. 2012; 2(33):192—198.

[13] Weisstein E.W. Multiple-angle formulas. Available at: https://mathworld.wolfram.com/ Multiple-AngleFormulas.html (accessed July 08, 2019).

[14] Даугавет И.К. Теория приближенных методов. Линейные уравнения. СПб.: БХВ-Петербург; 2006: 288.

Вычислительные технологии, 2022, том 27, № 1, с. 39-51. © ФИЦ ИВТ, 2022 ISSN 1560-7534

Computational Technologies, 2022, vol. 27, no. 1, pp. 39-51. © FRC ICT, 2022 elSSN 2313-691X

COMPUTATIONAL TECHNOLOGIES

DOI: 10.25743/ICT.2022.27.1.004

Finite element method: Proof of the existence for the accuracy limit when applying one-dimensional linear finite elements

VlNNIK Petr M.1'*, Vinnik Tatyana V.2, Khakimov Andrey A.3

1Baltic State Technical University "VOENMEH" named after D. F. Ustinov, 190005, Saint-Petersburg, Russia

2St. Petersburg State Technological Institute, 190013, Saint-Petersburg, Russia

3

*

Received August 16, 2021, revised November 9, 2021, accepted December 10, 2021.

Abstract

The finite element method is widely used to solve partial differential equations. The process of solving leads to a sparse system of linear algebraic equations with a singular matrix. It is known that the accuracy of solving systems of linear equations depends on the condition number of the system matrix. Previously, it was proposed to estimate the condition number for singular matrices using their non-zero eigenvalues, as for regular matrices. In the article, all eigenvalues of the stiffness matrix that arises when solving a one-dimensional problem using the same linear finite elements and all possible regularizations of the stiffness matrix — matrices obtained by deleting the row and column containing the selected diagonal element are analytically calculated. It is shown that for a stiffness matrix with an odd number of rows, the smallest of the condition numbers of regularization matrices is less than the number previously proposed as the condition number of a singular matrix. However, they are asymptotically equal. An estimate for the condition number of the matrix and its asvmptotics are given. It is shown that the generally accepted method to increase the accuracy of calculations by the finite element method, which consists in element refinement, has a theoretical accuracy limit.

Keywords: finite element method, eigenvalues of a stiffness matrix, condition number. Citation: Vinnik P.M., Vinnik T.V., Khakimov A.A. Finite element method: Proof of the existence for the accuracy limit when applying one-dimensional linear finite elements. Computational Technologies. 2022; 27(1):39-51. D01:10.25743/ICT.2022.27.1.004. (In Russ.)

References

1. Zienkiewicz O.C. The finite element method in engineering science, Second edition. McGraw-Hill-London; 1971: 521.

2. Horn R.A., Johnson C.R. Matrix analysis. Cambridge University Press; 2012: 662.

3. Faddeev D.K., Faddeeva V.N. Computational methods in linear algebra. Freeman, San Francisco, Calif.; 1963: 621.

4. Samarskii A.A., Vabishchevich P.N. Chislennye metody resheniya zadach konvektsii-diffuzii

[Numerical methods for solving convection-diffusion problems]. Moscow: Knizhnyy dom "LIBROKOM"; 2015: 248. (In Russ.)

5. Olshanskii M.A., Reusken A. Convergence analysis of a multigrid method for a convection-dominated model problem. SIAM Journal on Numerical Analysis. 2004; (42):1261-1291. D01:10.1137/S0036142902418679.

6. Segerlind L.J. Applied finite element analysis. N.Y., London, Sydney, Toronto: John Wiley & Sons, Inc.; 1985: 448.

7. Veerman J.J.P., Hammondy D.K., Baldiviesoz P.E. Spectra of tridiagonal matrices. Available at: https://arxiv.org/pdf/1801.04977.pdf (accessed November 02, 2021).

8. Noschese S., Pasquini L., Reichel L. Tridiagonal Toeplitz matrices: Properties and novel applications. Available at: https://www.math.kent.edu/~reichel/publications/toep3.pdf (accessed November 02, 2021).

9. Alvarez-Nodarse R., Petronilho J., Quinteroc N.R. Spectral properties of certain tridiagonal matrices. Linear Algebra and Its Applications. 2012; (436):682-698.

10. Godunov S.K., Malyshev A.N. On a special basis of approximate eigenvectors with local supports for an isolated narrow cluster of eigenvalues of a symmetric tridiagonal matrix. Computational Mathematics and Mathematical Physics. 2008; (48): 1089-1099. D0110.1134/S0965542508070026.

11. Batagaeva T.A., Kazakov A.L. Analytical and numerical study of condition number of quasi-Toeplitz tridiagonal matrix with unlimited dimension. Vestnik Irkutskogo Gosudarstvennogo Tekhnicheskogo Universiteta. 2015; 4(99):122-130. Available at: https://www.elibrary.ru/item. asp?id=23411978. (In Russ.)

12. Ivanov K.M., Vinnik P.M., Ivanov V.N. Numerical simulation of separation process in mechanical working. Vestnik of Samara University. Aerospace and Mechanical Engineering. 2012; 11(2):192—199. (In Russ.)

13. Weisstein E.W. Multiple-angle formulas. Available at: https://mathworld.wolfram.com /Multiple-AngleFormulas.html (accessed July 08, 2019).

14. Daugavet I.K. Teoriya priblizhennykh metodov. Lineynye uravneniya [Approximation Theory. Linear equations]. SPb.: BHV-Peterburg; 2006: 288. (In Russ.)