*2.2. Description of Structure Discretization*

The analyzed cross-sectional area of the MEMS structure was discretized using the mesh that can be described by the following equations [16]:

$$\mathfrak{q}\_k(t) = \mathfrak{q}(\mathbf{x}, \mathbf{y}, t), \mathbf{x} = \mathbf{i} \cdot \Delta \mathbf{x}, \mathbf{y} = \mathbf{j} \cdot \Delta \mathbf{y} \tag{7}$$

$$T\_k(t) = T(\mathbf{x}, y, t), \mathbf{x} = \mathbf{i} \cdot \Delta \mathbf{x}, y = j \cdot \Delta y \tag{8}$$

$$i \in \{1, 2, \ldots, n\_{\overline{x}}\}, j \in \{1, 2, \ldots, n\_{\overline{y}}\}, k \in \{1, 2, \ldots, n\_{\overline{x}} \cdot n\_{\overline{y}}\}$$

In the Equations above, *nx* and *ny* reflect the number of discretization nodes in both axes. It means that the product *nx* × *ny* describes the number of nodes used in discretization mesh. Nodes are numbered from the left side to the right side. First, the first bottom row is considered. When all nodes in this row have already had their numbers, the procedure is repeated in the second row from the bottom side. This process is continued until all nodes in a top row are numbered. The graphical interpretation of the applied discretization mesh for the investigated cross-section of the MEMS structure, as well as nodes numbering approach, are demonstrated in Figure 2.

**Figure 2.** The graphical interpretation of the proposed discretization mesh nodes numbering [16,21].

### *2.3. DPL Model Approximation Scheme for the Investigated MEMS Structure*

Considering the imposed FDM assumptions, the DPL model can be expressed in the following way [22]:

$$\begin{cases} \mathcal{M} \cdot \ddot{T}(t) + \mathcal{D} \cdot \dot{T}(t) + \mathcal{K} \cdot T(t) = \mathcal{b} \cdot u(t) \\ \qquad y(t) = \mathbf{c}^T \cdot T(t) \end{cases} \quad t \in \mathbb{R}\_+ \cup \ [0\} \tag{9}$$

where superscript *T* means a transposition. Moreover, *T* is the temperature function, while *T* and *T* are its first and second time derivative, respectively. All of them are *nx* · *ny* × 1 vectors. Other vectors and matrices can be described as follows [23]:

$$\begin{aligned} &I, M\_{\mathrm{FDM}}, M, D, K, \mathfrak{c}^T \in R^{n\_{\mathcal{X}} \cdot n\_{\mathcal{Y}} \times n\_{\mathcal{X}} \times n\_{\mathcal{Y}}}, \\ &\mathfrak{c}\_{\mathcal{V}}, k, \mathfrak{r}\_{\mathfrak{q}}, \mathfrak{r}\_{\mathcal{T}}, b, y \in R^{n\_{\mathcal{X}} \cdot n\_{\mathcal{Y}} \times 1}, \\ &u \in R \end{aligned}$$

$$\mathcal{M} = \operatorname{diag} \left( \mathfrak{r}\_q \right) \circ \operatorname{diag} \left( \mathfrak{c}\_v \right) \circ I \tag{10}$$

$$D = \operatorname{diag} \left( \mathbf{c}\_{\overline{\mathbf{v}}} \right) \circ \mathbf{I} - \frac{1}{\left( \Delta \mathbf{x} \right)^{2}} \cdot \operatorname{repmat} \left( \mathbf{r}\_{\overline{\mathbf{r}}} \right) \circ \operatorname{repmat} \left( \mathbf{k} \right) \circ \mathbf{M}\_{\mathrm{FDM}} \tag{11}$$

$$\mathbf{K} = -\frac{1}{\left(\Delta x\right)^2} \text{repmat}\left(\mathbf{k}\right) \circ \mathbf{M}\_{\text{FDM}}\tag{12}$$

$$b = a \cdot \begin{bmatrix} 0 & \cdots & 0 & 1 & \cdots & 1 & 0 & \cdots & 0 \end{bmatrix}^T, a \in \mathbb{R}\_+ \tag{13}$$

$$\mathbf{c} = I \tag{14}$$

.

..

where *I* is identity matrix, an operator ◦ means a Hadamard product, function diag(·) generates a diagonal matrix based on a given vector, and function repmat(·) replicates a considered vector and generates a matrix of desired dimensions. In addition, non-zero elements of vector *b* indicate nodes with internal heat generation. Moreover, *MFDM* is a matrix including FDM coefficients wgugc us determined based on Equation (6) and the considered boundary conditions. It is also worthwhile to note that the *MFDM* coefficients describing the air area between platinum resistors' surfaces, as well as the contact between them and the oxide, have been calculated considering a fractional order of the temperature function space derivative based on the Grünwald–Letnikov theory [22,24]. Thus, in this particular case, each one-dimensional part of Equation (6) is replaced by the following Equation [22]:

