Министерство науки и высшего образования Российской Федерации
Балтийский государственный технический университет «Военмех»

А.В. СОЛДАТКИН, Е.С. БАРАНОВА

ВВЕДЕНИЕ В МЕТОД
КОНЕЧНЫХ ЭЛЕМЕНТОВ
Учебное пособие

Санкт-Петербург
2020

УДК 517.962.1(075.8)
С 60
С 60

Солдаткин, А.В.
Введение в метод конечных элементов:
учебное пособие / А.В. Солдаткин, Е.С. Баранова; Балт. гос. техн. ун-т. – СПб., 2020. –
123 с.
ISBN 978-5-907324-05-3
Пособие посвящено введению в метод конечных
элементов и его приложению к решению краевых
задач для дифференциальных уравнений и уравнений
с частными производными. Рассмотрены схемы получения уравнений при помощи метода Галеркина и ряд
других вопросов. Теоретические основы проиллюстрированы примерами и задачами. Приведены примеры
компьютерного моделирования в среде Matlab.
Для студентов, магистров, аспирантов и инженерно-технических работников, специализирующихся в
применении численных методов.

УДК 517.962.1(075.8)

Р е ц е н з е н т вед. науч. сотрудник Ин-та вычислит. технологий
СО РАН, д-р физ.-мат. наук, проф. НГУ Г.С. Хакимзянов

Утверждено
редакционно-издательским
советом университета

ISBN 978-5-907324-05-3

© Авторы, 2020
© БГТУ, 2020

ВВЕДЕНИЕ
Метод конечных элементов (МКЭ) на сегодняшний день является
основным методом структурного анализа в целом ряде областей
науки и техники. Созданные на его основе комплексы программ для
мощных ЭВМ широко используются в строительстве, кораблестроении, аэрокосмической и автомобильной промышленности, акустике, оптике, медико-биологических исследованиях, разведке полезных ископаемых, при решении краевых задач механики сплошных
сред и электромагнетизма.
Широкое использование этого метода в значительной мере
объясняется простой физической интерпретацией его вычислительных операций, большой геометрической гибкостью и применимостью к широкому классу уравнений в частных производных. Он
позволяет достаточно точно описать криволинейные границы области
определения решения и краевые условия.
К настоящему времени существует большое количество работ,
посвященных теории и практике применения МКЭ. Рассматриваемый
круг вопросов необычайно широк, от теоретических основ метода
[2, 4, 5, 7] до различных аспектов, направленных на расширение области применения метода и повышения его эффективности [9].
Изложение основ метода часто зависит от специализации автора
(строительная механика, акустика и т.д.). Начинающему исследователю сложно получить информацию, на которую в первую очередь
следует обратить внимание. Настоящее пособие предназначено для
ознакомления с основами МКЭ и практического овладения этим
методом при решении некоторых математических задач, связанных с
решением дифференциальных уравнений в условиях различной
геометрии.
В пособии приводятся небольшие программные модули, созданные в Matlab для решения тех или иных задач, связанных непосредственно с методом или подготовкой данных или графического
отображения результатов. Такой подход способствует получению
практических навыков в процессе выполнения заданий, приведенных
в пособии.
Для усвоения курса не требуется глубоких знаний математики. Вполне достаточно тех навыков дифференцирования и интегрирования несложных функций, которыми обладает студент 2-3-го
курсов, а также элементарных сведений из алгебры матриц.
3

1. КЛАССИЧЕСКИЕ МЕТОДЫ ВЗВЕШЕННЫХ НЕВЯЗОК
1.1. Вводные замечания
Рассмотрим некоторые основные определения и свойства
системы функций
(1.1)
ϕ1 ( x), ϕ2 ( x),  , ϕn ( x) .
Если можно составить линейную комбинацию функций,
например
и для них справедливы свойства:
ϕ1 ( x) + ϕ2 ( x) = ϕ2 ( x) + ϕ1 ( x) ,

(α + β) ⋅ ϕ( x) = α ⋅ ϕ( x) + β ⋅ ϕ( x) ,
α ⋅ (ϕ1 ( x) + ϕ2 ( x)) = α ⋅ϕ1 ( x) + α ⋅ ϕ2 ( x) ,
где α и β ‒ числовые коэффициенты, то такие функции называют
элементами линейного пространства R.
Скалярное произведение функций ϕ1 ( x) и ϕ2 ( x) представляет
собой операцию
b

( ϕ1, ϕ2 )= ∫ ϕ1 ( x) ϕ2 ( x)dx.
a

Мера (норма) функции ϕ определяется следующим образом:

( ϕ, ϕ) .

ϕ=

Система функций (1.1) называется линейно-независимой, если
равенство с1 ⋅ϕ1 ( x) + с2 ⋅ ϕ2 ( x) +  + сn ⋅ ϕn ( x) =0 возможно только когда
все коэффициенты сi равны нулю.
Система линейно-независимых функций считается полной, если
для любой сколь угодно малой положительной величины ε можно
найти такое число n и набор постоянных сi , при которых для
произвольной допустимой функции u того же пространства R
n

(n)
справедливо неравенство u − u ( n ) < ε, где u=
∑ c i ⋅ϕi ( x) . Иначе

это условие можно представить в виде

u (n) → u
При этом функции ϕi ( x)
коэффициентами Фурье.

i =1

при n → ∞.
называются базисными, а ci

4

‒

Если для функций

(

)

=
ϕi , ϕ j 0 при i ≠ j и

ϕ i (x )

выполнены условия

ϕ i (x ) = 1 ,

( ϕi , ϕi )=1 , то базисные функции называются

ортогональными нормированными.
Норма функции u ( n ) для взаимно ортогональной последовательности {ϕn ( x)} определяется как

=
u (n)

n

∑ ci 2 ( ϕi , ϕi ) .

i =1

Следует отметить, что с помощью процесса ортогонализации
любая линейно-независимая система неортогональных функций
может быть заменена системой ортогональных нормированных
функций.
В теории рядов Фурье имеет место следующая теорема, которая
является основанием для проекционно-сеточных методов (метод
Галёркина, методы взвешенных невязок).
Теорема. Пусть {ϕn ( x)} – полная система функций с ненулевой
нормой, ортогональных на отрезке [a, b] . Если непрерывная функция
f ( x) ортогональна на отрезке [a, b] ко всем функциям ϕn ( x) , т.е.
b

∫ f ( x) ϕn ( x)dx =0 (n =1, 2,),

a

то f ( x) ≡ 0 при a ≤ x ≤ b.
Например, если вектор f ( f1, f 2 , f3 ) ортогонален к трем линейно-

независимым (базисным) векторам a(a1, a2 , a3 ), b(b1, b2 , b3 ), c(c1, c2 , c3 ) ,
т.е.
0  f1 = 0f ,
( f , a ) = 0  f1a1 + f 2 a2 + f3a3 =



(1.2)
f0 ,
0
0 ⇔  f2 =
( f , b) =⇔
 f1b1 + f 2b2 + f3b3 =
( f , c ) = 0
f c + f c + f c =
 f = 0f ,

 3
 11 2 2 3 3 0
то он тождественно равен нулю.
Определим оператор L( ) как действие, которое при применении
его к данной функции из некоторого класса функций приводит к
появлению другой функции p .
Множество функций, на которые действует оператор, называется
областью определения оператора.
5

Определение 1. Оператор называется линейным, если он обладает свойством
L(α1u1 + α 2u2 ) = α1L(u1 ) + α 2 L(u2 ) .
Определение 2. Оператор называется самосопряженным, если

( y, Lx )= ( x, Ly )

для любых элементов x, y .
Нетрудно понять, что если x, y ‒ векторы, а L ‒ матрица, то это
требование равносильно свойству симметрии матриц:
( y, Lx )= y т Lx
⇔L =
Lт .

т т
( Ly, x )= y L x
Определение 3. Оператор называется положительно определённым, если
( Lx, x ) ≥ 0 ∀x и ( Lx, x ) = 0 ⇔ x = 0 .
Приведенные определения 1 – 3 являются общими как для
алгебраических, так и для дифференциальных операторов.
Для проверки условий симметрии оператора можно использовать
правило интегрирования по частям (аналогично операции
транспонирования матрицы).
В качестве примера рассмотрим дифференциальный оператор
d  du 
Lu =−  p  + q ⋅ u , где p( x) ≥ c0 > 0; q( x) ≥ 0; ∀ x ∈ [0,1] ,
dx  dx 
действующий на множестве функций U , дважды дифференцируемых
в=
Ω {x | 0 < x < 1} и равных нулю на границах интервала.
Пусть u и v ‒ две произвольные функции класса U . Составим
произведение:





( Lu, v ) = ∫  −  p  ⋅ v + q ⋅ u ⋅ v dx .
dx
dx
1

0

d

du







Последовательно интегрируя первое слагаемое по частям,
получим
1

∫−

0

=x 1 =
x 1
1

d  du 
du
−vp
 p  ⋅ vdx =
dx=
dx
dx x



du dv
dv
dx =
up
dx
dx
dx x
0=
0
+∫ p

1
d  dv 
d  dv 
p
⋅
udx
=
−
∫ dx  p dx  ⋅ udx.


dx
dx




0
0

1

−∫

6

0

−

Далее убеждаемся, что это самосопряжённый и положительно
определённый оператор:





( Lu, v ) = ∫  −  p  ⋅ u + q ⋅ v ⋅ u dx = ( u, Lv ) ,
dx
dx
1

0

d

dv







2
1
1

 d  du 
 du 
2
u
,
Lu
p
u
q
u
dx
p
=
−
⋅
+
⋅
=
   + q ⋅ u 2 dx ≥ 0.
(
) ∫ 
∫




0  dx  dx 
0   dx 

du
= 0 .и u = 0 . Следовательно, условие
dx
( u, Lu ) = 0 равносильно тому, что u = 0 .
Отсюда получим

1.2. Классические и «слабые» формулировки краевых задач
для линейных дифференциальных уравнений
В отличие от метода конечных разностей в проекционных
методах основанием для нахождения приближенного решения служит
не исходная, а вариационная постановка задачи. Это позволяет
расширить класс функций, в котором ищется решение (слабые
формулировки), и приводит к дополнительному удобству в реализации краевых условий (естественные краевые условия).
Чтобы охарактеризовать «степень непрерывности», введем
следующую классификацию.
Рассмотрим функцию u ( x) , заданную в области Ω , и обозначим
через L2 (Ω) множество функций u ( x) , квадрат которых будет
интегрируемым в области Ω , т.е.


L2 (Ω) u ( x) ∫ u 2 dx < ∞
=

Ω


.


(1.3)

Наложение ограничений на непрерывность производных приводит к множеству функций, которые называют пространствами
Соболева:
2

2

 d ku  

2  du 

W
=
u ( x) ∫ u +   +  +  k  dx < ∞
 dx 
Ω

 dx  


(k )

7







(1.4)

Схематично примеры функций, принадлежащих различным
пространствам, приведены на рис.
1-3.
Обозначим через С k [0,1]

Рис. 1. Функция, интегрируемая
с квадратом u ∈ L2

множество функций, k раз непрерывно дифференцируемых на отрезке [0,1] , через С0k [0,1] – множество функций из С k [0,1] и
равных нулю при х = 0 .

Рис. 2. Функция, интегрируемая с квадратом первой производной u ∈ W (1)

Рис. 3. Функция, интегрируемая с квадратом второй производной

u ∈W

(2)

Задача 1. Классическая формулировка состоит в нахождении
функции u ( x) ∈ С 2 [0,1] , удовлетворяющей уравнению
−

d  du 
p =
 + q ⋅ u f , где 0 < x < 1,
dx  dx 

du
(1) = α. .
dx
Будем предполагать, что, p ∈ С1[0,1] , q, f ∈ C [0,1] , причем

и граничным условиям u (0) = 0,

p (1)

8

(1.5)

p( x) ≥ c0 > 0; q( x) ≥ 0; ∀ x ∈ [0,1].
Пусть v(x) – некоторая произвольная функция из класса С01[0,1] .
Умножая обе части уравнения (1.5) на эту функцию и интегрируя по
частям, приходим к следующей задаче.
Задача 2. Требуется найти функцию u ( x) , удовлетворяющую
условию
1
1
1
 du dv

(1.6)
∫  p dx dx + quv dx= ∫ fvdx + v(1)α для ∀v ∈ С0 [0,1]


0
0
Теорема 1. Задачи 1 и 2 эквивалентны.
Доказательство. Пусть u ( x) удовлетворяет задаче 1 и v( x) –
некоторая произвольная функция из класса С01[0,1] . Умножая
уравнение (1.5) на функцию v( x) и интегрируя по промежутку [0,1],
получим
1
 d  du 

0.
∫  − dx  p dx  + qu − f  vdx =



0
Интегрируя первое слагаемое по частям, с учетом граничных
условий получим
1
 du dv

∫  p dx dx + quv dx= ∫ fvdx + v(1)α для ∀v ∈ С0 [0,1].

0
0

1

1

Значит, из задачи 1 следует задача 2.
Обратно, пусть ∀v ∈ С01[0,1] , выполняется задача 2. Интегрируя
(1.6) по частям, получим









∫  − dx  p dx  + qu − f  vdx + v(1)  p dx − α  = 0 .

1

0

d

du





du







(1.7)

В силу произвольности в выборе v( x) ∈ С01[0,1] примем

 d  du 

0.
v( x) =
ϕ( x)  −  p  + qu − f  , где ϕ( x) > 0, ϕ(0) =
ϕ(1) =
 dx  dx 

Например, ϕ( x) = x(1 − x) . Тогда из равенства (1.7) получим

 d  du 
∫ ϕ( x)  − dx  p dx  + qu −



0

1

9

2


f  dx =
0.


Следовательно,

−

d  du 
0,
 p  + qu − f =
dx  dx 

а из (1.7)

 du

v(1)  p − α  = 0 .
dx


Выбрав теперь v( x) ∈ С01[0,1], v(1) ≠ 0 , получим

du
 du

 p − α  = 0 ⇔ p (1) = α .
dx
dx


Значит, из задачи 2 следует задача 1.
du
Граничное условие p (1) (1) = α не входит в область опредеdx
ления оператора, но входит в формулировку задачи 2. Такие условия
называют естественными. Напротив условие u (0) = 0 входит в
область определения функции u ( x) ∈ С01[0,1] и называется главным.
Для определения главных и естественных условий приведем (без
обоснования) правило: если дифференциальный оператор краевой
задачи имеет порядок 2m, то граничные условия, содержащие производные порядка меньше m считаются главными.
Формулировка задачи 2 называется слабой.
Введем обозначения:
1
 du dv

=
a (u , v) ∫  p
+ quv dx ,
dx
dx


0

(1.8)

1

ϕ(v) =
∫ fvdx .
0

В терминах этих обозначений уравнение (1.6) примет вид

a (u , v) =
ϕ(v) + v(1)α .
Здесь a (u , v) и ϕ(v) – соответственно билинейная и линейная
формы. Пусть u , v, w ‒ некоторые функции, а с1, с2 ‒ константы.
Пользуясь определением, нетрудно проверить, что они линейны по
каждому аргументу:
10

a (с1u + с2v, w)= с1a (u , w) + с2 a (v, w);
ϕ(с1u + с2v) = с1ϕ(u ) + с2ϕ(v).
Имеет место свойство симметрии:

a (u , v) = a (v, u );
(u , v) = (v, u ).
Проверим, что в нашем случае форма a (u , v) является
положительно определённой, т.е. a (u , u ) > 0 при u ≠ 0 , а условие
a (u , u ) = 0 равносильно условию u = 0. Действительно,
1

2

 du 
2
a (u , u=
) ∫  p
 + qu dx > 0,

0
  dx 

du
= u= 0.
dx
Следует отметить, что в формулировке (1.6) решение принадлежит классу С10 [0,1] (входят первые производные), в то время как в
классической формулировке (1.5) решение относится к классу
С02 [0,1] .
а условие a (u , u ) = 0 равносильно

1.3. Вариационная формулировка краевых задач для линейных
дифференциальных уравнений
Под функционалом будем понимать некоторое отображение
элементов функционального пространства на множество вещественных чисел.
Приведем теорему, которая при определённых условиях
позволяет представить решение краевой задачи как функцию, которая
доставляет минимум некоторому функционалу.
Теорема. Если a (u , v) – симметричная и положительно
1
определённая форма, ϕ(v) – линейная форма,=
F (u )
a (u , u ) − ϕ(u ) –
2
некоторый функционал, заданный на функциональном пространстве
U , то имеет место равносильность следующих двух задач:
1) найти такое u ∈ U , чтобы
a (u , v) =ϕ(v) для ∀v ∈ U;

2) найти такое u ∈ U , чтобы F (u ) ≤ F (v) для ∀v ∈ U .
11

(1.9)

Доказательство. Покажем, что задача 1 имеет единственное
решение. Предположим, что u1 и u2 ‒ два решения, введем в
рассмотрение вспомогательную функцию w= u1 − u2 , для которой
a ( w, w) = 0 и, следовательно, w = 0.
Предположим, что u – решение задачи 2, и покажем, что u также
является решением задачи 1.
Для v ∈ U и ∀ λ ∈ R

1
λ2
a (u , u ) + λa (u , v) +
a (v, v) − ϕ(u ) − λϕ(v) ≥
2
2
1
=
≥ F (u )
a (u , u ) − ϕ(u )
2

F (=
u + λv )

λ2
a (v, v) + λ(a (u , v) − ϕ(v)) ≥ 0 .
2
Из последнего неравенства, которое должно выполняться для
всякого λ ∈ R , следует, что a (u , v) − ϕ(v) = 0 .
Предположим, что u – решение задачи 1. Имеем
1
1
F (v) − F (u ) =F (u + (v − u )) − F (u ) = a (u , u ) + a (v − u , v − u ) − ϕ(u ) +
2
2
+ [ a (u , v − u ) − ϕ(v − u ) ] − F (u );
v − u ∈ U , так что, в силу предположения, a (u , v − u ) − ϕ(v − u ) = 0.

и

После

сокращений

получим

F (v) − F (u=
)

1
a (v − u , v − u ) ≥ 0.
2

Отсюда следует, что u – решение задачи 2.
Таким образом, краевая задача сводится к исследованию условий
экстремума некоторого функционала.
Рассмотрим функционал, соответствующий краевой задаче (1.5):

=
F (u )

2

1 1   du 
p
   + qu 2 − 2 fu dx − u (1)α .
∫
2 0   dx 


Пусть u ( x) ∈ C01[0,1] – решение, доставляющее минимум
функционалу, v( x) – некоторая произвольная функция класса

С01[0,1] . Рассмотрим
F=
(u + λv)

2

1 1   du
dv 
2
.
p
+
λ

∫

 + q (u + λv) − 2 f (u + λv) dx − (u (1) + λv(1))α
2 0   dx
dx 


12

Тем самым вопрос сводится к выполнению условий экстремума
функции одной переменной λ при λ= 0.
Величина λ ⋅

dF (u + λv)
называется первой вариацией функdλ
λ=0

ционала и обозначается δF .
Выполнив необходимое условие экстремума δF =
0 , получим
1
dF 1  du dv

= ∫p
+ quv dx − ∫ fvdx
=
− v(1)α 0 .
d λ 0  dx dx

0

Интегрируем первое слагаемое по частям:
1

vp

1
du
 d  du 

+ ∫  −  p  + qu − fv  vdx − v(1)α =0 .
dx 0 0  dx  dx 


Учитывая, что v(0) = 0 , получим









∫  − dx  p dx  + qu − f  vdx + v(1)  p dx − α 





1

d

du

0

du



=0 .
x =1

В силу произвольности функции v( x)

−

d  du 
 p  + qu= f , x ∈ [0,1] ,
dx  dx 
du
u (0) = 0, p (1) (1) = α.
dx
1.4. Метод Галёркина

Рассмотрим численную процедуру нахождения приближенного
решения краевой задачи
(1.10)
Lu = f , в области Ω,
(1.11)
F (u ) = 0, на границе ∂Ω;.
Приближенное решение аппроксимируется набором функций
ϕk ( x) , т.е.
n

u ( n )= ∑ α k ϕk ( x)
k =1

13

(1.12)

где α k – неизвестные коэффициенты, ϕk ( x) – линейно-независимые
функции, принадлежащие к полной последовательности. Будем
считать, что (1.12) удовлетворяет условиям (1.11) на границе области.
Подстановка (1.12) в (1.10) дает функцию ε , которую называют
невязкой:

=
ε Lu ( n) − f ≠ 0 .

(1.13)

Отметим, что при точном решении невязка равна нулю.
Стремятся, чтобы ошибка была равна нулю в среднем, полагая
равными нулю интегралы, взятые от невязки с некоторыми весовыми
функциями:

