Работа некоторых методов решения СЛАУ в МКЭ

[vitaly333]

Решил протестировать работу различных методов решения систем линейных уравнений.
Занимаюсь МКЭ , поэтому матрица всегда получается симметричной и положительно определенной (стержневые системы). Поэтому для решения СЛАУ чаще всего приходится пользоваться прямым методом Холесского(квадратных корней) и итерационным методом
Сопряженных Градиентов. На некоторых задачах различные методы дают различное друг от друга решение. Приведу вам небольшой пример:

Тип конструкции: плоская стержневая рама.
см.изображение
В узлах 0 и 2 есть граничные условия 1-ого рода - запрет перемещений по оси Y. Т.е
в узле 0 и 2 перемещения по оси Y равны 0. В узле 1 приложена нагрузка 12 Кн.

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

K1 = 
1.306165061   0   0.238778392   -1.306165061   1.112148398   0.238778392   0   0   0
0   1   0   0   0   0   0   0   0
0.238778392   0   0.955462467   -0.238778392   -0.245320266   0.477731234   0   0   0
-1.306165061   0   -0.238778392   1.981800485   -0.253351604   -0.422206912   -0.675635424   0   -0.18342852
1.112148398   0   -0.245320266   -0.253351604   2.581973396   -0.119320301   -0.858796794   0   0.125999965
0.238778392   0   0.477731234   -0.422206912   -0.119320301   1.72580511   0.18342852   0   0.385171321
0   0   0   -0.675635424   -0.858796794   0.18342852   0.675635424   0   0.18342852
0   0   0   0   0   0   0   1   0
0   0   0   -0.18342852   0.125999965   0.385171321   0.18342852   0   0.770342643

Вектор узловых нагрузок:

b1 = 
0
0
0
0
12
0
0
0
0

Решаем эту систему, воспользовавшись функцией lsolve из MatCAD:

x11 = lsolve(K1,b1)

x11 =
-81.084
40.16
-1.362
84.392
2.863
116.913
-43.398

Проверим получившийся результат – вычислим вектор невязки:

r11 = K1*x11 - b1
r11 =
-1.745E-15
4.943E-15
1.844E-15
-1.243E-14
0
0
-2.089E-15

Близко к нулю, но в идеале хотелось бы получить чтобы все элементы вектора невязки были равны точно нулю.

Решаем эту систему методом Холесского:

x12 = 

7.574879788886595
40.160100751029
87.29667559086803
84.39245474872334
2.8625983863313795
205.57242051314614
-43.39780790330944

r12 = K1*x12 - b1

r12 = 
-7.373E-15
5.223E-15
3.634E-15
-2.309E-14
3.178E-15
-6.035E-15
-7.248E-15

Решаем эту систему методом LU - разложения:

x16 = 
-21.787510139058092
40.16010075102904
57.93428566292343
84.39245474872344
2.8625983863313835
176.21003058520168
-43.397807903309484

r16 = K1*x16 - b1
r16 = 
6.508E-15
3.741E-15
-1.048E-14
3.553E-15
9.307E-15
-2.537E-15
0

Решаем эту систему методом сопряженных градиентов с точностью 1e -10:

x13 = 
-92.57311217541384
40.16010075102908
-12.851316373432214
84.39245474872354
2.8625983863313724
105.42442854884612
-43.39780790330953

r13 = K1*x13 - b1

Решаем эту систему методом сопряженных градиентов с точностью 1e -17:

-124.79332268534779
40.160100751029056
-45.07152688336625
84.39245474872347
2.862598386331376
73.20421803891205
-43.3978079033095

r13 = 
-1.392E-14
-1.705E-14
5.428E-14
5.507E-14
-2.826E-14
-3.902E-14
-8.573E-15

Воспользуемся функцией pcg(Точность 1e-10) из MatLab:

x15 = 
-92.5731
40.1601
-12.8513
84.3925
2.8626
105.4244
-43.3978

r15 = K1*x15 - b1

r15 = 
-92.5731
40.1601
-12.8513
84.3925
2.8626
105.4244
-43.3978 

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