1 (<sup>Δ</sup>*s*)<sup>α</sup>*<sup>x</sup>* · ) 2 *k*=0 (−<sup>1</sup>)*<sup>k</sup>* <sup>Γ</sup>(<sup>α</sup>*x*+<sup>1</sup>) <sup>Γ</sup>(*k*+<sup>1</sup>)·<sup>Γ</sup>(<sup>α</sup>*x*<sup>−</sup>*k*+<sup>1</sup>)*<sup>T</sup><sup>x</sup>* − *k* · Δ*x* + <sup>α</sup>*x*·Δ*<sup>x</sup>* 2 , *t* · 1 = = ( α*x*2 <sup>−</sup><sup>1</sup>)·*<sup>T</sup>*(*x*+2·Δ*x*,*y*,*<sup>t</sup>*)+(<sup>2</sup>−α*x*<sup>2</sup> <sup>−</sup>α*x*·( α*x*2 <sup>−</sup><sup>1</sup>))·*<sup>T</sup>*(*x*+Δ*x*,*y*,*<sup>t</sup>*) (<sup>Δ</sup>*x*)<sup>α</sup>*<sup>x</sup>* + + <sup>α</sup>*x*·(<sup>α</sup>*x*−<sup>1</sup>) 2 ·( α*x*2 <sup>−</sup><sup>1</sup>)−α*x*·(<sup>2</sup>−α*x*2 ) ·*<sup>T</sup>*(*<sup>x</sup>*,*y*,*<sup>t</sup>*)+(<sup>2</sup>−α*x*<sup>2</sup> )· <sup>α</sup>*x*·(<sup>α</sup>*x*−<sup>1</sup>) 2 ·*<sup>T</sup>*(*<sup>x</sup>*−Δ*x*,*y*,*<sup>t</sup>*) (<sup>Δ</sup>*x*)<sup>α</sup>*<sup>x</sup>* +( α*x*2 <sup>−</sup><sup>1</sup>)·*<sup>T</sup>*(*<sup>x</sup>*,*y*+2·Δ*x*,*<sup>t</sup>*)+(<sup>2</sup>−α*x*<sup>2</sup> <sup>−</sup>α*x*·( α*x*2 <sup>−</sup><sup>1</sup>))·*<sup>T</sup>*(*<sup>x</sup>*,*y*+Δ*x*,*<sup>t</sup>*) (<sup>Δ</sup>*x*)<sup>α</sup>*<sup>x</sup>* + + <sup>α</sup>*x*·(<sup>α</sup>*x*−<sup>1</sup>) 2 ·( α*x*2 <sup>−</sup><sup>1</sup>)−α*x*·(<sup>2</sup>−α*x*2 ) ·*<sup>T</sup>*(*<sup>x</sup>*,*y*,*<sup>t</sup>*)+(<sup>2</sup>−α*x*<sup>2</sup> )· <sup>α</sup>*x*·(<sup>α</sup>*x*−<sup>1</sup>) 2 ·*<sup>T</sup>*(*<sup>x</sup>*,*y*−Δ*x*,*<sup>t</sup>*) (<sup>Δ</sup>*x*)<sup>α</sup>*<sup>x</sup>* for α*x* ∈ (2, 2.5), Δ*x* → 0 (15)

In order to make the analysis more clear, the system Equations (9) of the second-order equations has been transformed into the first-order Equations [21]:

$$\begin{cases} \begin{aligned} \dot{E} \cdot \dot{\overline{T}}(t) &= \mathbf{A} \cdot \overline{T}(t) + \mathbf{B} \cdot u(t) \\ y(t) &= \mathbf{C}^T \cdot \overline{T}(t) \end{aligned} \quad t \in \mathbb{R}\_+ \cup \ \{0\} \end{aligned} \tag{16}$$

·

where − *T* and *T* are 2 · *nx* · *ny* × 1 vectors. First *nx* · *ny* elements of − *T* are the same, similar to the case of *T* vector, while its second part coincides with · *T*. The first half of · *T* includes elements of . *T*, whereas its second half states .. *T*. In addition, the remaining matrices and vectors visible in Equation (16) present as follows [21,23,25]:

$$E = \begin{bmatrix} I & \Theta \\ \Theta & M \end{bmatrix} A = \begin{bmatrix} \Theta & I \\ -\mathcal{K} & -\mathcal{D} \end{bmatrix} \mathcal{B} = \begin{bmatrix} \Theta\_1 \\ b \end{bmatrix}, \mathcal{C} = \begin{bmatrix} \mathcal{C} & \Theta \\ \Theta & \mathcal{c} \end{bmatrix} \tag{17}$$

where Θ and Θ1 are null matrices. The final solution is determined using the BDF method with a variable order between 1 and 5.