wi )
(ε ,=

0,=
i 1, 2,  , n ,

(1.14)

где wi – набор весовых функций.
Существует несколько разновидностей метода взвешенных
невязок: метод коолокаций, метод наименьших квадратов, метод
моментов, метод Галёркина. Все они различаются выбором весовых
функций.
Рассмотрим более подробно метод Галёркина, поскольку метод
конечных элементов является его дальнейшим развитием, связанным
с выбором базисных функций.
В основе метода лежит требование ортогональности базисных
функций ϕi к невязке ε :

 n

(1.15)
i 1, 2,  , n
i  0,=
 L( ∑ α k ϕk ) − f , ϕ=
 k =1

В случае линейного оператора (1.15) приводит к линейной
системе алгебраических уравнений для нахождения неизвестных αi .
Следует отметить, что метод Галёркина в равной степени применим и к нелинейным задачам.
Метод взвешенных невязок можно интерпретировать как процедуру для получения «слабых» решений.
Чтобы показать, как могут быть ослаблены требования к
непрерывности, и смягчить требование о том, чтобы базисные функции удовлетворяли всем краевым условиям, оставив среди них только
главные, рассмотрим уравнение второго порядка:

( ε, ϕ=
i)

L(u ) − f=

d 2u
dx 2

+u − =
x 0.
14

(1.16)

с главным граничным условием
при x =
0
u=
β

(1.17)

и естественным граничным условием
du
(1.18)
при x =
1.
=
α
dx
Потребуем, чтобы решение u удовлетворяло главному граничному условию, а весовая функция v – однородной форме главного условия, т.е.
(1.19)
u=
β
при x =
0, v =
0
при x =
0.
Можно записать:
1


d 2u
du

∫  2 + u − x  vdx + v(1) α − dx
0  dx


= 0,

(1.20)

x =1

(0.1)
где v( x) ∈ L2 , а u ( x) ∈ W (2) , т.е. функция и её первая производная
непрерывны.
Интегрируя (0.1) по частям, получаем

0

1

du dv 

∫ (u − x)v − dx dx dx + αv x=1 =0,


(1.21)

где u ( x) и v( x) ∈ W (1) , т.е. непрерывна только сама функция u ( x) .
Интегрируя уравнение (1.21) по частям, получаем
1

d 2v 
dv
−
+
=0,
u
x
v
u
dx + αv − u
(
)

∫
2
dx 0
dx 

0
1

(1.22)

где u ( x) ∈ L2 , а v( x) ∈ W (2) . Наиболее распространена формулировка
(1.21).
Пример. Найти функцию, удовлетворяющую уравнению

L(u ) − f=

d 2u
dx 2

+u + =
x 0

с граничными условиями
=
u 0=
при x 0;
du
= 0=
при x 1.
dx
15

(1.23)

(1.24)

Решение. Удовлетворяя главному граничному условию, представим искомую функцию u ( n) ( x) в виде
n

u ( n) ( x)= ∑ α k ϕk ( x),

(1.25)

k =1

где ϕk (=
x) x k +1 , =
k 1, 2, , n , удовлетворяя тем самым главному граничному условию u (0) = 0 , естественное граничное условие
du
(1) = 0 выполним в рамках формулировки (0.1), где α =0 .
dx
Примем
u (2) = α1ϕ1 ( x) + α 2ϕ2 ( x) = α1 x + α 2 x 2 .

(1.26)

Следовательно,

Вычислим невязку ε :

du (2)
=α1 + 2α 2 x .
dx

ε = 2α 2 + α1 x + α 2 x 2 + x

(1.27)

du
εb =
dx

(1.28)

и невязку на границе:
x =1

=α1 + 2α 2

и проведем ортогонализацию по отношению к ϕ1 ( x) и ϕ2 ( x) :
1

d 2u (2)

0

dx 2

1

2 (2)

∫ 
∫ 

0

d u

dx

2


du (2)
+ u (2) + x  ϕ1dx − ϕ1 (1)

dx


du
+ u (2) + x  ϕ2 dx − ϕ2 (1)

dx


1

x =1

= ∫ εxdx − εb x x =1 = 0,

(2)

0

1

= ∫ εx 2 dx − εb x 2
x =1

0

x =1

(1.29)

= 0.

После вычисления интегралов получим
2
3

3

4

3
1 
3

4   α1  
  =  .
17  α 2   1 

4

15 


Отсюда α1 =137 , α 2 =− 120 . , т.е.
139
278
16

(1.30)

=
u (2)

137
120 2
x−
x .
139
278

(1.31)

Точное решение:
sin x
− x.
cos1
Результаты вычислений приведены в табл. 1.
=
u

(1.32)

Таблица 1
x

u

0.5
1

0.3873
0.5574

du
dx
0.6242
0

u (2)

du (2)
dx

u (5)

0.3995
0.5540

0. 5540
0.1223

0.3873
0.5574

du (5)
dx
0.6242

8.25 ⋅10−6

Вычисление u (5) осуществлялось с помощью программы Matlab,
представленной ниже. Очевидно, что с увеличением номера приближения естественное граничное условие выполняется более точно.
Приведем программу и примеры реализации вычислений по
методу Галёркина в Matlab.
Листинг программы для решения дифференциального уравнения a2u '' + a1u ' + a0u =
f с краевыми условиями u ( x1 ) = u1, u ' ( x2 ) = u2 .
методом Галёркина:
syms L x P y(x)
syms C1 C2 C3 C4 C5 C6
syms x y(x);
poly=[-1 0 1 1];
bc=[0 0; 1 2];
n=5;
% размерность
a2=poly(1);
% коэффициенты уравнения
a1=poly(2);
a0=poly(3);
g=poly(4);
% правая часть
bc=bc.';
x1=bc(1);
% начало интервала
phi1=bc(2); % u(x1 ) = u1 .
x2=bc(3);
% конец интервала
L=x2
P=bc(4)
% u2 .
f1(x)=(x-x1);
% базисные функции
17

f2(x)=(x-x1)*x;
f3(x)=(x-x1)*x^2;
f4(x)=(x-x1)*x^3;
f4(x)=(x-x1)*x^4;
f5(x)=(x-x1)*x^5;
u(x)=phi1+C1*f1(x)+C2*f2(x)+C3*f3(x)+C4*f4(x)+C5*f5(x)
du(x)=diff(u(x));
d2u(x)=diff(du(x));
Eq=a2*d2u(x)+a1*du(x)+a0*u(x)-g;
eqn1=int(Eq*f1(x),x,x1,x2)+(P-du(x2))*f1(x2)
eqn2=int(Eq*f2(x),x,x1,x2)+(P-du(x2))*f2(x2)
eqn3=int(Eq*f3(x),x,x1,x2)+(P-du(x2))*f3(x2)
eqn4=int(Eq*f4(x),x,x1,x2)+(P-du(x2))*f4(x2)
eqn5=int(Eq*f5(x),x,x1,x2)+(P-du(x2))*f5(x2)
[C1,C2,C3,C4,C5]=solve(eqn1 == 0,eqn2 == 0,eqn3 == 0,eqn4 ==
0,eqn5 == 0)
x=(x_inf):.2:(x_sup);
y_apr=subs(u(x));
plot(x,y_apr,'b--o','Color','black'),grid
%yx = dsolve('-x*D2y+4*Dy+3*x*y = 8*sin(x)','y(2) = 0','Dy(5) = 2','x');
% Аналитическое решение (лучше прописывать явно)
yx = dsolve('-D2y+y = 1','y(0) = 0','Dy(1) = 2','x');
yx=subs(yx)
hold on
plot(x,yx,'b'),grid
hold on
e=exp(1);
e1=exp(-1);
c1=-e*(e+2)/(1+e*e);
c2=(-1+2*e)/(1+e*e);
plot(x,1+c1*exp(-x)+c2*exp(x),'r'),grid
hold off
Пример 1. Рассмотрим следующую краевую задачу:
 d 2u
+ u 1, где x ∈ [0,1],
− =
dx 2

du
u (0) 0,=
=
(1) 2.

dx
18

Её аналитическое решение:

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

и

аналити-

Метод Галеркина
2
Приближенное решение
Точное решение

1.8
1.6
1.4
1.2
1
0.8
0.6
0.4
0.2
0

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Рис. 4

Пример 2. Рассмотрим далее дифференциальное уравнение с
переменными коэффициентами и краевыми условиями смешанного
типа:
 d 2u
du
8sin x, где
− 2 + 4 + 3 xu =
dx
dx

u (2) = 0, du (5) = −2.

dx

x ∈[0,5],

Результаты сопоставления конечноэлементного и аналитического решения, полученного с помощью Matlab, приведены на
рис. 5.
19

Метод Галеркина
0.4
Приближенное решение
Точное решение

0.2
0
-0.2
-0.4
-0.6
-0.8
-1
-1.2
-1.4
-1.6

2

2.5

3

3.5

4

4.5

5

Рис. 5

Свойства системы и построенного с её помощью приближенного
решения зависят от выбора базисных функций ϕ1 ( x), ϕ2 ( x),, ϕn ( x) .
Если координатная система выбрана удачно, то с увеличением n
точность метода улучшается. Можно показать, что если координатная
система линейно-независима, то матрица К системы положительно
определена. Отсюда, в частности, вытекает однозначная разрешимость системы метода Галёркина.
Приведем примеры координатных систем метода Галёркина:
а) полиномиальная система:
ϕ1 ( x) = x(1 − x) , ϕ1 ( x) = x 2 (1 − x),  , ϕn ( x) = x n (1 − x).
б) тригонометрическая система:
ϕk ( x) = sin k πx , k = 1, 2,, n .
1.5. Метод Ритца
При определённых условиях краевая задача сводится к задаче
поиска минимума некоторого функционала F (u ) на некотором
классе допустимых функций. Этот класс определяется граничными
условиями, а также непрерывностью искомой функции и её производных.
20

Приближенное решение аппроксимируется набором функций
ϕk ( x) , таких, что
n

u ( n )= ∑ α k ϕk ( x) ,

(1.33)

k =1

где α k – неизвестные коэффициенты, ϕk ( x) – линейно-независимые
функции, принадлежащие к полной последовательности. Будем считать, что (1.33) удовлетворяет условиям на границе области при
любых значениях коэффициентов α k . Функционал принимает вид
n

n)
F (u (=
) F ( ∑ α k ϕk ( x))

(1.34)

k =1

Выберем постоянные α k так, чтобы
∂F
(1.35)
k 1, 2,, n.
= 0,=
∂α k
Если функционал F (u ) ‒ квадратичная функция u и u x то эти
уравнения линейны по отношению к неизвестным α k .
Обозначим через M n = F (u ( n ) ) . Очевидно, что M n ≥ M =
min F (u ).
Более того, M1 ≥ M 2 ≥ M 3 ≥  ≥ M n , так как любая линейная комбинация
функций ϕ1 ( x), ϕ2 ( x),, ϕn ( x) включена в линейные комбинации
функций ϕ1 ( x), ϕ2 ( x),, ϕn ( x), ϕn+1 ( x) .
Рассмотрим применение метода на простом примере нахождения
решения краевой задачи
 d 2u
2
0,
 2 +u + x =
 dx
u =
(1) 0.
 (o) u=

(1.36)

Сформулируем её как задачу нахождения минимума функционала
1

F (u )= ∫ (u x2 − u 2 − 2 x 2u )dx,

(1.37)

0

в классе функций, удовлетворяющих условиям u=
(0) u=
(1) 0 .
Представим искомую функцию u ( x) в виде ряда:
∞

u ( x) =x(1 − x) ∑ α k x k ;
k =0

21

(1.38)

в области 0 ≤ x ≤ 1 . Степенные функции в рассматриваемом интервале
образуют полную систему:
n

u ( n) ( x) =
x(1 − x) ∑ α k x k

(1.39)

k =0

и удовлетворяют краевым условиям (1.36) при любых значениях
коэффициентов α k .
Подставив u (0) ( x) =
α 0 x(1 − x) в равенство (1.37) и дифференцируя по α0 , получим
∂F 1 
= ∫ 2α 0 (1 − 4 x − 4 x 2 ) − 2α 0 x 2 (1 − 2 x + x 2 ) − 2( x3 − x 4 ) dx. (1.40)

∂α 0 0 

Вычисляя интеграл в (1.40) и учитывая (1.35), получаем α 0 =1 ,
6
поэтому первым приближением к u ( x) будет
1
x(1 − x)
6

u (0)=
( x)

(1.41)

Теперь аппроксимируем u ( x) с помощью u (1) ( x) . Подставляя
u (1) ( x) в интеграл (1.37) и вычисляя частные производные, получаем
∂F 1
= ∫ {2(1 − 2 x) α 0 (1 − 2 x) + α1 (2 x − 3 x 2  −


∂α 0 0

(1.42)

−2 α 0 x(1 − x) + α1 x 2 (1 − x)  x 2 (1 − x) − 2 x3 (1 − x)}dx =
0



и
∂F 1
= ∫ {2(2 x − 3 x 2 ) α 0 (1 − 2 x) + α1 (2 x − 3 x 2  −


∂α1 0

(1.43)

−2 α 0 x(1 − x) + α1x 2 (1 − x)  x 2 (1 − x) − 2 x 4 (1 − x)}dx =
0.



После вычисления интегралов
3
1
3
α1 + α 0 = ,

10
5
10

26
3
1

α1 + α 0 = .

10
15
105

Решая эту систему алгебраических уравнений, находим
22

(1.44)

=
α1

7
10
α0
, =
.
41
123

(1.45)

Точное решение (1.36):

 cos1 1 
2
=
−
u ( x) 2  2
 sin x − 2cos x + 2 − x .
 sin1 sin1 
В табл. 2 точное решение сравнивается с приближённым.

(1.46)

Таблица 2
x
0.25
0.50
0.75

u(x)
0.023
0.041
0.039

u0(x)
0.031
0.042
0.031

u1(x)
0.023
0.042
0.039

Листинг программы, решающей эту задачу в Matlab,
и результаты решения (рис. 6)
syms C1 C2 C3 C4 C5 L x
L=1
u(x)=x*(1-x)*(C1+C2*x+C3*x^2+C4*x^3+C5*x^4)
I=int(diff(u(x),x)^2-u(x)^2-2*(x^2)*u(x),x,0,L)
eqn1=diff(I,C1)
eqn2=diff(I,C2)
eqn3=diff(I,C3)
eqn4=diff(I,C4)
eqn5=diff(I,C5)
[C1,C2,C3,C4,C5]=solve(eqn1 == 0,eqn2 == 0,eqn3 == 0,eqn4 ==
0,eqn5 == 0)
x=(0):.1:(1);
y_apr=subs(u(x));
plot(x,y_apr,'b--o','Color','black'),grid
title('\fontsize{14} Метод Ритца')
exac=dsolve('D2u+u=-x^2','u(0)=0','u(1)=0','x');
exac=subs(exac);
y_exac=subs(exac);
hold on
plot(x,y_exac,'b')
grid on;
hold off
Результаты сопоставления аналитического и приближенного
решения задачи по методу Ритца приведены на рис. 6.
23

Метод Ритца
0.045
Приближенное решение
Точное решение

0.04
0.035
0.03
0.025
0.02
0.015
0.01
0.005
0

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Рис. 6

Следует отметить, что при использовании метода Галёркина
невязкой может быть
n

v = δu = ∑ ϕk ( x)δα k

(1.47)

k =1

Тогда уравнение метода взвешенных невязок приводит к системе
уравнений Галёркина в силу произвольности α k :
0
( Lu − f , δu ) =

(1.48)

Вместе с тем уравнение (1.48) может быть преобразовано к виду

( Lu − f , δu ) =

0 ⇔ δF = 0

(1.49)

и использовано для получения функционала задачи.
Рассмотрим, например, задачу

d  du 
 Lu =−  p  + q ⋅ u =f ,
dx  dx 

u=
(1) 0
 (0) u=

(1.50)

где p(x) ≥ c0 > 0; q(x) ≥ 0; ∀ x ∈ [0,1].
За счет выбора базисных функций обеспечим u=
(0) v=
(0) 0 .
Составим уравнения (1.48) с учетом (1.47) и (1.50):
24





( Lu − f , δu ) = ∫  −  p  + q ⋅ u − f  δudx =
dx
dx
1

d

0

du







0.

(1.51)

Интегрируя первое слагаемое по частям, с учетом u (0) = δ
=
δu (0) =
0 получим

( Lu − f , δu ) =δF (u ) =∫  p
1

0

du d δu

− quδu − f δu  dx =0 .
dx dx


(1.52)

Учитывая, что

(1.53)

Получим
2
1

1  du  1
δF (u ) = ∫  pδ   − qδ(u 2 ) − f δu  dx =


 dx  2
0 2

1  1  du 2 1
 
=
δ  ∫  p   − qu 2 − fu  dx  =
0.
 
0  2  dx  2

(1.54)

Следовательно,
1
1
=
F (u ) ∫ 
0
2

2

 du  1
p   − qu 2 − fu  dx .
 dx  2


(1.55)

Этот пример показывает, что метод Ритца является частным
случаем метода Галёркина.
Общим для обоих методов оказывается сложность выбора
системы базисных функций в области со сложной геометрией.
Необходимо также отметить тот факт, что при использовании
указанных методов матрица системы оказывается сплошь заполненной, т.е. все её элементы отличны от нуля. В этом принципиальное
отличие от разностного метода, при использовании которого матрица
системы линейных уравнений оказывается разряженной.
В 1943 г. Р. Курантом было замечено, что при специальном выборе базисных функций метод Галёркина приводит к системам линейных уравнений, по свойствам весьма близким к разностным.
25

В следующем разделе опишем простейший вариант метода конечных
элементов для одномерной задачи.
2. БАЗИСНЫЕ ФУНКЦИИ ДЛЯ ЭЛЕМЕНТОВ
2.1. Одномерный элемент
Рассмотрим основные типы элементов в одномерном случае, впоследствии полученные результаты распространим на двумерный.
Создание базисных функций базируется на использовании
интерполяционного многочлена Лагранжа. Базисные функции конструируются на каждом элементе. Рассмотрим одномерную дискретизацию области Ω на пять элементов (рис. 7).

Рис. 7

Для

элемента

[ x j , x j +1 ]
Ω( j ) =

неизвестную

функцию

u

посредством интерполяции будем аппроксимировать функцией
u ( h) ( x), x ∈ Ω( j ) :
u ( h) ( x=
) ax + b, x ∈ Ω( j ) .

(2.1)

Потребуем, чтобы в двух граничных точках Ω( j ) приближенная
функция u ( h) ( x) принимала значения неизвестной функции u:
(h)
=
u ( h ) ( x j ) u=
( x j +1 ) u j +1 .
j; u

(2.2)

Из этих условий можно определить a и b и переписать (2.1) в
виде
=
u ( h ) ( x) u j N j ( x) + u j +1 N j +1 ( x), ∀x ∈ Ω( j ) .

Функции N j ( x) и N j +1 ( x) называются функциями формы:
26

(2.3)

=
N j ( x)

x j +1 − x
x j +1 − x j

N j +1 ( x)
и=

x − xj
x j +1 − x j

, ∀x ∈ Ω( j ) .

(2.4)

Будем считать, что они равны нулю вне элемента
Ω =
[ x j , x j +1 ] . Из (2.4) видно, что они обладают свойствами:
( j)

1,
N j ( x) + N j +1 ( x) =
=
N j ( x j ) 1,=
N j ( x j +1 ) 0 и N=
N j +1 ( x j ) 0.
j +1 ( x j +1 ) 1 и =

Следует отметить, что очень удобно рассматривать локальную
координатную систему, поскольку область сводится к стандартной и
полученные результаты можно многократно использовать. Введем
параметр ξ ∈ [−1;1] следующим образом:
ξ=

2 x − x j +1 − x j
x j +1 − x j

(2.5)

.

Используя преобразование (2.5), можно записать (2.3) в виде
)
u ( h=
(ξ) u1N1 (ξ) + u2 N 2 (ξ),

∀ξ ∈ [−1,1],

(2.6)

1
1
где N1 (ξ=
)
(1 − ξ) и N 2 (ξ=
)
(1 + ξ) . Эти функции изображены
2
2
на рис. 8.
Функции формы более высокого порядка также могут быть
выражены в координатах ξ ∈ [−1;1] . Например, для квадратичного
элемента можно определить три функции формы (рис. 9):

Рис. 8

Рис. 9
27

1
N1 (ξ) = − ξ(1 − ξ);
2

N 2 (ξ)= (1 + ξ)(1 − ξ);

(2.7)

1
ξ(1 + ξ) .
2
Для того чтобы построить эти функции, можно рассмотреть
интерполяционный многочлен Лагранжа:
N3 (ξ)=

( x − x0 )  ( x − xk −1 )( x − xk +1 )  ( x − xn )
ψ (kn) =
,
( xk − x0 )  ( xk − xk −1 )( xk − xk +1 )  ( xk − xn )

(2.8)

где нижний индекс обозначает интерполяционную точку, верхний
n+1 – полное число точек интерполяции.
Например, чтобы получить N 2 (ξ) , положим x = ξ , возьмем три
интерполяционные точки x0 =ξ0 =−1, x1 =ξ1 = 0, x2 =ξ2 =1 . Из уравнения
(2.8) получим
( x − x1 )( x − x2 )
(ξ − 0)(ξ − 1)
1
=
N1 (ξ) =ψ 02 ( x) =
=
ξ(1 − ξ) .
( x0 − x1 )( x0 − x2 ) (−1 − 0)(−1 − 1) 2

(2.9)

Другие функции получены аналогично. Таким же образом можно
построить функции формы для кубического элемента (рис. 10):
9
1
1
27
1
(1 + ξ)( − ξ)(1 − ξ);
N1 (ξ) = − (1 − ξ)( + ξ)( − ξ); N 2 (ξ) =
16
3
3
16
3
(2.10)
27 1
9 1
1
( + ξ)(1 + ξ)(1 − ξ); N 4 (ξ) = − ( + ξ)( − ξ)(1 + ξ).
N 3 ( ξ) =
16 3
16 3
3

Рис. 10
28

2.2. Треугольный элемент
При решении двумерных задач часто используются прямоугольные и треугольные элементы. Последние получили наибольшее
распространение, так как если область аппроксимируется многоугольником, то его всегда можно разбить на треугольники.
Для линейного треугольного элемента неизвестная функция u в
пределах каждого элемента представляется в виде
u ( h) ( x, y ) =a + bx + cy ,

(2.11)

где a, b, c – некоторые постоянные, которые надо найти. Линейный
треугольный элемент имеет три узла, расположенные в его вершинах
(рис. 11). Узлы нумеруются против часовой стрелки числами 1,2,3 со
значениями u, обозначаемыми u1, u2, u3 соответственно.

Рис. 11

Уравнения, соответствующие (2.2):
u1 =u ( h ) ( x1 , y1 ) =a + bx1 + cy1 , u2 =u ( h ) ( x2 , y2 ) =a + bx2 + cy2 ,

(2.12)

u3 =
u ( h ) ( x3 , y3 ) =
a + bx3 + cy3 .

Решая эти уравнения относительно a, b, с и подставляя обратно в
(2.11), получаем
3

u ( h ) ( x, y ) = ∑ N j ( x, y )u j ,

(2.13)

j =1

где N j ( x, y ) ‒ функции формы, которые определяются выражением
)
N j ( x, y=

1
(a j + b j x + c j y ),
2A

в котором
29

где =
j 1, 2, 3 ,

(2.14)

1
1
и
=
A
1
2
1
В этих

x1 y1
1
x2=
y2
(b1c2 − b2c1 ) – площадь элемента.
2
x3 y3

уравнениях x j и y j ( j = 1, 2,3) обозначают координаты

вершин элемента.
Из (2.14) видно, что функции формы обладают следующим
важным свойством:
(2.15)
Другое важное свойство заключается в том, что N j ( x, y ) равно
нулю, если точка ( x, y ) находится на стороне элемента,
противоположной j-му узлу. Следовательно, величина u на стороне
элемента не связана со значением u в противоположном узле и двумя
точками рассматриваемой стороны. Это гарантирует непрерывность
решения при переходе через сторону элемента.
На рис. 11, б показана функция формы для треугольного элемента, которая также может быть выражена в локальных координатах, связанных с треугольником. Соединив произвольную точку P
с вершинами треугольника, получим три треугольника с площадями
A1, A2 , A3 (рис. 12).

Рис. 12
30

Тогда функции формы могут быть записаны в терминах координат, связанных с площадями:
N j ( x, y=
) N j ( L1, L2 , L3=
) L=
j 1, 2,3
j A j / A, =

Функции

(2.16)

L j удовлетворяют соотношению
L1 + L2 + L3 =
1.

(2.17)

Функции формы высшего порядка могут быть построены
подобным же образом. С другой стороны, они могут быть
сконструированы прямо в координатах, связанных с площадями:
N j ( L1 , L2 , L3) =
ψ II ( L1 )ψ JJ ( L2 )ψ K
K ( L3 ),

(2.18)

где ψ nk ( L j ) ( j = 1, 2,3) ‒ заданные полиномы Лагранжа, I + J + K =
M ‒
порядок полинома. Это соотношение может быть использовано для
создания функции формы любого порядка.
Пример. Получить функций формы из (2.18) для трех и шести
узловых треугольных элементов.
Решение. Для линейного треугольного элемента высший порядок полинома M=1. Из уравнения (2.8) в соответствии с рис. 12
получим
L1 − ( L1 )0
L −0
=1
=L
( L1 )1 − ( L1 )0 1 − 0 1

N1 ( L1 , L2 , L3 ) =N1 (1, 0, 0) =ψ11 =

Аналогично,

(2.19)

L2 − ( L2 )0
L2 − 0
ψ11 ( L2 ) =
=
=
N 2 ( L1 , L2 , L3 ) =
N 2 (0,1, 0) =
L2 ,
( L2 )1 − ( L2 )0 1 − 0
L3 − ( L3 )0
L3 − 0
N3 ( L1 , L2 , L3 ) =
N3 (0, 0,1) =
ψ11 ( L3 ) =
=
=
L .
( L3 )1 − ( L3 )0 1 − 0 3

Для квадратического треугольника с шестью узлами полином
имеет порядок М=2 (рис. 13). Это просто другой подход к той же самой проблеме.
Рассмотрим узел 1:
[ L − ( L1 )0 ][ L1 − ( L1 )1 ]
ψ (2, 0, 0) =
ψ 22 ( L1 ) = 1
=
N1 ( L1 , L2 , L3 ) =
[( L1 ) 2 − ( L1 )0 ][( L1 ) 2 − ( L1 )1 ]
( L1 − 0)( L1 − 0.5)
=
= L1 (2 L1 − 1).
(1 − 0)(1 − 0.5)
31

(2.20)

Для среднего узла, например 4:
[ L − ( L1 )0 ][ L2 − ( L2 )0 ]
ψ11 ( L1 )ψ11 ( L2 ) = 1
=
N 4 ( L1 , L2 , L3 ) =
N 4 (1,1, 0) =
[( L1 )1 − ( L1 )0 ][( L2 )1 − ( L2 )0 ] (2.21)
( L1 − 0)( L2 − 0)
= 4 L1L2 .
(1 − 0.5)(1 − 0.5)

Рис. 13

2.3. Четырехугольный элемент
Прямоугольный элемент имеет четыре узла. Неизвестная
функция u ( h) в пределах каждого элемента представляется в виде
u ( h) ( x, y ) =a + bx + cy + dxy ,

(2.22)

где a, b, c, d – некоторые постоянные, которые надо найти.
Чтобы получить функции формы N j , j = 1, 2,3, 4 в локальных
координатах, связанных с элементом, рассмотрим 2 × 2 квадрат
{(ξ, η) | −1 ≤ ξ, η ≤ 1} (рис. 14, а).

Рис. 14
32

Для простоты вершины элемента пронумерованы двойными
индексами (0,0), (1,0), (0,1) и (1,1), которые соответствуют узлам
1,2,3,4 (рис. 14, б). Функции формы могут быть получены в виде
n
,
N ( k ,l ) (ξ, η) =ψ m
k (ξ)ψ l (η), k , l = 0,1

(2.23)

n
где ψ m
k (ξ)ψ l (η) ‒ интерполяционные полиномы Лагранжа, m и n –
порядок полинома. Функции формы для четырехузлового элемента:

1
(1 − ξ)(1 − η),
4
1
N 2 (ξ, η) =ψ11 (ξ)ψ10 (η) = (1 + ξ)(1 − η),
(2.24)
4
1
N3 (ξ, η) = ψ11 (ξ)ψ11 (η) = (1 + ξ)(1 + η),
4
1
N 4 (ξ, η) =ψ10 (ξ)ψ11 (η) = (1 − ξ)(1 + η).
4
Этот прием может быть использован для построения функций
формы высшего порядка. Рассмотрим, например, 9-й узловой элемент
(рис. 15, б).
N1 (ξ, η) = ψ10 (ξ)ψ10 (η) =

Рис. 15

В каждом измерении имеется три узла, поэтому функции формы
будут иметь квадратичный порядок по ξ и по η :
N ( k ,l ) (ξ, η) =ψ km (ξ)ψ ln (η),

где k , l = 0,1, 2; m, n = 2.

Используя соотношение (2.8) и рис. 15, a, получим
33

(2.25)

(ξ − 0)(ξ − 1) (η − 0)(η − 1)
1
(η − 1)(ξ − 1)ηξ,
=
(−1 − 0)(−1 − 1) (−1 − 0)(−1 − 1) 4
(ξ − (ξ)0 )(ξ − (ξ) 2 )
(η − (η)1 ))(η − (η) 2 )
N5 (ξ, η)
=
(ξ)1 − (ξ)0 )(ξ)1 − (ξ) 2 ) (η)0 − (η)1 )(η)0 − (η) 2 )
N1 (ξ, η=
)

