Автоматихация производства ЭВА, страница 6

Рассмотрим систему из 4-х линейных алгебраических уравнений вида:

а⁺¹₁₁х₁ + а⁺¹₁₂ х₂ + а⁺¹₁₃ х₃ + а⁺¹₁₄ х₄ = b⁺¹₁

а⁺¹₂₁ х₁ + а⁺¹₂₂ х₂ + а⁺¹₂₃ х₃ + а⁺¹₂₄ х₄ = b⁺¹₂

а⁺¹₃₁ х₁ + а⁺¹₃₂ х₂ + а⁺¹₃₃ х₃ + а⁺¹₃₄ х₄ = b⁺¹₃

а⁺¹₄₁ х₁ + а⁺¹₄₂ х₂ + а⁺¹₄₃ х₃ + а⁺¹₄₄ х₄ = b⁺¹₄

Выразив из первого уравнения переменную x₁, имеем:

Подставляя полученное выражение для х₁ во 2-е, 3-е и 4-е уравнения и приводя подобные члены, приходим к системе:

а⁺¹₁₁х₁ + а⁺¹₁₂ х₂ + а⁺¹₁₃ х₃ + а⁺¹₁₄ х₄ = b⁺¹₁

0 + (а₂₂–а₂₁[а₁₂/а₁₁])х₂+(а₂₃–а₂₁[а₁₃/а₁₁])х₃ +(а₂₄–а₂₁[а₁₄/а₁₁])х₄ = b₂–b₁(а₂₁/а₁₁)

0 + (а₃₂–а₃₁[а₁₂/а₁₁])х₂+(а₃₃–а₃₁[а₁₃/а₁₁])х₃ +(а₃₄–а₃₁[а₁₄/а₁₁])х₄ = b₃–b₁(а₃₁/а₁₁)

0 + (а₄₂–а₄₁[а₁₂/а₁₁])х₂+(а₄₃–а₃₁[а₁₄/а₁₁])х₃ +(а₄₄–а₄₁[а₁₄/а₁₁])х₄ = b₄–b₁(а₄₁/а₁₁)

В трех последних уравнениях все коэффициенты а_pq и b_p должны иметь верхний индекс (+1), поскольку эти коэффициенты взяты из исходной системы. Далее указанный индекс будет использован для обозначения номера итерации решения исходной системы. Введем следующие обозначения:

a_pq^+(k+1)= a_pq^+k- a_pk^+k( a_kq^+k/ a_kk^+k) ; b_p^+(k+1)= b_p^+k- b_k^+k( a_pk^+k/ a_kk^+k) (14.1)

Тогда последнюю систему можно переписать в виде:

а⁺¹₁₁х₁ + а⁺¹₁₂ х₂ + а⁺¹₁₃ х₃ + а⁺¹₁₄ х₄ = b⁺¹₁

0 + а⁺²₂₂ х₂ + а⁺²₂₃ х₃ + а⁺²₂₄ х₄ = b⁺²₂

0 + а⁺²₃₂ х₂ + а⁺²₃₃ х₃ + а⁺²₃₄ х₄ = b⁺²₃

0 + а⁺²₄₂ х₂ + а⁺²₄₃ х₃ + а⁺²₄₄ х₄ = b⁺²₄

Выразив из второго уравнения переменную x₂, имеем:

Подставляя полученное выражение для х₂ в 3-е и 4-е уравнения и приводя подобные члены, приходим к системе:

а⁺¹₁₁х₁ + а⁺¹₁₂ х₂ + а⁺¹₁₃ х₃ + а⁺¹₁₄ х₄ = b⁺¹₁

а⁺²₂₂ х₂ + а⁺²₂₃ х₃ + а⁺²₂₄ х₄ = b⁺²₂

(а₃₃⁺²–а₃₂⁺² [а₂₃⁺²/а₂₂⁺²])х₃ + (а₃₄⁺²–а₃₂⁺² [а₂₄⁺²/а₂₂⁺²])х₄ = b⁺²₃–b₂⁺²(а₃₂⁺²/а₂₂⁺²)

