SUBROUTINE BANSOL{R,X,B,N,INT)

USES L.U FACTORISATION TO SOLVE FOR X, GIVEN R

DIMENSION R(65),X(65),B(5,65) IFdNT .EQ. 2)G0T0 3

INT = 1, TRIDIAGONAL SYSTEM

NP N - 1 DO 1 J = 1,NP JP J + 1

1 R(JP) R(JP) - B(2,JP)*R(J) X(N) = R(N)/B(3,N)

DO 2 J = 1,NP JA = N - J

X(JA) = (R{JA) - B(4,JA)*X(JA+1))/B(3,JA)

2 CONTINUE RETURN

INT 2, PENTADIAGONAL SYSTEM

ASSUMES FIRST EQUATION FORMED AT MIDSIDE NODE

SET IBC 0 IF LAST EQUATION FORKED AT MIDSIDE NODE

SET IBC = 1 IF LAST EQUATION FORMED AT CORNER NODE

3 IBC 1

NH = N/2

IF(2*NH .EQ. N)R(N+1) 0.

IF{2*NH .EQ. N)B(2,N+1) 0.

DO 6 I =1,2

DO 5 J = 1,NH

JA = 2*J

DO 4 К = 1,1

JB = JA - 1 + К

4 R(JB) * R(JB) - B(I,0B)*R(JB-1)

5 CONTINUE

6 CONTINUE

NEN = NH - IBC X(N) = R(N)/B(3,N)

IFdBC .EQ. l)X(N-l) = {R(N-l) - В (4,N-1) *X{N))/B(3.N-

DO 7 J 1,NEN

JA N - 2*J + 1 - IBC

X(JA) (R{JA) - B(4,JA)*X(JA+1) - В{5,JA)*X(JA+2))/В( X(JA-l) = (R(JA-l) - B(4,JA-1)*X(JA))/B(3,JA-1)

7 CONTINUE RETURN END

Рис. 6.19. Распечатка подпрограммы BANSOL.

нальных/пятидиагональных матриц является решение задачи Штурма - Лиувилля (п. 5.4.2).

Полезная проверка правильности составления подпрограмм решения уравнений (6.23) состоит в том, чтобы после определе-

ния V заполнить левые части уравнений (6.23) и убедиться в равенстве левых частей правым. Это, однако, не гарантирует точности построенного решения V, если матрица А плохо обусловленная, иначе говоря, если уравнения, составляющие систему (6.23), близки к линейным. Плохо обусловленные системы и их влияние на точность решения обсуждаются в книгах [Gerald, 1987; Dahlquist, Bjorck, 1974].

6.2.4. Обобщенный алгоритм Томаса

Применение конечно-разностных или конечно-элементных схем повышенного порядка приводит к матрицам с более широкими лентами, чем трехдиагональная матрица. Однако нетрудно провести обобщение алгоритма Томаса. Допустим, что нам необходимо решить следующую пятидиагональную систему:

bi Ci h

(h b2 C2 /2

3 3 bz /3

et at bi

ejsji ai &ЛГ-1 yv-i % bf

Первый этап состоит в исключении всех элементов ei, что дает

bi с, /,

2 *2 С2 h

3 К

3 С,

а; й; < п

N

Следующий этап по существу совпадает с первым этапом обычного алгоритма Томаса. Иначе говоря, все элементы а\ исключаются, что дает

1 </.

1 </г

JV-1

(6.33)

(6.34)

После этого для решения уравнений (6.33) требуется обратная подстановка следующего общего вида:

(6.35)

Очевидно, что для решения пятидиагональной системы потребуются две прямые и одна обратная прогонки.

Различные этапы обобщенного алгоритма Томаса можно интерпретировать как серии операций (прямых прогонок), требуемых для приведения матрицы А к верхнетреугольному виду.