(ξ + 1)(ξ − 1) (η − 0)(η − 1)
1
=
=
− (1 − ξ2 )η(1 − η).
(0 + 1)(0 − 1) (−1 − 0)(−1 − 1)
2

Аналогично
1
(1 − η)(1 + ξ)ηξ, N3 (ξ, η=
)
4
1
)
(1 − ξ)(1 + η)ηξ, N 6 (ξ, η
)
N 4 (ξ, η
=
=
4
1
N 7 (ξ, η) = (1 − ξ2 )(1 + η)η, N8 (ξ, η) =
2
)
N 2 (ξ, η=

Эти функции изображены на рис. 16.

Рис. 16
34

1
(1 + η)(1 + ξ)ηξ,
4
1
(1 + ξ)(1 − η2 )ξ,
2
1
− (1 − ξ)(1 − η2 )ξ.
2

(2.26)

3. МЕТОД КОНЕЧНЫХ ЭЛЕМЕНТОВ ДЛЯ РЕШЕНИЯ
ОБЫКНОВЕННЫХ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ
3.1. Классическая и интегральная формулировка задачи
Рассмотрим применение метода на примере простейшей краевой
задачи:
−

d  du 
p =
 + q ⋅ u f , где 0 < x < 1,
dx  dx 

=
u (0) 0,=
u (1) 0.

(3.1)

Предположим, что p ∈ С1[0,1], q, f ∈ C [0,1] , причем
p ( x) ≥ c0 > 0; q ( x) ≥ 0; ∀x ∈ [0,1].

Обозначим через С0k [0,1] – множество функций из С k [0,1] и
равных нулю в граничных точках отрезка.
Умножим (3.1) на произвольную функцию v ∈ C01[0,1] , проинтегрируем полученное равенство по отрезку [0,1] и используем формулу
интегрирования по частям для преобразования интеграла, содержащего производные функции u. В результате получим


quv  dx ∫ fvdx
∫  p dx dx + =

1

1

du dv

0



∀v ∈ C01[0,1].

(3.2)

0

Используем тождество (3.2) для построения приближенного
решения задачи (3.1) методом Галёркина. Введем в рассмотрение
базисные или координатные функции N1 ( x), N 2 ( x), , N n ( x) ∈ С01[0,1] .
Приближенное решение задачи (3.1) будем искать в виде
n

(n)
u=
∑ c i ⋅Ni ( x),

(3.3)

i =1

определяя коэффициенты с1 , с2 ,  сn из системы линейных алгебраических уравнений:
1



du ( n ) dN

1

(n)
i
∫  p dx dx + qu Ni  dx = ∫ fNi dx,
0
0


i = 1, 2  , n

(3.4)

Учитывая (3.3), запишем систему более подробно:
n

kij c j b=
∑=
i , i 1, 2,  n,
j =1

35

(3.5)

где
1

dN dN j
=
+ qNi N j dx ,
kij ∫  p i
dx dx

0

1

bi = ∫ fNi dx .

Или в матричной форме:
,
где K
KC B=
=

(3.6)
(3.7)

0

kij )
B
(=
i , j =1
n

(b1 , b2 ,  bn ) т . .

(3.8)

Все действия, описанные выше, представляют собой процедуру
метода Галёркина. Метод конечных элементов предполагает
использование специальных базисных функций.
3.2. Дискретизация области, выбор базисных функций
Построим на отрезке [0,1] сетку ω= {0= x1 < x2 <  < xn= 1}
(вообще говоря, неравномерную). Ячейки сетки ei = [ xi ; xi +1 ],
=
i 1, 2,  , n − 1 , назовем конечными элементами.
Положим hi =xi +1 − xi , i =1, 2 , n − 1, h = max hi . Будем рас1≤i ≤ n −1

сматривать приближенное решение в виде линейной на каждом
конечном элементе функции, равной нулю при x = 0 и x = 1 . Тогда
имеет место выражение (3.3), где Ni – линейные на каждом элементе,
непрерывные на отрезке [0,1] функции, удовлетворяющие условию
1, i = j ,
Ni ( x j ) = 
1, 2, , n.
0, i ≠ j , i =

(3.9)

Рассмотрим, как определяются кусочно-линейные базисные
функции. Каждому внутреннему узлу x ставится в соответствие
кусочно-линейная функция
формы:

i

N i ( x) , называемая также функцией

x ≤ xi −1 ,
0,
x− x
i −1

, xi −1 ≤ x ≤ xi , =
i 2,  , n − 1,
 hi −1
N i ( x) = 
 xi +1 − x , x ≤ x ≤ x ,
i
i +1
 hi

x ≥ xi +1.
0,
36

(3.10)

Для граничных узлов x и xn
1
 x − xn−1
 x2 − x
, xn−1 ≤ x ≤ xn ,
, x1 ≤ x ≤ x2 ,


N1 ( x) =  h1
N n ( x) =  hn−1
0,
0,
x ≥ x2 ;
x ≤ xn−1.



(3.11)

На рис. 17 представлены базисные функции для области,
разбитой на три элемента.

Рис. 17

Пользуясь свойством функций Ni , легко показать, что они
линейно независимы.
Предположим обратное, т.е. существует нетривиальная линейная
комбинация этих функций, равная нулю:
(3.12)
с1 N1 ( x) + с2 N 2 ( x) +  сn N n ( x=
) 0, ∀x ∈ [0,1].
Подставив в равенство (3.12)
=
x x=
i , i 1, 2,  , n , в силу свойств
функций формы, получим
(3.13)
=
сi 0,=
i 1, 2,  , n,
что противоречит предположению.
Функции формы Ni ( x) иногда называют базисными функциями с
локальным носителем, поскольку множество точек, где Ni ( x) ≠ 0
(носитель функции Ni ( x) ), есть малый отрезок, представляющий
37

собой объединение соседних элементов. Как следствие этого, Ni −1 ( x)
и Ni +1 ( x) ортогональны. В целом базис не является ортогональным,

так как функции Ni ( x) и N j ( x) неортогональны, если i − j ≤ 1
(рис. 18). Поэтому матрица системы является трехдиагональной.

Рис. 18

Рассмотрим более подробно формирование разрешающей
системы линейных алгебраических уравнений на примере системы из
трех конечных элементов [0,1] =[0, 1 ] ∪  1 , 2  ∪ [ 2 ,1] . Будем считать,
3  3 3  3
что
=
p ( x) 1,=
q ( x) 0,=
f ( x) 1. Обозначим
1

dNi dN j
dx,
0 dx dx

( Ni′, N ′j ) = ∫

(3.14)

1

( f , Ni ) = ∫ fNi dx.
0

В силу свойств функций формы ci = u ( n ) ( xi ) разрешающая
система (3.8) имеет вид
0
0
 ( N1′, N1′ ) ( N1′, N 2′ )
  u1   ( f , N1 ) 
 ′ ′
  

′
′
′
′
(
N
,
N
)
(
N
,
N
)
(
N
,
N
)
0
2
2
2
3
 2 1
  u2  =  ( f , N 2 )  (3.15)

0
( N3′ , N 2′ ) ( N3′ , N3′ ) ( N3′ , N 4′ )   u3   ( f , N3 ) 

  

0
0
( N 4′ , N3′ ) ( N 4′ , N 4′ )   u4   ( f , N 4 ) 

38

В этом выражении частично
учтена область определения
функций Ni . Например, ( N1′, N3′ ) = 0 , так как области определения
этих функций не пересекаются D( N1 ) ∩ D( N3 ) =
∅ . Вместе с тем
некоторые интегралы имеют промежуток интегрирования, включающий соседние элементы, например:
2
3

1

2

3 dN dN
3 dN dN
dN 2 dN 2
2 dx
2 dx .
( N 2′ , N 2′ ) ∫ =
dx ∫ 2
=
+∫ 2
dx
dx
dx
dx
dx
dx
1
0
0

(3.16)

3

Если область интегрирования относится к элементу с номером
е = 1,2,3, обозначим его соответствующим индексом. Введем
обозначения:
1
3

1
3

0

0

(1)
N1′, N1′ ) ∫=
N1′N1′dx k11
N1′, N 2′ ) ∫ =
N1′N 2′ dx
(=
; (=
1
3

1
3

( N=
2′ , N1′ ) ∫ N=
2′ N1′dx

(1)
k21
f , N1 ) ∫=
fN1dx
; (=

0

(1)
k12
;

. (3.17)

b1(1) ;

0

1
3

2
3

0

1
3

(1)
(2)
+ k11
k22
( N 2′ , N 2′ ) =
;
∫ N 2′ N 2′ dx + ∫ N 2′ N 2′ dx =

1
3

2
3

0

1
3

( f , N 2 ) = ∫ fN 2 dx + ∫ fN 2 = b2(1) + b1(2) ;
2
3

=
( N3′ , N 2′ ) ∫=
N3′ N 2′ dx
1
3

(2)
k21
;

2
3

1

1
3

2
3

2
3

=
( N 2′ , N3′ ) ∫=
N3′ N 2′ dx

( N3′ , N3′ ) =
∫ N3′ N3′ dx + ∫ N3′ N3′ dx

1
3

(2)
=
k22

(3)
+ k11
;

(2)
k12
;
2
3

1

1
3

2
3

( f , N3 ) =
∫ fN3dx + ∫ fN3 =

В новых обозначениях система (3.15) имеет вид
39

 k (1)
 11
(1)
 k21

 0

 0

(1)
k12

0

(1)
(2)
k22
+ k11

(2)
k12

(2)
k21

0

(3)
+ k11
(3)
k21

(2)
k22

0   u   b1(1) ) 
 1


0   u2   b2(1) + b1(2) 
=


(3)  u 
(2)
(3)
k12
 3
 b2 + b1 


(3) 
  u4   b(3) 
k22



2

(3.18)

или
KU = B,

(3.19)

где
 k (1)
 11
 (1)
K = k21
 0
 0


(1)
k12
(1)
k22

0
0

0 0 0 0
 
(2)
0 0  +  0 k11
 
(2)
0 0   0 k21
0 0   0 0

0
(2)
k12
(2)
k22

0

0 0
 
0 0
+
0 0
0   0

0

0

0

0

(3)
0 k11
(3)
0 k21

0 

0 
(3)  ,
k12

(3) 
k22


(3.20)

 b(1)   0   0 
 1   (2)   0 

 (1)   b  
=
F  b2  +  1  +  (1)  .
(2)
b
 0   b2   3 
 0   0   b(3) 

 
  4 

Матрица K (e)

 k (e) k (e) 
 b(e) 
11
12
=
 и вектор  1  , e 1, 2,3, (0.1)
 k (e) k (e) 
 b(e) 
 21
 2 
22 

соответственно называются матрицей жесткости и вектором нагрузки
элемента. Сам процесс сборки глобальной матрицы системы и
вектора нагрузки называется аccемблированием.
3.3. Формирование дискретной конечноэлементной схемы
Составим матрицу жесткости и вектор нагрузки для элемента
[ xi , x j ] . В соответствии с (3.9)
=
Ni

xj − x
x − xi
,Nj
.
=
h
h

Вычисляя составляющие матрицы жесткости, получаем
40

(3.22)

xj

(e)
= ∫ N1′N1′dx =
k11
xi

x

j
1
1
(e)
(e)
= k12
= ∫ N1′N 2′ dx = − ,
, k21
h
h
xi

xj

(3.23)

1
,
N 2′ N 2′ dx
∫=
h
xi

(e)
=
k22

1  1 −1
K (e) = 
.
h  −1 1 

Составляющие вектора нагрузки:
b1(e)

xj

xj − x
( x j − x)2
h
=
N
dx
=
dx
=
−
=
,
∫ i
∫ h
2h
2
xi
xi
xj

xj

xi

xj

xj

xi

xi

=
b2(e) ∫=
N j dx ∫

(3.24)

xj

x − xi
( x − xi ) 2
h
=
dx
=
.
h
2h
2
xi

Таким образом,
 1 −1
  u1 

u 
1  −1 2 −1
 2  =
−1 2 −1  u3 
h

 
−1 1   u 4 


1
 
h  2 .
2  2
 
1

(3.25)

Далее учтем, что u1 = u4 = 0:
1
  u1 
0

u 
 
2
−
1

  2  = h2  1  .

  u3 
1
−1 2

 
 
1   u4 

 0

(3.26)

В координатном виде
u1

2u2


−u2




=0
−u3

=
h2 .

+2u3

=
h
=0
2

u4

Отсюда
41

(3.27)

=
u1 0,=
u2 h 2=
, u3 h 2=
, u4 0.

(3.28)

Точное решение задачи:
x x2
− .
2 2
Результаты вычислений представлены в табл. 3
u ( x)=

(3.29)

Таблица 3
x

u ( x)

u

0

0

0

1
3
2
3

1
9
1
9

1
9
1
9

1

0

0

Отметим некоторые свойства матрицы K согласно (3.6):
1. Матрица K обладает свойством симметрии. Действительно,
1

dN dN j
=
+ qNi N j dx
kij ∫  p i
dx dx

0

⇔

=
kij k ji. .

2. Согласно правилу (3.9) функция N j ( xi=
) 0, ∀i ≠ j и, как
следствие этого,
1

dN dN j
=
+=
kij ∫  p i
qNi N j dx 0,
dx dx

0

если j > i + 1

Поэтому матрица является разряженной, т.е. содержит большое
число нулей и имеет ленточную структуру.
3. Матрица K обладает свойством положительной определённости, что гарантирует однозначную разрешимость системы
линейных алгебраических уравнений.
Определение. Квадратная матрица K размера n × n называется
положительно определённой, если C T KC ≥ 0, для любого вектора С
размера n ×1 , C T KC = 0 ⇔ С = 0 .
Теорема Матрица K из соотношения (3.6) является положительно определённой.
Пусть С ( с1 , с1 , , сn ) – произвольный вектор размера n ×1 .
Составим функцию
42

n

u = ∑ ci Ni .

(3.30)

i =1

Далее
n n

n n

n

n

j= 1 i= 1

j= 1 i= 1

i= 1

j= 1

=
С т КС ∑
=
( Ni , N j )c j a ( ∑ ci N=
∑ ci Kij c j ∑ ∑ ci a=
i , ∑ c j N j ) a (u , u ),

где
1
 du dv

(3.31)
=
+ quv dx
a (u , v) ∫  p
dx
dx


0
‒ cимметричная и положительно определённая билинейная форма
(см. (1.8)). Отсюда
т
(3.32)
С=
КС a (u , u ) ≥ 0.,

а

из

условия

a (u , u ) = 0

u =0,

т.е.

для

любого

x ∈ [0,1]

n

u ( x) ∑
c j N j ( x) 0 . Это равносильно
=
=
j =1

с=
= с=
1 с=
2 
n 0,

(3.33)

так как функции N j ( x) линейно независимы.
Пользуясь свойствами матрицы K системы линейных алгебраических уравнений (3.8), докажем, что она имеет единственное решение.
Предположим противное, т.е. ∃u1 , u2 , так что Ku1 = B и Ku2 = B.
Пусть C= u1 − u2 , тогда КС =0 ⇒ С T KC =0 ⇒ C =0 ⇔ u1 =u2 . Это
обстоятельство при решении линейных алгебраических систем позволяет хранить в памяти компьютера лишь половину ленты (с учетом
симметрии матрицы).
Ширина ленты матрицы зависит от способа нумерации узлов сетки, и методы решения линейных алгебраических уравнений, учитывающие ленточную структуру, должны учитывать это обстоятельство.
Рассмотрим для примера случаи а и б с разными способами
нумерации (рис. 19).

Рис. 19
43

Если J > I + 1 , ненулевые части Ni и N j не пересекаются.
В случае нумерации а получим систему с трехдиагональной
матрицей (рис. 20). В случае нумерации б, которая отличается от а
тем, что второй узел имеет номер n − 1 , а предпоследний ‒ номер 2,
ширина ленты матрицы равна n − 1 . Этот случай изображен на рис. 21.

Рис. 20

В общем случае расстояние от главной диагонали до последней
ненулевой диагонали, параллельной главной, можно определить по
формуле
(3.34)
W= R + 1 ,
где R – максимальная по элементам величина разности между
номерами узлов в элементе.
44

Рис. 21

Как следует из рассмотренных выше примеров, основными
моментами метода конечных элементов являются:
1) дискретизация области, т.е. представление её в виде объединения конечных элементов;
2) вычисление матриц элементов и вектора нагрузки;
3) ассемблирование этих матриц и вектора нагрузки в глобальную матрицу и вектор;
4) учет граничных условий;
5) решение системы линейных алгебраических уравнений;
6) графическое отображение результатов решения задачи;
Эти этапы всегда присутствуют в методе конечных элементов
(хотя в одномерном случае дискретизация может быть элементарной).
Рассмотрим подробнее эти шаги и приведем программу Matlab
для решения основных рассматриваемых задач. Будем считать, что
область Ω =[a, b] уже разбита на конечные элементы, а их граничные
точки (узлы сетки) каким-либо образом пронумерованы.
Рассмотрим, например, элемент, лежащий между узлами i и j ,
т.е. [ xi , x j ] . Считая, что приближенное решение есть линейная функция от х на элементе u= c1 + c2 x , запишем это выражение в иной
45

форме, выразив c1 и c2 через значения функции u в граничных точках
элемента:
u ( xi ) =
c1 + c2 xi =
ui ,
(3.35)
u ( x j ) =+
c1 c2 x j =
u j.
Из этих условий получим
с1 =
с2 =

u j − ui
x j − xi

,

ui x j − u j xi
x j − xi

(3.36)
.

Тогда искомое решение можно записать в виде
=
u ( x) Ni ( x)ui + N j ( x)u j ,

(3.37)

где
Ni ( x) =

xj − x

,
h
x − xi
N j ( x) =
,
h
h x j − xi .
=

Рис. 22

(3.38)

Функции формы элемента представлены
на рис. 22.
Приступая к выводу основных уравнений, в
качестве исходной будем считать интегральную
формулировку (3.2). Представляя интегралы,
входящие в (3.2), в виде суммы интегралов по
элементам ( N e – число элементов), получим
(e)
Ne x j

 du dv



+ quv − vf  dx =
0.
∑ ∫ p

e =1 x( e )  dx dx

(3.39)

i

Представим функции u и ν в виде разложения по базисным функциям:
n

n

=
u ( x) ∑
=
ui Ni ( x), v( x) ∑ vi Ni ( x).
=i 1 =i 1

(3.40)

В соответствии с определением базисных функций в пределах
элемента с номером е
u ( x) =
ui Ni ( x) + u j N j ( x), v( x) =+
vi Ni ( x) v j N j ( x).
46

(3.41)

Введем матричные обозначения:
du
=
dx

т

N ]{U }, v( x)
[=

=
u ( x)

{V }т [ N ] ,
т

(3.42)

dN j
dNi

 0],
dx
dx

(3.43)

dv
dx

