значения которого совпадают со значениями функции в этих точках, и точно интегрируется (фиг. 8.9).
Так как п значений функции определяют полином степени п- 1, ошибка имеет порядок 0(Д)", где Д - расстояние между точками. В результате получаем известные квадратурные формулы Ньютона - Котеса, в соответствии с которыми интеграл можно записать в виде
~ (8.24)
при интегрировании в пределах от -1 до -fl (фиг. 8.9, а). Например, если п = 2, то
/ = f(-l) + f(l)- (8.25)
При /1 = 3 получается известная формула трапеций
/ = [f(-l) + 4f(0)-f f(l)], (8.26)
при /1 = 4 - формула одной трети Симпсона
/ = [/(-l) + 3f(-)-f 3f(l)-ff (1)]. (8.27)
Формулы интегрирования для различных п (до /г = 20 включительно) приведены у Копаля [8].
Квадратура Гаусса. Если значения функции вычисляются не в априорно заданных точках, а так, чтобы достигалась наилучшая возможная для данного количества точек точность, интегрирование выполняется точнее. Если положить, что
I п.
1= \ f(s)rf = EHy.
(8.28)
и снова представить подынтегральную функцию в виде полинома, то нетрудно увидеть, что для п точек интегрирования получим 2п неизвестных (/,• и г), и, следовательно, можно построить и точно проинтегрировать полино.м степени 2п - 1 (фиг. 8.9,6). В результате ошибка будет иметь порядок 0(Д)2п.
Решить полученную систему уравнений трудно, но с помощью некоторых математических приемов [7] решение можно получить в полиномах Лежандра. Поэтому этот метод часто называется квадратурой Гаусса - Лежандра.
В табл. 8.1 приведены координаты точек и весовые коэффициенты для интегрирования по Гауссу.
Таблица 8.1
Абсциссы в весовые коэффициенты квадратурной формулы Гаусса
f{x)dx= 2 Hjf(aj)
±а
1 = 2
0,57735
02 691
89 626
„ О
1,00000
00 000
00 000
0,77459
66692
41 483
п = о
0,55555
55 555
55 556
0,00000
00 000
00 000
„ а
0,88888
88 888
88 889
0,86113
63115
94 053
п -
0,34785
48 451
37 454
0,33998
10 435
84 856
0,65214
51 548
62 546
п = 5
0,90617
98 459
38 664
0,23692
68 850
56 189
0,53846
93101
05 683
0,47862
86 704
99 366
0,00000
00 ООО
00 000
0,56888
88 888
88 889
п = 6
0,93246
95 142
03 152
0,17132
44 923
79 170
0,66120
93 864
66265
0,36076
15 730
48 139
0,23861
91 860
83197
0,46791
39 345
72 691
п = 7
0,94910
79 123
42 759
0,12948
49 661
68870
0,74153
II 855
99 394
0,27970
53914
89 277
0,40584
51 513
77 397
0,38183
00 505
05 119
0,00000
00 ООО
00 ООО
п = 8
0,41795
91 836
73 469
0,96028
98 564
97 536
0.10122
85 362
90 376
0,79666
64 774
13 627
0,22238
10 344
53 374
0,52553
24 099
16 329
0,31370
66 458
77 887
0,18343
46 424
95 650
0,36268
37 833
78 362
« = 9
0,96816
02 395
07 626
0,08127
43 883
61 574
0,83603
И 073
26 636
0,18064
8! 606
94 857
0,61337
14 327
00 590
0,26061
06 964
02 935
0,32425
34234
03 809
0,31234
70 770
40 003
0,00000
00 000
00 000
0,33023
93 550
01 260
п= 10
0,97390
65 285
17 172
0,06667
13 443
08 688
0,86506
33 666
88 985
0,14945
13 491
50 581
0,67940
95 682
99 024
0,21908
М 625
15 982
0,43339
53 941
29247
0,26926
67 193
09 996
0,14887
43 389
81 631
0,29552
42 247
14753
6 Зак. 613
В методе конечных элементов сложность вычислений заключается в определении значений интегрируемой функции f, поэтому в дальнейшем будет применяться только формула Гаусса, использующая минимальное число значений функции.
Можно получить формулы для интегрирования с заданной степенью точности выражений вида
(8.29)
при известной функции w{l), если опять заменить f() полиномом [7].
8.9. Прямоугольник или прямая призма
Проще всего взять интеграл
/=5 5 f{l,r))dldri.
(8.30)
вычисляя сначала значение внутреннего интеграла в предположении, что Т1 - постоянная, т. е. используя формулу
5 f (I, Ti) d = Hjf {l„ r\)= (Ti). (8.31)
-I /=1
Вычисляя затем внешний интеграл, получаем
I п п п
7 = 5 Ф (Т1) dTi = Ягф (Т1,) = X Я, Яf (, ,1,) =
1=1 1=
1-1 i-i
Аналогично для прямой призмы имеем i I I
/= S 5 S fil, г\, l)dldr]d-. -1 -1 -I
п п п
(8.32)
(8.33)
Здесь предполагалось, что число точек иитегрировання в каждом направлении одинаково.
Интересно отметить, что двойное суммирование легко заменить простым по (п X п) точкам для прямоугольника (или по п? точкам для куба). На фиг. 8.10 показано девять точек, необходимых для точного интегрирования полинома пятого порядка по каждой из переменных.
Однако к решению этой задачи можно было бы подойти с другой стороны и интегрировать точно полином пятого порядка по двум переменным. В этом случае в каждой точке пришлось бы определить две координаты и значение подынтегральной функции f, входящей в весовую формулу типа
I I т
/=5 5 f(g, Ti)drfTi = 2eD,f(,,7i,).
(8.34)
Фиг. 8.10. Точки интегрирования для л-= 3 в квадратной области. (Точно интегрируется по-.тином пятой степени по каждой переменной.)
Оказывается, что при этом для достижения той же точности достаточно было бы использовать только семь точек. Формула такого типа для трехмерного «кирпичика» получена и с успехом использована Айронсом [9].
8.10. Треугольник или тетраэдр
Для треугольника интегралы по L-координатам имеют вид
/ = 5 5 /(L„ Ц, L,)dL,dL,.
(8.35)
m=I /-I (=-1
Оиять можно было бы использовать п гауссовых точек и получить сумму типа, рассмотренного в предыдущем разделе. Однако теперь пределы интегрирования сами содержат переменные, поэтому для второго интегрирования удобно использовать выражения вида (8.29) квадратуры Гаусса, где w - лнией-ная функция. Это было сделано Радо [10, 11]. В табл. 8.2 приведены весовые коэффициенты, входящие в выражение
I=tiwf{Li, L„ Ц),
Таблица 8.2
Коэффициенты для интегрирования по формулам Гаусса и Радо
Количество
точек интегрирования в кажлом направлении
J = h п
AI \1] /=1, п
.4SI71 1 = 1, п
п = 1
0,3333333333
0,75
(1.0)
(0,25)
П = 2
0,2113248654
0,1550510257
0,3764030627
0,7886751346
0,6449489743
0,5124858262
(1,0)
(0,1111111111)
п = 3
0,1(27016654
0,2777777778
9,0885879595
0,2204622112
0,4444444444
0,4094668644
0,3881934688
0,8872983346
0,2777777778
0,7876594618
0,3288443200
(1,0)
(0,0.625)
И = 4
0,0694318442
0,1739274226
0,0571041961
0,1437135608
0,3300094782
0,3260725774
0,2768430136
0,2813560151
0,6699905218
0,3260725774
0,5835904324
0,3118265230
0.9305681558
0,1739274226
0,8602401357
0,223103901!
(1,0)
(0,04)
и = 5
0,0469100770
0,1184634425
0,0398098571
0,1007941926
0,2307653449
0,2393143353
0,1980134179
0,2084506672
0,2844444444
0,4379748102
0,2604633916
0,7692346551
0,2393143353
0,6954642734
0,2426935942
0,9530899230
0,1184634425
0,9014649142
0,1598203766
(1,0)
(0,277777778)
L2 = AJ{j){\-L,),
W = AS(i)H(f)(\~L,).
(8,36)
Аналогичные соотношения можно получить и для тетраэдра.
На фиг. 8.11 показано расположение точек интегрирования в треугольниках при п = 1, 2, 3. Видно, что они расположены неравномерно и несимметрично. Кроме того, в направлениях и, 1г, Ьз точность интегрирования различна. Другое (более изящное) расположение точек, предложенное Хаммером и др. [12], позволяет существенно упростить расчет; весовые коэффи-
Фиг. 8.11. Точки интегрирования для треугольника при использовании метода Гаусса - Радо.
циенты для выражений, аналогичных (8.34), приведены в табл. 8.3 [13].
Можно убедиться, что точек всегда столько или немного больше, чем требуется для получения полных полиномов заданного порядка.
Очевидно, что эти результаты можно обобщить и на тетраэдры. В табл. 8.4 представлены некоторые формулы из работы [12].
8.11. Заключительные замечания
В этой главе показано, каким образом можно построить большое количество криволинейных элементов. Необходимость использования методов численного интегрирования потребовала описания некоторых из них. Дальнейшие подробности можно найти в различных учебниках по численному анализу.
Очевидно, что численное интегрирование является приближенным. Вопрос о том, какая степень точности необходима в практических задачах, будет рассмотрен в следующей главе. Там же будут приведены основные принципы организации программ при использовании численного интегрирования.
