ORMAT Tridiagonal in a, b, c 4´4[1] for\tbtri.wpj
1)
Note that a1 and c4 are not in equation . Write the set of equations
eqn # / col1 / col2 / col3 / Col4 / Rhs11 / b1 u1 / c1 u2 / 0 / 0 / = / r1
2 / a2 u1 / b2 u2 / c2 u3 / 0 / = / r2
3 / 0 / a3 u2 / b3 u3 / c3 u4 / = / r3
4 / 0 / 0 / a4 u3 / b4 u4 / = / r4
Divide equation 1 by b1
eqn # / col1 / col2 / col3 / Col4 / Rhs11 / u1 / c1/b1 u2 / 0 / 0 / = / r1/b1
2 / a2 u1 / b2 u2 / c2 u3 / 0 / = / r2
3 / 0 / a3 u2 / b3 u3 / c3 u4 / = / r3
4 / 0 / 0 / a4 u3 / b4 u4 / = / r4
Define U1 =r1/bt and Define g1 = c1 /b1.
eqn # / col1 / col2 / col3 / Col4 / Rhs11 / u1 / g1u2 / 0 / 0 / = / U1
2 / a2 u1 / b2 u2 / c2 u3 / 0 / = / r2
3 / 0 / a3 u2 / b3 u3 / c3 u4 / = / r3
4 / 0 / 0 / a4 u3 / b4 u4 / = / r4
Multiply the first equation by a2 and subtract it from the second equation.
eqn # / col1 / col2 / col3 / Col4 / Rhs11 / u1 / g1u2 / 0 / 0 / = / U1
2 / 0 / (b2-a2g1)u2 / c2 u3 / 0 / = / r2-a2U1
3 / 0 / a3 u2 / b3 u3 / c3 u4 / = / r3
4 / 0 / 0 / a4 u3 / b4 u4 / = / r4
Define new bt = (b2-a2g1) and divide equation 2 by this
eqn # / col1 / col2 / col3 / Col4 / Rhs11 / u1 / g1u2 / 0 / 0 / = / U1
2 / 0 / u2 / (c2/bt) u3 / 0 / = / (r2-a2U1)/bt
3 / 0 / a3 u2 / b3 u3 / c3 u4 / = / r3
4 / 0 / 0 / a4 u3 / b4 u4 / = / r4
Define U2 =(r2-a2U1)/bt and Define g2 = c2 /bt
eqn # / col1 / col2 / col3 / Col4 / Rhs11 / u1 / g1u2 / 0 / 0 / = / U1
2 / 0 / u2 / g2u3 / 0 / = / U2
3 / 0 / a3 u2 / b3 u3 / c3 u4 / = / r3
4 / 0 / 0 / a4 u3 / b4 u4 / = / r4
Multiply the second equation by a3 and subtract it from the third equation.
eqn # / col1 / col2 / col3 / Col4 / Rhs11 / u1 / g1u2 / 0 / 0 / = / U1
2 / 0 / u2 / g2u3 / 0 / = / U2
3 / 0 / 0 / (b3-a3g2) u3 / c3 u4 / = / r3-a3U2
4 / 0 / 0 / a4 u3 / b4 u4 / = / r4
Define new bt = (b3-a3g2) and divide equation 3 by this
eqn # / col1 / col2 / col3 / Col4 / Rhs11 / u1 / g1u2 / 0 / 0 / = / U1
2 / 0 / u2 / g2u3 / 0 / = / U2
3 / 0 / 0 / u3 / (c3/bt) u4 / = / (r3-a3U2)/bt
4 / 0 / 0 / a4 u3 / b4 u4 / = / r4
Define U3 =(r3-a3U2)/bt and Define g3 = c3 /bt
eqn # / col1 / col2 / col3 / Col4 / Rhs11 / u1 / g1u2 / 0 / 0 / = / U1
2 / 0 / u2 / g2u3 / 0 / = / U2
3 / 0 / 0 / u3 / g3u4 / = / U3
4 / 0 / 0 / a4 u3 / b4 u4 / = / r4
Define U4 =(r4-a4U3)/bt and Define g4 = c4 /bt – skip a few steps
eqn # / col1 / col2 / col3 / Col4 / Rhs11 / u1 / g1u2 / 0 / 0 / = / U1
2 / 0 / u2 / g2u3 / 0 / = / U2
3 / 0 / 0 / u3 / g3u4 / = / U3
4 / 0 / 0 / u4 / = / U4
This last is the N’th term and we have found u4. Then u3 = U3-g3u4 and so on back to u1. Note that U and u can share the same storage locations. This is the Press code from section 2.6 with the addition of real*8 and some comments TRIDAG.FOR. Note that gamma needs to be stored in a work space dimensioned 500.
In Weighted Fourier Series.doc the limitation on the size of an array is determined by the storage. If the H’ terms are all one or zero the equation set involves the same a, b, and c for each term. Thus for\tbtri.wpj uses inputs of a, b, and c as numbers and creates a direct access temporary file for gamma.
Tridiagonal Matrices[2]
To solve Ax=b, for i = 2,…,N
Then back substitute for i = N-1, … ,1
“You might expect that we would next assemble a program to compute the inverse of a tridiagonal matrix. There is only one problem: The inverse of a tridiagonal matrix is not necessarily sparse.”
Diagonal in a, b, c,d,e 8´8
eqn # / col1 / col2 / col3 / Col4 / C5 / C6 / C7 / C8 / Rhs11 / c A11 / dA21 / eA31 / 0 / 0 / 0 / 0 / 0 / = / I11
2 / b A11 / c A21 / dA31 / eA41 / 0 / 0 / 0 / 0 / = / I12
3 / a A11 / b A21 / c A31 / dA41 / eA51 / 0 / 0 / 0 / = / I13
4 / 0 / a A21 / b A31 / c A41 / dA51 / eA61 / 0 / 0 / = / I14
5 / 0 / 0 / a A21 / b A41 / cA51 / dA61 / eA71 / 0 / = / I15
6 / 0 / 0 / 0 / a A41 / b A51 / cA61 / dA71 / eA81 / = / I16
7 / 0 / 0 / 0 / 0 / a A51 / b A61 / cA71 / dA81 / = / I17
8 / 0 / 0 / 0 / 0 / 0 / a A61 / b A71 / cA81 / = / I18
Multiply 1 by b/c and subtract from 2, then by a/c and subtract from 3
eqn # / col1 / col2 / col3 / Col4 / C5 / C6 / C7 / C8 / Rhs11 / c A11 / dA21 / eA31 / 0 / 0 / 0 / 0 / 0 / = / I11
2 / 0 / (c-d b/c) A21 / (d-e b/c)A31 / eA41 / 0 / 0 / 0 / 0 / = / I12-b/c I11
3 / 0 / (b-d b/a) A21 / (c –e b/a)A31 / dA41 / eA51 / 0 / 0 / 0 / = / I13-b/a I11
4 / 0 / a A21 / b A31 / c A41 / dA51 / eA61 / 0 / 0 / = / I14
5 / 0 / 0 / a A21 / b A41 / cA51 / dA61 / eA71 / 0 / = / I15
6 / 0 / 0 / 0 / a A41 / b A51 / cA61 / dA71 / eA81 / = / I16
7 / 0 / 0 / 0 / 0 / a A51 / b A61 / cA71 / dA81 / = / I17
8 / 0 / 0 / 0 / 0 / 0 / a A61 / b A71 / cA81 / = / I18
Define, ctt=(c –e b/a), bt = b – d (b/a), btt = d – d (b/a), It=I12-I11 (b/c) and Itt = I13 – I11 (b/a) – drop the first row and column
1. ct = c – d (b/c)
2. bt = b – d (b/a)
3. dt = d – e (b/c)
4. ctt = c – e (b/a)
eqn # / col2 / col3 / Col4 / C5 / C6 / C7 / C8 / Rhs12 / ct A21 / dt A31 / eA41 / 0 / 0 / 0 / 0 / = / It
3 / bt A21 / ctt A31 / dA41 / eA51 / 0 / 0 / 0 / = / Itt
4 / a A21 / b A31 / c A41 / dA51 / eA61 / 0 / 0 / = / I14
5 / 0 / a A21 / b A41 / cA51 / dA61 / eA71 / 0 / = / I15
6 / 0 / 0 / a A41 / b A51 / cA61 / dA71 / eA81 / = / I16
7 / 0 / 0 / 0 / a A51 / b A61 / cA71 / dA81 / = / I17
8 / 0 / 0 / 0 / 0 / a A61 / b A71 / cA81 / = / I18
Again make the subtractions, watch out for the ctt. Eventually these reduce to 3 equations in 3 unknowns. While one can continue, this is Gauss Jordan.
Band Diagonal Systems[3]
Band diagonal systems have m1 elements below the diagonal and m2 above the diagonal
3 / 1 / 0 / 0 / 0 / 0 / 04 / 1 / 5 / 0 / 0 / 0 / 0
9 / 2 / 6 / 5 / 0 / 0 / 0
0 / 3 / 5 / 8 / 9 / 0 / 0
0 / 0 / 7 / 9 / 3 / 2 / 0
0 / 0 / 0 / 3 / 8 / 4 / 6
0 / 0 / 0 / 0 / 2 / 4 / 4
has N (number of rows 7, m1 = 2, and m2 = 1. It is stored as
X / X / 3 / 1X / 4 / 1 / 5
9 / 3 / 6 / 5
3 / 5 / 8 / 9
7 / 9 / 3 / 2
3 / 8 / 4 / 6
2 / 4 / 4 / X
The routine bandec calculates the upper triangle in the space of A above. This routine is used with banbks to solve band-diagonal sets ofequations.
The routines bandec and banbks are based on the Handbook routines bandet1 and bansol1 in Wilkinson and Reinsh[4]
[1] Press – section 2.6
[2] A.L. Gasrcia, Numerical Methods for Physics, Prentice Hall (1994), pp 253-255
[3] Press, Numerical Recipes in C, Cambridge University Press 1988,1992 section 2.4
[4] Wilkinson, J. H., and Reinsch, C., 1971, Linear Algebra, vol II of Handbook for Automatic Computation (New York: Springer-Verlag, Chapter 1/6.