// solves triadiagonal linear systems // using LU decomposition. // // Data: // // N - size of the system (NxN) // DM - (i, i-1) diagonal // D - (i, i) diagonal // DP - (i, i+1) diagonal // // Author: Alexander Bogomolny, CTK Software, Inc. // URL: http://www.cut-the-knot.com // Date: November 30, 2000 // Copyright: A. Bogomolny // Permission to use and modify the file is therefore granted // as long as this comment remains unchanged. Do this at your // own risk. // class Tridiagonal { int N; double[] DM, D, DP; // three diagonals // LU decomposition double[] A, C; // L has two non zero diagonals A & DM, U has units on main diagonal and C public Tridiagonal(int n, double[] m, double[] d, double[] p) { N = n; DM = m; D = d; DP = p; FindLU(); } // the same object may be used repeatedly // public void SetMatrix(int n, double[] m, double[] d, double[] p) { N = n; DM = m; D = d; DP = p; FindLU(); } public void SetMatrix(double[] m, double[] d, double[] p) { DM = m; D = d; DP = p; FindLU(); } // standard LU decomposition. // See K. Atkinson "An Introduction // to Numerical Analysis", p 455 (Chapter 8) // public void FindLU() { A = new double[N]; C = new double[N-1]; A[0] = D[0]; C[0] = DP[0]/D[0]; for (int i = 1; i < N-1; i++) { A[i] = D[i] - DM[i-1]*C[i-1]; C[i] = DP[i]/A[i]; } A[N-1] = D[N-1] - DM[N-2]*C[N-2]; } // solve the system for a given // right-hand side r // public double[] Solve(double[] r) { double[] z = new double[N]; z[0] = r[0]/A[0]; for (int i = 1; i < N; i++) z[i] = (r[i] - DM[i-1]*z[i-1])/A[i]; double[] x = new double[N]; x[N-1] = z[N-1]; for (int i = N-2; i >= 0; i--) x[i] = z[i] - C[i]*x[i+1]; return x; } }