Mathematical Physics - Volume II - Numerical Methods

3.6 Program for solving of elliptical problems

103

&

nband COMMON /cnode / x(2,maxn) , u(maxn) COMMON /cmatrx/ gk(maxn,maxband) , gf(maxn) PRINT * , ’solving equations ...’ CALL trib(gk,nnode,nband,maxn,maxband) CALL rhsb(gk,u,gf,nnode,nband,maxn,maxband) RETURN END

SUBROUTINE trib(a,n,ib,l1,l2) DIMENSION a(l1,l2) c this subroutine triangularizes a banded and symmetric matrix a c only the upper half-band of matrix is stored c storage is in the form of a rectangular array l1 x l2

c the half-band width is ib c the number of equations is n DO 100 i = 2 , n m1 = min0(ib-1,n-i+1) DO 50 j = 1 , m1 sum = 0.0 k1 = min0(i-1,ib-j) DO 20 k = 1 , k1

sum = sum + a(i-k,k+1)*a(i-k,j+k)/a(i-k,1)

20

CONTINUE a(i,j) = a(i,j) - sum

50

CONTINUE

100 CONTINUE RETURN END

SUBROUTINE rhsb(a,x,b,n,ib,l1,l2) DIMENSION a(l1,l2) , x(1) , b(1) c for the linear sytem a*x=b with the matrix a c triangularized by routine trib c this routine performs the forward substitution c into b and back substitution into x

c the half-band width of a is ib c the number of equations is n np1 = n + 1 DO 100 i = 2 , n sum = 0.0 k1 = min0(ib-1,i-1) DO 50 k = 1 , k1

sum = sum + a(i-k,k+1)/a(i-k,1)*b(i-k)

50

CONTINUE b(i) = b(i) - sum

100 CONTINUE c begin back-substitution x(n) = b(n)/a(n,1)