N x′ ]{U },
[=

{V }т [ N x′ ] ,

где
[ N ] = [0  Ni
[ N x′ ] = [0 

 N j  0],

и
 u1 
 v1 
 
 
u
v
2
=
{U } =
,
{V }  2  .
(3.44)
  
  
 
 
 un 
 vn 
Тогда подынтегральную функцию в (3.39) можно записать в виде
Ne x j



+ quv − fv dx =
∑ ∫p

e =1 xi  dx dx
du dv

Ne  x j



т
т
= {V }т ∑  ∫ p [ N x′ ] [ N x′ ] + q [ N ] [ N ] dx {U } − {
e =1 
 xi



)

(

Ne  x j


т
} − {V }т ∑  ∫ f [ N ] dx .
e =1 
 xi


(3.45)

т,

Отсюда, в силу произвольности {V }

K {U } = {B} ,

(3.46)

где
Ne

Ne

K ∑
K (e) ;
B ∑ B (e) ;
=
=
=e 1 =e 1
xj
т
(e)
K
=
∫ p N x′ N x′ + q
xi

(

xj

[ ][ ]
т

B (e) = ∫ f [ N ] dx.
xi

47

[ N ]т [ N ]) dx;

(3.47)

Матрица K (e) и столбец B (e) представляют собой матрицы
жесткости и вектор нагрузки элемента.
Выполняя действия с матрицами в (3.47), можно получить

K (e)

  
0 
 0


(e)
(e)
kii
 kij
 
 


B (e)

  
 ,
=
 
k (jie)  k (jje)
 


 0



0 


  
 (e) 
 bi 
  


 b(je) 


  



(3.48)

где
kii(e)

 dN 2

2
i
=
∫ p  dx  + qNi dx,


xi 

xj

k (jje)

 dN 2

j
2
=
∫ p  dx  + qN j dx,


xi 


xj

xj
xj
(3.49)
 dNi dN j

 dN j dNi

(e)
,
kij(e) =
p
qN
N
dx
k
+
=
∫  dx dx
∫ p  dx dx + qN j Ni dx,
i j
ji


xi 
xi 

b ( e )i

xj

xj

xi

xi

fNi dx,
b(e) j ∫ fN j dx.
∫=

3.4. Примеры решения задач
Рассмотрим задачу
 d 2u
u 1, 0 < x < 1 ,
 2 +=
 dx
u (0) 0,=
=
u (1) 0.


(3.50)

Прежде всего запишем интегральную формулировку для данного
уравнения с помощью метода взвешенных невязок. Пусть функция
v( x) удовлетворяет однородным условиям
=
v(0) 0=
и v(1) 0 , u ( x)
– приближенное решение.
2
Умножая невязку 1 − d u − u на функцию v( x) и интегрируя по
dx 2
промежутку [0,1], получаем
48


d 2u 
v
1
−
u
−
0.
(3.51)

 dx =
∫ 
dx 2 
0 
Интегрируя первое слагаемое по частям, с учетом того, что
=
v(0) 0=
и v(1) 0 , приходим к уравнению
1

1

du dv

1

1

0

0

0.
∫ dx dx dx − ∫ vudx + ∫ vdx =

0

(3.52)

Составим матрицу и вектор элемента. Для этого воспользуемся
формулами (3.49), приняв p = 1, q = −1 :
xj
 dNi dN j

(e)
(e)
k=
k
=
∫  dx dx − Ni N j dx,
i, j
j ,i

xi 
2
2
 dNi 

∫   dx  − ( Ni ) dx,

xi  

xj − x
dN j 1
x − xi
dNi
1
,
,Nj =
,
Ni =
=
− ,
=
hi
hi
dx
hi
dx
hi
(e)
k=
i ,i

(e)
k=
j, j

xj 

xj

x
xj
2
 dNi dN j 
1 j
1
1
 dNi 
,
.
dx
dx
=
−
=
−
∫  dx dx 
∫  dx  dx =
2 ∫
h
h
hi xi

i
i

xi 
xi 

xj
xj
xj
x j − x x − xi
hi
hi
2
,
N
N
dx
dx
=
=
=
∫ i j
∫ h
∫ ( Ni ) dx 3 .,
6
h
i
i
xi
xi
xi

ki(,ei ) =
k (je, )j =
hi −

hi
,
3

ki(,ej) =
k (je,i) =
−

1 hi
− .
hi 6

Вычислим вклад элемента V e = ( xi , x j ) в вектор правых частей
(3.46):
xj
xj
( x − x j ) 2 x j hi
xj − x
,
bi(e) =
N
dx
==
dx
−
=
∫ i
∫ h
2
hi
xi
i
xi
xi
xj

xj

xi

xi

=
b(je) ∫=
N j dx ∫

( x − xi ) 2 x j hi
x − xi
.
=
dx
=
2
hi
hi
xi

Следует отметить, что в точке 0 и 1 не равны нулю только
базисные функции,=
т.е. N1 (0) 1=
и N n (1) 1 . Таким образом, вектор
и матрица для элемента V e = ( xi , x j ) имеют вид
49

Ke

0

0


ki ,i
ki , j


0
0
=

0
0


k j ,i
k j, j


0
0 



0

1 hi

−

hi 3

0
0

1 h

− − i

hi 6

0

0




0
, B (e)
=
0

1 hi

−

hi 3

0 



1 h
− − i
hi 6

0

 − hi
 2
0

0
 h
− i
 2
0







.







Вычисляя компоненты матрицы элемента, можно простым
суммированием получить матрицу всей системы. Процесс формирования матрицы системы и глобального вектора правых частей в
МКЭ назывется аcсемблированием (или сборкой). Матрица системы
называется матрицей жесткости.
Запишем, например, систему для трех элементов и четырех узлов.
Положим, что длина каждого элемента V e = ( xi , x j ) равна h:

K

1 h
 1 h

0
0 
− −
 h−3
h 6


 − 1 − h 2 1 − h  − 1 − h

0 


 h 6
h 6
h 3
=
 , B (e)
1 h
1 h

1 h
2 −  − − 
− −
 0
h 6
h 6
h 3


1 h
1 h
1 h 

− −
− −
− 
 0
h 6
h 6
h 3 


h 

 −2 


− h − h 
 2 2


 − h − h .
 2 2




h 

 −

2 


Если отрезок [0,1], где ищется решениe, разбит всего на три
элемента, то h = 1/ 3 :

 26
 9
55
 −18





−55
18
52
9
−55
18

−55
18
52
9
−55
18

  u   −61 
 1   
  u2   −31 
=
−55   u   −1 
18   3   3 
26   u4   −1 
9 
6

(3.53)

Учитывая граничные условия, получаем
1







52
9
−55
18

−55
18
52
9

  u1   0 
    −1 
  u2  =  3  .
  u   −1 
 3   3 
1   u4   0 
50

(3.54)

Запишем эту систему в координатном виде:
6
−104u2 + 55u3 =
6

−
=
−0,1224, u1 =
u3 =
u4 =
0.
55u2 − 104u3 =6 ⇔ u2 =
49
u= u= 0
4
 1
Точное решение задачи:

(3.55)

cos1 − 1
1
2
sin x, u ( ) =
u( ) =
-0.1237.
sin1
3
3
Далее рассмотрим формирование матриц элементов в процессе
решения задач.
Для примера решим следующую задачу:
u ( x) =
1 − cos x +

 2 d 2u
du
− 2x =
− 4u x 2 , 10 < x < 20
x
dx
 dx 2
u (10) 0,=
=
u (20) 100.


(3.56)

Умножая уравнение (3.56) на некоторую функцию ν, удовлетворяющую условиям
=
v(10) 0,=
v(20) 0 , и интегрируя, получим
слабую формулировку задачи:
x = 20

20
du
du 
 2 dv du


.
− ∫ x 2vdx +  x 2v 
∫  x dx dx + 4 xv dx + 4uv dx =
dx  x =10


10 
10
20

(3.57)

В силу условий на функцию v на границе интервала последнее
слагаемое в правой части (3.57) несущественно.
Для вычисления матриц элементов можно использовать символьные вычисления. Следующий фрагмент программы вычисляет
матрицы элемента [ xi , x j ] :
function main
clear
syms i j he x xi xj K k11 k12 k21 k22
phi1=(xj - x )/he; %функция формы N i
phi2=(x - xi)/he; %функция формы N j
he=sym('he','real'); % шаг сетки
xi=sym('xi','real');
xj=sym('xj','real');
x=sym('x','real');
p=sym('x^2');

%Определение функции
51

q=sym('4*x');
b=sym('4');
c=sym('-x^2');
K=[k11 k12; k21 k22];
%матрица элемента
H=[phi1 phi2];
%вектор нагрузки
H';
H1=[diff(phi1) diff(phi2)];
H1';
K1=int(p*H1'*H1,x,xi,xj); %Вычисление интеграла
K2=int(q*H'*H1,x,xi,xj); %Вычисление интеграла
K3=int(b*H'*H,x,xi,xj); %Вычисление интеграла
K=(K1+K2+K3);
K=simplify(simplify(K))
F=int(c*H',x,xi,xj);
F=simplify(simplify(F));
end
Результаты её выполнения:
2
3
3
2
1  x j − 4 x j xi + 6 x j xi + xi
h 2  − xi 3 + 2 xi 2 x j − x j 3 )

4
3
4
1  −3 xi + 4 xi x j − x j 
=
.
12h  − xi 4 + 4 xi x j 3 − 3 x j 4 




,
− xi 3 + 4 xi 2 x j − 6 xi x j 2 + 3 x j 3  (3.58)
xi 3 − 2 xi x j 2 + x j 3

K (e) =

B (e)

Листинг программы, реализующей алгоритм
конечных элементов для решения этой задачи
function f2em1D()
% Пример решения задачи
 2 d 2u
du
− 2x =
− 4u x 2 , 10 < x < 20
% x
2
dx
 dx
u (10) 0,=
=
u (20) 100.

% используется 10 элементов
% Описание переменных
% k = матрица элемента
% f = вектор элемента
% kk = глобальная матрица
% ff = глобальный вектор
% index = номера узлов сетки связанные с элементом
52

метода

% bcdof = номера узлов связанных с граничными условиями
% bcval = значения величин граничных условий
%связанных с узлами сетки в bcdof.
%Вычисление матрицы и вектора элемента
clear
syms i j he x xi xj K k11 k12 k21 k22
phi1=(xj - x )/he;
phi2=(x - xi)/he;
he=sym('he','real');
xi=sym('xi','real');
xj=sym('xj','real');
x=sym('x','real');
p=sym('x^2'); %Определение функции
q=sym('4*x');
b=sym('4');
c=sym('-x^2');
K=[k11 k12; k21 k22];
H=[phi1 phi2];
H';
H1=[diff(phi1) diff(phi2)];
H1';
K1=int(p*H1'*H1,x,xi,xj); %Вычисление интеграла
K2=int(q*H'*H1,x,xi,xj); %Вычисление интеграла
K3=int(b*H'*H,x,xi,xj); %Вычисление интеграла
KE=(K1+K2+K3);
KE=simplify(simplify(KE))%матрица элемента в символьном виде
FE=int(c*H',x,xi,xj); %вектор элемента в символьном виде
FE=simplify(simplify(FE));
% Ввод исходных данных
ndof=1;
nnode=11;% число узлов
nel=10; % число элементов
nnel=2; % число узлов на элемент
ndof=1; % number of неизвестных на узел
sdof=nnode*ndof; % полное число узлов
% ввод координат узлов сетки
gcoord(1)=10; gcoord(2)=11; gcoord(3)=12; gcoord( 4)=13;
gcoord(5)=14;gcoord(6)=15; gcoord(7)=16; gcoord(8)=17;
gcoord(9)=18; gcoord(10)=19; gcoord(11)=20;
% ввод данных о связи между элементами
nodes(1,1)=1; nodes(1,2)=2; nodes(2,1)=2; nodes(2,2)=3;
53

nodes(3,1)=3;nodes(3,2)=4; nodes(4,1)=4; nodes(4,2)=5;
nodes(5,1)=5; nodes(5,2)=6;nodes(6,1)=6; nodes(6,2)=7;
nodes(7,1)=7; nodes(7,2)=8;nodes(8,1)=8;nodes(8,2)=9;
nodes(9,1)=9; nodes(9,2)=10; nodes(10,1)=10;
nodes(10,2)=11;
% ввод данных для граничных условий
bcdof(1)=1; bcval(1)=0;
bcdof(2)=11; bcval(2)=100;
% инициализация матриц и векторов
ff=zeros( sdof, 1);
kk=zeros( sdof,sdof );
index=zeros( nnel *ndof, 1);
% вычисление матрицы и вектора элемента и их сборка
for iel=1:nel % цикл по всем элементам
nl=nodes(iel,1); nr=nodes(iel,2);
xl=gcoord(nl); xr=gcoord(nr);
eleng=xr-xl;
index=feeldofl (iel,nnel,ndof);
xi=xl;
xj=xr;
he=eleng;
k=feode21(subs(KE)) % вычисление матрицы элемента
f=fefll(subs(FE)) % вычисление вектора элемента
[kk,ff]=feasmbl2(kk,ff,k,f,index); % сборка матриц и векторов
элементов
end % конец цикла по элемента
% применение главных граничных условий
[kk,ff]=feaplyc2 (kk,ff, bcdof, bcval );
full (kk)
full(ff)
% решение матричного уравнения
fsol=kk\ff;
% аналитическое решение
for i=1:nnode
x=gcoord(i);
esol(i)=0.00102*x^4-0.16667*x^2+64.5187 /x;
end
%plot(x,esol,'r')
figure(1);
plot(gcoord,esol,'b--o',gcoord,fsol,'r'),grid;
%
54

title('\fontsize{14} Метод конечных элементов');
hleg1=legend('Приближенное решение','Точное решение')
set(hleg1,'Location','NorthWest')
set(hleg1,'Interpreter','none')
% печать точного и конечноэлементного решения
num=1:1:sdof;
results=[num' fsol esol']
function [kk,ff]=feaplyc2(kk,ff,bcdof,bcval)
%--------------------------% Цель:
% Применение граничных условий в матричном уравнении
[kk]x=ff
%
% [kk,ff]=feaplybc(kk,ff,bcdof,bcval)
%
% Описание переменных:
% kk – матрица системы перед применением условий
% ff – вектор перед применением условий
% bcdof – вектор, содержащий номера граничных узлов
% bcval – вектор, содержащий значения величин заданных на
границе
%
% Для примера граничные узлы 1 and 11
% и граничные значения функции О.О и 100,
% соответственно.Тогда, bcdof(1)=1 and bcdof(2)=11; и
% bcval(l)=0.O and bcval(2)=100.
n=length(bcdof);
sdof=size(kk);
%
for i=1:n
c=bcdof(i);
for j=1:sdof
kk(c,j)=0;
end
%
kk(c,c)=1;
ff(c)=bcval(i);
end
function [kk,ff]=feasmbl2 (kk,ff,k,f,index)
%------------------------55

% Цель:
% Сборка матриц элементов в глобальную матрицу системы
% Сборка векторов элементов в глобальный вектор системы
%
% [kk,ff] =feasmbl2 (kk,ff ,k,f,index)
%
% Описание переменных::
% kk – матрица системы
% ff – вектор системы
% k – матрица элемента
% f – вектор элемента
% index – вектор, содержащий узлы связанные с элементами
%------------------------%
edof = length(index);
for i=1:edof
ii=index(i);
ff(ii)=ff(ii)+f(i);
for j=1:edof
jj=index(j);
kk(ii,jj)=kk(ii,jj)+k(i,j);
end
end
function [index]=feeldofl(iel,nnel,ndof)
% Цель:
% Вычисление узлов связанных с каждым элементом
%
% [index]=feeldofl(iel,nnel,ndof)
%
% Описание переменных:
% index - глобальный вектор узлов, связанных с элементом
% iel - element number whose system dofs are to Ье determined
% nnel – число узлов на элемент
% ndof - number of dofs per node
%,-------------------------%
edof = nnel*ndof;
start = (iel-1 )*(nnel-1 )*ndof;
%
for i=1:edof
index(i )=start+i;
56

end
function [f]=fefll(FE)
%
% Цель:
% вычисление вектора элемента
%
% [f]=fefll(xl,xr)
%
% Описание переменных:
% f – матрица элемента
% FЕ –вектор полученный замещением символов значениями
для элемента
f=FE;
function [k]=feode21( KE)
%,------------------------% Цель:
% вычисление матрицы элемента
%
% [k]=feode2l( КЕ)
% Описание переменных:
% k – матрица элемента
% КЕ –матрица полученная замещением символов значениями
для элемента
k=KE ;
На рис. 23 представлены данные, полученные в результате
решения задачи. Сопоставление точного и конечноэлементного
решения приведено в табл. 4.
Таблица 4
№ узлов
1
2
3
4
5
6
7
8
9
10
11

Конечноэлементное
решение
0.0000
0.6046
2.4650
5.8421
11.0278
18.3432
28.1371
40.7850
56.6888
76.2761
100.0000
57

Аналитическое решение
0.0151
0.6321
2.5268
5.9280
11.1255
18.4380
28.2116
40.8190
56.6588
76.1553
99.7579

Метод конечных элементов
100
Приближенное решение
Точное решение
80

60

40

20

0

-20
10

11

12

13

14

15

16

17

18

19

20

Рис. 23

Далее рассмотрим
граничными условиями:

задачу

d 2u

с

главными

и

естественными

du
+ 3 xu =
8sin x;
dx
dx
du
u (2) = 0,
(5) = −2.
dx
Слабая формулировка задачи (3.59):
−

2

+4

(3.59)



∫  x dx dx + 5v dx + 3xuv=
dx 8∫ v sin xdx + 5 (5).
dx

5

2

dv du

1

du



du

(3.60)

0

Необходимая модификация программы заключается в изменении
данных, относящихся к коэффициентам уравнения, и учете граничных условий. Использовалось разбиение области на десять конечных
элементов. Матрица и вектор элемента, вычисленные с помощью
программы:
2
2
2
3
2
2
2
3
1  −10h + 2( xi − x j ) − h (3 xi + x j ) 10h + 2( xi − x j ) + ( xi + x j )h ) 
,
4h 2  2( xi 2 − x j 2 ) − 10h 2 − ( xi + x j )h3 10h 2 − 2( xi 2 − x j 2 ) + h3 ( xi + 3x j )) 


8  sin(xi ) - sin(x j ) - xi cos(xi ) + x jcos(xi ) 
= 
.
h  -sin(xi ) + sin(x j ) + xi cos(x j ) - x jcos(x j ) 
58

K (e) =

B (e)

В результате сборки получена следующая матрица K и вектор B .
Последний элемент вектора B содержит слагаемое 5u ′(5) . Для
выполнения граничного условия подставляются данные u ′(5) = −2 :
1


 -9.3442   16.7133 -5.2992

-10.2992 18.8933

-11.2542













-6.2542
24.3461
15.8507

 u1 
0
 u  

 2  

1.7763
 u  

 3  

1.2280
 u  

-10.8507
0.5481
 4  



 u5 

26.2146 -7.9095
-0.0104
 


-12.9095 24.8520 -8.8645
-0.6546
 u 6 =






-13.8645 27.0320 -9.8195
-1.3022
u
 7  

 u  

-14.8195 29.2120 -10.7745
-1.8335
 8  

-15.7745 28.7196 -8.5115
-2.5156


 u


9 

-13.5115 31.2175 -12.9392 
-2.6917
 u  

-17.9392 20.1442  10   -1.1760+5u′(5) 
u 
 11 

Результаты сопоставления
решения приведены на рис. 24.

точного

и

конечноэлементного

Метод конечных элементов
0.4
Приближенное решение
Точное решение

0.2
0
-0.2
-0.4
-0.6
-0.8
-1
-1.2
-1.4

2

2.5

3

3.5

4

4.5

5

Рис. 8

Если дифференциальный оператор краевой задачи является
самосопряженным, то система KU = B., получаемая в MKЭ, является
системой с симметричной матрицей. В рассмотренных выше примерах для учета главного граничного условия, например u(0) = 1,
первое уравнение в системе заменялось на соответствующее условие
59

в узле u1 = 1. Это приводит к нарушению симметрии матрицы. Для
сохранения симметрии используется прием, который рассмотрим на
следующем примере.
Пусть требуется учесть что u1 = 150 и u5 = 40:
−46u2
 55u1
−46u +140u
1
2

4
u
46
u
−
 1
2




+4u3
−46u3
+110u3
−46u3
4u3

−46u4
+142u4
−46u4

+4u5
−46u5
+65u5

= 500
= 2000
.
= 1000
= 2000
= 900

На первом этапе приравняем нулю все внедиагональные
элементы в первом и пятом уравнении, правые части уравнений в
этих строках заменяем соответственно на 55u1 и 65u5 . Тогда получим
 55u1
−46u +140u
1
2

4
u
46
u
−
 1
2




−46u3
+110u3
−46u3

−46u4
+142u4

+4u5
−46u5
65u5

= 8250
= 2000
.
= 1000
= 2000
= 2600

Второй этап состоит в исключении столбцов, которые
умножаются на u1 или u5.
Завершая этот шаг, получим систему с симметричной матрицей:
55u1

140u2

−46u2





−46u3
+110u3
−46u3

= 8250
= 8900
= 240
= 3840

−46u4
+142u4
65u5

= 2600

3.5. Сходимость метода
Рассмотрим далее вопрос о сходимости метода. Предварительно
получим некоторые вспомогательные результаты.
Лемма 1. Пусть u ( x) – непрерывная и кусочно-непрерывно
дифференцируемая на отрезке [0,1] функция, u (0) = 0 . Тогда
60

2
C

u

где u

С

1

≤ ∫ u ′x2 dx,

(3.61)

0

= max | u ( x) | .
0≤ x≤1

Доказательство. По формуле Ньютона‒Лейбница получим
1

u ( x=
) ∫ u ′(ξ)d ξ, для ∀x ∈ [0,1].

(3.62)

0

Далее используем известное неравенство Коши‒Буняковского:
b

∫ f ⋅ g dx

≤

b

b

a

a

2
2
∫ f dx ∫ g dx ,

a

(3.63)

2

x
x
1
x

u 2 ( x)=  ∫ 1 ⋅ u ′(ξ)d ξ  ≤ ∫ 12 d ξ ∫ u ′2 d ξ ≤ ∫ u ′2 d ξ.
0
0
0
0

Следовательно,

u

2
C

(3.64)

1

≤ ∫ u ′2 (ξ)d ξ .

(3.65)

0

Лемма 2. Пусть u ( x) ∈ C02 [0,1] ,
n

=
u I ( x) ∑ u (xi )ϕi ( x)
i =1

(3.66)

Если кусочно-линейная функция совпадает в узлах сетки с
функцией u(x), то
u − uI

C

≤

h2
u ′′ C
8

(3.67)

где h max( xi +1 − xi ) .
=
Доказательство. Пусть x ∈ [ xi , xi +1 ], для i =
1, 2,  , n − 1 . Используя представление для остаточного члена интерполяционного
многочлена в форме Лагранжа, получаем
u ′′(ξ)
u ( x) −=
uI ( x)
( x − xi )( x − xi +1 ), ξ ∈ [ xi , xi +1 ].
(3.68)
2
Прибавляя и отнимая в каждой скобке выражения ( x − xi )( xi +1 − x)
величину xi + xi +1 , после несложных преобразований получаем
2
61

2

2

x +x  h
h 
( x − xi )( xi +1 − x=
)   −  x − i i +1  ≤  
2
2 
 2

2

(3.69)

n

v = ∑ vi Ni ( x) ,

(3.70)

i =1

откуда и следует утверждение леммы.
Теорема. Пусть u – решение задачи (3.1), u – приближённое
решение, полученное по методу конечных элементов. Тогда
u − u

≤ c h 2 u ′′ C ,

C

(3.71)

где с – положительная постоянная, зависящая только от коэффициентов (3.1).
Доказательство. Обозначим через W линейное множество всех
функций:
n

где x ∈ [0,1] .

=
v ∑ vi Ni ( x),
i =1

(3.72)

Покажем, что если u ( x) ‒ решение задачи (3.1), то
1

1

0

0

 ,
a (u , v ) ≡ ∫ ( pu ′v′ +=
quv )dx ∫ fvdx

∀v ∈ V .

(3.73)

Действительно, интегрируя первое слагаемое по частям и используя
v (0)
= v=
(1) 0 , получим
1

1

0

0


 ,
a (u , v ) ≡ ∫ [−( pu ′)′ +=
qu ]vdx
∫ fvdx

Значит,

1

=
a (u , v )

 ,
∫ fvdx

∀v ∈ V .

∀v ∈ V .

(3.74)

0

Из определения функции u вытекает тождество
=
a (u , v )

1

 ,
∫ fvdx

для ∀v ∈ V .

0

Таким образом, a (u , v ) = a (u, v ) .
Используя тождества (3.74) и (3.75), нетрудно получить
62

(3.75)

a (u − u I , u − u I )= a (u , u − u I ) − a (u I , u − u I )=
1

= ∫ f (u − u I )dx −a (u I , u − u I ) = a (u − u I , u − u I ).

.

(3.76)

0

Поскольку функция (u − uI )′ постоянна на каждом элементе, то,
применяя формулу интегрирования по частям, получаем
xi +1

xi +1

xi

xi

− ∫ p′(u − uI )(u − uI )′dx +
∫ p(u − uI )′(u − uI )′dx =
x
p (u − u I )(u − u I )′ xi +1
i

+

(3.77)

xi +1

=
− ∫ p′(u − u I )(u − u I )′dx.
xi

Следовательно, применяя лемму 1, получаем
1

a (u − u I , u − u I ) =
∫ [− p′(u − uI )(u − uI )′ + q(u − uI )(u − uI )]dx ≤
0

(3.78)

1

≤ u − u I ∫ [− p′(u − u I )′ + q (u − u I )]dx.
0

Далее
xi +1

n −1 
xi +1
′
′
′



−
−
=
−
−
+
p
u
u
dx
p
u
u
(
)
(
)
∑
∫
∫ p′′(u − uI )dx  =
I
I xi
i =1 
xi
0


1

(3.79)

1

= ∫ p′′(u − u I )dx.
0

Тогда
a (u − u I , u − u I ) ≤ u − u I
≤ u − uI

С

≤ с u − uI

u − u I
С

1
С

∫ [ p′′(u − uI ) + q(u − uI )]dx ≤

0

1
С

u − u I

∫ ( p′′ + q)dx ≤

(3.80)

0
1

, где с = | ∫ ( p′′ + q )dx | .
С
0

Применяя лемму 1, получаем
a (u − u I , u − u I ) ≤ с u − u I

Далее
63

1
С

2
∫ (u′ − u′I ) dx .

0

(3.81)

1

a (u − u I , u − u I =
) ∫ [ p (u ′ − u ′I ) 2 + q (u − u I ) 2 ]dx ≥
0
1

≥ c0 ∫ (u ′ − u ′I ) dx ≥c0 u − u I
0

(3.82)

1

2

C

∫ (u′ − uI′ ) dx .
2

0

Таким образом,
c0 u − u I

1
C

2
∫ (u′ − u′I ) dx ≤ a(u − uI , u − uI ) =a(u − uI , u − uI ) ≤

0

(3.83)

1

≤ с u − uI

С

∫ (u′ − u′I ) dx .
2

0

Отсюда
u − u I

Следовательно,
u − u

C

C

≤ c u − uI

= u − u I + u I − u

C

C

≤ u − uI

.

C

(3.84)
+ u − uI

C

≤

(3.85)
h2
′′
.
u
C
C
8
Таким образом, отклонение точного решения от приближенного
представляет собой величину порядка h 2 .
≤ (1 + с)

≤ (1 + c) u − u I

4. МЕТОД КОНЕЧНЫХ ЭЛЕМЕНТОВ ДЛЯ РЕШЕНИЯ
УРАВНЕНИЙ С ЧАСТНЫМИ ПРОИЗВОДНЫМИ
4.1. Классическая и интегральная формулировка задачи
Рассмотрим обобщение метода на двумерные задачи.
Пусть Ω некоторая область в R 2 , ∂Ω – её граница, n –
единичный вектор внешней нормали (рис. 25).
Определим скалярное произведение:

( ϕ1, ϕ2 )= ∫ ϕ1 ( x) ϕ2 ( x)dx.
Ω

(4.1)

При изучении свойств дифференциального оператора для одномерных задач применялось правило интегрирования по частям. В двумерном случае для этого используются формулы Грина. Сформулируем некоторые результаты, необходимые для дальнейшего.
64

Рис. 25

Будем считать, что область Ω удовлетворяет условиям формулы



Остроградского‒Гаусса. Пусть =
F ( x, y ) iF1 ( x, y ) + jF2 ( x, y ), ( x, y ) ∈ Ω,
где F1 ( x, y ), F2 ( x, y ) – некоторые дифференцируемые функции. Кроме
того, будем считать, что в каждой точке ( x, y ) ∈ ∂Ω, границы области
 
Ω задан единичный вектор внешней нормали n = n( x, y ) . Обозначим
 ∂F ∂F
1+
2 .
=
divF
∂x
∂y
Формула Остроградского‒Гаусса, которую иногда называют
теоремой о дивергенции:

 
dy ∫ F ⋅ n ds .
∫∫ divF dx=
(4.2)
Ω

∂Ω

Предположим, что граница области ∂Ω состоит из двух
непустых непересекающихся частей, ∂Ω1 и ∂Ω 2 , т. е.
∂Ω = ∂Ω1 ∪ ∂Ω 2 ,

где

∂Ω1 ∩ ∂Ω 2 = ∅.

Будем считать, что в области Ω и на её границе заданы функции
α, p, q удовлетворяющие условиям:
α : ∂Ω 2 → R ‒ непрерывная функция α > 0 ;
p : Ω → R, некоторая непрерывно дифференцируемая функция;
p ∈ C1[Ω],
p ≥ c0 > 0 ( x, y ) ∈ Ω ;
q ∈ C 0 (Ω),
q ≥ 0 ‒ непрерывная функция, заданная в области Ω .
Определим класс U дважды непрерывно дифференцируемых в
области Ω функций и удовлетворяющих условиям
65

∂u
p=
+ αu 0 на ∂Ω 2 ,
∂n
∂u


=
U u ∈ C 2=
+ αu 0 на ∂Ω 2  ,
(Ω) : u 0 на ∂Ω1 , p =
∂n


где ∂u ‒ нормальная производная функции u на границе ∂Ω :
∂n

∂u
=∇u ⋅ n ,
∂n
∂u
∂u – градиент функции u.
∇u= i + j
∂x
∂y
Лемма. В классе функций u ∈ U оператор
=
u 0 на ∂Ω1 ,

Lu =
− div( p∇u ) =
−

∂  ∂u  ∂  ∂u  .
 p −  p 
∂x  ∂x  ∂y  ∂y 

положителен и обладает свойством симметрии, т.е.
=
( Lu, v ) ( u, Lv ) , ∀u, v ∈U .

( u, Lu ) > 0,

( u, Lu ) = 0

(4.3)

(4.4)

(4.5)

⇔ u = 0.

Доказательство. Действительно, пусть

F =−vp∇u ,
∀u , v ∈ U .

(4.6)

Тогда

 ∂  ∂u  ∂  ∂u  
∂  ∂u  ∂  ∂u 
−  vp  −  vp  =
−v   p  +  p   −
divF =
∂x  ∂x  ∂y  ∂y 
 ∂x  ∂x  ∂y  ∂y  
∂v ∂u
∂v ∂u
−p
−p
.
∂x ∂x
∂y ∂y

Или, в векторных обозначениях,

divF = −vdiv( p∇u ) − p∇u ⋅∇v .

(4.7)

(4.8)

Подставим это выражение в (4.2):
∂u

∫∫ ( −vdiv( p∇u ) − p∇u ⋅∇v ) dx dy = − ∫ vp ∂n ds,
Ω
∂Ω2
Поменяв местами роли функций u и ν, получим
∂v
∫∫ ( −udiv( p∇v) − p∇v ⋅∇u ) dx dy = − ∫ up ∂n ds,
Ω
∂Ω 2
Вычитая (4.10) из (4.9), получаем
66

∀u , v ∈ U .

∀u , v ∈ U .

(4.9)

(4.10)


∂Ω 2 

∂u

∂v 

− ∫  vp − up  ds.
∫∫ ( −vdiv( p∇u ) + udiv( p∇v) ) dx dy =
∂n
∂n
Ω



(4.11)

Формулы (4.9) и (4.11) называют соответственно первой и второй
формулой Грина.
В силу определения класса U
∂u
∂v
p
= −αu ,
p
= −αv,
( x, y ) ∈ Ω 2 .
(4.12)
∂n
∂n
Таким образом,
∂v 
 ∂u
− ∫  vp − up =
ds
ds 0.
∫ ( vαu − uαv )=
(4.13)
∂n
∂n 
∂Ω 2 
∂Ω 2
Следовательно,
dy ∫∫ ( vLu − uLv ) dx=
dy 0,
∫∫ ( −vdiv( p∇u ) + udiv( p∇v) ) dx=
Ω

Ω

(4.14)

т.е. оператор L симметричен. Подставив в (4.9) вместо функции ν
функцию u, получим
∂u
2
∀u ∈ U .
∫∫ −udiv( p∇u ) − p(∇u ) dx dy =− ∫ up ∂n ds,
(4.15)
Ω
∂Ω2

(

)

Так как p ∂u = −αu, ( x, y ) ∈ ∂Ω 2 , то отсюда
∂n
2
2
∫∫ uLu dx dy= ∫∫ p(∇u ) dx dy + ∫ αu ds ≥ 0,
Ω

Ω

Далее

∀u ∈ U .

∂Ω 2

=
∇u 0, ( x, y ) ∈ Ω,

0⇔
∫∫ uLu dx dy= =
u 0, ( x, y ) ∈ ∂Ω 2 .
Ω

(4.16)

(4.17)

Из первого уравнения в (4.17) следует, что u = c, где c ‒
постоянная, а из второго с = 0, что равносильно u = 0.
Рассмотрим следующие две задачи, аналогичные рассмотренным
в одномерном случае.
∂u

Задача 1. Пусть U= u ∈ C 2 (Ω) : =
u 0 на ∂Ω1 , p + α=
u 0
∂n


на ∂Ω2  , где n ‒ внешняя нормаль. Требуется найти такую функцию


u ∈ U , что

Lu = −

∂  ∂u  ∂  ∂u 
 p  −  p  + qu = f ,
∂x  ∂x  ∂y  ∂y 
67

( x, y ) ∈ Ω .

(4.18)

Задача 2. Пусть U=

{u ∈ C1(Ω) : u=

}

0 на ∂Ω1 . Требуется най-

ти такую функцию u ∈ U , что
 ∂u ∂v

∂u ∂v 

∫∫ p  ∂x ∂x + ∂y ∂y  dxdy +=
∫ α u vds ∫∫ fv dxdy
Ω





∂Ω2

∀v ∈ U .

Ω

(4.19)

Теорема. Задачи 1 и 2 эквивалентны.
Доказательство. Покажем сначала, что задача 1 имеет не более
одного решения. Предположим, что u1 , u 2 ∈ U ‒ два различных
решения задачи.
Введем в рассмотрение вспомогательную функцию =
v u1 − u 2 .
Тогда Lv = 0 в области Ω и v ∈ U . В этом случае ( Lv, v ) = 0 ,
следовательно, v = 0 , что равносильно u1 = u 2 .
Предположим, что u является решением задачи 1, т.е.
−

∂  ∂u  ∂  ∂u 
 + qu − f 0,
 p −  p =
∂x  ∂x  ∂y  ∂y 

( x, y ) ∈ Ω. .

(4.20)

Умножим (4.20) на произвольную функцию v ∈ U и проинтегрируем
по области Ω :
0.
∫ ( −vdiv( p∇u ) − vf ) dxdy =

(4.21)

Ω

Преобразовав первое слагаемое по формуле (4.9), получим
∂u

= ∫∫ p∇u ⋅∇v dx dy − ∫ vp ds, ∀u , v ∈ U .
∫∫ −vdiv( p∇u ) dx dy
∂n
Ω
Ω
∂Ω 2

(4.22)

Поскольку p ∂u = −αu, ( x, y ) ∈ ∂Ω 2 , то
∂n
 ∂u ∂v

∂u ∂v 

∫∫ p  ∂x ∂x + ∂y ∂y  dxdy +=
∫ α u vds ∫∫ fv dxdy
Ω





Ω2

∀v ∈ U

Ω

(4.23)

Будем считать что u явлется решением задачи 2, т.е. для ∀v ∈ U
функция u = 0 на границе ∂Ω1 и удовлетворяет уравнению (4.19).
Преобразуем первое слагаемое по формуле (4.9):
∂u
∫∫ p∇u ⋅∇v dx dy= ∫∫ −vdiv( p∇u ) dx dy + ∫ vp ∂n ds, ∀u, v ∈ U . (4.24)
Ω
Ω
∂Ω 2
После подстановки этого выражения в (4.23) получим
68

∂u



= 0, ∀u , v ∈ U .
∫∫ v ( −div( p∇u ) − f ) dx dy + ∫ v  p ∂n + αu  ds

Ω
∂Ω 2 

(4.25)

Обозначим
γ = −div( p∇u ) − f , ( x, y ) ∈ Ω,
(4.26)
∂u
( x, y ) ∈ Ω 2 .
=
β p + αu ,
∂n
Чтобы показать, что уравнения (4.18) выполняются, нужно доказать,
что
=
γ 0, ( x, y ) ∈ Ω,
(4.27)
=
β 0, ( x, y ) ∈ ∂Ω 2 .

Выберем сначала
v = γϕ ,
( x,=
y ) ∈ Ω , ϕ 0,
где ϕ > 0,
нение (4.25) примет вид

(4.28)

( x, y ) ∈ ∂Ω , ϕ∈ C (Ω) . Тогда урав0

2
∫ γ ϕdxdy =0 ,

(4.29)

Ω

что равносильно
=
γ 0, ( x, y ) ∈ Ω . Далее выберем
v = βψ,

(4.30)

где ψ > 0, ( x, y ) ∈ ∂Ω 2 , ψ =0, ∈ ∂Ω1 , ψ ∈ C (∂Ω) .
Тогда (4.25) с учетом того,
что γ 0, ( x, y ) ∈ Ω , примет вид
=
0

2
∫ β ψdxdy =0 ,

∂Ω 2

(4.31)

откуда
=
β 0, ( x, y ) ∈ ∂Ω 2 . Таким образом теорема доказана.
Следует отметить, что условие
=
u 0 на ∂Ω1 , является
u
главным, а условие p ∂=
+ αu 0 на ∂Ω 2 – естественным.
∂n
Введем обозначения
 ∂u ∂v ∂u ∂v 
+
a=
(u , v) ∫∫ p 
 dxdy + ∫ α u vds;
Ω  ∂x ∂x ∂y ∂y 
∂Ω 2
ϕ(v) =
∫∫ fv dxdy .
Ω

В терминах этих обозначений уравнение примет вид
69

a (u , v) = ϕ(v) .

Билинейная форма a (u , v) , так же как и в предыдущем примере,
является симметричной и положительно определённой.
Согласно теореме (1.9) краевую задачу (4.18) можно
сформулировать как задачу определения экстремума функционала
1
=
F (u )
a (u , u ) − ϕ(u )
(4.32)
2
в классе функций U=
F (u )=

{u ∈ C1(Ω) : u=

}

0 на ∂Ω1 , или

1
2
2
∫ p(∇u ) dxdy + ∫ α u vds − ∫ fudxdy .
2Ω
∂Ω 2
Ω

Приравняв нулю первую вариацию функционала, получим
 ∂u ∂ (δu ) ∂u ∂ (δu )

=
δF (u ) ∫∫ 
+
− f δu  dxdy + ∫ =
α uδuds
∂y ∂y

Ω  ∂x ∂x
∂Ω 2

(4.33)

0.

Применим к первым двум слагаемым формулу (4.9) v = δu :
∂u
∫∫ p∇u∇(δu )dxdy= ∫∫ −δudiv( p∇u )dxdy + ∫ p ∂n (δu )ds .
Ω
Ω
∂Ω 2
Таким образом,

 ∂u

δF (u=
) ∫∫ [ − div( p∇u ) − f ] δudxdy + ∫  p + α u  δuds= 0 .

Ω
∂Ω2  ∂n

В силу произвольности δu
−div
=
( p∇u ) f , ( x, y ) ∈ Ω ,

(4.34)

∂u
p=
+ αu 0, ( x, y ) ∈ ∂Ω 2 .
∂n
Отметим частный случай: если p = 1 то уравнение (4.34) принимает
вид
=
u (0) 0, ( x, y ) ∈ ∂Ω1 ,

 ∂ 2u ∂ 2u 
−div( p∇u ) = −  2 + 2  = −∆u ,
 ∂x
∂y 


(4.35)

2
2
где ∆u= ∂ u + ∂ u . Уравнение
∂x 2 ∂y 2

=
−∆u f , ( x, y ) ∈ Ω

(4.36)
называется уравнением Пуассона. Иногда его рассматривают
совместно с условием одного типа, заданным на всей границе:
70

=
u 0,

( x, y ) ∈ ∂Ω

(4.37)
– задача Дирихле или первая краевая задача для уравнения Пуассона.
Решение (4.36) с условием
∂u
= 0,
( x, y ) ∈ ∂Ω
(4.38)
∂n
называется задачей Неймана или второй краевой задачей для
уравнения Пуассона.
Третьей краевой задачей для уравнения Пуассона называется
краевая задача с условием
∂u
=
+ αu 0 ( x, y ) ∈ ∂Ω .
(4.39)
∂n
Рассмотренная краевая задача (4.18), в которой на разных частях
границы задаются различные краевые условия, называется
смешанной краевой задачей.
4.2. Дискретизация области, выбор базисных функций
Рассмотрим применение метода конечных элементов для
решения простейших задач. При разбиении области используем
треугольные элементы с базисными функциями линейного типа.
Путем обхода всех узлов сетки сформируем систему линейных
алгебраических уравнений и учтем граничные условия.
Метод сборки системы путем обхода узлов практически не
применяется при машинной реализации и служит здесь для учебных
целей. Формирование матрицы системы из матриц элементов
рассматривается в следующем разделе.
Рассмотрим первую краевую задачу (задачу Дирихле) для уравнения Лапласа в случае, когда область =
Ω [0,1] × [0,1]. представляет
собой квадрат:
 ∂ 2u ∂ 2u
+ 2 0, ( x, y ) ∈ Ω,
 =
2
∂y
 ∂x
=
u (0, y ) 0,
y ∈ [0;1],

u ( x,1) 0,
x ∈ [0;1],
=
=
u ( x, 0) 0,
x ∈ [0;1],

u
(1,
y
)
=
y
(1
−
y
),
y ∈ [0;1].



71

(4.40)

Умножим (4.40) на некоторую функцию v , равную нулю на
границе области ∂Ω , проинтегрируем по области Ω с использованием формулы Грина.
В результате задача (4.40) сводится к нахождению функции u,
принимающей на границе области ∂Ω заданные значения и
удовлетворяющей уравнению
 ∂u ∂v

∂u ∂v 

0
∫∫  ∂x ∂x + ∂y ∂y  dxdy =

(4.41)

для любой функции, v принимающей нулевые значения на границе
области ∂Ω .
Введем для этого случая на области Ω квадратную сетку с шагом
h, h = 1/ n :
Ω

xi , yi ) (k h, l h),=
k , l 0,1,  , n} .
=
ω {(=

(4.42)

Узел сетки с номером i определяется как пересечение строки с
номером k и столбца с номером l (рис. 26). Соответственно значение
функции u ( h ) и базисной функции Ni ( x, y ) в узле с номером i
обозначим
u ( h ) ( xi , yi ) ≡ uk(,hl) ,
(4.43)
Ni ( x, y ) ≡ N k,l ( x, y ).

Рис. 9

Каждую ячейку сетки разделим на два треугольника диагональю,
параллельной биссектриссе первого координатного угла. Получим
72

разбиение области Ω на треугольные элементы Ω(e) , e =
1, 2,  , 2n 2 ,
т.е. триангуляцию области Ω (рис. 27).

Рис. 10

Функция Ni ( x, y ) («шапочка Куранта» ‒ рис. 27, а), как нетрудно
проверить, тождественно равна нулю вне области Ω k ,l , представляющей собой объединение треугольных элементов, имеющих
общую вершину.
Таким образом, решение задачи сводится к нахождению
неизвестных значений функции в узлах конечноэлементной сетки.
Функцию u ( h ) ( x, y ) можно представить в виде
n −1

u ( h ) ( x, y ) = ∑ uk( h,l) N k ,l ( x, y ) ,
k ,l =1

(4.44)

где принято следующее соответствие:
(k ,l ) ⇔ i

и (k ′,l ′) ⇔ j .

(4.45)

Примем v = N k,l , подставим (4.44) в (4.41). После несложных
преобразований получим
n −1

 ∂N k ,l ∂N k ′,l ′ ∂N k ,l ∂N k ′,l ′ 
+
∫∫ N k ,l f dxdx . (4.46)
dxdy =
∂x
∂y
∂y 
Ω  ∂x
Ω

∑ ukh′,l ′ ∫∫ 

k , l=1

Понятно, что интеграл может быть отличен от нуля только в
случае, когда области определения функций N k ,l и N k ′,l ′
пересекаются. Для фиксированных k,l таких областей Ω k′,l ′ может
быть только шесть (рис. 27, б), а именно: Ω k −1,l , Ω k −1,l-1 , Ω k ,l-1 ,
73

Ω k +1,l , Ω k +1,l+1 , Ω k ,l+1 . При этом в левой части (4.46) остается семь
слагаемых:
 ∂N 2  ∂N 2 
ukh,l ∫∫  k ,l  +  k ,l  dxdy +
∂x   ∂y  
Ω k ,l 


 ∂N k ,l ∂N k −1,l ∂N k ,l ∂N k −1,l
+ukh−1,l
+
∫∫  ∂x
∂x
∂y
∂y
Ω k ,l ∩Ω k −1,l 
+ukh−1,l-1

 ∂N k ,l ∂N k −1,l-1 ∂N k ,l ∂N k −1,l-1
+

∂x
∂y
∂y
Ω k ,l ∩Ω k −1,l-1  ∂x

∫∫


dxdy +


 ∂N k ,l ∂N k +1,l ∂N k ,l ∂N k +1,l 
+

dxdy +
∂x
∂y
∂y 
Ω k ,l ∩Ω k +1,l  ∂x

+ukh+1,l

∫∫

+ukh+1,l+1
+ukh,l-1


dxdy +


 ∂N k ,l ∂N k +1,l+1 ∂N k ,l ∂N k +1,l+1 
+

dxdy +
∂x
∂y
∂y

Ω k ,l ∩Ω k +1,l+1  ∂x

∫∫

 ∂N k ,l ∂N k ,l-1 ∂N k ,l ∂N k ,l-1
+

∂x
∂y
∂y
Ω k ,l ∩Ω k ,l-1  ∂x

+ukh,l+1

∫∫


dxdy +


 ∂N k ,l ∂N k ,l+1 ∂N k ,l ∂N k ,l+1
+

∂x
∂y
∂y
Ω k ,l ∩Ω k ,l+1  ∂x

∫∫


∫∫ N k ,l f dxdx. (4.47)
dxdy =

Ω k ,l

Первое слагаемое
 ∂N 2  ∂N 2 
k ,l
k ,l
4.
∫∫  ∂x  +  ∂y  dxdy =

 
 
Ω kl 

Результаты вычислений остальных слагаемых в (4.47) сведены в
табл. 5 – 10.
Таблица 5
n

∂N k −1,l-1
∂x

∂N k −1,l-1
∂y

∂N k ,l
∂x

∂N k ,l
∂y

1
2
3
4
5
6

-1/h
0
0
0
0
0

0
0
0
0
0
-1/h

0
-1/h
-1/h
0
1/h
1/h

1/h
1/h
0
-1/h
-1/h
0

Таким образом, получим
74

 ∂N k,l ∂N k −1,l−1 ∂N k,l ∂N k −1,l−1 
0.
+

dxdy =
∂x
∂y
∂y

Ω k ,l ∩Ω k −1,l−1  ∂x

∫∫

Аналогично

Таблица 6

n

∂N k ,l-1
∂x

∂N k ,l-1
∂y

∂N k ,l
∂x

∂N k ,l
∂y

1
2
3
4
5
6

0
0
0
0
-1/h
-1/h

0
0
0
0
0
1/h

0
-1/h
-1/h
0
1/h
1/h

1/h
1/h
0
-1/h
-1/h
0

 ∂N k ,l ∂N k ,l-1 ∂N k ,l ∂N k ,l-1
+

∂x
∂y
∂y
Ω k , l ∩Ω k , l −1  ∂x

∫∫


−1.
dxdy =


Таблица 7
n

∂N k +1,l
∂x

∂N k +1,l
∂y

∂N k ,l
∂x

∂N k ,l
∂y

1
2
3
4
5
6

0
0
0
-1/h
0
0

0
0
0
1/h
1/h
0

0
-1/h
-1/h
0
1/h
1/h

1/h
1/h
0
-1/h
-1/h
0

 ∂N k ,l ∂N k +1,l ∂N k ,l ∂N k +1,l
+

∂x
∂y
∂y
Ω k , l ∩Ω k +1, l  ∂x

∫∫


−1.
dxdy =


Таблица 8
n

∂N k +1,l+1
∂x

∂N k +1,l+1
∂y

∂N k ,l
∂x

∂N k ,l
∂y

1
2
3
4
5
6

0
0
0
1/h
0
0

0
0
1/h
0
0
0

0
-1/h
-1/h
0
1/h
1/h

1/h
1/h
0
-1/h
-1/h
0

75

 ∂N k ,l ∂N k +1,l+1 ∂N k ,l ∂N k +1,l+1
+

∂x
∂y
∂y
Ω k , l ∩Ω k +1, l +1  ∂x

∫∫


0.
dxdy =


Таблица 9
n

∂N k ,l+1
∂x

∂N k ,l+1
∂y

∂N k ,l
∂x

∂N k ,l
∂y

1

0

0

0

1/h

2

1/h

0

-1/h

1/h

3

1/h

-1/h

-1/h

0

4

0

0

0

-1/h

5

0

0

1/h

-1/h

6

0

0

1/h

0

 ∂N k ,l ∂N k ,l+1 ∂N k ,l ∂N k ,l+1
+

∂x
∂y
∂y
Ω k , l ∩Ω k , l +1  ∂x

∫∫


−1.
dxdy =


Т а б л и ц а 10
n

∂N k −1,l
∂x

∂N k −1,l
∂y

∂N k ,l
∂x

∂N k ,l
∂y

1

1/h

-1/h

0

1/h

2

0

-1/h

-1/h

1/h

3

0

0

-1/h

0

4

0

0

0

-1/h

5

0

0

1/h

-1/h

6

0

0

1/h

0

 ∂N k ,l ∂N k −1,l ∂N k ,l ∂N k −1,l
+

∂x
∂y
∂y
Ω k , l ∩Ω k −1, l  ∂x

∫∫


−1 .
dxdy =


Таким образом, уравнение (4.47), записанное для узла (k.l),
примет вид
−ukh−1,l − ukh+1,l + 4ukh,l − ukh,l-1 − ukh,l+1 =
0.
76

(4.48)

Точно такое же уравнение получается при замене в (4.40) производных
разделёнными разностями:
∂ 2u
∂x

2

∂ 2u
∂y 2

=

uk −1,l − 2uk ,l + uk +1,l

=

uk ,l-1

h2
− 2uk ,l + uk ,l+1
h2

,

(4.49)
.

Пронумеруем узлы сетки, как показано на рис. 28. Далее, в соответствии со сказанным выше, составим
Рис. 28
систему линейных алгебраических
уравнений.
Для узлов, находящихся в углах 1, 5, 21, 25, строка матрицы K
содержит три ненулевых элемента, для узлов на границе 2 – четыре, а
для внутренних узлов, например 7, – пять.
В узле 1 функция формы определена в треугольнике с узлами 1,
6, 2 и имеет производные по х и у:
 ∂N 2  ∂N 2 
∂N1
1 ∂N1 1
2 h2
=
− ;
=
⇒ ∫∫  1  +  1   dxdy =
=
1.
2 2
∂x
∂у h Ω(1)  ∂x   ∂y  
h
h



Найдем производные по у от функции формы, относящейся к узлу 2:
 ∂N ∂N ∂N ∂N 
∂N 2 1 ∂N 2
1 h2
1
=
=
− 2
=
− .
;
0 ⇒ ∫∫  2 1 + 2 1  dxdy =
h
∂x
∂у
∂x ∂x
∂y ∂y 
2
h 2
Ω(1) 

Аналогично вычисляем коэффициент для узла 6:
 ∂N ∂N ∂N ∂N 
∂N 6
∂N 6
1
1 h2
1
=
0;
=
− ⇒ ∫∫  6 1 + 6 1  dxdy =
− 2
=
− .
∂x
∂у
h Ω(1)  ∂x ∂x
∂y ∂y 
2
h 2

Таким образом, уравнение системы, относящееся к узлу 1, имеет вид
u1 − 0.5u2 − 0.5u6 =
0.

Результирующая система линейных алгебраических уравнений,
полученная в результате обхода всех узлов, имеет следующий вид
(для сокращения обе части уравнений системы умножены на два):
77

2 -1 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0   u1   0 
 
1 4 -1 0 0 0 - 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0   u2   0 
0 -1 4 -1 0 0 0 - 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0   u3   0 
   
0 0 -1 4 -1 0 0 0 - 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0   u4   0 
0 0 0 -1 2 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0   u5   0 
   
1 0 0 0 0 4 - 2 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0   u6   0 
0 - 2 0 0 0 - 2 8 - 2 0 0 0 - 2 0 0 0 0 0 0 0 0 0 0 0 0 0   u7   0 
   
0 0 - 2 0 0 0 - 2 8 - 2 0 0 0 - 2 0 0 0 0 0 0 0 0 0 0 0 0   u8   0 
0 0 0 - 2 0 0 0 - 2 8 - 2 0 0 0 - 2 0 0 0 0 0 0 0 0 0 0 0   u9   0 
 
0 0 0 0 -1 0 0 0 - 2 4 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0   u10   0 
   
0 0 0 0 0 -1 0 0 0 0 4 - 2 0 0 0 -1 0 0 0 0 0 0 0 0 0   u11   0 
0 0 0 0 0 0 - 2 0 0 0 - 2 8 - 2 0 0 0 - 2 0 0 0 0 0 0 0 0   u12   0 
   
0 0 0 0 0 0 0 - 2 0 0 0 - 2 8 - 2 0 0 0 - 2 0 0 0 0 0 0 0   u13  =  0 
0 0 0 0 0 0 0 0 - 2 0 0 0 - 2 8 - 2 0 0 0 - 2 0 0 0 0 0 0   u14   0 
   
0 0 0 0 0 0 0 0 0 -1 0 0 0 - 2 4 0 0 0 0 -1 0 0 0 0 0   u15   0 
0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 4 - 2 0 0 0 -1 0 0 0 0   u16   0 
   
0 0 0 0 0 0 0 0 0 0 0 - 2 0 0 0 - 2 8 - 2 0 0 0 - 2 0 0 0   u17   0 
   
0 0 0 0 0 0 0 0 0 0 0 0 - 2 0 0 0 - 2 8 - 2 0 0 0 - 2 0 0   u18   0 
0 0 0 0 0 0 0 0 0 0 0 0 0 - 2 0 0 0 - 2 8 - 2 0 0 0 - 2 0   u19   0  (4.50)
   
0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 - 2 4 0 0 0 0 -1   u20   0 
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 2 -1 0 0 0   u21   0 
   
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 - 2 0 0 0 -1 4 -1 0 0   u22   0 
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 - 2 0 0 0 -1 4 -1 0   u23   0 
   
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 - 2 0 0 0 -1 4 -1   u24   0 
 
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 -1 2   u25   0 

зом:

Выполнив граничные условия, поступаем следующим обра-

1) обнуляем строку, соответствующую граничному узлу,
2) на главной диагонали ставим единицу,
3) в вектор-столбец правой части заносим значение функции в
данном узле.
В результате получим
78

1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0   u1   0 
 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0  u  


 2   0 
 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0   u3   0 

  

 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0   u4   0 
 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0   u5   0 

  

 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0   u6   0 
 0 -1 0 0 0 -1 4 -1 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0   u   0 

 7  

 0 0 -1 0 0 0 -1 4 -1 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0   u8   0 
 0 0 0 -1 0 0 0 -1 4 -1 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0   u   0 

 9  

 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0   u10   0.1875 
 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0  u  


  11   0 
 0 0 0 0 0 0 -1 0 0 0 -1 4 -1 0 0 0 -1 0 0 0 0 0 0 0 0   u12   0 

  

 0 0 0 0 0 0 0 -1 0 0 0 -1 4 -1 0 0 0 -1 0 0 0 0 0 0 0   u13  =  0 
 0 0 0 0 0 0 0 0 -1 0 0 0 -1 4 -1 0 0 0 -1 0 0 0 0 0 0   u14   0 

  

 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0   u15   0.2500 
 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 1-1 0 0 0 -1 0 0 0 0   u   0 

  16  

0
0
0
0
0
0
0
0
0
0
0
-1
0
0
0
-1
4
-1
0
0
0
-1
0
0
0

  u17   0 
 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 -1 4 -1 0 0 0 -1 0 0   u   0 

  18  

 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 -1 4 -1 0 0 0 -1 0   u19   0 
 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0  u  


  20   0.1875 
 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0   u21   0 

  

 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0   u22   0 
 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0   u23   0 

  

 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0   u24   0 
 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1  u   0 
 (4.51)

 25  

Полученное решение изображено графически на рис. 29 в виде
изолиний.
На рис. 30 представлена картина, соответсвующая аналитическому решению этой задачи:
∞

sh(π(2m − 1) x) sin(π(2m − 1) y )

m =1

(π(2m − 1))3 sh(π(2m − 1))

u ( x, y ) = 8 ∑

79

.

(4.52)

Метод конечных элементов (сетка 5Х5)
1

0.05

0.9

0.1

0.05

0.8

5
0.1
0.7

0.05

0.1

0.6
0.5

5
0.1

0.4

0.2

0.3

1
0.

0.0
5

0.2
0.1

0.05
0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Рис. 29

Аналитическое решение
1

0.05

0.9

0.1

5
0.0

0.8

15
0.

0.
2

0.7

0.05

0.1

0.6
0.5

0.15

0.4

0.2

0.3

0.
1

0.0
5

0.2
0.1

0.05
0.1

0.2

0.3

0.4

0.5

Рис. 30
80

0.6

0.7

0.8

0.9

1

4.3. Конечноэлементная формулировка задачи.
Матрицы элементов
Рассмотрим общий случай решения уравнения
−

∂  ∂u  ∂  ∂u 
p  + qu f ,
 p  − =
∂x  ∂x  ∂y  ∂y 

( x, y ) ∈ Ω

(4.53)

с граничными условиями
∂u
p=
+ αu 0 ( x, y )∂Ω 2 .
(4.54)
∂n
Пусть v – некоторая функция, удовлетворяющая главным краевым условиям
=
v 0 ( x, y ) ∈ ∂Ω1 .
(4.55)
=
u 0 ( x, y ) ∈ ∂Ω1 ,

Применяя метод Галёркина, умножим обе части на функцию v и,
следуя процедуре, описанной в подразд. 4.1, приходим к задаче
определения функции u из интегрального тождества:
 ∂u ∂v

∂u ∂v 

∫∫ p  ∂x ∂x + ∂y ∂y  dxdy + ∫ α u vds =
∫∫ fv dxdy
Ω





∂Ω2

Ω

(4.56)

для любой функции v , удовлетворяющей условию (4.55).
Следуя методу конечных элементов, представим область Ω в
виде объединения N el попарно непересекающихся конечных
элементов Ω(e) с границей ∂Ω(e) (рис. 31):
Ω = ∪Ω(e)

Ω(i ) ∩ Ω( j ) = ∅, ∀i ≠ j.

Рис. 31
81

(4.57)

Будем считать, что пространство конечных элементов U ( h )
состоит из функций, непрерывных в области Ω , линейных на каждом
элементе и равных нулю на границе ∂Ω1 , т.е.
u | u − C 0 [Ω], u − линейная функция на
U (h) = 
(e)
и u 0 на ∂Ω1
 треугольнике Ω =

каждом 

.


Ясно, что U ( h ) ⊂ U .
Базис пространства U ( h ) составляет набор линейно-независимых
базисных функций { Ni ( x, y )}n , где Ni ( x, y ) – функции, непрерывные
i =1

в области Ω , равные нулю на границе ∂Ω1 , линейные на каждом
элементе Ω(e) , удовлетворяющие условиям
0, ( x, y ) ≠ ( xi , yi ),
N i ( x, y ) = 
(4.58)
1, ( x, y ) = ( xi , yi ).
Приближенное решение u(h) задачи (4.56) будем искать в виде
функций, непрерывных в области Ω , линейных на каждом элементе и
равных нулю на границе ∂Ω1 , т.е. u (h ) ∈ U (h ) .

v(h)

Представим приближенное решение u ( h ) и тестовую функцию
в виде разложения по базисным функциям Ni ( x, y ) :
n

n

=
u (h) ∑
=
ui( h ) Ni ( x, y ),
v ( h ) ∑ ci Ni ( x, y ),

(4.59)

=i 1 =i 1
где ui( h ) в силу (4.58) ‒ значения неизвестной функции ui( h ) в узлах
сетки ui( h ) = ui( h ) ( xi , yi ) , ci – некоторые произвольные коэффициенты.

Подставив соотношения (4.59) в (4.56) получим

(h)
(h)
(h) (h)
(h)
0.
∫∫ p∇u ∇v dxdy + ∫ α u v ds − ∫∫ fv dxdy =
Ω

∂Ω2

Ω

(4.60)

Выражения (4.59) для функций u ( h ) и v ( h ) c использованием
матричных обозначений примут вид
=
u (h) [ =
N ]{U },
v ( h ) [ N ]{V },
(4.61)
∇u ( h ) =
[∇N ]{U },
∇v ( h ) =
[∇N ]{V },
где
т
=
{U } [u=
{V } [v1v2  vn ]т ,
1u2  un ] ,

[ N ] =[ N1 N 2  N n ],

[∇N ] =∇
[ N1∇N 2  ∇N n ].
82

(4.62)

Тогда (4.60) можно записать в виде

т
т
{V }т  ∫∫ p [∇N ] [∇N ] + q [ N ] [ N ] dxdy +
∂Ω
 Ω

т
т
0
+ ∫ α [ N ] [ N ] ds − ∫∫ f [ N ] dxdy {U } =
∂Ω 2
Ω


)

(

(4.63)

Отсюда, в силу произвольности {V }т ,
K {U } = {B} ,

где

(

т

)

т

(4.64)
т

K= ∫∫ p [∇N ] [∇N ] + q [ N ] [ N ] dxdy + ∫ α [ N ] [ N ] ds
Ω

∂Ω 2

т

B = ∫∫ f [ N ] dxdy.
Ω

Представив интегралы в (4.64) в виде суммы интегралов по
элементам, получим
n

n

K= ∑ K (pe) + K q(e) + K α(e) ,

B= ∑ B (e) ,

=e 1 =e 1

(4.65)

где
K (pe) =

т
∫∫ p [∇N ] [∇N ]dxdy,

Ω

т

K α(e) =
∫ α [ N ] [ N ] ds,
∂Ω(2e )

где

т
∫∫ q [ N ] [ N ]dxdy,

K q(e) =

Ω( e )

(e)

т

B (e) =
∫∫ f [ N ] dxdy,
Ω( e )

– сторона треугольника, принадлежащая границе ∂Ω 2 .
В пределах треугольного элемента в строке базисных функций
отличны от нуля всего три элемента с номерами i, j , k , которые
соответствуют глобальным номерам узлов элемента, N1 , N 2 , N3 ‒
функции формы треугольного элемента (2.14):
(e)

∂Ω 2

[ N ] = [0  N1  N 2  N3  0] .



i

j

k

(4.66)

Поэтому в матрицах K (pe) = ( pij ), K q(e) = (qij ), K α(e) = (αij ) размером
n × n будут отличны от нуля всего лишь девять элементов:
83

=
pij

 ∂N ∂N

∂N ∂N 

j
j
i
i
∫∫ p  ∂x ∂x + ∂y ∂y  dxdy,

Ω( e ) 

qij = ∫∫ qNi N j dxdy,

(4.67)

Ω( e )

αij=

∫ αNi N j ds,

∂Ω(2e )

где i, j принимают значения глобальных номеров узлов элемента:
=
Ni N=
N=
N3 .
1, N j
2 , Nk
Выражение для матрицы K α(e) зависит от того, по какой из сторон
производится интегрирование. Например, если интегрирование
производится по стороне между узлами i и j , то N k = 0 , если вдоль
стороны с узлами j и k , то Ni = 0 . Аналогично при интегрировании
по стороне между узлами k и i N j = 0 .
Обозначив стороны треугольного элемента через
ветственно, получим выражения для матрицы K α(e) :

K α(e)

K α(e)

K α(e)

 N2
Ni N j
 Ni 
 i
 
N 2j
=
α (e) ∫  N j  [ Ni N j 0]ds =
α (e) ∫  Ni N j
Lij 
Lij 

0 
 0
0


Lij , L jk , L ki

соот-

0 

0 ds,

0 


 N2
0
Ni N k 
 Ni 
(4.68)
 i

 
0
0 ds,
=
α (e) ∫ 0  [ Ni 0 N k ]ds =
α(e) ∫  0


L ki 
L ki

2 
N N
 Nk 
0
Nk
 i k

 0
0
0 
0 


 

(e)
(e)
2
Nj
N j N k ds.
=
α ∫  N j  [0 N j N k ]ds =
α ∫  0
L jk 
L jk 


 Nk 
 0
Nk N j
N k2 



Матрица

K (e) = K (pe) + K q(e) + K α(e)

кости элемента:
84

называется матрицей жест-

 p11

=
K
 p21
p
 31
(e)

p12
p22
p32

p13   q11
 
p23  +  q21
p33   q31

q12
q22
q32

q13   α11 α12
 
q23  +  α 21 α 22
q33   α31 α32

α13 

α 23  .
α33 

Аналогично вектор B (e) называется вектором нагрузки:
b1 
 
=
B
b2  , b1
=
b 
 3
(e)

=
=
, b2 ∫∫ fN j dxdy
, b3 ∫∫ fN k dxdy. (4.69)
∫∫ fNi dxdy

Ω( e )

Ω( e )

Ω( e )

Применение формул (4.65) приводит к методу сборки «по
элементам»: сначала вычисляются матрицы элементов, а затем путем
суммирования ‒ матрица глобальной системы. Этот метод отличается
от метода сборки «по узлам», когда рассматриваются лишь элементы,
примыкающие к данному узлу.
Для вычисления матрицы жесткости элемента и вектора нагрузки
используются функции формы треугольного элемента (2.14). Для
определения места матрицы элемента, при суммировании в глобальной матрице используется так называемая матрица инцидентности, которая ставит в соответствие каждому элементу глобальные номера узлов.
Функции p, q, α могут определяться экспериментально, например, в узлах конечноэлементной сетки или таблицей значений в узлах
некоторой прямоугольной сетки.
Для вычисления матрицы жесткости элемента может потребоваться численное интегрирование, в основном используются процедуры, основанные на квадратурных формулах Гаусса. Эти вопросы
выходят за рамки данного пособия.
Представим расчетные формулы для вычисления матрицы элемента, считая, что в пределах каждого элемента эти функции
постоянны:
p ( x, y ) =
p ( e ) , q ( x, y ) =
q ( e ) , α ( x, y ) =
α(e) ,

В соответствии с (2.14)
)
N j ( x, y=

1
(a j + b j x + c j y ), где =
j 1, 2,3,
2A

a1 =
x2 y3 − y2 x3 ; b1 =
y2 − y3 ; c1 =
x3 − x2 ;
a2 =−
x3 y1 y3 x1; b2 =
y3 − y1; c2 =
x1 − x3 ;
a3 =−
x1 y2 y1 x2 ; b3 =
y1 − y2 ; c3 =
x2 − x1
85

f ( x, y ) =
f (e) .

(4.70)

и
1
1
=
A
1
2
1

x1

y1

1
x2=
y2
(b1c2 − b2c1 ) площадь элемента.
2
x3 y3

и y j ( j = 1, 2,3) ‒ координаты вершин

В этих уравнениях x j

элемента.
Вычислим сначала матрицу K (pe) = ( pij ) . В пределах элемента
подынтегральная функция постоянна:
 ∂N ∂N j ∂Ni ∂N j  p (e)
p i
bi b j + ci c j
+
=
∂y ∂y  4 A2
 ∂x ∂x

(

)

(4.71)

Поэтому интегрирование сводится к умножению подынтегральной функции на площадь элемента A:
p (e)
pij =
bi b j + ci c j , i, j =1, 2,3 .
4A

(

)

(4.72)

Для вычисления интегралов ∫∫ qN N dxdy удобно использовать
i j
Ω( e )

барицентрические координаты L1 , L2 , L3 (2.19):
=
L1 N=
N=
1 , L2
2 , L3 N 3 .

(4.73)

L1 + L2 + L3 =
1 . Линейные уравнения Li ( x=
i 1, 2,3 пред, y ) 0,=
ставляют собой стороны треугольника П 2 П3 , П3П1 , П1П 2 ,
Li ( x=
, y ) 1,=
i 1, 2,3 в вершинах П1, П 2 , П3 .

Иначе говоря, треугольник П1П 2 П 3 в плоскости ( x, y )
преобразуется в стандартный треугольник P1P2P3 в плоскости ( L1 , L2 ) ,
причем P1 (1, 0) , P2 (0,1) , P3 (0, 0) (рис. 32). Обратное преобразование
плоскости ( L1 , L2 ) в плоскость ( x, y ) выражается формулами
x = x3 + ( x1 − x3 ) L1 + ( x2 − x3 ) L2 ,
y = y3 + ( y1 − y3 ) L1 + ( y2 − y3 ) L2 .
86

(4.74)

Рис. 32

Вычисляя якобиан преобразования, получим
∂x
∂x
∂L1 ∂L2
x1 − x3 x2 − x3
=
=
J ( L1 , L2 ) =
y1 − y3 y2 − y3
∂y
∂y
∂L1 ∂L2

(4.75)

= ( x1 − x3 )( y2 − y3 ) − ( x2 − x3 )( y1 − y3 ) = 2 A.

Таким образом,
(e)
=
∫∫ q Ni N j dxdy
П1П 2 П3

= 2 Aq

(e)
=
∫∫ q Li L j J ( L1, L2 ) dL1dL2

P1P2P3

(e)

∫∫ Li L j dL1dL2 .

(4.76)

P1P2P3

Выведем формулу для вычисления интегралов:
I (=
i, j , k )

i j k
i j
k
=
=
∫∫ L1L2 L3 dL
∫∫ L1L2 (1 − L1 − L2 ) dL
1dL2
1dL2
PP
1 2P3

PP
1 2P3
1

1− L1

0

0

(4.77)

= ∫ L1i dL1 ∫ L2j (1 − L1 − L2 ) k dL2 .

Интегрируем по частям:
1

1− L1

0

0

=
I (i, j , k ) ∫ L1i dL1 ∫ L2j +1 (1 − L1 − L2 ) k −1dL2 .

Следовательно,

I (i,=
j, k )

k
I (i, j + 1, k − 1).
j +1

Применяем последующее интегрирование по частям:
87

(4.78)

j !k !
j !k ! 1 i 1− L1 j + k
I (i, j=
dL2
+ k , 0)
∫ L1dL1 ∫ L2=
( j + k )!
( j + k )! 0
0

I (=
i, j , k )

=

j !k ! 1 i
j + k +1
dL1.,
∫ L1 (1 − L1 )
( j + k + 1)! 0

(4.79)

Дальнейшее интегрирование по частям:
I (i, j , k )

1
j !k !( j + k + 1)( j + k ) 1
i + j + k +1
(1 − L1 )0 dL1. (4.80)
∫L
( j + k + 1)!(i + 1)(i + 2)  (i + j + k + 1) 0 1

После элементарного интегрирования получим
=
I (i, j , k )

i j k
=
∫∫ L1L2 L3 dL1dL2

P1P2P3

i ! j !k !
.
(i + j + k + 2)!

(4.81)

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

Aq (e)
2!0!0!
)
2 0 0
(e)
=
L
L
L
dL
dL
Aq
2 Aq (e=
2
,
∫∫ 1 2 3 1 2
(2 + 0 + 0 + 2)!
6
P1P2P3

(e)
1 1 0
(e)
q=
∫∫ L1L2 L3dL1dL=
12 q=
21 2 Aq
2 2 Aq
P1P2P3

q=
13 q=
31 2 Aq

(e)

∫∫

P1P2P3

=
L11L02 L13dL1dL
2

2 Aq

(e)

Aq (e)
1!1!0!
=
,
(1 + 1 + 0 + 2)!
12
Aq (e)
1!0!1!
=
,
(1 + 0 + 1 + 2)!
12

q22

0!2!0!
Aq (e)
)
0 2 0
2 Aq (e=
,
Aq (e)
∫∫ L1 L2 L3dL1dL2 2=
(0 + 2 + 0 + 2)!
6
P1P2P3

q33

0!0!2!
Aq (e)
)
0 0 2
(e)
=
2 Aq (e=
2
.
L
L
L
dL
dL
Aq
∫∫ 1 2 3 1 2
(0 + 0 + 2 + 2)!
6
P1P2P3

(4.82)

Приведем вычисление матрицы K α(e) для случая, когда интегрирование осуществляется вдоль стороны между узлами i и j
(сторона П1П 2 ).
При использовании L координат эта сторона переходит в сторону
треугольника P1P2 , уравнение которой в координатах L1 , L2 имеет вид
1 − L1 − L2 =
0 (см. рис. 32). С учетом этого (4.74) принимают вид
x =x1 + ( x1 − x2 ) L1 ,
y =y1 + ( y1 − y2 ) L1.
88

(4.83)

Используя эти выражения для вычисления элемента длины дуги ds,
получим
ds =

dx 2 + dy 2 =

( x1 − x2 ) 2 + ( y1 − y2 ) 2 dL1 = Lij dL1 ,

(4.84)

где Lij – длина рассматриваемой стороны треугольника. Таким
образом,

 L12
L1L2
0 
 Ni 
1



α (e) ∫  N j  [ Ni N j 0]ds =
α (e)Lij ∫  L2 L1 L22
K α(e) =
0 dL1 .

0
Lij 

 0
0
0 
0 


Вычислив интегралы, получим
1

(4.85)

1

∫ L1 dL1 = 3 ,
2

0
1

1

1

2
2
,
∫ L2 dL1 =
∫ (1 − L1 ) dL1 =
3

0
1

(4.86)

0

1

1

∫ L1L2 dL1 =∫ L1 (1 − L1 )dL1 =6.

0

Матрица

K α(e)

0

имеет вид

2 1 0
α (e)Lij 

=
1 2 0 .

6 

0 0 0
Аналогично для других случаев
2 0 1
α (e)Lki 

0 0 0,
K α(e) =
6 

1 0 2
K α(e)

α L jk
(e)

K α(e) =

6

0 0 0


 0 2 1 .
 0 1 2



(4.87)

(4.88)

К этим же результатам можно прийти, применяя формулы
a b c
=
∫ L1 L2 L3 d γ
γ

a !b !c !
γ,
(a + b + c + 1)!

a !b !c !
∫∫ L1 L2 L3 dA = (a + b + c + 2)! A.
A
a b c

89

(4.89)

Эти формулы были предложены для вычисления интегралов по
стороне и по площади элемента в 1973 г. в статье Martin A. Eisenberg,
Lawrence E. Malvern «On finite element integration in natural coordinate»
в журнале «International Journal for Numerical Methods in Engineering»
(Volume 7, Issue 4, 1973).
Вычислим составляющие вектора нагрузки для элемента (4.69)
B (e) :
b1
=

Ω( e )

b2
=

=
∫∫ fN j dxdy f

Ω( e )

b3
=

1!0!0!
A
f (e) A f (e) ,
=
(1 + 0 + 0 + 2)!
6

(e)
L11L02 L03dA
=
=
∫∫ fNi dxdy f ∫∫
A
(e)

0 1 0
=
∫∫ L1 L2 L3dA
A

(e)
0 0 1
=
∫∫ fN k dxdy f =
∫∫ L1 L2 L3dA

Ω( e )

A

0!1!0!
A (4.90)
f (e) A f (e) ,
=
(0 + 1 + 0 + 2)!
6
0!1!0!
A
f (e) A f (e) .
=
(0 + 1 + 0 + 2)!
6

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

Рис. 33
90

В настоящем пособии для триангуляции Делоне-области использована имеющаяся в свободном доступе программа triangle J. R.
Shewchuk [8].
Основные параметры программы и примеры. Исходные данные
программы записываются в виде текстового файла с расширением
*.poly.
Сама программа работает из командной строки, где в качестве
параметров указываются различные ключи, определяющие режим
работы программы (с полным списком параметров можно ознакомиться в приложении или в первоисточнике).
В результате работы программы создаются файлы, которые можно использовать для последующего разбиения. Например, эта серия создана путем последовательного выполнения команд
triangle -pqc box.poly, triangle -rpa0.2 box.1, triangle -rpa.05 box.2,
triangle -rpa.0125 box.3
Графический вывод триангуляции области (рис. 34) можно
осуществить с помощью программы Matlab show.m.

Рис. 34
91

При этом необходимо указать имя файла с данными области.
В качестве дополнительных опций можно выводить номера узлов, а
также номера элементов.
Прокомментируем содержание исходного файла box.poly с описанием планарного графа. Этот файл содержит разделы с описанием
вершин, их координат, сегментов:
# прямоугольник с 8 вершинами, один маркер.
8201
# Внешний прямоугольник состоит из узлов:
1 0 0 0
2 0 3 0
3 3 0 0
4 3 3 33
# Внутренний квадрат состоит из узлов:
5 11 0
6 12 0
7 21 0
8 22 0
# Пять сегментов с граничными маркерами.
5 1
1 1 2 5 # Левая сторона внешнего прямоугольника.
# Квадратное отверстие состоит из сегментов:
2 57 0
3 78 0
4 8 6 10
5 65 0
# одно отверстие в середине.
1
1 1.5 1.5
В результате выполнения программы (triangle -rpa0.2 box.1)
создаются текстовые файлы с данными о координатах узлов
и принадлежности узлов к той или иной границе (файл box.2.node)
и файл box.2.ele c данными о номерах узлов для каждого элемента.
Если граница области является сложной, то для получения достоверной картины приходится вводить довольно большой объем
данных.
Этот процесс отчасти можно автоматизировать, представляя границу области в виде совокупности участков, описываемых линейными или квадратичными сплайнами.
Представим границу области в виде системы из трех участков,
описываемых линейными или квадратичными сплайнами (рис. 35).
92

Рис. 35

Программа обрабатывает эти данные и создает два текстовых
файла с описаниями узлов и сегментов, принадлежащих границе
области. По этим файлам программа triangle строит триангуляцию
области (рис. 36).
dant.1

5

4

--Y axis--

3

2

1

0

-1

1

2

3

4

6
5
--X axis--

7

8

9

10

Рис. 36

4.5. Примеры решения задач
Далее на той же сетке (см. рис. 28) рассмотрим задачу для
уравнения Лапласа c условиями смешанного типа:
 ∂ 2u ∂ 2u
Ω [0;1] × [0;1],
0,
( x, y ) ∈=
 2+ =
∂y 2
 ∂x
=
u (0, y ) 0,
y ∈ [0;1],

 ∂u
=
y ) 0, y ∈ [0;1],
 (1,
 ∂x
x, 0) 0, x ∈ [0;1],
u (=

u
x) =100sin(πx / 2), x ∈ [0;1].
(1,



93

(4.91)

Аналитическое решение этой задачи:
u ( x, y ) =

100 sh(πy / 2) sin(πx / 2)
.
sh(π / 2)

(4.92)

Программа решения этой задачи приведена ниже. Исходные данные в примере вводятся в тексте программы. К ним относятся данные
о координатах узлов конечноэлементной сетки, номерах узлов для
каждого элемента, перечисление граничных узлов и значений искомой функции в узлах в сответствии с граничными условиями. Следует
отметить, что в данном случае выполняются лишь главные граничные
условия.
В отличие от примера, рассмотренного в подразд. 4.2, сборка
системы осуществлялась «по элементам» с использованием формул,
представленных выше. Матрицы систем линейных алгебраических
уравнений до и после выполнения граничных условий представлены
формулами (4.93) и (4.94).
 2 -1 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0   u1   0 

   
1 4 -1 0 0 0 - 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0   u2   0 
 0 -1 4 -1 0 0 0 - 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0   u3   0 

   
 0 0 -1 4 -1 0 0 0 - 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0   u4   0 
 0 0 0 -1 2 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0   u5   0 

   
1 0 0 0 0 4 - 2 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0   u6   0 
0 - 2 0 0 0 - 2 8- 2 0 0 0 - 2 0 0 0 0 0 0 0 0 0 0 0 0 0 u  0

 7   
 0 0 - 2 0 0 0 - 2 8 - 2 0 0 0 - 2 0 0 0 0 0 0 0 0 0 0 0 0   u8   0 
0 0 0 - 2 0 0 0 - 2 8- 2 0 0 0 - 2 0 0 0 0 0 0 0 0 0 0 0 u  0

 9   
 0 0 0 0 -1 0 0 0 - 2 4 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0   u10   0 

   
 0 0 0 0 0 -1 0 0 0 0 4 - 2 0 0 0 -1 0 0 0 0 0 0 0 0 0   u11   0 
 0 0 0 0 0 0 - 2 0 0 0 - 2 8 - 2 0 0 0 - 2 0 0 0 0 0 0 0 0   u12   0 

   
 0 0 0 0 0 0 0 - 2 0 0 0 - 2 8 - 2 0 0 0 - 2 0 0 0 0 0 0 0   u13  =  0 
0 0 0 0 0 0 0 0 - 2 0 0 0 - 2 8- 2 0 0 0 - 2 0 0 0 0 0 0 u  0

  14   
 0 0 0 0 0 0 0 0 0 -1 0 0 0 - 2 4 0 0 0 0 -1 0 0 0 0 0   u15   0 
 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 4 - 2 0 0 0 -1 0 0 0 0   u   0 

  16   
 0 0 0 0 0 0 0 0 0 0 0 - 2 0 0 0 - 2 8 - 2 0 0 0 - 2 0 0 0   u17   0 

   
 0 0 0 0 0 0 0 0 0 0 0 0 - 2 0 0 0 - 2 8 - 2 0 0 0 - 2 0 0   u18   0 
 0 0 0 0 0 0 0 0 0 0 0 0 0 - 2 0 0 0 - 2 8 - 2 0 0 0 - 2 0   u19   0 

   
 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 - 2 4 0 0 0 0 -1   u20   0 
 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 2 -1 0 0 0   u   0 

  21   
 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 - 2 0 0 0 -1 4 -1 0 0   u22   0 
 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 - 2 0 0 0 -1 4 -1 0   u   0 

  23   
 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 - 2 0 0 0 -1 4 -1   u24   0 

   
 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 -1 2   u25   0 

94

(4.93)

0 
1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0   u1  
 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0  u  

38,2683
2

  

 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0   u3   70.7107 

  

 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0   u4   92.3880 
 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0   u5   1000 

  

0 
 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0   u6  
 0 -1 0 0 0 -1 4 -1 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0   u  
0 

 7  

0 
 0 0 -1 0 0 0 -1 4 -1 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0   u8  
 0 0 0 -1 0 0 0 -1 4 -1 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0   u  
0 

 9  

0 
 0 0 0 0 -1 0 0 0 -1 2 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0   u10  
 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0  u  
0 

  11  
 0 0 0 0 0 0 -1 0 0 0 -1 4 -1 0 0 0 -1 0 0 0 0 0 0 0 0   u12  
0 

  

0 
 0 0 0 0 0 0 0 -1 0 0 0 -1 4 -1 0 0 0 -1 0 0 0 0 0 0 0   u13  = 
 0 0 0 0 0 0 0 0 -1 0 0 0 -1 4 -1 0 0 0 -1 0 0 0 0 0 0   u14  
0 

  

0 
 0 0 0 0 0 0 0 0 0 -1 0 0 0 -1 2 0 0 0 0 -1 0 0 0 0 0   u15  
 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0  u  
0 

  16  

u
0
0
0
0
0
0
0
0
0
0
0
-1
0
0
0
-1
4
-1
0
0
0
-1
0
0
0
0 

  17  
 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 -1 4 -1 0 0 0 -1 0 0   u  
0 

  18  

0 
 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 -1 4 -1 0 0 0 -1 0   u19  
 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 -1 2 0 0 0 0 -1   u  
0 

  20  
 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0   u21  
0 

  

u
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
1
0
0
0
0 

  22  
 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0   u23  
0 

  

u
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
1
0
0 

  24  
 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1  u  
0 

 25  

(4.94)

Результаты приближенного и аналитического решения
приведены графически в виде изолиний на рис. 37 и 38 соответственно. Следует отметить неплохое количественное соответствие,
несмотря на сравнительно грубую сетку.
Листинг программы, реализующей решение задачи Дирихле
для уравнения Лапласа методом конечных элементов
function fem2D()
% Описание переменных
% k = матрица жесткости элемента
% f = вектор нагрузки элемента
% kk = результирующая матрица
% ff = результирующий вектор
% gcoord = значения координат для каждого узла
% nodes = номера узлов для каждого элемента
95

% bcdof = вектор, содержащий номера узлов связанных с
%граничными условиями
% bcval = вектор, содержащий величины значений
%функции для узлов связанных с граничными условиями
%----------------------------nel=32; % число элементов
nnel=3;% число узлов для одного элемента
ndof=1;
nnode=25; % число узлов
sdof=nnode*ndof;
% Ввод координат для всех узлов
gcoord(1,1)=0.0;
gcoord(1,2)= 1.00;
gcoord(2,1)=0.25;
gcoord(2,2)= 1.00;
gcoord(3,1)=0.5;
gcoord(3,2)= 1.00;
gcoord( 4,1 )=0.75;
gcoord( 4,2)= 1.00;
gcoord(5,1)=1.0;
gcoord(5,2)= 1.00;
gcoord(6,1)=0.0;
gcoord(6,2)= 0.75;
gcoord(7,1 )=0.25;
gcoord(7,2)= 0.75;
gcoord(8,1)=0.5;
gcoord(8,2)= 0.75;
gcoord(9,1)=0.75;
gcoord(9,2)= 0.75;
gcoord(10,1)=1.0;
gcoord(10,2)= 0.75;
gcoord(11,1)=0.0;
gcoord(11,2 )=0.50;
gcoord(12,1)=0.25;
gcoord(12,2)=0.50;
gcoord(13,1 )=0.5;
gcoord(13,2)=0.50;
gcoord(14,1)=0.75;
gcoord(14,2)=0.50;
gcoord(15,1 )=1.0;
gcoord(15,2)=0.50;
gcoord(16,1)=0.0;
gcoord(16,2)=0.25;
gcoord(17,1)=0.25;
gcoord(17,2)=0.25;
gcoord(18,1)=0.5;
gcoord(18,2)=0.25;
gcoord(19,1)=0.75;
gcoord(19,2)=0.25;
gcoord(20,1)=1.0;
gcoord(20,2)=0.25;
gcoord(21,1)=0.0;
gcoord(21,2)=0.00;
gcoord(22,1)=0.25;
gcoord(22,2)=0.00;
gcoord(23,1)=0.5;
gcoord(23,2)=0.00;
gcoord(24,1)=0.75;
gcoord(24,2)=0.00;
gcoord(25,1)=1.0;
gcoord(25,2)=0.;
% Ввод данных о номерах узлов для элемента
nodes(1,1)=1;
nodes(1,2)=6; nodes(1,3)=2;
nodes(2,1)=6;
nodes(2,2)=7; nodes(2,3)=2;
96

nodes(3,1)=7;
nodes(3,2)=3; nodes(3,3)=2;
nodes(4,1)=7;
nodes(4,2)=8; nodes(4,3)=3;
nodes(5,1)=8;
nodes(5,2)=4; nodes(5,3)=3;
nodes(6,1)=8;
nodes(6,2)=9; nodes(6,3)=4;
nodes(7,1)=9;
nodes(7,2)=5; nodes(7,3)=4;
nodes(8,1)=9;
nodes(8,2)=10; nodes(8,3)=5;
nodes(9,1)=11;
nodes(9,2)=7; nodes(9,3)=6;
nodes(10,1)=11; nodes(10,2)=12; nodes(10,3)=7;
nodes(11,1)=12; nodes(11,2)=8; nodes(11,3)=7;
nodes(12,1)=12; nodes(12,2)=13; nodes(12,3)=8;
nodes(13,1 )=13; nodes(13,2)=9; nodes(13,3)=8;
nodes(14,1)=13; nodes(14,2)=14; nodes(14,3)=9;
nodes(15,1)=14; nodes(15,2)=10; nodes(15,3)=9;
nodes(16,1)=14; nodes(16,2)=15; nodes(16,3)=10;
nodes(17,1)=16; nodes(17,2)=12; nodes(17,3)=11;
nodes(18,1)=16; nodes(18,2)=17; nodes(18,3)=12;
nodes(19,1)=17; nodes(19,2)=13; nodes(19,3)=12;
nodes(20,1)=17; nodes(20,2)=18; nodes(20,3)=13;
nodes(21,1)=18; nodes(21,2)=14; nodes(21,3)=13;
nodes(22,1)=18; nodes(22,2)=19; nodes(22,3)=14;
nodes(23,1)=19; nodes(23,2)=15; nodes(23,3 )=14;
nodes(24,1)=19; nodes(24,2)=20; nodes(24,3)=15;
nodes(25,1 )=21; nodes(25,2)=17; nodes(25,3)=16;
nodes(26,1 )=21; nodes(26,2)=22; nodes(26,3)=17;
nodes(27,1)=22; nodes(27,2)=18; nodes(27,3)=17;
nodes(28,1)=22; nodes(28,2)=23; nodes(28,3)=18;
nodes(29,1)=23; nodes(29,2)=19; nodes(29,3)=18;
nodes(30,1)=23; nodes(30,2)=24; nodes(30,3)=19;
nodes(31,1)=24; nodes(31,2)=20; nodes(31,3)=19;
nodes(32,1)=24; nodes(32,2)=25; nodes(32,3)=20;
% Ввод данных для граничных условий
bcdof(1)=1;
bcval(1)=0;
bcdof(2)=22; % номер узла на границе
bcval(2)=0; % значение функции в этом узле
bcdof(3)=23;
bcval(3)=0;
bcdof(4)=24;
bcval(4)=0;
97

bcdof(5)=5;
bcval(5)=100;
bcdof(6)=6;
bcval(6)=0;
bcdof(7)=11;
bcval(7)=0;
bcdof(8)=16;
bcval(8)=0;
bcdof(9)=21;
bcval(9)=0;
bcdof(10)=2;
bcval(10)=bf(0.25,1.);
bcdof(11)=3;
bcval(11)=bf(0.5,1.);
bcdof(12)=4;
bcval(12)=bf(0.75,1.);
bcdof(13)=25;
bcval(13)=0;
%
%--------------% инициализация матриц и векторов
%,---------------ff =zeros( sdof, 1);
kk=zeros( sdof,sdof);
%
%----------------------% Цикл по элементам вычисление матриц и векторов
элементов и их сборка
%1--------------------for iel=1:nel % цикл по полному числу элементов
iel
%
nd(1)=nodes(iel,1); % первый узел для (iel)-го элемента
nd(2)=nodes(iel,2); % второй узел для (iel)-го элемента
nd(3)=nodes(iel,3); % третьий узел для (iel)-го элемента
x1=gcoord(nd(1),1); y1=gcoord(nd(1),2); % координаты
x2=gcoord(nd(2),1); y2=gcoord(nd(2),2); % вершин
x3=gcoord(nd(3),1); y3=gcoord(nd(3),2); % треугольника
%Формирование матрицы и вектора элемента
k=zeros(3,3);
98

A=0.5*(x2*y3+x1*y2+x3*y1-x2*y1-x1*y3-x3*y2);
% А-площадь треугольника
k(1,1)=((x3-x2)*(x3-x2)+(y2-y3)*(y2-y3))/(4*A)
k(1,2)=((x3-x2)*(x1-x3)+(y2-y3)*(y3-y1) )/(4*A);
k(1,3)=((x3-x2)*(x2-x1)+(y2-y3)*(y1-y2))/(4*A);
k(2,1)=k(1,2);
k(2,2)=( (x1-x3)*(x1-x3)+(y3-y1)*(y3-y1) )/(4*A);
k(2,3)=((x1-x3)*(x2-x1)+(y3-y1)*(y1-y2))/(4*A);
k(3,1)=k(1,3);
k(3,2)=k(2,3);
k(3,3)=((x2-x1)*(x2-x1)+(y1-y2)*(y1-y2))/(4*A);
% Сборка матрицы системы
for i=1:3
ii=nd(i);
for j=1:3
jj=nd(j);
kk(ii,jj )=kk(ii,jj) + k(i,j);
end
end
end
% Выполнение граничных условий
[kk, ff]=feaplyc2 (kk,ff, bcdof, bcval );
%Решение системы
fsol=kk\ff;
% вычисление аналитического решения в узлах
isolines(fsol);
for i=1:nnode
x=gcoord(i,1); y=gcoord(i,2);
esol(i)=bf(x,y);
end
% Печать точного и конечноэлементного решения
num=1:1:sdof;
store=[num' fsol esol']
function bb=bf(x,y)
bb=100*sinh(3.1415927*y/2)*sin(3.1415927*x/2)/sinh(3.1415927/2);
%--------------------------function [kk,ff]=feaplyc2(kk,ff,bcdof,bcval)
% Цель:
% Применение граничных условий
99

n=length(bcdof)
sdof=size(kk)
%
for i=1:n
c=bcdof(i);
for j=1:sdof
kk(c,j)=0;
end
%
kk(c,c)=1;
ff(c)=bcval(i);
end
%---------------------------function [iso]= isolines (B)
n1=25;
n2=25;
H=repmat(0,n1,n2);
XMIN=0.;
XMAX=1.;
YMIN=0.;
YMAX=1.;
DX=(XMAX-XMIN)/(n1);
DY=(YMAX-YMIN)/(n2);
XI=zeros(n1,n2);
YI=zeros(n1,n2);
for I=1:n1
for J=1:n2
YI(I,J)=YMIN+J*DY;
XI(I,J)=XMIN+I*DX;
H(I,J)=bf(XI(I,J),YI(I,J));
end
end
Z=zeros(5,5);
for I=1:5
for J=1:5
Z(I,J)=B(J+(I-1)*5);
end
end
X1=zeros(5,5);
100

Y1=zeros(5,5);
for I=1:5
for J=1:5
Y1(I,J)=1-(I-1)*0.25;
X1(I,J)=(J-1)*0.25;
end
end
ZI = interp2(X1, Y1, Z, XI, YI);
figure(1);
[c,h] = contour(XI,YI,H);
grid on;
clabel(c,h);
title('\fontsize{14} Аналитическое решение');
hold off;
figure(2);
[c,h] = contour(XI,YI,ZI);
grid on
clabel(c,h);
title('\fontsize{14} Метод конечных элементов (сетка 5Х5)');
hold off;
Метод конечных элементов (сетка 5Х5)
70

80

90

60

40

50

30

20

0.9

10

1

80
70

0.8

30

40
50

20

10

0.6

60

50

0.7

40
0.5

30
30

0.4

20

0.3

10

20

0.2

10

10

0.1
0.1

0.2

0.3

0.4

0.5

Рис. 37
101

0.6

0.7

0.8

0.9

1

Аналитическое решение
40

50

30

20

0.9

10

1

60

70

80

90
80

70

0.8

30

40
50

20

10

0.6

60

50

0.7

40
0.5

30

0.4

30

20

0.3

10

20

0.2

10

10

0.1
0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Рис. 38

Далее рассмотрим уравнение Лапласа с краевыми условиями
смешанного типа:
 ∂ 2u ∂ 2u
+ 2 0,
( x=
, y ) ∈ Ω {( x, y ) | 0 < x < 1; 1 − x < y < 1},
 =
2
∂y
 ∂x
u (1, y ) =
y (1 − y ),
y ∈ [0;1],

 ∂u
0, y =
1 − x, x ∈ [0;1]
 ∂n =

x(1 x), x ∈ [0;1].
u (1, x) =−

(4.95)

Аналитическое решение этой задачи:
∞

sh(π(2m − 1) / 2) ch(π(2m − 1) / 2)(1 − 2 x) / 2) sin(π(2m − 1) y )

m =1

(π(2m − 1))3 sh(π(2m − 1))

u 16 ∑

+

∞

sh(π(2m − 1) / 2) ch(π(2m − 1) / 2)(1 − 2 y ) / 2) sin(π(2m − 1) x)

m =1

(π(2m − 1))3 sh(π(2m − 1))

+16 ∑

(4.96)
.

Для решения задачи использовалась конечноэлементная сетка
(рис. 39), подготовленная с помощью программы triangle.
102

1

1

49
42

0.9

Конечно-элементная сетка
39
2
3
4
26
46

43

0.8

12

38

31

23

20

45

11

0.4

6

18

41

29

24

40

0.5

5

33
13

15
0.6

21
25

22

37

0.7

--ось Y --

28

44
16

27

36
30

14

32

7

34
35

19

0.3

17

10

47

50

8

48

51

0.2
0.1
9

0
0

0.2

0.4
0.6
--ось X --

0.8

1

Рис. 39

Результаты конечноэлементного и аналитического решения
задачи графически, в виде линий равных значений, приведены на
рис. 40 и 41 соответственно.
Конечно-элементное решение
1

0.2
07
0.20 69
0 35

0.8

000..0.1.1
1144434
881963
9951

7
69
12 32
0. 563 66
0...1116711 4
0
0
0...11778 3 8
00...11.1187855668
00.1
00.19
011
9330

55
033
0000011
0..22336
88
.10.118899556
4
000..1
7 811334
0
.1
000.1
.17 66
0.163

0. 0 00.2
2000.0..222.222229977
07..6221155023366
209 033
76
9

0.9

0.7

0.
2
00..22 15 0
0
0..2 22223 3
29 7 66
7

355
00.2.200003 7 69
0
00..220 03
15
2
.
0

0.6
0.5

0.2
007.
6290
7
00.2.200609
03355
0.1930
1
0.18568
0.17834
0.171

0.4
0.3
0.2
0.1
0
0

0.2

0.4

Рис. 40
103

0.6

0.8

1

Аналитическое решение
1

0.
20
0.2 9 82
03 4
1
0.19
7

0.8

0.197
9 82
0.203 410.20

0.
0.2 216
0.22 22 622
9 033

0.7

-- ось Y --

00. .1
133
9239
65

67
5271 7
145 8 598
0 .115 4 38
00...1167717 799
00.1 84 1
1 90 6
00..1

2
98
20 41
00. .29073
0.1 6
0.190

0.235 44
0.2
03
00.2.22229
16 63
22

0.9

0.6
0.5

0.2
09
0.2082
3 41
0.1
7
0.19906
0.184 19

0.4
0.3
0.2
0.1
0
0

0.2

0.4
0.6
--ось X --

0.8

1

Рис. 41

Сопоставление точного и приближенного решения в узлах
конечноэлементной сетки приведено в табл. 11.
Т а б л и ц а 11
Номер узла
1
2
3
4
5
6
7
9
10
11
12
13
14
15
16
17
18
19

Точное решение
0
0.1875
0.2500
0.1875
0
0.1875
0.2500
0
0.1954
0.2034
0.1795
0.2022
0.2077
0.1997
0.2047
0.2110
0.1963
0.2049
104

Приближенное решение
0.0000
0.1875
0.2500
0.1875
0.0000
0.1875
0.2500
0.0000
0.1939
0.2053
0.1939
0.2041
0.2079
0.2046
0.2106
0.2109
0.1945
0.2046

Номер узла
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51

Точное решение
0.1489
0.1094
0.1886
0.1972
0.2033
0.1697
0.2344
0.2016
0.2166
0.2038
0.2182
0.2025
0.2047
0.2067
0.2141
0.2064
0.2344
0.2005
0.1928
0.2344
0.2019
0.2033
0.0845
0.2209
0.1950
0.1094
0.1557
0.2344
0.1475
0
0.1896
0.1094

Приближенное решение
0.1455
0.1094
0.1887
0.1972
0.2042
0.1691
0.2344
0.2021
0.2188
0.2039
0.2177
0.2057
0.2055
0.2101
0.2144
0.2058
0.2344
0.2059
0.2023
0.2344
0.2053
0.2018
0.1459
0.2243
0.2061
0.1094
0.1807
0.2344
0.1459
0.1094
0.1882
0.1094

В качестве следующего примера рассмотрим задачу восстановления вещественной части аналитической функции u = z2 по её
значениям на границе области:
 ∂ 2u ∂ 2u
+ 2 0;
( x, y ) ∈ Ω,
=
2
∂y
 ∂x

u x2 − y 2 ,
( x, y ) ∈ ∂Ω.
=

(4.97)

Для решения задачи использовалась разбиение нерегулярной
двусвязной области программой triangle (рис. 42). Исходные данные
105

для триангуляции были получены путем представления границ в виде
совокупности линейных и квадратических сплайнов.
Конечно-элементная сетка
16

14

ось --Y --

12

10

8

6

4

2
0

2

4

6

8
10
ось--X --

12

14

16

Рис. 42
0

100

200

300

400

500

600
0

100

200

300
nz = 3351

Рис. 43
106

400

500

600

18

Матрица системы является разреженной и имеет ленточную
структуру. Из рис. 43 видно, что при выбранном способе нумерации
узлов сетки ширина ленты не является оптимальной. Вопросы
оптимизации ширины ленты за счет изменения нумерации узлов
выходят за рамки данного пособия и рассмотрены в работе [3].
Результаты решения на рис. 44 соответствуют аналитическому
решению u = x2 ‒ y2.
Метод конечных-элементов
16

-2 00

0
-1 0

0
-1 5

0
-5

0

14

50
0
-1 0

ось--Y --

12

-5 0

10

0

10
0

-1 50

50

8

0

-5 0

50

6

0

4

2
0

2

4

6

10
8
ось--X --

12

14

16

18

Рис. 44

4.6. Задания для практических занятий
Билинейная интерполяция. Представление данных
в многосвязной области
Пусть в узлах некоторой прямоугольной
сетки
заданы
значения функции H(x, y) (рис. 45) и
[ xi , xi +1 ] × [ y j , y j +1 ] ‒ типичный элемент, в
вершинах которого заданы функции
=
Hij H(
=
xi , y j ],
Hi +1. j H( xi +1 , y j ],
H(
=
xi , y j +1 ] , Hi +1. j +1 H( xi +1 , y j +1 ] .
107

Рис. 45

Требуется

определить
( x, y ) ∈ [ xi , xi +1 ] × [ y j , y j +1 ] .

значение

в

H( x, y )

Рассмотрим сначала значение функции
( xi , y j ) и ( xi +1 , y j ) .

H( x, y )

точке

в точках

Интерполируя, получим
xi +1 − x
x − xi
Н( x, y j )
Н( xi , y j ) +
Н( xi +1 , y j ) ,
=
xi +1 − xi
xi +1 − xi
xi +1 − x
x − xi
Н( x, y j +1 )
Н( xi , y j +1 ) +
Н( xi +1 , y j +1 ) .
=
xi +1 − xi
xi +1 − xi
Затем для точки ( x, y ) запишем:
y j +1 − y
y − yj
Н( x, y )
Н( x, y j ) +
Н( x, y j +1 ) .
=
y j +1 − y j
y j +1 − y j
Таким образом, для определения функции
воспользоваться выражением

( xi +1 − x)( y j +1 − y )Hij + ( xi +1 − x)( y − y j )Hij +1

H( x, y )

( xi +1 − xi )( y j +1 − y j )
+

можно

H( x, y )

+

( x − xi )( y j +1 − y )Hi +1 j + ( x − xi )( y j − y )Hi +1 j +1
( xi +1 − xi )( y j +1 − y j )

.

Перейдем здесь к координатам
x=

xi + xi +1
+ x′,
2

y=

y j + y j +1
2

+ y′.

Тогда


x −x
  y j +1 − y j
( xi +1 − xi )( y j +1 − y j )H( x′, y′=
)  i +1 i − x′  
− y′  Hij +
2
2



y j +1 − y j
x −x

+  i +1 i − x′   y′ +
2
2




 Hij +1 +



x −x
  y j +1 − y j
+  i +1 i + x′  
− y′  Hi +1 j +
2
2




x −x
  y j +1 − y j
+  i +1 i + x′  
+ y′  Hi +1 j +1.
2
2




Введем обозначения:
hx = xi+1 – xi,

hy = yj+1 – yj,
108

hy

h
  hy
h

hx hy H ( x′, y′) =  x − x′   − y′  H ij +  x − x′   y′ +
2
 2
 2
 2




 H ij +1 +



h
  hy
h
  hy
+  x + x′   − y′  H i +1 j +  x + x′   + y′  H i +1 j +1.
 2
 2
 2
 2


′
′
2x
2 y . В результате
Перейдем теперь к координатам
=
ξ =
, η
hx
hy

приходим к функциям формы четырехугольного элемента в естественной системе координат (2.24), полученной ранее:

(1 − ξ )(1 − η ) H + (1 − ξ )(1 + η ) H + (1 + ξ )(1 − η ) H +
H(ξ ,η ) =
ij
ij +1
i +1 j
4
4
4
(1 + ξ )(1 + η ) H
.
+
i +1 j +1
4
Задание 1. Пусть заданы =
A [a, b] × [c, d ],
=

{( x, y) ∈ A : ( x − x0 )2 + ( y − y0 )2 < R20 } ,
В

узлах

прямоугольной

xi = a + (i − 1)dx,

сетки

R=

n.
1 ≤ i ≤ n, 1 ≤ j ≤ n

где dx =

y j = c + ( j − 1)dy,

R0 , x0 , y0 ,

[ xi , y j ],

b−a
d − c , зада, dy =
n −1
n −1

ны значения функции H(x, y).
Требуется с помощью билинейной интерполяции определить
функцию H(x, y) и нарисовать линии постоянного уровня в двусвязной области Ω = A \ R .
Исходные данные для различных вариантов представлены ниже.
=
a 2,=
b 16,=
c 2,=
d 16, =
R0 2, =
n 10 .
№

x0

y0

H(x, y)

1

5

11

x2 + y 2

2

4

10

x2 − y 2

3

8

12

xy

109

№

x0

y0

H(x, y)

4

5

7

x3 + y 3

5

5

4

x2 − y 2

6

8

10

x3 − y 3

7

5

10

x2 y 2

8

7

10

x4 + y 4

9

5

11

x2 + y 2

10

8

10

x2 + y 2

11

7

9

x3 − y 3

12

11

10

x3 y 3

13

8

8

x2 + y 2

14

8

10

x2 − y 2

15

5

11

x2 y 2

16

5

12

x3 + y 3

17

11

10

x2 + y 2

18

5

8

x2 + y 2

Пример. Рассмотрим случай
=
a 2,=
b 16,=
c 2,=
d 16,
R=
0 2, x=
0 5, y=
0 10,

=
n 20, H( x, y=
) x2 + y 2 .

Для решения этой задачи используем Matlab. Результаты
приведены на рис. 46.
110

16

30 0

35
0

25 0

14

20 0
25
0

12 15 0

45
0

40
0
30
0

35
0

20
0

10
0
30

6

0
25

10
0

0
20

15
0

8
50

6

8

10

12

20 0

4

15 0

10 0

50

2
2

14

25 0

4

16

Рис. 46

Листинг программы
Вычисление интерполирующей функции выполнено в виде
функции HIJ(x, y).
function H2=HIJ(XT,YT)
n1=20;
n2=20;
H=zeros(n1,n2);
XMIN=-2.;
XMAX=20.;
YMIN=-2.;
YMAX=20.;
X=linspace(XMIN,XMAX,n1);
Y=linspace(YMIN,YMAX,n2);
for I=1:n1
for J=1:n2
H(I,J)=X(I)^2+Y(J)^2;
end
end
for I=1:n1-1
111

for J=1:n2-1
if(((XT-X(I))*(X(I+1)-XT)>=0)&((YT-Y(J))*(Y(J+1)-YT)>=0))
H2=((X(I+1)-XT)*(Y(J+1)-YT)*H(I,J)+(XT-X(I))*(Y(J+1)YT)*H(I+1,J)+(X(I+1)-XT)*(YT-Y(J))*H(I,J+1)+(X(I)-XT)*(Y(J)YT)*H(I+1,J+1))/((X(I+1)-X(I))*(Y(J+1)-Y(J)));
end
end
end
Главная программа, вычерчивающая линии:
n1=200;
n2=200;
H=repmat(0,n1,n2);
XMIN=2.;
XMAX=16.;
YMIN=2.;
YMAX=16.;
DX=(XMAX-XMIN)/n1;
DY=(YMAX-YMIN)/n2;
XI=zeros(n1,n2);
YI=zeros(n1,n2);
for I=1:n1
for J=1:n2
YI(I,J)=YMIN+J*DY;
XI(I,J)=XMIN+I*DX;
H(I,J)=HIJ(XI(I,J),YI(I,J));
if((XI(I,J)-5)^2+(YI(I,J)-10)^2<4)
H(I,J)=NaN;
end
end
end
figure(1);
t = 0:0.1:2*pi;
x=5+2*cos(t);
y=10+2*sin(t);
plot(x,y,'r');
grid on
hold on
[c,h] = contour(XI,YI,H);
clabel(c,h);
112

Задание 2. Аппроксимация границ области сплайнфункциями. Для аппроксимации нерегулярных участков границы
можно использовать сплайн-функции.
Рассмотрим процесс построения квадратичного сплайна. Пусть в
узлах x1 , x2 , xn известны значения функции f1 , f 2 ,  f n .
Содержание метода сводится к представлению на каждом из
отрезков x ∈ [ xi , xi +1 ], i =
1, 2  n − 1 аппроксимирующей функции
многочленом второй степени Si=
( x) , i 1, 2  n − 1 . При этом
необходимо
обеспечить
выполнение
условий
интерполяции, а на внутренних границах ‒ условия сопряжения производных.
Удовлетворяя условиям интерполяции и сопряжения, положим
f −f
Si ( x) = fi + i +1 i ( x − xi ) + ci ( x − xi )( x − xi +1 ) , i =1,  n − 1.
xi +1 − xi
Из

условия

равенства

'
производных =
,
S 'i ( xi ) S=
i −1 ( xi ) , i 1

= 1, 2,  , n − 1 получим реккурентное соотношение для вычисления
коэффициентов сплайна сi :
сi =−сi −1

( xi − xi −1 )  fi − fi −1 fi +1 − fi
+−
+
( xi +1 − xi )  xi − xi −1 xi +1 − xi


1
, i =2,  , n − 1.

 ( xi +1 − xi )
дважды,
получим
S ''i ( x) = 2ci .

Дифференцируя
сплайн
Следовательно, используя интерполяционный многочлен Лагранжа
по узлам x1 , x2 , x3 , можно взять
с1 =

f1 ( x3 − x2 ) + f 2 ( x3 − x1 ) + f3 ( x2 − x1 ) .
( x3 − x1 )( x3 − x2 )( x2 − x1 )

Пусть в области Ω из задания 1 задан ряд точек ( xi , yi ), 1 ≤ i ≤ n ,
=
xmin min(
=
xi ), xmax max(
=
xi ), ymin min(
=
yi ), ymax max( yi )
i

i

i

i

1 ≤ i ≤ n , примыкающих к границе области Ω . Введем в рассмотрение функцию y = g ( x) с помощью сплайна, проходящего через
заданные точки.
Обозначим через R1 множество точек:

=
R1

{( x, y ) : ( xmin ≤ x ≤ xmax ) ∧ ( y ≤ g ( x))} .

Включим построенную границу в область Ω1 , полагая Ω1 =Ω \ R1 ,
так что Ω1 ⊆ Ω .
113

Цель работы состоит в построении линий уровня функции
H( x, y ) в нерегулярной многосвязной области Ω1 .
Данные по включаемому криволинейному участку принять в
соответствии с номером варианта.
1)

2)
16

