17.6. Vectorized double numerical differentiation

Introduction

Section 17.5 demonstrated that both the strain–displacement matrix $\boldsymbol B$ and the tangent stiffness matrix can be obtained by numerical differentiation. Its direct implementation, however, is very expensive because it repeats the same sequence of calculations for every element, Gauss point and perturbation.

Section 17.6 retains, in principle, the same double numerical differentiation algorithm and reorganizes the calculations through vectorization and batched array operations. Central differences, the constitutive law, the multiplicative update of the deformation gradient and the Newton procedure are preserved. Both $\boldsymbol B$ and the tangent are still differentiated numerically; no symmetrization of the tangent is introduced. The difference is the computational organization, not a new finite-element formulation.

The vectorized implementation presented in this section was developed with the assistance of Chiara and GPT-6 Astra.

Computational organization

A Q4 element has eight nodal displacement degrees of freedom. The element displacement vector, the three in-plane engineering strain components, the strain–displacement matrix, the element internal force vector and the tangent stiffness matrix have the following dimensions:

\boldsymbol u_{\rm el}\in\mathbb R^{8\times1},\qquad \boldsymbol e\in\mathbb R^{3\times1},\qquad \boldsymbol B\in\mathbb R^{3\times8}

\boldsymbol f_{\rm el}\in\mathbb R^{8\times1},\qquad \boldsymbol k_{\rm el}\in\mathbb R^{8\times8}

Each column of $\boldsymbol B$ is calculated by central differences, exactly as in Section 17.5. For one Gauss point and one element state, the calculation includes one unperturbed strain-tensor evaluation, eight positive perturbations and eight negative perturbations:

1+2\times8=17

Each evaluation produces three strain components. Thus, 17 complete strain-tensor evaluations correspond to 51 scalar strain-component evaluations:

17\times3=51

The first level of grouping concerns the numerical tangent, which also requires 17 evaluations of the element internal force vector. Instead of calculating these outer displacement states successively, stiff.m uses repmat to construct a three-dimensional array of size $8\times n_{\rm el}\times17$ . Indexing then introduces the positive and negative ddK displacement perturbations: one state is unperturbed, eight have a positive perturbation and eight have a negative perturbation. All seventeen states are passed together to fel_el. They read the same preceding load-step deformation-gradient history; only the unperturbed evaluation is stored.

For the full mesh used in the test (the example from Section 17.3):

n_H=8,\qquad n_L=60,\qquad n_{\rm el}=480

The four Gauss points per element give 1920 Gauss points. With seventeen outer displacement states, each batch therefore contains:

4\times480=1920,\qquad17\times1920=32\,640

deformation-gradient states. These are not processed through 32 640 successive scalar calls: they are held in arrays and treated by page-wise and element-wise array operations.

The second level of grouping concerns the numerical evaluation of $\boldsymbol B$ . In fel_el.m, the element, Gauss-point and outer-perturbation indices are retained together in arrays; the Gauss-point and element indices are combined into a single batch dimension. A $3\times8$ matrix $\boldsymbol B$ is needed for each of the 32 640 outer states. For each of its eight columns, the positive strain perturbation is evaluated for all 480 elements, all four Gauss points and all seventeen outer states together; the negative perturbation is evaluated over the same complete batch.

Including the unperturbed strain evaluation, this gives 17 strain-tensor evaluations for each outer state. The equivalent count per element and Newton iteration is:

17\times17\times4=1156

For the complete mesh, the two counts must be distinguished:

17\times17\times4\times480=554\,880\quad\text{complete strain-tensor evaluations}

3\times554\,880=1\,664\,640\quad\text{scalar strain-component evaluations}

Vectorization does not reduce this mathematical count. It changes how the calculations are grouped and executed.

The two numerical differentiations have different roles. For a Q4 element, the first obtains the strain–displacement matrix from the strain vector:

\boldsymbol e\xrightarrow{\partial/\partial\boldsymbol u_{\rm el}}\boldsymbol B

