2014 dxdy logo

Научный форум dxdy

Математика, Физика, Computer Science, Machine Learning, LaTeX, Механика и Техника, Химия,
Биология и Медицина, Экономика и Финансовая Математика, Гуманитарные науки




Начать новую тему Ответить на тему
 
 Работа некоторых методов решения СЛАУ в МКЭ
Сообщение30.04.2007, 15:29 


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

Тип конструкции: плоская стержневая рама.
см.изображение
В узлах 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 вычеркиванием столбцов и строк и без него.
Для более наглядного представления см здесь

 Профиль  
                  
 
 
Сообщение02.05.2007, 08:55 
Заслуженный участник
Аватара пользователя


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

 Профиль  
                  
 
 
Сообщение02.05.2007, 14:59 


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

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

Цитата:
У Вас механизм

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

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

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

 Профиль  
                  
 
 
Сообщение02.05.2007, 15:07 
Заслуженный участник
Аватара пользователя


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

 Профиль  
                  
 
 
Сообщение03.05.2007, 14:00 


26/11/06
76
Цитата:
Нужно закрепить одну из двух точек в направлении X.


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

 Профиль  
                  
 
 
Сообщение03.05.2007, 17:50 


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

 Профиль  
                  
 
 
Сообщение03.05.2007, 19:51 
Заслуженный участник
Аватара пользователя


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

 Профиль  
                  
 
 
Сообщение04.05.2007, 00:14 


26/11/06
76
Цитата:
Считайте детерминант перед обращением матрицы жесткости. Он не должен быть равен нулю.


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

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

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


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

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

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

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

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

 Профиль  
                  
 
 
Сообщение04.05.2007, 07:45 
Заслуженный участник
Аватара пользователя


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

 Профиль  
                  
 
 
Сообщение04.05.2007, 23:37 


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

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

 Профиль  
                  
 
 
Сообщение15.05.2007, 07:51 


02/05/06
56
Zai писал(а):
Считайте детерминант перед обращением матрицы жесткости. Он не должен быть равен нулю.


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

 Профиль  
                  
 
 Re: Работа некоторых методов решения СЛАУ в МКЭ
Сообщение06.05.2011, 15:08 


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

 Профиль  
                  
Показать сообщения за:  Поле сортировки  
Начать новую тему Ответить на тему  [ Сообщений: 12 ] 

Модераторы: Karan, Toucan, PAV, maxal, Супермодераторы



Кто сейчас на конференции

Сейчас этот форум просматривают: нет зарегистрированных пользователей


Вы не можете начинать темы
Вы не можете отвечать на сообщения
Вы не можете редактировать свои сообщения
Вы не можете удалять свои сообщения
Вы не можете добавлять вложения

Найти:
Powered by phpBB © 2000, 2002, 2005, 2007 phpBB Group