(а₄₃⁺²–а₄₂⁺² [а₂₃⁺²/а₂₂⁺²])х₃ + (а₄₄⁺²–а₄₂⁺² [а₂₄⁺²/а₂₂⁺²])х₄ = b⁺²₄–b₂⁺²(а₄₂⁺²/а₂₂⁺²)

Коэффициент при неизвестной х3 в третьем уравнении логично было бы обозначить как а₃₃⁺³. Попробуем получить его формально, используя первую формулу в выражении (14.1). С этой целью обозначим p=3 (№ строки) , q=3 (№ столбца), k=2 (номер текущей итерации) и подставим эти индексы в (14.1):

a₃₃⁺⁽²⁺¹⁾= a₃₃⁺²- a₃₂⁺²( a₂₃⁺²/ a₂₂⁺²) ;

Получили очевидное совпадение результатов. Вычислим аналогично остальные коэффициенты при неизвестных в третьем и четвертом уравнениях:

a₃₄⁺⁽²⁺¹⁾= a₃₄⁺²- a₃₂⁺²( a₂₄⁺²/ a₂₂⁺²) =a₃₄⁺³

a₄₃⁺⁽²⁺¹⁾= a₄₃⁺²- a₄₂⁺²( a₂₃⁺²/ a₂₂⁺²) =a₄₃⁺³

a₄₄⁺⁽²⁺¹⁾= a₄₃⁺²- a₄₂⁺²( a₂₃⁺²/ a₂₂⁺²) =a₄₄⁺³

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

где: b₃⁺³= b⁺²₃–b₂⁺²(а₃₂⁺²/а₂₂⁺²) и b₄⁺³= b₄⁺² –b₂⁺²(а₄₂⁺²/а₂₂⁺²).

После третьей итерации система примет вид:

а₁₁⁺¹х₁+ а₁₂⁺¹х₂+ а₁₃⁺¹х₃+ а₁₄⁺¹х₄ = b₁⁺¹

0+ а₂₂⁺²х₂+ а₂₃⁺²х₃+ а₂₄⁺²х₄ = b₂⁺²

0+ 0+ а₃₃⁺³х₃+ а₃₄⁺³х₄ = b₃⁺³

0+ 0+ 0+ а₄₄⁺⁴х₄ = b₄⁺⁴

где: a₄₄⁺⁴= a₄₄⁺³–a₄₃⁺³(а₃₄⁺³/а₃₃⁺³) и b₄⁺⁴= b₄⁺³ –b₃⁺³(а₄₃⁺³/а₃₃⁺³).

Решение полученной системы выполняем методом обратной прогонки. Из четвертого уравнения вычисляем неизвестную Х4:

х₄ = b₄⁺⁴ /а₄₄⁺⁴

Из третьего уравнения вычисляем неизвестную Х3:

х₃ = [b₃⁺³ -(а₃₄⁺³х₄)]/а₃₃⁺³

Из второго уравнения вычисляем неизвестную Х2:

х₂+ = [b₂⁺² – (а₂₄⁺²х₄ + а₂₃⁺²х3)]/а₂₂⁺²

Наконец, из первого уравнения вычисляем неизвестную Х1:

х₁ = [b₁⁺¹ – (а₁₄⁺¹х₄+а₁₃⁺¹х₃+а₁₂⁺¹х₂)]/а₁₁⁺¹

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

Однако, непосредственно перед решением система должна быть преобразована, поскольку, как правило, некоторые компоненты неизвестного вектора Ф (в программе – это вектор Xr) узловых значений известны. Так, в большинстве задач теории поля некоторые граничные значения искомой величины заданы; во всех задачах теории упругости должны быть фиксированы некоторые перемещения с тем, чтобы исключить перемещение среды жесткого тела.

То есть, элементы матриц [K] и {F} системы: [K]{Ф}={F} необходимо преобразовать, чтобы ее решение после преобразование давало правильный результат. При этом желательно не изменять программу решения самой системы, поскольку это повлечет за собой трудности при программировании.

Пусть в исходной системе задано фиксированное значение р-й переменной (Х_р=Q). Преобразование системы проводим по шагам:

• коэффициенты р-й строки, кроме диагонального коэффициента, равного а_pp⁺¹, приравниваем нулю;

• свободный член в р-й строке заменяем произведением: (а_pp⁺¹Q);

