Computational methods Jacobian matrix

Two kinds of method can be used to solve this system Newton type methods, in which it is necessary to compute a Jacobian matrix, and direct iterative methods, in which it is necessary to compute only the functions fj, f2,. . ., f . 

By hand, compute the Jacobian matrix of this system, and generate the refined estimate starting from an initial guess ofxl = [2 1]. In the reduced-step Newton method, would this estimate be accepted ... 

Steady-state solutions are found by iterative solution of the nonlinear residual equations R(a,P) = 0 using Newton s methods, as described elsewhere (28). Contributions to the Jacobian matrix are formed explicitly in terms of the finite element coefficients for the interface shape and the field variables. Special matrix software (31) is used for Gaussian elimination of the linear equation sets which result at each Newton iteration. This software accounts for the special "arrow structure of the Jacobian matrix and computes an LU-decomposition of the matrix so that qu2usi-Newton iteration schemes can be used for additional savings. 

If the Newton-Raphson method is used to solve Eq. (1), the Jacobian matrix (df/3x)u is already available. The computation of the sensitivity matrix amounts to solving the same Eq. (59) with m different right-hand side vectors which form the columns — (3f/<5u)x. Notice that only the partial derivatives with respect to those external variables subject to actual changes in values need be included in the m right-hand sides. 

A more efficient way of solving the DFT equations is via a Newton-Raphson (NR) procedure as outlined here for a fluid between two surfaces. In this case one starts with an initial guess for the density profile. The self-consistent fields are then calculated and the next guess for density profile is obtained through a single-chain simulation. The difference from the Picard iteration method is that an NR procedure is used to estimate the new guess from the density profile from the old one and the one monitored in the single-chain simulation. This requires the computation of a Jacobian matrix in the course of the simulation, as described below. 

In this section we consider how Newton-Raphson iteration can be applied to solve the governing equations listed in Section 4.1. There are three steps to setting up the iteration (1) reducing the complexity of the problem by reserving the equations that can be solved linearly, (2) computing the residuals, and (3) calculating the Jacobian matrix. Because reserving the equations with linear solutions reduces the number of basis entries carried in the iteration, the solution technique described here is known as the reduced basis method. ... 

Fig. 4.4. Comparison of the computing effort, expressed in thousands of floating point operations (Aflop), required to factor the Jacobian matrix for a 20-component system (Nc = 20) during a Newton-Raphson iteration. For a technique that carries a nonlinear variable for each chemical component and each mineral in the system (top line), the computing effort increases as the number of minerals increases. For the reduced basis method (bottom line), however, less computing effort is required as the number of minerals increases.

Such a scheme is sometimes called a soft Newton-Raphson formulation because the partial derivatives in the Jacobian matrix are incomplete. We could, in principle, use a hard formulation in which the Jacobian accounts for the devia-tives dy/dm,i and daw/dm,i. The hard formulation sometimes converges in fewer iterations, but in tests, the advantage was more than offset by the extra effort in computing the Jacobian. The soft method also allows us to keep the method for calculating activity coefficients (see Chapter 8) separate from the Newton-Raphson formulation, which simplifies programming. 

Note that the Jacobian matrix dh/dx on the left-hand side of Equation (A.26) is analogous to A in Equation (A.20), and Ax is analogous to x. To compute the correction vector Ax, dh/dx must be nonsingular. However, there is no guarantee even then that Newton s method will converge to an x that satisfies h(x) = 0. 

ODE solver. Relative to non-stiff ODE solvers, stiff ODE solvers typically use implicit methods, which require the numerical inversion of an Ns x Ns Jacobian matrix, and thus are considerably more expensive. In a transported PDF simulation lasting T time units, the composition variables must be updated /Vsm, = T/At 106 times for each notional particle. Since the number of notional particles will be of the order of A p 106, the total number of times that (6.245) must be solved during a transported PDF simulation can be as high as A p x A sim 1012. Thus, the computational cost associated with treating the chemical source term becomes the critical issue when dealing with detailed chemistry. 

On the contrary, the Newton method gives robusmess and quadratic rate of convergence. The main drawback of this algorithm is the computation of a Jacobian matrix and the resolution of a large linear system at each step. This method has been largely described and used successfully. The reader is refered to [9] for more details. 

To solve the equations of problem P ), the optimization algorithms used are Levenberg-Marquardt and Trust-Region procedures. These methods enable computation of the solution by using the Jacobian matrix and the Hessian matrix (or its approximation) related to the objective function E(Y) [57]. 

