# Vectorized Numerical Algorithms to Solve Internal Problems of Computational Fluid Dynamics
## 5 Discretization of Derivatives  
Consider a finite-volume discretization of the derivative of a mesh function **U** defined on a sequential set of \(N\) equally spaced nodes. This dataset is represented using a linear array containing \(N\) elements. To calculate the derivative, the function y = diff(x,n) is defined, which returns finite difference of order \(n\) that satisfies the recurrent relation diff(x,n) = diff(x,n-1). If x is one-dimensional array  
x = [x(1) x(2)... x(n)],  
Figure 5: Vectorized flux computation (arrow shows direction of flux). Figure 6: Example of the link between fluxes and vector structure. Figure 7: Calculation of fluxes through cell faces (arrows show direction of flux). then diff(x) is a vector of finite differences between neighboring elements  
diff(x) = [x(2)-x(1) x(3)-x(2)... x(n)-x(n-1)]. The number of elements of vector x is one less than the number of elements of vector diff(x). The derivative approximation of order \(n\) is diff(y,n)./diff(x,n), where \"./\" operation denotes element-wise division of vectors. The derivative is calculated in another way, using the method of addressing mesh data structures. The integration interval \(x_{0}\leqslant x\leqslant x_{N}\) is divided into \(N\) equal segments, and denote the length of each of them as \(\Delta x=x_{i}-x_{i-1}\) (Figure 8a). A number of all intervals varies from \(1\) to \(N-1\). An array of indices \(\mathbf{I}\), in which the sequential increase in indices corresponds to the order of its elements, is  
\[\mathbf{I}=\{1,2,\ldots,P-1,P,P+1,\ldots,N-1,N\}^{\prime}=1:N.\]  
The finite-difference template is shown in Figure 8b, where \(P\) is the midpoint of each segment, and the left and right boundaries of the interval are \(L\) and \(R\), respectively (there are \(N-1\) points \(P\), \(L\), and \(R\) in total). The arrays of indices for the left and right elements of mesh structure are introduced  
\[\mathbf{L} = \{1,2,\ldots,P-1,P,P+1,\ldots,N-2,N-1\}^{\prime}=1:N-1;\] \[\mathbf{R} = \{2,3,\ldots,P-1,P,P+1,\ldots,N-1,N\}^{\prime}=2:N.\]  
The vector subtraction operation is defined as  
\[\Delta\mathbf{U}=\mathbf{U}(\mathbf{R})-\mathbf{U}(\mathbf{L}),\]  
where  
\[ \begin{pmatrix}\Delta\mathbf{U}(1)\\ \Delta\mathbf{U}(2)\\ \vdots\\ \Delta\mathbf{U}(N-1)\end{pmatrix} = \begin{pmatrix}\mathbf{U}(\mathbf{R}(1))- \mathbf{U}(\mathbf{L}(1))\\ \mathbf{U}(\mathbf{R}(2))-\mathbf{U}(\mathbf{L}(2))\\ \vdots\\ \mathbf{U}(\mathbf{R}(N-1))-\mathbf{U}(\mathbf{L}(N-1))\end{pmatrix} .\]  
Then, the formula for vectorized calculation of derivative has the form  
\[\frac{d\mathbf{U}(x)}{dx}=\frac{\mathbf{U}(\mathbf{R})-\mathbf{U}(\mathbf{L}) }{\Delta x}+O(\Delta x^{2}).\]  
In a pseudo programming language that supports vectorized operations, the algorithm of derivative calculation is represented as the following sequence of steps. % indexes of intervals and boundary nodes iI=1:N-1; iL=iI; iR=iI+1; % mesh step dx=1/(N-1); % derivative y'(x) y=(y(iR)-y(iL))/dx;  
Figure 8: Partition of the integration interval in \(x\) direction (**a**) and finite-difference template (**b**). The presented method of addressing mesh data structures underlies the calculation of derivatives of mesh functions and the construction of vectorized algorithms for solving boundary value problems. For the left and right boundaries, index vectors \(\mathbf{I}_{L}\) and \(\mathbf{I}_{R}\) are introduced, having dimension \(M\), whose components are equal to \(1\) if the left/right boundary condition is specified for the component \(k\) of the solution vector, and \(0\) if there is no left/right boundary condition for the \(k\) component of the solution vector. Two vectors \(\mathbf{L}\) and \(\mathbf{R}\) are formed. Their components contain, at the places determined by the index vectors \(\mathbf{I}_{L}\) and \(\mathbf{I}_{R}\), left and right boundary values, respectively, and the remaining components have undefined (arbitrary) values. Arbitrarily specified components of the vectors \(\mathbf{L}\) and \(\mathbf{R}\) do not take part in the calculations.