Solution Using Matrix Inversion

A set of simultaneous linear equations can also be solved by using matrices, as shown in Chapter 9. The solution matrix is obtained by multiplying the matrix of constants by the inverse of the matrix of coefficients. Applying this simple solution to the spectrophotometric data used above, the inverted matrix is obtained by selecting a 3R x 3C array of cells, entering the array formula 

The solution matrix is obtained by selecting a 3R x 1C array, then entering the array formula 

---

This allows us to represent partial differential equations as found in the balance equations using the collocation method. Equation (11.47) is a solution to a partial differential equation represented by a system of linear algebraic equations, formed by the interpolation coefficients, oij, and the operated radial functions. The interpolation coefficients are solved for using matrix inversion techniques to approximately satisfy the partial differential equation... 

Regardless of the type of model used, a method must be chosen for the self-consistent solution of the polarizable degrees of freedom. Direct solution via matrix inversion is nearly always avoided by most researchers in the field, because of the prohibitive O(N ) scaling with system size, N. Both iterative and predictive methods reduce the scaling to match that of the potential evaluation [O(N ) for direct summation 0(N In N) for Ewald-based meth-ods ° " 0 N) if interactions are neglected beyond some distance cutoff], but the cost of the iterations means that the predictive methods are always more efficient. Extended Lagrangian methods have been implemented for all four types of polarizable... 

There are various ways to obtain the solutions to this problem. The most straightforward method is to solve the full problem by first computing the Lagrange multipliers from the time-differentiated constraint equations and then using the values obtained to solve the equations of motion [7,8,37]. This method, however, is not computationally cheap because it requires a matrix inversion at every iteration. In practice, therefore, the problem is solved by a simple iterative scheme to satisfy the constraints. This scheme is called SHAKE [6,14] (see Section V.B). Note that the computational advantage has to be balanced against the additional work required to solve the constraint equations. This approach allows a modest increase in speed by a factor of 2 or 3 if all bonds are constrained. 

The above model was solved numerically by writing finite difference approximations for each term. The equations were decoupled by writing the reaction terms on the previous time steps where the concentrations are known. Similarly the equations were linearized by writing the diffusivities on the previous time step also. The model was solved numerically using a linear matrix inversion routine, updating the solution matrix between iterations to include the proper concentration dependent diffusivities and reactions. 

Gaussian elimination is a very efficient method for solving n equations in n unknowns, and this algorithm is readily available in many software packages. For solution of linear equations, this method is preferred computationally over the use of the matrix inverse. For hand calculations, Cramer s rule is also popular. 

Equation (15.23) displays the feature of locality that the blending functions should possess in order to be computationally advantageous that is, during the process of matrix inversion, one wishes the calculation to proceed quickly. As mentioned earlier, the use of linear approximation functions results in at most five terms on the left side of the equation analogous to (15.23), yielding a much crader approximation, but one more easily calculated. The current choice of Bezier functions, on the other hand, is rapidly convergent for methods such as relaxation, possesses excellent continuity properties (the solution is guaranteed to look and behave reasonably), and does not require substantial computation. 

We thus can obtain the solution vector Y from (6.107) explicitly. Then we can use Y in (6.104) to compute the vector X via another matrix inversion using formula (6.104). 

Matrix inversion is not widely used in practice, but from a theoretical point of view is extremely useful, because it allows us to calculate the minimum number of projections that are required for a complete reconstruction. If we have p projections of a structure, and each projection contains r rays, a reconstruction procedure amounts to solving a system of p-r equations in n2 unknowns, and algebra tells us that a solution exists only if the number of linearly independent equations is equal to the number of the unknowns. 

The standard method of solving the problem by the equilibrium constant approach is to use linearized matrix inversion. Convergence assumes, of course, that the solution not only exists but that it is unique. If a system can have several thermodynamically metastable states (local minima in the Gibbs function) then several nonunique solutions are possible. 

Matrix A is a 3N x 3N dense matrix. For a small number of unknowns, direct solvers are practical, especially in the case of multiple sources. One can use different types of iterative methods, discussed in Chapter 4, for the solution of this problem. However, if N is large, the storage of A is extremely memory consuming, not to mention the complexity of direct matrix inversion. 

Note that when the problem is stiff (Pe>1000, or O > 20), N > 20 node points might be needed for accurate solutions. Consequently, inverting the A matrix symbolically involves lot of computational effort as the order of the matrix increases with N. It is recommended that you specify the values for the parameters and convert the entries of the A matrix to decimal points using the following Maple command before the matrix inversion ... 

---