The two points retained for the next step are x and either Xq or x, the choice being made so that the pair of values/(Jc), and either/(jc ) or/(jc ), have opposite signs to maintain the bracket on x. (This variation is called regula falsi or the method of felse position.) In Figure L.7, for the k + l)st stage, x and Xq would be selected as the end points of the secant line. Secant methods may seem crude, but they work well in practice. The details of the computational aspects of a sound algorithm to solve multiple equations by the secant method are too lengthy to outline here (particularly the calculation of a new Jacobian matrix from the former one instead refer to Dennis and Schnabel"). 

In the classical Newton-Raphson technique, the Jacobian matrix is inverted every iteration in order to compute the corrections AT] and Al]. The method of Tomich, however, uses the Broyden procedure (Broyden, 1965) in subsequent iterations for updating the inverted Jacobian matrix. 

This form of the general Jacobian element allows for the straightforward solution of the SSOZ equation for molecules of arbitrary symmetry. However, in the numerical solution using Gillan s methods, most of the computation time is involved in calculating the elements of the Jacobian matrix, rather than in the calculation of its inverse or in the calculation of transforms. Indeed, as the forward and backward Fourier transforms can be carried out using a fast Fourier transform routine, the time-limiting step is the double summation over / and j in Eq. (4.3.36). With this restriction in mind, it is... 

C.M. Bethke (52) has shown that significant numerical advantages in such calculations can be realized by switching into the basis set a mineral species that is in partial equilibrium with the aqueous phase. This avoids expansion of the size of the Jacobian matrix and reduces computation time. A method based on this concept is being developed for use in the 3270 version of EQ3/6. The concept appears to show promise for improvement of the "optimizer" algorithm as well as the Newton-Raphson one. 

The residues (Nb-kNp at each node) are reduced to zero (a small positive number fixed by specifying an error tolerance at input) iteratively by computing corrections to current values of the unknowns using the Newton-Raphson method (14). Elements of the Jacobian matrix required by this method are computed from analytical expressions. The system of equations to be solved for the corrections has block tridiagonal form and is solved by use of a published software routine (1.5b... 

Tomich25 was the first to apply Broyden s method (developed in chap. 15) to the solution of distillation problems. Broyden s method is based on the use of numerical approximations of the partial derivatives appearing in the jacobian matrix. The approach proposed by Broyden permits the inverse of the jacobian matrix to be updated each trial after the first through the use of Householder s formula.6 Thus, it is necessary to invert the jacobian matrix only once. Since approximate values for the partial derivatives are used, procedure 2 generally requires more trials than does procedure 1. However, since the evaluation of the partial derivatives and the inversion of the jacobian matrix are not generally required after the first trial of procedure 2, it requires less computer time per trial than does procedure 1. 

In the class of methods proposed by Broyden, the partial derivatives df/dxj in the jacobian matrix Jk of Eq. (4-29) are generally evaluated only once. In each successive trial, the elements of the inverse of the jacobian matrix are corrected by use of the computed values of the functions. An algebraic example will be given after the calculational procedure proposed by Broyden has been presented. 

After the Broyden correction for the independent variables has been computed, Broyden proposed that the inverse of the jacobian matrix of the Newton-Raphson equations be updated by use of Householder s formula. Herein lies the difficulty with Broyden s method. For Newton-Raphson formulations such as the Almost Band Algorithm for problems involving highly nonideal solutions, the corresponding jacobian matrices are exceedingly sparse, and the inverse of a sparse matrix is not necessarily sparse. The sparse characteristic of these jacobian matrices makes the application of Broyden s method (wherein the inverse of the jacobian matrix is updated by use of Householder s formula) impractical. 

Numerical methods used to solve a system of ODEs are widely available in computational libraries and through texts such as Numerical Recipes. Certain considerations arise in the use of these standard techniques for nonlinear systems, particularly in models of chemical systems, which often consist of systems of stiff equations that require special care. Stiff equations are characterized by the presence of widely differing time scales, which leads to eigenvalues of the Jacobian matrix differing by many orders of magnitude. 

This requires n + 1 passes through the recycle loop to complete the Jacobian matrix for just one iteration of the Newton-Raphson method that is, for n = 10, eleven passes are necessary, usually involving far too many computations to be competitive. 

Since the system is three-dimensional in nature, the VdelR condition may be used to compute the critical a policy for the system. Since the system is no longer expressed in terms of concentrations, however, the specific critical a policy given by Equation 7.3 cannot be used. Rather, the Jacobian matrix expressed in terms of mass fractions dr(z) must be computed and then used to compute (p(z), as in Chapter 7. The computations are somewhat lengthy, and hence this approach will not be adopted here. Instead, the parallel complement automated AR construction method discussed in Chapter 8 shall be employed, providing a quick means to compute the AR for use in comparisons. 

Methods III and IV have reduced computational complexities of 0(N ). Method III is the most efficient tqrproach for the joint space inertia matrix of all those listed for N < 6, and it computes the manipulator Jacobian matrix as well. Method IV, although the simultaneous Jacobian computation has been eliminated, is the most efficient tqrproach for computing the joint space inertia matrix for N >6. Thus, these two methods together provide the most efficient calculation of H for all N. 

---