## 6 Numerical Methods and Algorithms
### Models at the Macro-Scale
#### 6.1.2 Solution of the \"Discrete in Space\"
Only the equation system at the line 7 of Algorithm 2 has to be reformulated by introducing the BFGS-update of the inverse of the secant matrix. See Algorithm 3 for details on the BFGS reformulation of Algorithm 2. Algorithms 2 and 3 require the solution of a linear problem (see, respectively, lines 6 and 7) that can be solved by methods described in the next point 3. Details about convergence and other numerical issues are available in the literature (see [145]). ```
1:procedureBFGSNewtonSolver(\(\mathbf{G}\), \(\mathbf{u}_{0}\), \(N\), \(\mathbf{u}_{N}\))
2:Input:\(\mathbf{G}\), \(\mathbf{u}_{0}\), \(N\)
3:Output:\(\mathbf{u}_{N}\)
4:\(\mathbf{g}_{1}\leftarrow\)Compute \(\mathbf{g}_{1}=-\mathbf{G}(\mathbf{u}_{0})\)
5:\(\mathbf{J}_{0}^{\text{ON}}\leftarrow\)Compute \(\mathbf{J}_{0}^{\text{ON}}\)
6:\(\mathbf{K}_{0}^{\text{ON}}\leftarrow\)Solve \(\mathbf{J}_{0}^{\text{ON}}\mathbf{K}_{0}^{\text{ON}}=\mathbf{I}\)
7:\(\mathbf{w}_{1}\leftarrow\)Compute \(\mathbf{w}_{1}=\mathbf{K}_{0}^{\text{ON}}\mathbf{g}_{1}\)
8:for\(n=1\) to \(N-1\)do
9:\(\mathbf{u}_{n}\leftarrow\)Compute \(\mathbf{u}_{n}=\mathbf{u}_{n-1}+\mathbf{w}_{n}\)
10:\(\mathbf{G}_{n}\leftarrow\)Compute \(\mathbf{G}_{n}=\mathbf{G}(\mathbf{u}_{n-1})\)
11:\(\mathbf{a}_{n}\leftarrow\)Compute \(\mathbf{a}_{n}=\frac{1}{\mathbf{g}_{1}^{\text{IN}}\mathbf{w}_{n}}\mathbf{w}_{n}\)
12:\(\mathbf{b}_{n}\leftarrow\)Compute \(\mathbf{b}_{n}=-\Bigg{\{}\mathbf{g}_{n}-\Big{[}-\frac{\mathbf{w}_{n}^{\text{IN}}\mathbf{g}_{ n}}{\mathbf{w}_{n}^{\text{IN}}\mathbf{G}_{n}}\Big{]}^{-\frac{1}{2}}\mathbf{G}_{n}\Bigg{\}}\). 13:\(\mathbf{K}_{n}^{\text{ON}}\leftarrow\)Compute \(\mathbf{K}_{n}^{\text{ON}}=\big{(}\mathbf{1}+\mathbf{a}_{n}\mathbf{b}_{n}^{\text{IN}}\big{)} \mathbf{K}_{n-1}^{\text{ON}}\big{(}\mathbf{1}+\mathbf{b}_{n}\mathbf{a}_{n}^{\text{IN}}\big{)}\),
14:\(\mathbf{w}_{n+1}\leftarrow\)Compute \(\mathbf{w}_{n+1}=\mathbf{K}_{n}^{\text{ON}}\mathbf{g}_{n+1}\)
15:endfor
16:\(\mathbf{u}_{N}\leftarrow\)Compute \(\mathbf{u}_{N}=\mathbf{u}_{N-1}+\mathbf{w}_{N}\)
17:endprocedure
```  
**Algorithm 3** The algorithm implementing the _\"BFGS\"_ reformulation of the Newton-Raphson method. 3. **The linear problem solver**: Let us consider the linear problem of the form \[u_{n+1} =f(u_{n},\dot{u}_{n},) \text{explicit schema}\]9 (125) where \[u_{n+1} =f(u_{n},\dot{u}_{n},u_{n+1}) \text{implicit schema}\]0 Since the matrices coming from the discretization of problems of interest for this review are sparse if not even structured--for example, see matrices generated for the solution of Navier-Stokes equations for viscoelastic fluids in Section 6.1.3, which have a lot of zero elements (or have zero blocks)--particular attention will be spent in this work to describe methods to be used for the solution of (125) when \(\mathbf{A}\) is sparse (or structured) (see [146] for a more precise definition of the sparsity concept). There are three fundamental classes of methods that are used to solve Equation (125): **Direct methods** The direct solution of (125) is generally based on a technique called \"_LU decomposition_\" [147]. This technique consists of factoring the matrix \(\mathbf{A}\) as in Equation (126) (by algorithms whose computational complexity is generally on the order of \(M^{3}\)), and it is based on the concept that triangular systems of equations are \"easy\" to solve: \[u_{n+1} =f(u_{n},\dot{u}_{n},u_{n+1}) \text{implicit schema}\]1 (126)