\boldsymbol e=\begin{Bmatrix}e_x\\e_y\\\gamma_{xy}\end{Bmatrix},\qquad\boldsymbol B\in\mathbb R^{3\times8}.

Each column of $\boldsymbol B$ is obtained by central finite differences. At each Gauss point, the calculation includes $1+2\times8=17$ strain-vector evaluations: one unperturbed, eight with positive perturbations and eight with negative perturbations. The second differentiation obtains the element tangent stiffness matrix from the element internal force vector:

\boldsymbol f_{\rm el}\xrightarrow{\partial/\partial\boldsymbol u_{\rm el}}\boldsymbol k_{\rm el},\qquad\boldsymbol f_{\rm el}\in\mathbb R^{8\times1},\qquad\boldsymbol k_{\rm el}\in\mathbb R^{8\times8}.

This calculation again requires $1+2\times8=17$ internal-force evaluations: one unperturbed, eight with positive perturbations and eight with negative perturbations. The arrows represent numerical differentiation.

In the product below, the first factor 17 counts the outer states needed to differentiate $\boldsymbol f_{\rm el}$ and obtain $\boldsymbol k_{\rm el}$ . The second counts the inner strain evaluations needed to obtain $\boldsymbol B$ for each outer state. Thus, for each element and Gauss point:

17\times17=289=1+16+16+256.

The four terms count one evaluation with neither perturbation, sixteen with only the inner perturbation $\Delta_B$ , sixteen with only the outer perturbation $\Delta_K$ , and 256 with both. Most strain evaluations therefore contain two superposed perturbations: one for the numerical evaluation of $\boldsymbol B$ , the other for the numerical evaluation of $\boldsymbol k_{\rm el}$ .

For 480 elements and four Gauss points per element, the complete count per Newton iteration is:

17\times17\times4\times480=554\,880

complete strain-vector evaluations. Since each vector contains three components, this corresponds to:

3\times554\,880=1\,664\,640

scalar strain-component evaluations.

The program does not execute all 554 880 evaluations in a single operation. In fel_el.m, they are organized into seventeen calls to strain_el: one without an inner perturbation, eight with positive inner perturbations and eight with negative inner perturbations. Each call includes all seventeen outer displacement states, all four Gauss points and all 480 elements, so each batch contains:

17\times4\times480=32\,640

strain-vector evaluations. The seventeen batches therefore correspond to the inner differentiation, while the seventeen outer states are present together in every batch. Equivalently, each outer state receives one unperturbed and sixteen inner-perturbed strain evaluations at every Gauss point.

The 554 880 strain-vector evaluations supply 32 640 matrices $\boldsymbol B$ , one for each element, Gauss point and outer state. Integration over the four Gauss points produces $32\,640/4=8\,160$ element internal force vectors $\boldsymbol f_{\rm el}$ . Finally, the seventeen force vectors for each element provide $8\,160/17=480$ element tangent stiffness matrices $\boldsymbol k_{\rm el}$ .

The main array operations follow this organization. pagemtimes applies matrix multiplication to many matrix pages at once, while pagemldivide solves the families of small systems associated with the Jacobians. Implicit array expansion applies common geometry and perturbation data across whole batches of elements, Gauss points and outer states, without a separate scalar operation sequence for each state.

The strain_el.m routine receives arrays of deformation gradients and evaluates Hencky or Biot strain for the entire batch. Repeated scalar eig calls are avoided: the spectral functions of the symmetric $2\times2$ matrices are evaluated using algebraically equivalent expressions applied directly to arrays.

The global stiffness matrix is also assembled in a fast sparse form, following the same general philosophy introduced in Section 9.5. Vectorized indexing supplies the global row and column indices, and sparse assembles the element tangent contributions. No parallel workers, GPU or Parallel Computing Toolbox are used.

Thus, the principal acceleration is obtained not by reducing the number of finite-difference evaluations, but by changing their computational organization: operations that were previously executed one after another are now performed on large arrays of element, Gauss-point and perturbation data.