16

14

14

12

12

10

10

8

8

6

6

4

4

2
2

10

8

6

4

12

14

2
2

16

3)

4

6

8

10

12

14

16

4

6

8

10

12

14

16

4)
16
14
12
10
8
6
4
2
2

4

6

8

10

12

14

16

5)

6)
16

16

14

14

12

12

10

10

8

8

6

6

4

4

2
2

4

6

8

10

12

14

2
2

16

114

7)

8)
16

16

14

14

12

12

10

10

8

8

6

6

4

4

2
2

4

6

8

10

12

14

2
2

16

9)

4

6

8

10

12

14

16

4

6

8

10

12

14

16

4

6

8

10

12

14

16

10)
16

16

14

14

12

12

10

10

8

8

6

6

4

4

2
2

4

6

8

10

12

14

2
2

16

11)

12)
16

16

14

14

12

12

10

10

8

8

6

6

4

4

2
2

4

6

8

10

12

14

2
2

16

115

13)

14)
16

16

14

14

12

12

10

10

8

8

6

6

4

4

2
2

4

6

8

10

12

14

2
2

16

15)

4

6

8

10

12

14

16

4

6

8

10

12

14

16

4

6

8

10

12

14

16

16)
16

16

14

14

12

12

10

10

8

8

6

6

