Программирование МКЭ в МATLAB. Даутов P.З. - 57 стр.

UptoLike

Составители: 

2.4. Учет краевых условий. Формирование системы МКЭ. 57
мы рассмотрели функцию assemba вычисления K и F . Необходима
аналогичная функция вычисления H и G; в этой же функции есте-
ственно вычислить вектор столбец u
D
размера n
p
такой, что все его
элементы равны нулю, кроме элементов с номерами j i
d
, равными
u
D j
.
Матрица K
0
получается из матрицы A вычеркиванием его строк
и столбцов с номерами из множества i
d
. Чтобы реализовать эту опе-
рацию, удобно ввести матрицу N размера n
p
× m, m = numel(i
n
),
которую удобнее определить через его транспонированную N
T
ра-
венством
N
T
(u
1
, u
2
, . . . , u
n
p
)
T
= (u
i
1
, u
i
2
, . . . , u
i
n
)
T
.
Таким образом, матрица N
T
реализует требуемую операцию вычер-
кивания элементов u с номерами из i
d
. Легко проверить, что
N=s p a r s e ( [ i_1 , i_2 , \ ld o t s , i_n ] , 1 :m, 1 , np ,m) ;
Следовательно, если вычислить H, G, N и u
D
, то K
0
и F
0
легко по-
лучить по формулам
K
0
= N
T
(K + H)N, F
0
= N
T
((F + G) (K + H)u
D
).
Отметим, что после решения системы K
0
u
0
= F
0
формула u =
Nu
0
+ u
D
позволяет получить вектор узловых значений решения u
h
схемы МКЭ (дополняя вектор u
0
граничными значениями u
D
). Имея
в виду расчетные формулы для H и G, полученные ранее, приходим
к следующей функции.
f u n c t i o n [N,H, uD,G]= assembe ( bc , p , e )
% ASSEMBE: Assembles boundary c o n d i t i o n s co n t r i b u t i o n s i n a PDE problem .
%
% The f o l l o w i n g c a l l i s a l s o al lo we d :
% N=ASSEMBE( bc , p , e )
% [N,H]=ASSEMBE( bc , p , e )
% [N,H, uD]=ASSEMBE( bc , p , e )
np=s i z e ( p , 2 ) ;
H=s p a r s e (np , np ) ; G=s p a r s e ( np , 1 ) ; N=speye (np , np ) ;
uD=s p a r s e ( np , 1 ) ;
i f a l l ( bc.bsN==i n f ) && ( numel ( b c . s g )==0) && ( numel ( bc.mu)==0)
re t u rn % on a l l boundary homogeneous Neumann b . c .
end
2.4. Учет краевых условий. Формирование системы МКЭ.                                               57


мы рассмотрели функцию assemba вычисления K и F . Необходима
аналогичная функция вычисления H и G; в этой же функции есте-
ственно вычислить вектор столбец uD размера np такой, что все его
элементы равны нулю, кроме элементов с номерами j ∈ id , равными
uD j .
     Матрица K0 получается из матрицы A вычеркиванием его строк
и столбцов с номерами из множества id . Чтобы реализовать эту опе-
рацию, удобно ввести матрицу N размера np × m, m = numel(in ),
которую удобнее определить через его транспонированную N T ра-
венством
              N T (u1 , u2 , . . . , unp )T = (ui1 , ui2 , . . . , uin )T .
Таким образом, матрица N T реализует требуемую операцию вычер-
кивания элементов u с номерами из id . Легко проверить, что
N=s p a r s e ( [ i_1 , i_2 , \ l d o t s , i_n ] , 1 : m, 1 , np ,m) ;

Следовательно, если вычислить H, G, N и uD , то K0 и F0 легко по-
лучить по формулам

            K0 = N T (K + H)N,                      F0 = N T ((F + G) − (K + H)uD ).

Отметим, что после решения системы K0 u0 = F0 формула u =
N u0 + uD позволяет получить вектор узловых значений решения uh
схемы МКЭ (дополняя вектор u0 граничными значениями uD ). Имея
в виду расчетные формулы для H и G, полученные ранее, приходим
к следующей функции.

f u n c t i o n [ N, H, uD ,G]= assembe ( bc , p , e )
% ASSEMBE: Assembles boundary c o n d i t i o n s c o n t r i b u t i o n s i n a PDE p r o b l e m .
%
% The f o l l o w i n g c a l l i s a l s o a l l o w e d :
% N=ASSEMBE( bc , p , e )
% [ N,H]=ASSEMBE( bc , p , e )
% [ N, H, uD]=ASSEMBE( bc , p , e )
np=s i z e ( p , 2 ) ;

H=s p a r s e ( np , np ) ;   G=s p a r s e ( np , 1 ) ; N=s p e y e ( np , np ) ;
uD=s p a r s e ( np , 1 ) ;

i f a l l ( bc.bsN==i n f ) && ( numel ( b c . s g )==0) && ( numel ( bc.mu)==0)
   return       % on a l l boundary homogeneous Neumann b . c .
end