• уравнения, содержащие переменную Х_р, преобразуем, вычитая из обеих частей каждого из них произведение (а_qp⁺¹Q), где q – номер строки (qp).

Проиллюстрируем это на примере системы уравнений:

Здесь, согласно условию задачи, фиксирована одна степень свободы узлового параметра {Т₁=150}. Преобразование системы проводим по шагам:

• коэффициенты 1-й строки, кроме диагонального коэффициента, равного К₁₁=46.6, приравниваем нулю:

• 46,6T₁ + 0 + 0 + 0 = 1000

• свободный член в 1-й строке заменяем произведением: (К₁₁Т₁)=6990:

• 46,6T₁ + 0 + 0 + 0 = 6990

• переменная Т₁ входит еще во второе уравнение, поэтому вычитаем из левой и правой части 2-го уравнения произведение К₂₁Т₁=(-21,7150):

0 + 93,2T₂ – 21,7T₃ + 0 = 2000–(-3255)= 5255

Таким образом, искомая система для решения примет вид:

46,6T₁ + 0 + 0 + 0 = 6990

0 + 93,2T₂ – 21,7T₃ + 0 = 5255

0 - 21,7T₂ + 93,2T₃ – 21,7T₄ = 2000

0 + 0 – 21,7T₃ + 56,6T₄ = 1400

что совпадает с системой из раздела 12.

В программе решения системы уравнений, приводимой на стр.27, преобразование выполняется оператором:

For i:=1 to n do If defX[i]=1 Then UppCase(i);

Собственно преобразование выполняется подпрограммой UppCase, которая в качестве параметра принимает номер фиксированной степени свободы. Последний выбирается из исходного линейного массива defX, в котором каждый фиксированный параметр должен быть заранее помечен единицей. Учебная pascal-программа преобразования и решения системы из n ЛАУ приведена в приложении:

5. Метод конечных элементов (МКЭ)

5.1. Типы конечных элементов.

Используемые в настоящее время численные методы рассматривают ДУ непосредственно в той форме, в которой в которой они были выведены (без каких-либо дальнейших математических преобразований и манипуляций) при помощи: (а) аппроксимации дифференциальных операторов конечно-разностными алгебраическими операторами, действующими в последовательности узлов, находящихся в области; (б) при помощи представления самой области элементами среды, не являющимися бесконечно малыми (то есть конечными элементами), которые в совокупности аппроксимируют реальную систему. Наиболее громоздкой и трудно программируемой операцией в МКЭ является учет граничных условий задачи, причем точность полученного численного решения полностью зависит от степени измельченности сетки, определяющей узловые точки. Отсюда следует, что в процессе решения задачи программисту приходится иметь дело с системами алгебраических уравнений очень высокого порядка.

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

К недостаткам метода следует отнести: (а) в основе МКЭ лежит дискретизация всего ИТО, что неизбежно ведет к очень большому количеству конечных элементов (особенно в трехмерных задачах); (б) МКЭ часто приводит к нереальным разрывам значений физических величин между смежными элементами.

При решении задач методом КЭ используются одномерные, 2 и 3-мерные КЭ. Одномерный КЭ показан на рисунке 9.1. Площадь поперечного сечения одномерного КЭ может изменяться по длине, но в большинстве практических задач ее считают постоянной. Наиболее часто такой элемент используется в одномерных задачах распространения тепла и в задачах строительной механики.

Для построения дискретной модели двумерной области используют треугольники и 4-хугольники, которые могут иметь как прямолинейные, так и криволинейные стороны.

Собственно процесс дискретизации ИТО может быть разделен на два этапа: (а) разбиение ИТО на КЭ и (б) нумерация элементов и узлов. Последний этап может существенно повлиять на эффективность вычислений.

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

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

На практике используются три типа элементов: симплекс-элемент, комплекс-элемент и мультиплекс-элемент.

Симплекс - элементам соответствуют полиномы, содержащие константу и линейные члены. Число коэффициентов в таком полиноме на 1 больше размерности координатного пространства. Например, полином:  = ₁ +₂ x +₃ y представляет собой симплексную функцию для двумерного треугольного элемента. Этот полином линеен по X и Y и содержит три коэффициента, потому что треугольник имеет три узла.