4

4

2
2

4

6

8

10

12

14

2
2

16

17)

18)
16

16

14

14

12

12

10

10

8

8

6

6

4

4

2
2

4

6

8

10

12

14

2
2

16

Пример. Пусть область Ω соответствует иллюстрации программы, приведенной в задании 1.
116

Учтем, что
xi

5

6

7

7.8

10

12

16

yi

2

3.2

3.9

4

5

8

8

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

30 0

25 0

35
0

25 0

14 20 0

40
0
30
0

35
0

20 0

12 15 0

45
0

25
0
10

15
0

6

0
30

8

20
0

10
0
50

4

2
2

4

6

10

8

12

14

16

Рис. 47

Листинг программы
Вычисление интеполирующей функции выполнено в виде функции spl(x, y, tt) .
function X = spl(xy,tt,j)
[n,m]=size(xy);
[k,h]=size(tt);
t=zeros(m);
t(1)=0.;
for i=2:m
117

dt=(xy(1,i)-xy(1,i-1))^2+(xy(2,i)-xy(2,i-1))^2;
t(i)=t(i-1)+sqrt(dt);
end
dt=1./t(m);
for i=1:m
t(i)=t(i)*dt;
end
kp=m-1
c=zeros(2,kp);
for jp=1:2
for i=1:kp
if(i==1)
c(jp,1)=(xy(jp,1)*(t(3)-t(2))-xy(jp,2)*(t(3)-t(1))+
xy(jp,3)*(t(2)-t(1)))/((t(2)-t(1))*(t(3)-t(2))*(t(3)-t(1)));
else
c(jp,i)=-c(jp,i-1)*(t(i)-t(i-1))/(t(i+1)-t(i))+
(-(xy(jp,i)-xy(jp,i-1))/(t(i)-t(i-1))+(xy(jp,i+1)
-xy(jp,i))/(t(i+1)-t(i)))/(t(i+1)-t(i));
end
end
end
X=zeros(size(tt));
for s=1:h
for i=1:m-1
if((tt(s)-t(i))*(t(i+1)-tt(s))>0)
X(s)=xy(j,i)+(xy(j,i+1)
-xy(j,i))*(tt(s)-t(i))/(t(i+1)-t(i))+c(j,i)*(tt(s)-t(i))*(tt(s)-t(i+1));
end
end
end
return
end
Основной модуль программы
n1=200;
n2=200;
H=repmat(0,n1,n2);
XMIN=2.;
XMAX=16.;
YMIN=2.;
YMAX=16.;
DX=(XMAX-XMIN)/n1;
118