Кроме того решения получаются также различными при выборе разных вариантов учёта граничных условий. Так как в выше приведенном примере перемещения по направлению закрепленных связей заведомо равны нулю(x[2]=0,x[8]=0; i=2,i=8) то учёт граничных условий производится путём вычеркивания i-ой строки и i – ого столбца матрицы жесткости (а также i-тых элементов векторов перемещений и нагрузок ).
Поэтому удаляем 2 и 8 – ую строку и 2 и 8 –ой столбец из матрицы жесткости а также
2 –ой и 8-ой элементы из векторов нагрузки и перемещений. Порядок матрицы уменьшается с 9 – ти до 7-ми Такой способ учета граничных условий подходит только для ручных расчетов небольших матриц а удалять строки и столбцы из статических больших матриц, изменяя этим самым их размерность очень дорогостоящая операция на ЭВМ , поэтому на практике i- ую строку и i-ый столбец делают нулевыми , а диагональный коэффициент K[i][i] = 1. Также в i-ую позицию вектора нагрузки засылают 0. При таком способе учёта граничных условий не изменяется размерность матрицы , никакие строки и столбцы не исключаются.

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

Матрица жесткости 9 x 9:

K = 
1.306165061   0   0.238778392   -1.306165061   1.112148398   0.238778392   0   0   0
0   1   0   0   0   0   0   0   0
0.238778392   0   0.955462467   -0.238778392   -0.245320266   0.477731234   0   0   0
-1.306165061   0   -0.238778392   1.981800485   -0.253351604   -0.422206912   -0.675635424   0   -0.18342852
1.112148398   0   -0.245320266   -0.253351604   2.581973396   -0.119320301   -0.858796794   0   0.125999965
0.238778392   0   0.477731234   -0.422206912   -0.119320301   1.72580511   0.18342852   0   0.385171321
0   0   0   -0.675635424   -0.858796794   0.18342852   0.675635424   0   0.18342852
0   0   0   0   0   0   0   1   0
0   0   0   -0.18342852   0.125999965   0.385171321   0.18342852   0   0.770342643

Вектор нагрузок:

b = 
0
0
0
0
12
0
0
0
0

Решаем функцией lsolve из MatCad:

x = lsolve(K,b)

x = 
1.26E+17
0
47.923
1.26E+17
98.239
3.348E-15
1.26E+17
0
-47.678

r = K*x - b

r = 
-4.697
0
-1.235
0.519
-6.546
-0.204
4.184
0
2.063

Как вы видите здесь вектор невязки получается далеко не нулевой и это говорит о том что это решение никуда не годиться.

Решаем методом Холесского:

x2 = 
-49.766265919852245
0.0
40.160100751028956
29.955529882129103
84.39245474872327
2.862598386331361
148.23127480440712
0.0
-43.39780790330938

r2 = K*x2 - b

r2 =
-5.126E-9
0
-5.82E-8
1.229E-9
-2.656E-8
-5.497E-8
3.896E-9
0
-5.109E-8

Решаем эту систему методом LU - разложения:

x6 = 
-12.612695546824247
0.0
40.160100751029006
67.10910025515724
84.39245474872341
2.862598386331378
185.38484517743547
0.0
-43.39780790330947

r6 = K*x6 -b

r6 = 
-5.126E-9
0
-5.82E-8
1.229E-9
-2.656E-8
-5.497E-8
3.896E-9
0
-5.109E-8

Решаем эту систему методом сопряженных градиентов с точностью 1e -10:

x3 = 
-92.57311217541357
0.0
40.16010075102897
-12.85131637343223
84.39245474872334
2.8625983863313635
105.42442854884582
0.0
-43.397807903309406

r3 = K*x3-b

r3 = 
-5.126E-9
0
-5.82E-8
1.229E-9
-2.656E-8
-5.497E-8
3.896E-9
0
-5.109E-8

Решаем эту систему методом сопряженных градиентов с точностью 1e -17:

x4 = 
-732.6769555414033
0.0
40.160100751029034
-652.9551597394219
84.39245474872342
2.862598386331401
-534.6794148171435
0.0
-43.39780790330951

r4 = K*x4 - b

r4 = 
-5.126E-9
0
-5.82E-8
1.229E-9
-2.656E-8
-5.497E-8
3.896E-9
0
-5.109E-8

PCG из Matlab (Точность 1e-10):

x5 = 
-92.5731
0
40.1601
-12.8513
84.3925
2.8626
105.4244
0
-43.3978

