## 6 Numerical Methods and Algorithms
### Models at the Macro-Scale
#### 6.1.2 Solution of the \"Discrete in Space\"
\begin{table}
\begin{tabular}{c c} \hline \hline  & The matrix \(J_{n}\) is approximated by \(\mathbf{J}_{n}^{DN}\) where \\  & \(\mathbf{J}_{n}^{DN}=\left[n\beta_{k}^{DN}\right]_{k=1,\ldots,M^{\prime}}\) \\
**Discrete Newton** & and where \\  & \({}^{n}\beta_{k}^{DN}=\frac{1}{h_{k}}(G(\mathbf{u}_{n}+e_{k}\delta_{k})-G(\mathbf{u}_{n}))\). \\  & \({}^{n}\beta_{k}^{DN}\) should be considered an approximation of the \(k\)-th column of \(J_{n}\) obtained by a difference finite method with step size \(e_{k}\). The value of \(e_{k}\) must be selected so that the tangent matrix approximation is as accurate as feasible. Its value should be very tiny. \\ \hline
**Modified Newton** & The matrix \(\mathbf{J}^{MN}_{\phantom{MN}n}\) substitutes \(\mathbf{I}_{n}\). It is considered to be equal to \(J_{0}\) for all steps \(n\), i.e., \\  & \(\mathbf{J}^{MN}_{\phantom{MN}n}=\mathbf{J}_{0}\), \(\forall n=1,\ldots,N\) \\ \hline  & The tangent matrix \(\mathbf{I}_{n}\) is approximated by the “_secant matrix_” \(\mathbf{I}_{n}^{DN}\) for which the following relation is valid: \\  & Letting \(\mathbf{K}_{n}^{DN}=\left(\mathbf{I}_{n}^{DN}\right)^{-1}\), then \(\mathbf{K}_{n}^{DN}\) can be considered an approximation of the inverse of the tangent matrix \(J_{n}\) for which the following relation is valid : \\  & \(\mathbf{K}_{n}^{DN}\mathbf{g}_{n}=\mathbf{w}_{n}\) \\
**Quasi-Newton** & where \(\mathbf{w}_{n}=\mathbf{u}_{n}-\mathbf{u}_{n-1}\) and \(\mathbf{g}_{n}=-(G(\mathbf{u}_{n})-G(\mathbf{u}_{n-1}))\). A number of update algorithms exist for the explicit determination of \(\mathbf{K}_{n}^{DN}\). Here, the BFGS method will be presented since it has been observed to be successfully applied in the solution of problems from Finite Element Methods. It was based on the following update for \(\mathbf{K}_{n}^{DN}\): \\  & \(\mathbf{K}_{n}^{DN}=\left(\mathbf{1}+\mathbf{a}_{n}\mathbf{b}_{n}^{T}\right)\mathbf{K}_{n-1}^{DN} \left(\mathbf{1}+\mathbf{b}_{n}\mathbf{a}_{n}^{T}\right)\), \\  & where \\  & \(a_{n}=\frac{1}{g_{n}^{T}\mathbf{w}_{n}}\mathbf{w}_{n}\), \\  & \(b_{n}=-\left\{\mathbf{g}_{n}-\left[-\frac{\mathbf{w}_{n}^{T}\mathbf{g}_{n}}{\mathbf{w}_{n}^{T} \mathbf{G}(\mathbf{u}_{n-1})}\right]^{-\frac{1}{2}}G(\mathbf{u}_{n-1})\right\}\). \\ \hline \hline \end{tabular}
\end{table} 
Table 8: The Newton–Raphson-like schemes obtained from the original formulation by using different approaches to approximate/compute the Jacobi-matrix of \(\mathbf{G}\). where \(\mathbf{L}\) and \(\mathbf{U}\) are, respectively, lower triangular and upper triangular matrices and where \(\mathbf{P}\) and \(\mathbf{Q}\) are permutation matrices. A permutation matrix is used to represent row or column interchanges in matrix \(\mathbf{A}\), and it is obtained from the identity matrix \(\mathbf{I}\) by applying on it the same sequence of row or column interchanges. The row or column permutations on \(\mathbf{A}\) are often required for both numerical and sparsity issues [146]. The \(LU\) decomposition facilitates the solution of (125) effectively, particularly when solving several systems that share the same matrix \(\mathbf{A}\). In fact, if a matrix \(\mathbf{A}\) already has a \(LU\) decomposition available, the linear system in (125) can be solved by (1) the solution of \(L\mathbf{z}=\mathbf{P}\mathbf{y}\) (forward substitution), followed by (2) the solution of \(\mathbf{U}\mathbf{w}=\mathbf{z}\) (backward substitution) where \(\mathbf{w}=\mathbf{Q}^{T}\mathbf{x}\). Then, the values of all the unknowns \(x_{i},\ \forall i=1,\ldots,M\) can be obtained by permuting the rows of \(\mathbf{w}\) based on the permutation matrix \(\mathbf{Q}\). The execution of both a forward and a backward substitution has a general computational complexity of \(O\big{(}2M^{2}\big{)}\)). If the matrix \(\mathbf{A}\) is symmetric and positive definite (that is, if \(\mathrm{i}x^{T}\mathbf{A}\mathbf{x}>0\) for all nonzero vectors \(x\)), then the \(LU\) decomposition of \(\mathbf{A}\) becomes \(\mathbf{P}\mathbf{A}\mathbf{P}^{T}=\mathbf{L}\mathbf{L}^{T}\) and is called \"_Cholesky_\" factorization.