DY=(YMAX-YMIN)/n2;
n=7;
xy = zeros (2,n);
xy(1,1)= 5.00; xy(2,1)=2.000;
xy(1,2)= 6.00; xy(2,2)= 3.2 ;
xy(1,3)= 7.00; xy(2,3)= 3.90;
xy(1,4)= 7.80; xy(2,4)= 4.00;
xy(1,5)= 10.0; xy(2,5)= 5.00;
xy(1,6)= 12.0; xy(2,6)= 8.00;
xy(1,7)= 16.0; xy(2,7)= 8.00;
x4=xy(1,:);
y4=xy(2,:);
tt=linspace(0.001,0.9995,25);
xx=spl(xy,tt,1);
yy=spl(xy,tt,2);
hold on
plot(x4,y4,'bo');
XI=zeros(n1,n2);
YI=zeros(n1,n2);
for I=1:n1
for J=1:n2
YI(I,J)=YMIN+J*DY;
XI(I,J)=XMIN+I*DX;
H(I,J)=HIJ(XI(I,J),YI(I,J));
if((XI(I,J)-5)^2+(YI(I,J)-10)^2<4)
H(I,J)=NaN;
end
end
end
for I=1:n1
for J=1:n2
YI(I,J)=YMIN+J*DY;
XI(I,J)=XMIN+I*DX;
for K=2:25
if((xx(K)-XI(I,J))*(XI(I,J)-xx(K-1))>0)
if(YI(I,J)<yy(K))
H(I,J)=NaN;
end
end
end
119

