области: °/./-2 -бд; is

4 2Дд:,

-(д17+А3) + 1(Д/ + Д/) (5.84)

+ ЗДх,.+, + 15

-1 .

*/./+2- бДл:., 15-

/ Дл:.- \ / 2 Дл:, \ 4

Особая форма выражений (5.84) и (5.85) является следствием ограничений = 0.5 (л:/ 2 + Xj) и = 0.5 (х/ + jc/+2). Эти ограничения могут быть ослаблены, что, однако, приведет к усложнению алгебраических выражений.

В граничном узле m = J формулы (5.77) и (5.81) принимают вид

-1 Ал:,

Подставляя (5.74) и (5.69) в уравнение (5.72) и проводя не-терполяции, но центрированных в узлах в середине сторон, формулы (5.77) и (5.81) принимают вид 4 2 AjCy

16 Да:,

/./ = -1д17 + -1. (5-88)

4 2 AjCy

/ / + 1 = ЗА] + ~15

/2Ах,\ /16Ajc/\ Г2АХ:\ ,

Система уравнений (5.76) является трехдиагональной, если применяется линейная интерполяция, а в случае применения квадратичной интерполяции - попеременно трехдиагональной и пятидиагональной. Эффективные модификации метода исключения по Гауссу могут быть использованы для решения уравнений (5.76) в том или ином из указанных случаев и описываются в пп. 6.2.2 и 6.2.3.

Программа STURM применяет метод Галёркина с конечными элементами к уравнению Штурма - Лиувилля (5.66). Распечатка программы STURM приводится на рис. 5.19. Основные параметры, используемые в программе STURM, описываются в табл. 5.9.

Таблица 5.9. Параметры, используемые в программе STURM

Программа STURM считывает (и выводит на печать) исходные управляющие параметры (строки 11-21). Задаются значения х, точное решение, а также узловые значения функции F в правой части уравнения (5.66) (строки 28-42). Вычисляются вклады в матрицу В и вектор G в уравнениях (5.76) (строки 43-69) и проводится корректировка с целью учета граничных условий. Система уравнений (5.76) с ленточной матрицей решается путем обращения к подпрограмме BANFAC (строка 83), выполняющей факторизацию В в произведение L-U, и обращения к подпрограмме BANSOL (строка 86) для решения факторизованной системы. Распечатка и описание подпрограмм BANFAC и BANSOL дается в п. 6.2.3. Далее рассчи-

12 К. Флетчер, т. 1

178 Гл. 5. Методы взвешенных невязок 1

3 С STURM COMPUTES THE SOLUTION OF THE STURM-LIOUVILLE PROBLEM,

4C P Q=1,F -SIG(AL*SIN(L-0.5)*PI*X), USING THE FIHITE ELEMENT

5 с METHOD ON A UNIFORM GRID WITH EITHER LINEAR OR QUADRATIC

! 6 с INTERPOLATION.

7 с THE BANDED SYSTEM IS SOLVED WITH BANFAC, BANSOL (SECTION 6.2).

9 DIMENSION Y(65),YEX(65),X(65),B(5,65),G(65),F(65),A(5),YD{65)

10 с

11 OPEN(1.FILE= STURM.DAT)

12 OPEN(6, FILE- STURM.OUT *)

13 READ(1,1)INT,IPR,NX,A

14 1 FORMAT(3I5,5E10.3)

15 с

16 IFdNT .EQ. 1)VRITE(6,2)

17 IFdNT .EQ. 2)WRITE(6,3)

18 VRITE(6,4)NX,A

19 2 FORMATC STURM-LIOUVILLE PROBLEM, FEMtLINEAR INTERPOLATION)

20 3 FORMATC STURM-LIOUVILLE PROBLEM, FEMtQUADRATIC ХНТГГРОЬАТЮК)

21 4 FORMATC NX-M3,* A \5E10.3, )

22 С

23 PI 3.1415927

24 NXP = NX - 1

25 ANX NXP

26 DX l./ANX

27 DXS = DX*DX

28 С

29 с SET GRID, EXACT SOLUTION AND RHS(F)

30 с

31 DO 5 I 1,NX

32 AI I - 1

33 X(I) AI*DX

34 F(I) = 0.

35 YEX(I) = 0.

36 DO 5 L 1,5

37 AL = L

38 . AL (AL - 0.5)*PI

39 DUM = A(L)*SIN(AL*X(I))

40 F(I) F(I) - DUM

41 YEX(I) YEX(I) - DUM/(1. - AL*AL) 42. 5 CONTINUE

43.с

44 с SET UP ARRAYS, В AND G

45 с

46 DO 7 I 2,NX,INT

47 IM I - 1

48 IFdNT .EQ. 2)G0T0 6

49 B(1,IM) = 0.

50 B(2,IM) l./DXS 1./6.

51 B(3,IM) = -2./DXS 2./3.

52 B(4,IM) = B(2,IM)

53 B(5,IM) = 0.

54 G(IM) = (Fd-l) + 4.*F(I) + Fa-H))/6.

55 GOTO 7

Рис. 5.19. Распечатка программы STURM (начало).