r5 = K*x5 - b

r5 =
4.504E-5
0
-1.211E-5
2.242E-5
1.517E-4
-9.227E-6
-6.746E-5
0
4.12E-6

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

[Zai]

Изложеные Вами материалы немного непонятны. Если Вы решаете плоскую задачу, то почему все же три степени свободы в точках. МКЭ решает только статически определимые\неопределимые задачи. У Вас механизм. Кроме вертикального закрепления в точках 0 и 2 , Вам нужно закрепить их горизонтально.

[vitaly333]

Цитата:

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

Потому что перемещения в узлах у плоской рамы могут быть по трем направлениям по x - от действия продольной силы, по y - от действия поперечной силы и по z (поворот) от действия момента.

Цитата:

У Вас механизм

Что вы подрузамеваете под этим словом.

Цитата:

Кроме вертикального закрепления в точках 0 и 2 , Вам нужно закрепить их горизонтально

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

[Zai]

Если третья степень свободы угол поворота, то у Вас статически неопределимая задача а не механизм. Конечный элемент не стержень( сжатие растяжение), а балка(сжатие-растяжение-изгиб). Тем не менее Ваша модель может смещаться вдоль X без реакции. Нужно закрепить одну из двух точек в направлении X.

[vitaly333]

Цитата:

Нужно закрепить одну из двух точек в направлении X.

Закрепил. Всё теперь действительно все решения сходятся!. Спасибо вам. Тогда такой вопрос - почему МКЭ на статически неопределимых задачах дает неодинаковое решение при выборе различных методов решения и учета ГУ?

[VLarin]

В принципе, я примерно так и говорил ранее - изначальная матрица жесктости была кривой и решения не имела. Поэтому при решении получался полный бред, на разных решателях получались разные результаты:)))
Если задача поставлена правильно, то никаких трудностей не возникает и связка CG+cholesky хорошо справляется с матрицами МКЭ.

[Zai]

Считайте детерминант перед обращением матрицы жесткости. Он не должен быть равен нулю.

[vitaly333]

Цитата:

Считайте детерминант перед обращением матрицы жесткости. Он не должен быть равен нулю.

Помоему вы не в эту тему написали.

Добавлено спустя 23 минуты 29 секунд:

Цитата:

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

В принципе да , но как проверять перед решением СЛАУ матрицу на такую "кривизну" при минимальных затратах ресурсов ЭВМ. Как сделать так, чтобы пользователь, работающий с системой, не смог создать такую "плохую" матрицу.

Цитата:

CG+cholesky хорошо справляется с матрицами МКЭ.

Всё таки не очень хорошо справляется Не получилось у меня решить вышеприведенную задачу (уже с закреплением по оси х) с ITL - овским предобуславливателем. И ксати при тестировании различных задач (а протестировал я много) практически ни одну связка CG+cholesky(ITL) не решила. Связка Мой переобуславливатель + CG некоторые решила но не все. "Одинокий" CG как ни странно очень хорошо работает без предобуславливателей , практически со всеми задачами он справился на ура. От своего прекондитинера отказался из-за того что во 1-ых: на больших задачах очень долго выполняется.
во 2-ых: некоторые не решает , вследствии жесткого отбрасывания новых появившихся ненулевых элементов. в 3-их: всю работу , как вы и говорили, делает именно он а CG ничего не остается

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

Кстати у вас МКЭ как я понял применяется совсем в другой области , а у меня в строительстве поэтому возможно что у вас матрицы получаются совсем другие(немного другой структуры) и связка CG+cholesky у вас работает а у меня нет.

[Zai]

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

[vitaly333]

Цитата:

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

Строить обратную матрицу это слишком дорогого стоит (эквивалентно решению n - систем с n - неизвестными). Поэтому у себя в программе я использую треугольное разложение (прямые методы, а именно Холесский).

[Comga]

Zai писал(а):

Считайте детерминант перед обращением матрицы жесткости. Он не должен быть равен нулю.

В силу погрешностей вычисления детерминант вырожденной системы может быть не равен нулю. А абсолютное значение детерминанта мало что говорит о системе. Стандартный путь - вычисление обусловленности матрицы.

[metallica4ever]

Не подскажете, где можно найти код алгоритма CG (сопряженных градиентов), использующий предварительное разложение по Холецкому?