end
end
figure(1);
t = 0:0.1:2*pi;
x=5+2*cos(t);
y=10+2*sin(t);
plot(x,y,'r');
hold on
grid on
hold on
[c,h] = contour(XI,YI,H);
clabel(c,h);
hold on
plot(xx,yy,'r');
Задание 3. Восстановление аналитической функции по её
вещественной или мнимой части, заданной на границе области.
Эта задача сводится к решению задачи Дирихле для уравнения
Лапласа, которая подробно рассматривалась в подразд. 4.5:
∂ 2u ∂ 2u
=
+
0,
∂x 2 ∂y 2
=
u Re(ϕ( z )),

( x, y ) ∈ Ω
( x, y ) ∈ ∂Ω

где ϕ( z ) ‒ аналитическая в области Ω функция.
На эскизах принять одну клеточку за единицу. В качестве ϕ( z )
принять функции z2 и z3 для четного и нечетных вариантов соответственно.
Для решения задачи адаптировать программу, листинг которой
приведен в подразд. 4.5. Требуемая модификация заключается в
изменении координат узлов элементов, номеров узлов элементов,
граничных узлов и граничных условий.
В отчете представить интегральную формулировку задачи, варианты с исходной и уточненной сеткой. Для исходной сетки использовать треугольники с исходными данными, уточненную сетку построить с помощью программы triangle. Привести алгебраическую
систему до и после учета граничных условий. Данные представить в
табличном виде: номер узла, точное и конечноэлементное решение,
абсолютная погрешность решения.
120

1)

2)

3)

4)

5)

6)

7)

8)

9)

10)

11)

12)

13)

14)

15)

16)

17)

121

Библиографический список
1. Даутов Р.З., Корчевский М.М. Введение в теорию метода конечных элементов. Казань, 2004.
2. Деклу Ж. Метод конечных элементов. М.: Мир, 1976.
3. Джордж А., Лю Дж. Численное решение больших разреженных систем
уравнений. М.: Мир, 1984.
4. Зенкевич О.К. Метод конечных элементов в технике. М.: Мир, 1989.
5. Митчел Э., Уэйт Р. Метод конечных элементов для уравнений с частными производными. М.: Мир, 1981.
6. Михлин С.Г. Вариационные методы в математической физике. М.:
Гостехиздат, 1957.
7. Стренг Г., Фикс Дж. Теория метода конечных элементов. М.: Мир, 1977.
8. Jonathan Richard Shewchuk, Delaunay Refinement Algorithms for Triangular
Mesh Generation, Computational Geometry: Theory and Applications 22(1-3):21-74,
May 2002.
9. Шайдуров В.В. Многосеточные методы конечных элементов. М.: Наука,
1989.

122

ОГЛАВЛЕНИЕ
ВВЕДЕНИЕ ................................................................................................................... 3
1. КЛАССИЧЕСКИЕ МЕТОДЫ ВЗВЕШЕННЫХ НЕВЯЗОК .................................. 4
1.1. Вводные замечания ........................................................................................... 4
1.2. Классические и «слабые» формулировки краевых задач для линейных
дифференциальных уравнений .......................................................................... 7
1.3. Вариационная формулировка краевых задач для линейных дифференциальных уравнений ......................................................................................... 11
1.4. Метод Галёркина ............................................................................................. 13
1.5. Метод Ритца..................................................................................................... 20
2. БАЗИСНЫЕ ФУНКЦИИ ДЛЯ ЭЛЕМЕНТОВ ...................................................... 26
2.1. Одномерный элемент ...................................................................................... 26
2.2. Треугольный элемент ...................................................................................... 29
2.3. Четырехугольный элемент ............................................................................. 32
3. МЕТОД КОНЕЧНЫХ ЭЛЕМЕНТОВ ДЛЯ РЕШЕНИЯ ОБЫКНОВЕННЫХ
ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ ............................................................ 35
3.1. Классическая и интегральная формулировка задачи ................................... 35
3.2. Дискретизация области, выбор базисных функций...................................... 36
3.3. Формирование дискретной конечноэлементной схемы ............................... 40
3.4. Примеры решения задач ................................................................................. 48
3.5. Сходимость метода ......................................................................................... 60
4. МЕТОД КОНЕЧНЫХ ЭЛЕМЕНТОВ ДЛЯ РЕШЕНИЯ УРАВНЕНИЙ С
ЧАСТНЫМИ ПРОИЗВОДНЫМИ ........................................................................ 64
4.1 Классическая и интегральная формулировка задачи .................................... 64
4.2. Дискретизация области, выбор базисных функций...................................... 71
4.3. Конечно-элементная формулировка задачи. Матрицы элементов............. 81
4.4. Триангуляция области .................................................................................... 90
4.5. Примеры решения задач ................................................................................. 93
4.6. Задания для практических занятий .............................................................. 107
Библиографический список ....................................................................................... 122

Солдаткин Андрей Владимирович, Баранова Елена Семёновна
Введение в метод конечных элементов
Редактор Г.М. Звягина
Корректор Л.А. Петрова
Компьютерная верстка: А.В. Мещерякова
Подписано в печать 07.09.2020. Формат 60х84/16. Бумага документная.
Печать трафаретная. Усл.-п. л. 7,15. Тираж 200 экз. Заказ № 123.
Балтийский государственный технический университет
Типография БГТУ
190005, С.-Петербург, 1-я Красноармейская ул., д.1