Numerical comparison – reduced mesh

Reduced mesh: $5\times30$ , 186 nodes, 150 Q4 elements.

Load steps	Method	$ddB=ddK$	Total Newton iterations	$u$ [mm]	$v$ [mm]	$t$ [s]
5	17.4	—	54	142.8025	276.1211	1.78
5	17.5	$10^{-5}$	43	142.8081	276.1717	152.17
5	17.6	$10^{-5}$	43	142.8081	276.1717	0.915
10	17.4	—	104	142.8000	276.1126	3.49
10	17.5	$10^{-5}$	72	142.8054	276.1631	274.72
10	17.6	$10^{-5}$	71	142.8054	276.1631	1.375
20	17.4	—	197	142.7986	276.1086	6.080
20	17.5	$10^{-5}$	120	142.8039	276.1590	556.35
20	17.6	$10^{-5}$	112	142.8039	276.1590	2.299

The ANSYS reference displacements are:

u_{\mathrm{ANSYS}}=142.795\ \mathrm{mm},\qquad v_{\mathrm{ANSYS}}=276.148\ \mathrm{mm}

For this reduced mesh, Section 17.6 practically reproduces the displacement results of Section 17.5. With five load steps, the number of Newton iterations is also identical. With ten and twenty load steps, small differences in the total iteration count occur without significant changes in the reported displacements. These iteration differences should not be over-interpreted.

The computation time is reduced by approximately $166\times$ , $200\times$ , and $242\times$ for five, ten, and twenty load steps, respectively. The main benefit is the much lower cost of each Newton iteration.

Numerical comparison – full mesh

Full mesh: $8\times60$ , 549 nodes, 480 Q4 elements.

Method	Load steps	Newton iterations	$u$ [mm]	$v$ [mm]	CPU time [s]	Speed-up vs 17.5
17.4 analytical – `stiff`	20	189	155.8090	286.5564	15.9168	—
17.4 analytical – `stiff_mult`	20	207	155.6850	286.3191	19.0258	—
17.5 direct numerical	20	118	155.8118	286.5712	1207.9648	1.0
17.6 vectorized numerical	20	110	155.6870	286.3331	5.0377	239.8

The ANSYS reference displacements are:

u_{\mathrm{ANSYS}}=155.676\ \mathrm{mm},\qquad v_{\mathrm{ANSYS}}=286.321\ \mathrm{mm}

In this test, the vectorized numerical implementation is approximately $239.8$ times faster than the direct numerical implementation. It is also faster than the two analytical implementations of Section 17.4:

\frac{15.9168}{5.0377}\approx3.16\qquad\text{for stiff}

\frac{19.0258}{5.0377}\approx3.78\qquad\text{for stiff\_mult}

These ratios apply to the test, mesh and implementations considered here. They are not universal properties of the methods.

Conclusions

The vectorized implementation retains the numerical formulation of Section 17.5: the strain–displacement matrix and the tangent stiffness matrix are both obtained by central numerical differentiation.
The reduction in computation time is dramatic, without changing the constitutive model, perturbations or Newton algorithm.
The very high cost observed in Section 17.5 is therefore not an unavoidable consequence of double numerical differentiation; a substantial part of it results from the sequential implementation of a highly repetitive algorithm.

Section 17.5 concluded: “Greater analytical simplicity can be purchased at a very high numerical cost.” The present implementation shows that this cost can be reduced radically by a more efficient computational organization.

For the present example, careful vectorization makes the fully numerical implementation not only practical, but even faster than the analytical implementations considered in Section 17.4.

Download

The Q4 vectorized double numerical differentiation program can be downloaded below. Run main.m; the Hencky/Biot selector and the perturbations ddB=1e-5 and ddK=1e-5 are set in stiff.m.

MATLABQ4 vectorized double numerical differentiationEight MATLAB files and READMEDOWNLOAD ZIP