'************************************************************************* '* This program solves a tridiagonal linear set of equations M * U = R * '* --------------------------------------------------------------------- * '* SAMPLE RUN: * '* * '* Input data file (tridiag.dat): * '* * '* 5 * '* 0.0 1.0 -2.0 * '* 10.0 5.0 -7.5 * '* 2.22 -3.25 0.00456 * '* 2.0 4.0 -3.3 * '* -1.0 3.0 0.0 * '* * '* (The first line contains the size of system * '* Then the three columns give the three main diagonals of matrix M, the * '* first column begins with a zero and the third column ends with zero) * '* * '* Output file (tridiag.lst): * '* * '* -------------------------------------------------------------------- * '* Tridiagonal system of equations * '* -------------------------------------------------------------------- * '* upper subdiagonal: * '* -2.000000 -7.500000 0.004560 -3.300000 0.000000 * '* main diagonal: * '* 1.000000 5.000000 -3.250000 4.000000 3.000000 * '* lower subdiagonal: * '* 0.000000 10.00000 2.220000 2.000000 -1.000000 * '* right side member: * '* 1.000000 1.000000 1.000000 1.000000 1.000000 * '* -------------------------------------------------------------------- * '* System solution: * '* -0.136497 -0.568249 -0.694162 1.202870 0.734290 * '* -------------------------------------------------------------------- * '* * '* Reference: "Numerical Recipes By W.H. Press, B.P. Flannery, S.A. Teu- * '* kolsky, W.T. Vetterling, Cambridge University Press, 1987 * '* [BIBLI 08]. * '* * '* Basic version from FORTRAN by J-P Moreau, Paris * '* (www.jpmoreau.fr) * '************************************************************************* defint i-n defdbl a-h,o-z dim a(25), b(25), c(25), r(25), u(25) 'irc : INTEGER cls print print " Input datafile name (without .dat): "; : input nom\$ open nom\$+".dat" for input as #1 open nom\$+".lst" for output as #2 input #1, n for i=1 to n input #1, a(i),b(i),c(i) r(i)=1# next i close #1 print #2, print #2,"-------------------------------------------------------" print #2," Tridiagonal system of equations" print #2,"-------------------------------------------------------" print #2," upper subdiagonal:" for i=1 to n print #2,using " ##.######"; c(i); next i print #2, print #2," main diagonal:" for i=1 to n print #2,using " ##.######"; b(i); next i print #2, print #2," lower subdiagonal:" for i=1 to n print #2,using " ##.######"; a(i); next i print #2, print #2," right side member:" for i=1 to n print #2,using " ##.######"; r(i); next i print #2, print #2,"-------------------------------------------------------" 'TRIDAG(a,b,c,r,u,n,irc) gosub 1000 if irc<>0 then print #2, " ERROR IN TRIDAG : RETURN CODE="; irc else print #2, " System solution:" for i=1 to n print #2,using " ##.######"; u(i); next i end if print #2, print #2,"-------------------------------------------------------" print print " Results in ";nom\$;".lst." close #2 END 'of main program 1000 'SUBROUTINE TRIDAG '******************************************************************** '* Solves for a vector U of length N the tridiagonal linear set * '* M U = R, where A, B and C are the three main diagonals of matrix * '* M(N,N), the other terms are 0. R is the right side vector. * '******************************************************************** dim GAM(25) IF B(1)=0 THEN irc=1 return end if BET=B(1) U(1)=R(1)/BET for J=2 to N 'Decomposition and forward substitution GAM(J)=C(J-1)/BET BET=B(J)-A(J)*GAM(J) IF BET=0 THEN 'Algorithm fails irc=2 return end if U(J)=(R(J)-A(J)*U(J-1))/BET next J for J=N-1 to 1 step -1 'Back substitution U(J)=U(J)-GAM(J+1)*U(J+1) next J irc=0 Return 'End of file tridiag.bas