ORMAT
Notes on BAND and LISOLV
J. W. Evans and P. Kar
Dept. of Materials Science and Engineering
University of California, Berkeley, CA 94720
These are notes on software frequently used in the past for numerical solution of partial differential equations (pdes). BAND can be found in Appendix C of “Electrochemical Systems” by Prof. J. S. Newman, Chem. Eng. Dept., UC, Berkeley (Prentice-Hall, 1973 and later editions). The algorithm is a convenient one for solving sets of (more than one) pdes in one spatial dimension. Generally the equations would be coupled, as for example in some electrochemical problems where equations for potential and concentration must be solved simultaneously. Newman’s Appendix C contains a FORTRAN code for BAND and an associated matrix inversion subroutine MATINV.
LISOLV is an algorithm in the software TEACH that was developed at Imperial College, London in the 1970s. It was intended for numerically solving pdes in two (or more) spatial dimensions by solving along a line traversing the solution domain and then sweeping that line across the domain. LISOLV treats one equation at a time so that it must be applied iteratively if multiple coupled equations are to be solved. [This is not a concern in TEACH as TEACH is intended for solution of flow/transport in turbulent fluids where the inherent non-linearity requires an iterative solution anyway.]
These notes provide descriptions of the algorithms (to complement Newman in the case of BAND) and MATLAB versions of the software, together with examples of the testing of the MATLAB versions.
BAND
Newman starts with the set of coupled linear differential equations in the spatial dimension x:
1)
There are n equations of the form and n unknowns, c1, c2,……ck,……cn.
The summation sign in applies to the whole equation. k is the index indicating the dependent variable; i is the equation number. The summation arises because, for example, in the equation where c5 is the dependent variable the values of c1,c2……and their derivatives may appear. The coefficients ai,k etc. are permitted to vary from one equation to the next and to be functions of position (x).
Newman then provides the finite difference approximations to the derivatives appearing in :
2)
and
3)
These are the approximations to the derivatives at the point xj in the mesh used for the numerical solution. The mesh is uniformly spaced with a separation h between mesh points. These approximations are substituted in and the set of equations multiplied by h2 to yield the set of difference equations:
4)
evaluated at position xj-1 = xj – h. Bi,k(j) is the coefficient, in equation i, at the position xj, of the dependent variable ck evaluated at position xj.
The coefficients in are given by:
5)
6)
7)
and 8)
[A nomenclature slightly different from Newman’s has been used in eqn. , but not in most subsequent equations.]
The set of equations applies to all interior mesh points in the domain but other equations have to be invoked for the boundary points that reflect the boundary conditions. Newman provides for very general boundary conditions in the form:
9)
This equation would be applied (with different coefficients, perhaps) at both boundaries, e.g. at x = 0 and x = L if the domain runs from 0 to L. It allows the values of the dependent variables, or their derivatives, or the link between the variables and derivatives, to be specified at each of the two boundaries. Of course equation applies right up to the boundary of the domain; however it cannot be simply be used at the boundary point because it does not incorporate the boundary condition and, furthermore, in its finite difference form, equation, contains a point which would be outside the domain. Newman uses a method known as the fictitious point method (image point in his terminology) for overcoming these difficulties. This entails incorporating one additional fictitious point just outside each end of the domain. A finite difference equation, based on the boundary condition (BC) is then written for each fictitious point and equation , invoked for the boundary point itself. [For example, if the BC is that the gradient is zero then the value of the dependent variable at the fictitious point is set equal to the value one in from the boundary.] In setting up the mesh for the calculations using m mesh points, point 1 is the fictitious point at one end of the domain and point m that at the other (see Fig. C-1 in Newman). Equation is then used for points 2 through m-1. For point 1 the fictitious point method gives:
10)
with
11)
12)
and
13)
These coefficients are supplied by the user to BAND so that BAND solves equation in conjunction with solving equation . If no derivatives appear in the boundary condition then it is not necessary to invoke the fictitious point method and the mesh can start at the boundary. In that case we still use equation but put
14)
and
15)
with Gi(1) as before.
We can use similar equations to employ the fictitious point method at the other boundary:
16)
Newman does not give the coefficients in this equation. They are:
17)
18)
and
19)
Again, if no derivatives appear in the BCs then the mesh can end at the boundary and then:
20)
and
21)
with Gi(m) given by .
This completes the formulation of the finite difference equations and Newman describes the development of the algorithm for their solution, an extension of Thomas’ algorithm on pages 417 and 418. No elaboration of that development is needed here (except to note an apparent typographical error in the second equation (C-21) where ai,l should probably be Ai,l). Newman then provides the FORTRAN program ‘BAND’ implementing the algorithm, along with a FORTRAN subroutine ‘MATINV’ that is called by BAND. There are two things to note here.
- The variable name ‘SAVE’ is used (thrice) in MATINV and was an innocuous variable name at the time the book was written. However, recent versions of FORTRAN make SAVE an executable statement and modern FORTRAN compilers will balk at this variable name. Change it to something else (SAYVE, for example).
- When the book was written computer memory was at a premium and BAND was written to save memory by supplying the coefficients Ai,k etc. to BAND every time BAND is called. Because BAND is called for each mesh point in turn then this can lead to inefficiencies. This is true for the solution of equations where the coefficients are position dependent; they must then be assigned on each call to BAND. Future versions of BAND might store all necessary coefficients in arrays and this would obviate the need to assign them on each call to BAND.
MATLAB version of BAND
The following is a version of BAND (and MATINV) obtained by “translating” the FORTRAN version in Newman.
function [y]=band(J)
global A B C D E G X Y NPOINTS N NJ NP1 DETERMINANT
% ***********************
% Case where J negative
% ***********************
if (J-2) < 0
NP1=N+1;
for I=1:N
D(I,2*N+1)=G(I);
for L=1:N
LPN=L+N;
D(I,LPN)=X(I,L);
end
end
% ********************************************
% Calling Matinv
% ********************************************
matinv(N,2*N+1);
if (DETERMINANT == 0)
sprintf('DETERMINANT = 0 at J= %6.3f',J)
break
else
for K=1:N
E(K,NP1,1)=D(K,2*N+1);
for L=1:N
E(K,L,1) = -D(K,L);
LPN = L+N;
X(K,L) = -D(K,LPN);
end
end
return
end
end
% ***************************
% Case where J-2 equal to zero
% ***************************
if((J-2)==0)
for I=1:N
for K=1:N
for L=1:N
D(I,K)=D(I,K)+A(I,L)*X(L,K);
end
end
end
end
% ***********************************
% Case where J-2 is greater than zero
% **********************************
if (J-NJ) < 0
for I=1:N
D(I,NP1) = -G(I);
for L=1:N
D(I,NP1)=D(I,NP1)+A(I,L)*E(L,NP1,J-1);
for K=1:N
B(I,K)=B(I,K)+A(I,L)*E(L,K,J-1);
end
end
end
else
for I=1:N
for L=1:N
G(I)=G(I)-Y(I,L)*E(L,NP1,J-2);
for M=1:N
A(I,L)=A(I,L)+Y(I,M)*E(M,L,J-2);
end
end
end
for I=1:N
D(I,NP1) = -G(I);
for L=1:N
D(I,NP1)=D(I,NP1)+A(I,L)*E(L,NP1,J-1);
for K=1:N
B(I,K)=B(I,K)+A(I,L)*E(L,K,J-1);
end
end
end
end
% ********************************************
% Calling Matinv
% ********************************************
matinv(N,NP1);
if DETERMINANT ~= 0
for K=1:N
for M=1:NP1
E(K,M,J) = -D(K,M);
end
end
else
sprintf('DETERMINANT = 0 at J= %6.3f',J)
break
end
if (J-NJ)<0
return
else
for K = 1:N
C(K,J) = E(K,NP1,J);
end
for JJ = 2:NJ
M = NJ - JJ + 1;
for K=1:N
C(K,M) = E(K,NP1,M);
for L=1:N
C(K,M)=C(K,M)+E(K,L,M)*C(L,M+1);
end
end
end
for L=1:N
for K=1:N
C(K,1) = C(K,1) + X(K,L)*C(L,3);
end
end
end
%*********************************************************************
function [DETERMINANT] = matinv(N,M)
global A B C D E G X Y NPOINTS NJ DETERMINANT
DETERMINANT=1.0;
for I=1:N
ID(I)=0;
end
for NN=1:N
BMAX = 1.1;
for I=1:N
if(ID(I)~=0)
continue
else
BNEXT = 0;
BTRY = 0;
for J=1:N
if(ID(J)==0)
if(abs(B(I,J)) > BNEXT)
BNEXT=abs(B(I,J));
if(BNEXT > BTRY)
BNEXT = BTRY;
BTRY = abs(B(I,J));
JC = J;
end
end
end
end
end
if(BNEXT < (BMAX*BTRY))
BMAX=BNEXT/BTRY;
IROW=I;
JCOL=JC;
end
end
if(ID(JC)==0)
ID(JCOL)=1;
else
DETERMINANT=0;
return
end
if(JCOL==IROW)
F=1/(B(JCOL,JCOL));
for J=1:N
B(JCOL,J)=B(JCOL,J)*F;
end
for K=1:M
D(JCOL,K)=D(JCOL,K)*F;
end
else
for J=1:N
SAVE=B(IROW,J);
B(IROW,J)=B(JCOL,J);
B(JCOL,J)=SAVE;
end
for K=1:M
SAVE=D(IROW,K);
D(IROW,K)=D(JCOL,K);
D(JCOL,K)=SAVE;
end
end
for I=1:N
if(I ~= JCOL)
F=B(I,JCOL);
for J=1:N
B(I,J)=B(I,J)-F*B(JCOL,J);
end
for K=1:M
D(I,K)=D(I,K)-F*D(JCOL,K);
end
end
end
end
return
Testing the MATLAB version of BAND
The MATLAB software below tests BAND for two cases. The first is a simple ordinary differential equation with an analytical solution and is not discussed further here. The second test employs the coupled equations:
22)
and
23)
The analytical solutions are
24)
and
25)
The table below compares the analytical solution with that obtained numerically using BAND with a mesh of 101 grid points. Refer to the software for the values of the coefficients used in BAND.
global A B C D G X Y N NJ
NPOINTS=101;
N=1;
NJ=NPOINTS;
H=(3.14159/2)/(NPOINTS-1);
%FIRST TEST WITH THE ODE D2CBYDX=C WITH C(0)=0, C(1)=1
for I=1:NPOINTS
A(1,1)=1;
B(1,1)=-2+H^2;
D(1,1)=1;
G(1)=0;
if(I == 1)
A(1,1)=0;
X(1,1)=0;
B(1,1)=1;
D(1,1)=0;
G(1)=0;
end
if(I == NPOINTS)
A(1,1)=0;
Y(1,1)=0;
B(1,1)=1;
D(1,1)=0;
G(1)=1;
end
band(I);
end
% %CALCULATE ANALYTICAL SOLUTION AND WRITE COMPARISON WITH NUMERICAL VALUES
for I=1:NPOINTS
ANS(I)=sin((3.14159/2)*(I-1)/(NPOINTS-1));
end
for I=1:NPOINTS
sprintf('\t %2.4f \t %2.4f ',C(1,I),ANS(I));
end
% NOW A PAIR OF SIMULTANEOUS EQUATIONS
% EQUATIONS ARE:
% D2YBYDT2 = Y-X WITH Y(0) = 1, Y(1)=0
% D2YBYDT2 = X-Y WITH X(0) = 0, X(1)=1
H = 1/(NPOINTS-1);
N=2;
for I=1:NPOINTS
A(1,1)=1;
A(1,2)=0;
A(2,1)=0;
A(2,2)=1;
B(1,1)=-2-H^2;
B(1,2)=H^2;
B(2,1)=H^2;
B(2,2)=B(1,1);
D(1,1)=A(1,1);
D(1,2)=A(1,2);
D(2,1)=A(2,1);
D(2,2)=A(2,2);
G(1)=0;
G(2)=0;
if(I == 1)
A(1,1)=0;
A(1,2)=0;
A(2,1)=0;
A(2,2)=0;
X(1,1)=0;
X(1,2)=0;
X(2,1)=0;
X(2,2)=0;
B(1,1)=1;
B(1,2)=0;
B(2,1)=0;
B(2,2)=1;
D(1,1)=0;
D(1,2)=0;
D(2,1)=0;
D(2,2)=0;
G(1)=1;
G(2)=0;
end
if(I == NPOINTS)
A(1,1)=0;
A(1,2)=0;
A(2,1)=0;
A(2,2)=0;
Y(1,1)=0;
Y(1,2)=0;
Y(2,1)=0;
Y(2,2)=0;
B(1,1)=1;
B(1,2)=0;
B(2,1)=0;
B(2,2)=1;
D(1,1)=0;
D(1,2)=0;
D(2,1)=0;
D(2,2)=0;
G(1)=0;
G(2)=1;
end
band(I);
end
% NOW FOR ANALYTICAL SOLUTION
V = exp(sqrt(2));
for I=1:NPOINTS
T=(I-1)/(NPOINTS-1);
YANS(I)=-0.5*(-V-exp(-(T-1)*sqrt(2))+1+exp(T*sqrt(2)))/(V-1);
XANS(I)=+0.5*(+V-exp(-(T-1)*sqrt(2))-1+exp(T*sqrt(2)))/(V-1);
end
% WRITE COMPARISON OF NUMERICAL AND ANALYTICAL SOLUTION
% Y=C(1,*), X=C(2,*)
for I=1:NPOINTS
fprintf('\t %2.4f \t %2.4f \t %2.4f \t %2.4f',C(1,I),C(2,I),YANS(I),XANS(I))
fprintf('\n')
end
Matlab Output:
Test1
Numerical Analytical
0.0000 0.0000
0.0157 0.0157
0.0314 0.0314
0.0471 0.0471
0.0628 0.0628
0.0785 0.0785
0.0941 0.0941
0.1097 0.1097
0.1253 0.1253
0.1409 0.1409
0.1564 0.1564
0.1719 0.1719
0.1874 0.1874
0.2028 0.2028
0.2181 0.2181
0.2334 0.2334
0.2487 0.2487
0.2639 0.2639
0.2790 0.2790
0.2940 0.2940
0.3090 0.3090
0.3239 0.3239
0.3387 0.3387
0.3535 0.3535
0.3681 0.3681
0.3827 0.3827
0.3972 0.3971
0.4115 0.4115
0.4258 0.4258
0.4399 0.4399
0.4540 0.4540
0.4679 0.4679
0.4818 0.4818
0.4955 0.4955
0.5090 0.5090
0.5225 0.5225
0.5358 0.5358
0.5490 0.5490
0.5621 0.5621
0.5750 0.5750
0.5878 0.5878
0.6004 0.6004
0.6129 0.6129
0.6252 0.6252
0.6374 0.6374
0.6495 0.6494
0.6613 0.6613
0.6730 0.6730
0.6846 0.6845
0.6959 0.6959
0.7071 0.7071
0.7181 0.7181
0.7290 0.7290
0.7396 0.7396
0.7501 0.7501
0.7604 0.7604
0.7705 0.7705
0.7804 0.7804
0.7902 0.7902
0.7997 0.7997
0.8090 0.8090
0.8182 0.8181
0.8271 0.8271
0.8358 0.8358
0.8443 0.8443
0.8526 0.8526
0.8607 0.8607
0.8686 0.8686
0.8763 0.8763
0.8838 0.8838
0.8910 0.8910
0.8980 0.8980
0.9048 0.9048
0.9114 0.9114
0.9178 0.9178
0.9239 0.9239
0.9298 0.9298
0.9354 0.9354
0.9409 0.9409
0.9461 0.9461
0.9511 0.9511
0.9558 0.9558
0.9603 0.9603
0.9646 0.9646
0.9686 0.9686
0.9724 0.9724
0.9759 0.9759
0.9792 0.9792
0.9823 0.9823
0.9851 0.9851
0.9877 0.9877
0.9900 0.9900
0.9921 0.9921
0.9940 0.9940
0.9956 0.9956
0.9969 0.9969
0.9980 0.9980
0.9989 0.9989
0.9995 0.9995
0.9999 0.9999
1.0000 1.0000
Test2
Numerical Analytical
1.0000 0.0000 1.0000 0.0000
0.9884 0.0116 0.9884 0.0116
0.9770 0.0230 0.9770 0.0230
0.9656 0.0344 0.9656 0.0344
0.9543 0.0457 0.9543 0.0457
0.9431 0.0569 0.9431 0.0569
0.9320 0.0680 0.9320 0.0680
0.9210 0.0790 0.9210 0.0790
0.9101 0.0899 0.9101 0.0899
0.8993 0.1007 0.8993 0.1007
0.8885 0.1115 0.8885 0.1115
0.8778 0.1222 0.8778 0.1222
0.8672 0.1328 0.8672 0.1328
0.8566 0.1434 0.8566 0.1434
0.8462 0.1538 0.8462 0.1538
0.8358 0.1642 0.8358 0.1642
0.8254 0.1746 0.8254 0.1746
0.8152 0.1848 0.8152 0.1848
0.8050 0.1950 0.8050 0.1950
0.7948 0.2052 0.7948 0.2052
0.7848 0.2152 0.7848 0.2152
0.7747 0.2253 0.7747 0.2253
0.7648 0.2352 0.7648 0.2352
0.7548 0.2452 0.7548 0.2452
0.7450 0.2550 0.7450 0.2550
0.7351 0.2649 0.7351 0.2649
0.7254 0.2746 0.7254 0.2746
0.7157 0.2843 0.7157 0.2843
0.7060 0.2940 0.7060 0.2940
0.6963 0.3037 0.6963 0.3037
0.6867 0.3133 0.6867 0.3133
0.6772 0.3228 0.6772 0.3228
0.6676 0.3324 0.6676 0.3324
0.6581 0.3419 0.6581 0.3419
0.6487 0.3513 0.6487 0.3513
0.6392 0.3608 0.6392 0.3608
0.6298 0.3702 0.6298 0.3702
0.6204 0.3796 0.6204 0.3796
0.6111 0.3889 0.6111 0.3889
0.6018 0.3982 0.6018 0.3982
0.5924 0.4076 0.5924 0.4076
0.5831 0.4169 0.5831 0.4169
0.5739 0.4261 0.5739 0.4261
0.5646 0.4354 0.5646 0.4354
0.5553 0.4447 0.5553 0.4447
0.5461 0.4539 0.5461 0.4539
0.5369 0.4631 0.5369 0.4631
0.5276 0.4724 0.5276 0.4724
0.5184 0.4816 0.5184 0.4816
0.5092 0.4908 0.5092 0.4908
0.5000 0.5000 0.5000 0.5000
0.4908 0.5092 0.4908 0.5092
0.4816 0.5184 0.4816 0.5184
0.4724 0.5276 0.4724 0.5276
0.4631 0.5369 0.4631 0.5369