Numerical Analysis and Stability of the Moore-Gibson-Thompson-Fourier Model

Abstract.

This work is concerned the Moore-Gibson-Thompson-Fourier Model. Our contribution will consist in studying the numerical stability of the Moore-Gibson-Thompson-Fourier system. First we introduce a finite element approximation after the discretization, then we prove that the associated discrete energy decreases and later we establish a priori error estimates. Finally, we obtain some numerical simulations.

1. Introduction

In this paper, we consider a Moore-Gibson-Thompson (MGT) equation

(1) $$u_{ttt}+\alpha u_{tt}+\beta A{u}_t+\gamma Au=0,$$

which describes the evolution of the unknown function $u=u(x,t):\mathrm{\Omega }\times [0,\mathrm{\infty })\to ℝ$, where $\mathrm{\Omega }\subset {ℝ}^N$ is a bounded domain with a sufficiently smooth boundary $\partial \mathrm{\Omega }$. The equation includes various parameters such as $\alpha ,\beta ,\gamma >0$, which are fixed structural parameters.

Originally, for the Laplace-Dirichlet operator $A=-\mathrm{\Delta }$ this equation was introduced to model wave propagation in viscous thermally relaxing fluids [14, 17], with its first appearance dating back to a paper by Stokes [15]. Over time, researchers have discovered that the MGT equation finds applications in a wide range of physical phenomena, including viscoelasticity and thermal conduction. Notably, it has been interpreted as a model for vibrations in a standard linear viscoelastic solid [9, 10].

In the particular case where $A={\mathrm{\Delta }}^2$, with proper boundary conditions (see [12]), the MGT equation appears as a possible model for the vertical displacement in viscoelastic plates.

The mathematical analysis of the MGT equation has attracted significant attention, resulting in a vast literature with numerous studies and references available [4, 5, 11, 13, 2, 3, 16]. The main findings can be summarized as follows:

For any positive values of the parameters $\alpha ,\beta ,\gamma$, the MGT equation generates a strongly continuous semigroup of solutions. However, the behavior of these solutions depends significantly on the constant $\nu$, defined as follows:

$$\nu =\alpha \beta -\gamma .$$

For $A={\mathrm{\Delta }}^2$ the semigroup of MGT equation is analytic and exponentially stable the case of $\nu >0$.

In present paper, we consider the MGT-Fourier system

(2) $$\{\begin{array}{c}{u}_{ttt}+\alpha u_{tt}+\beta {\mathrm{\Delta }}^2{u}_t+\gamma {\mathrm{\Delta }}^{2}u=-\eta \mathrm{\Delta }\theta ,\hfill \\ {\theta }_t-\kappa \mathrm{\Delta }\theta =\eta \mathrm{\Delta }{u}_{tt}+\alpha \eta \mathrm{\Delta }{u}_t\hfill \end{array}$$

where, the unknown function $u=u(x,t)$ represents the vibration of flexible structures and $\theta =\theta (x,t)$ the difference of temperature between the actual state and a reference temperature where $x\in \mathrm{\Omega },t\in (0,\mathrm{\infty })$. The parameters $\alpha ,\beta ,\gamma$ and $\kappa$ are positive real numbers, $\eta \ne 0$ and $\mathrm{\Omega }=[0,1]$ is a bounded domain. We assume the initial conditions

(3) $$u(x,0)=u_0,u_t(x,0)=u_1,u_{tt}(x,0)=u_2,\theta (x,0)={\theta }_0$$

where $u_0,u_1,u_2,{\theta }_0:\mathrm{\Omega }\to ℝ$ are assigned initial data. The system is complemented with the boundary conditions

(4) $u(x,t)=\frac{\partial u}{\partial \nu }(x,t)=\theta (x,t)=0,x\in \partial \mathrm{\Omega }$

Now, we intoduce a new variable $z=u_t+\alpha u$ and using $\nu =\alpha \beta -\gamma$. Consequently, the system $(\text{<a href="\#S1.E2" title="Equation 2 ‣ 1. Introduction ‣ Numerical Analysis and Stability of the Moore-Gibson-Thompson-Fourier Model" class="ltx\_ref ltx\_font\_italic"><span class="ltx\_text ltx\_ref\_tag">2</span></a>})$ is equivalent to

(5) $$\{\begin{array}{c}{z}_{tt}+\frac{\gamma }{\alpha }{\mathrm{\Delta }}^{2}z+\frac{\nu }{\alpha }{\mathrm{\Delta }}^2{u}_t+\eta \mathrm{\Delta }\theta =0,\hfill \\ {\theta }_t-\kappa \mathrm{\Delta }\theta -\eta \mathrm{\Delta }{z}_t=0.\hfill \end{array}$$

Associated to (2)–(4), we consider the energy functional

(6) $E(t)$ $=\frac{1}{2}\left({\Vert z_t\Vert }^2+\frac{\gamma }{\alpha }{\Vert \mathrm{\Delta }z\Vert }^2+\frac{\nu }{\alpha }{\Vert \mathrm{\Delta }{u}_t\Vert }^2+{\Vert \theta \Vert }^2\right)$

$=\frac{1}{2}\left({\Vert u_{tt}+\alpha u_t\Vert }^2+\frac{\gamma }{\alpha }{\Vert \mathrm{\Delta }{u}_t+\alpha \mathrm{\Delta }u\Vert }^2+\frac{\nu }{\alpha }{\Vert \mathrm{\Delta }{u}_t\Vert }^2+{\Vert \theta \Vert }^2\right)$

Theorem 1 ([8]).

The semigroup associate to $(\text{<a href="\#S1.E2" title="Equation 2 ‣ 1. Introduction ‣ Numerical Analysis and Stability of the Moore-Gibson-Thompson-Fourier Model" class="ltx\_ref"><span class="ltx\_text ltx\_ref\_tag">2</span></a>})$ is analytic and exponentially stable for $\nu >0$.

As a results from [8] the energy $(\text{<a href="\#S3.EGx1" title="1. Introduction ‣ Numerical Analysis and Stability of the Moore-Gibson-Thompson-Fourier Model" class="ltx\_ref ltx\_font\_italic"><span class="ltx\_text ltx\_ref\_tag">1</span></a>})$ decays exponentially for $\nu >0$, that is, there exist two positive constants ${ϵ}_1,{ϵ}_2$ such that

$$E(t)\le {ϵ}_1{e}^{-{ϵ}_{2}t};\text{for all }t\ge 0.$$

and satisfies

(7) $$E^{\prime }(t)=-\nu {\Vert \mathrm{\Delta }{u}_t\Vert }^2-\kappa {\Vert \nabla \theta \Vert }^2\le 0,$$

For further details, refer to [8].

2. Numerical approximation

In this section, we propose a finite element approximation to system $(\text{<a href="\#S1.E2" title="Equation 2 ‣ 1. Introduction ‣ Numerical Analysis and Stability of the Moore-Gibson-Thompson-Fourier Model" class="ltx\_ref ltx\_font\_italic"><span class="ltx\_text ltx\_ref\_tag">2</span></a>})$ with boundary conditions $(\text{<a href="\#S1.E4" title="Equation 4 ‣ 1. Introduction ‣ Numerical Analysis and Stability of the Moore-Gibson-Thompson-Fourier Model" class="ltx\_ref ltx\_font\_italic"><span class="ltx\_text ltx\_ref\_tag">4</span></a>})$ and initial conditions $(\text{<a href="\#S1.E3" title="Equation 3 ‣ 1. Introduction ‣ Numerical Analysis and Stability of the Moore-Gibson-Thompson-Fourier Model" class="ltx\_ref ltx\_font\_italic"><span class="ltx\_text ltx\_ref\_tag">3</span></a>})$.

We introduce and study finite elements in space and an implicit Euler type scheme based on finite differences in time. We prove that the discrete energy decays.

Introducing new variables $y=z_t$,$v=u_t$,$\mathrm{\Phi }=-\mathrm{\Delta }z$ and $\mathrm{\Psi }=-\mathrm{\Delta }v$; we rewrite system $(\text{<a href="\#S1.E5" title="Equation 5 ‣ 1. Introduction ‣ Numerical Analysis and Stability of the Moore-Gibson-Thompson-Fourier Model" class="ltx\_ref ltx\_font\_italic"><span class="ltx\_text ltx\_ref\_tag">5</span></a>})$

(8)

In order to obtain the weak form associated with system $(\text{<a href="\#S2.E8" title="Equation 8 ‣ 2. Numerical approximation ‣ Numerical Analysis and Stability of the Moore-Gibson-Thompson-Fourier Model" class="ltx\_ref ltx\_font\_italic"><span class="ltx\_text ltx\_ref\_tag">8</span></a>})$, we multiply the equations by test functions $\chi ,\xi ,\omega ,\zeta \in H_0^1(0,1)$ and integrate by parts.

(9)

For our purposes, we considered $J$ a nonnegative integer and $h=\frac{1}{J}$a subdivision of the interval $(0,1)$ given by , such that $x_j=jh$, for all $j=0,\mathrm{\dots },J$. We take

(10) $S^h=\{g\in H^1(0,1)|g\in C\left([0,1]\right),$ ${g|}_{(x_j,x_{j+1})}\text{is a linear polynomial, with}$

$j=0,\mathrm{\dots },J-1\}$

and

$$S_0^h=\left\{g\in S^h\right|g(0)=g(1)=0\}.$$

For a given final time $T$ and a positive integer $N$, let $\mathrm{△}t=T/N$ be the time step and $t_n=n\mathrm{\Delta }t,n=0,\mathrm{\dots },N.$

The finite element method for $(\text{<a href="\#S2.E9" title="Equation 9 ‣ 2. Numerical approximation ‣ Numerical Analysis and Stability of the Moore-Gibson-Thompson-Fourier Model" class="ltx\_ref ltx\_font\_italic"><span class="ltx\_text ltx\_ref\_tag">9</span></a>})$ using the backward Euler scheme is to find $y_h^n,{\theta }_h^n,{\mathrm{\Phi }}_h^n,{\mathrm{\Psi }}_h^n\in S_0^h$ such that, for $n=1,\mathrm{\dots },N$ and for all ${\chi }_h,{\xi }_h,{\omega }_h,{\zeta }_h\in S_0^h$

(11)

where

(12) $$v_h^n=\frac{u_h^n-u_h^{n-1}}{\mathrm{△}t},z_h^n=v_h^n+\alpha u_h^n,\text{and}{y}_h^n=\frac{z_h^n-z_h^{n-1}}{\mathrm{△}t},$$

are approximations to $u_t\left(t_n\right),v\left(t_n\right)+\alpha u\left(t_n\right),z_t\left(t_n\right)$ respectively.

By leveraging the properties of inner products and norms, we derive the following identity, which will be frequently used:

(13) $$(a-b,a)=\frac{1}{2}\left({\Vert a-b\Vert }^2+{\Vert a\Vert }^2-{\Vert b\Vert }^2\right).$$

The next result is a discrete version of the energy decay property satisfied by the solution of system $(\text{<a href="\#S1.E2" title="Equation 2 ‣ 1. Introduction ‣ Numerical Analysis and Stability of the Moore-Gibson-Thompson-Fourier Model" class="ltx\_ref ltx\_font\_italic"><span class="ltx\_text ltx\_ref\_tag">2</span></a>})$.

We introduce the following discrete energy,

(14) ${ℰ}_h^n$ $=\frac{1}{2}\left({\Vert y_h^n\Vert }^2+\frac{\gamma }{\alpha }{\Vert {\mathrm{\Phi }}_h^n\Vert }^2+\frac{\nu }{\alpha }{\Vert {\mathrm{\Psi }}_h^n\Vert }^2+{\Vert {\theta }_h^n\Vert }^2\right).$

Theorem 2.

The discrete energy decay to zero, that is,

(15) $$\frac{{ℰ}_h^n-{ℰ}_h^{n-1}}{\mathrm{△}t}\le 0,$$

holds for $n=1,2,\mathrm{\dots },N$.

Proof.

Taking ${\chi }_h=y_h^n$ and ${\xi }_h={\theta }_h^n$ in $(\text{<a href="\#S2.E11" title="Equation 11 ‣ 2. Numerical approximation ‣ Numerical Analysis and Stability of the Moore-Gibson-Thompson-Fourier Model" class="ltx\_ref ltx\_font\_italic"><span class="ltx\_text ltx\_ref\_tag">11</span></a>})$.

(16)

Summing equations of system $(\text{<a href="\#S2.E16" title="Equation 16 ‣ 2. Numerical approximation ‣ Numerical Analysis and Stability of the Moore-Gibson-Thompson-Fourier Model" class="ltx\_ref ltx\_font\_italic"><span class="ltx\_text ltx\_ref\_tag">16</span></a>})$, we have

(17) $\frac{1}{\mathrm{\Delta }t}(y_h^n-y_h^{n-1},y_h^n)+\frac{\gamma }{\alpha }(\nabla {\mathrm{\Phi }}_h^n,\nabla y_h^n)+\frac{\nu }{\alpha }(\nabla {\mathrm{\Psi }}_h^n,\nabla y_h^n)$ $+\frac{1}{\mathrm{\Delta }t}({\theta }_h^n-{\theta }_h^{n-1},{\theta }_h^n)+$

$+\kappa (\nabla {\theta }_h^n,\nabla {\theta }_h^n)=0.$

Recalling $(\text{<a href="\#S2.E12" title="Equation 12 ‣ 2. Numerical approximation ‣ Numerical Analysis and Stability of the Moore-Gibson-Thompson-Fourier Model" class="ltx\_ref ltx\_font\_italic"><span class="ltx\_text ltx\_ref\_tag">12</span></a>})$ and $(\text{<a href="\#S2.E13" title="Equation 13 ‣ 2. Numerical approximation ‣ Numerical Analysis and Stability of the Moore-Gibson-Thompson-Fourier Model" class="ltx\_ref ltx\_font\_italic"><span class="ltx\_text ltx\_ref\_tag">13</span></a>})$, we have

(18) $$\frac{1}{\mathrm{\Delta }t}(y_h^n-y_h^{n-1},y_h^n)=\frac{1}{2\mathrm{\Delta }t}\left({\Vert y_h^n-y_h^{n-1}\Vert }^2+{\Vert y_h^n\Vert }^2-{\Vert y_h^{n-1}\Vert }^2\right).$$

Next,

(19) $\frac{\gamma }{\alpha }(\nabla {\mathrm{\Phi }}_h^n,\nabla y_h^n)$ $=-\frac{\gamma }{\alpha }({\mathrm{\Phi }}_h^n,\mathrm{\Delta }{y}_h^n)=$

$=-\frac{\gamma }{\alpha }({\mathrm{\Phi }}_h^n,\frac{\mathrm{\Delta }{z}_h^n-\mathrm{\Delta }{z}_h^{n-1}}{\mathrm{△}t})$

$=\frac{\gamma }{\alpha }({\mathrm{\Phi }}_h^n,\frac{{\mathrm{\Phi }}_h^n-{\mathrm{\Phi }}_h^{n-1}}{\mathrm{△}t})$

$=\frac{\gamma }{2\alpha \mathrm{\Delta }t}\left({\Vert {\mathrm{\Phi }}_h^n-{\mathrm{\Phi }}_h^{n-1}\Vert }^2+{\Vert {\mathrm{\Phi }}_h^n\Vert }^2-{\Vert {\mathrm{\Phi }}_h^{n-1}\Vert }^2\right).$

Similarly,

(20) $\frac{\nu }{\alpha }(\nabla {\mathrm{\Psi }}_h^n,\nabla y_h^n)$ $=-\frac{\nu }{\alpha }({\mathrm{\Psi }}_h^n,\mathrm{\Delta }{y}_h^n)$

$=-\frac{\nu }{\alpha }({\mathrm{\Psi }}_h^n,\frac{\mathrm{\Delta }{z}_h^n-\mathrm{\Delta }{z}_h^{n-1}}{\mathrm{△}t})$

$=-\frac{\nu }{\alpha }({\mathrm{\Psi }}_h^n,\frac{\mathrm{\Delta }(v_h^n+\alpha u_h^n)-\mathrm{\Delta }(v_h^{n-1}+\alpha u_h^{n-1})}{\mathrm{△}t})$

$=-\frac{\nu }{\alpha }({\mathrm{\Psi }}_h^n,\frac{\mathrm{\Delta }{v}_h^n-\mathrm{\Delta }{v}_h^{n-1}}{\mathrm{△}t})-\nu ({\mathrm{\Psi }}_h^n,\frac{\mathrm{\Delta }{u}_h^n-\mathrm{\Delta }{u}_h^{n-1}}{\mathrm{△}t})$

$=\frac{\nu }{\alpha }({\mathrm{\Psi }}_h^n,\frac{{\mathrm{\Psi }}_h^n-{\mathrm{\Psi }}_h^{n-1}}{\mathrm{△}t})+\nu ({\mathrm{\Psi }}_h^n,{\mathrm{\Psi }}_h^n)$

$=\frac{\nu }{2\alpha \mathrm{\Delta }t}\left({\Vert {\mathrm{\Psi }}_h^n-{\mathrm{\Psi }}_h^{n-1}\Vert }^2+{\Vert {\mathrm{\Psi }}_h^n\Vert }^2-{\Vert {\mathrm{\Psi }}_h^{n-1}\Vert }^2\right)+\nu {\Vert {\mathrm{\Psi }}_h^n\Vert }^2.$

Also,

(21) $$\frac{1}{\mathrm{\Delta }t}({\theta }_h^n-{\theta }_h^{n-1},{\theta }_h^n)=\frac{1}{2\mathrm{\Delta }t}\left({\Vert {\theta }_h^n-{\theta }_h^{n-1}\Vert }^2+{\Vert {\theta }_h^n\Vert }^2-{\Vert {\theta }_h^{n-1}\Vert }^2\right).$$

Thus,

$\frac{1}{2\mathrm{\Delta }t}\left({\Vert y_h^n-y_h^{n-1}\Vert }^2+{\Vert y_h^n\Vert }^2-{\Vert y_h^{n-1}\Vert }^2\right)+$

(22) $+\frac{\gamma }{2\alpha \mathrm{\Delta }t}\left({\Vert {\mathrm{\Phi }}_h^n-{\mathrm{\Phi }}_h^{n-1}\Vert }^2+{\Vert {\mathrm{\Phi }}_h^n\Vert }^2-{\Vert {\mathrm{\Phi }}_h^{n-1}\Vert }^2\right)$

$+\frac{\nu }{2\alpha \mathrm{\Delta }t}\left({\Vert {\mathrm{\Psi }}_h^n-{\mathrm{\Psi }}_h^{n-1}\Vert }^2+{\Vert {\mathrm{\Psi }}_h^n\Vert }^2-{\Vert {\mathrm{\Psi }}_h^{n-1}\Vert }^2\right)+\nu {\Vert {\mathrm{\Psi }}_h^n\Vert }^2$

$+\frac{1}{2\mathrm{\Delta }t}\left({\Vert {\theta }_h^n-{\theta }_h^{n-1}\Vert }^2+{\Vert {\theta }_h^n\Vert }^2-{\Vert {\theta }_h^{n-1}\Vert }^2\right)+\kappa {\Vert \nabla {\theta }_h^n\Vert }^2=0.$

We deduce that

$0$ $=\frac{1}{2\mathrm{\Delta }t}\left({\Vert y_h^n-y_h^{n-1}\Vert }^2+{\Vert y_h^n\Vert }^2-{\Vert y_h^{n-1}\Vert }^2\right)$

$+\frac{\gamma }{2\alpha \mathrm{\Delta }t}\left({\Vert {\mathrm{\Phi }}_h^n-{\mathrm{\Phi }}_h^{n-1}\Vert }^2+{\Vert {\mathrm{\Phi }}_h^n\Vert }^2-{\Vert {\mathrm{\Phi }}_h^{n-1}\Vert }^2\right)$

(23) $+\frac{\nu }{2\alpha \mathrm{\Delta }t}\left({\Vert {\mathrm{\Psi }}_h^n-{\mathrm{\Psi }}_h^{n-1}\Vert }^2+{\Vert {\mathrm{\Psi }}_h^n\Vert }^2-{\Vert {\mathrm{\Psi }}_h^{n-1}\Vert }^2\right)+\nu {\Vert {\mathrm{\Psi }}_h^n\Vert }^2$

$+\frac{1}{2\mathrm{\Delta }t}\left({\Vert {\theta }_h^n-{\theta }_h^{n-1}\Vert }^2+{\Vert {\theta }_h^n\Vert }^2-{\Vert {\theta }_h^{n-1}\Vert }^2\right)+\kappa {\Vert \nabla {\theta }_h^n\Vert }^2$

$\ge \frac{{ℰ}_h^n-{ℰ}_h^{n-1}}{\mathrm{△}t}.$

Which implies $\frac{{ℰ}_h^n-{ℰ}_h^{n-1}}{\mathrm{△}t}\le 0$ and this completes the proof. ∎

Now, we prove a main error estimates result.

Theorem 3.

There exists a positive constant $C$, independent of the discretization parameters $h$ and $\mathrm{\Delta }t$ such that for all ${\{{\chi }_h^i,{\xi }_h^i\}}_N^{i=0}\subset S_0^h$,

(24) $\underset{0\le n\le N}{\max }\left\{{\Vert y^n-y_h^n\Vert }^2+{\Vert {\mathrm{\Phi }}^n-{\mathrm{\Phi }}_h^n\Vert }^2+{\Vert {\mathrm{\Psi }}^n-{\mathrm{\Psi }}_h^n\Vert }^2+{\Vert {\theta }^n-{\theta }_h^n\Vert }^2\right\}\le$

$\le C\mathrm{\Delta }t{\displaystyle \sum _{i=1}^N}(\parallel y_t^i-\delta y^i{\parallel }^2+\parallel {\mathrm{\Phi }}_t^i-\delta {\mathrm{\Phi }}^i{\parallel }^2+\parallel {\mathrm{\Psi }}_t^i-\delta {\mathrm{\Psi }}^i{\parallel }^2+\parallel {\theta }_t^i-\delta {\theta }^i{\parallel }^2$

$+\parallel \nabla y^i-\nabla {\chi }_h^i{\parallel }^2+\parallel \nabla {\theta }^i-\nabla {\xi }_h^i{\parallel }^2)+C\max {}_{0\le n\le N}\{\parallel y^n-{\chi }_h^n{\parallel }^2+\parallel {\theta }^n-{\xi }_h^n{\parallel }^2\}$

$+\frac{C}{\mathrm{\Delta }t}{\displaystyle \sum _{i=1}^{N-1}}\left({\Vert y^i-{\chi }_h^i-\left(y^{i+1}-{\chi }_h^{i+1}\right)\Vert }^2+{\Vert {\theta }^i-{\xi }_h^i-\left({\theta }^{i+1}-{\xi }_h^{i+1}\right)\Vert }^2\right)$

$+C\left({\Vert y^0-y_h^0\Vert }^2+{\Vert {\mathrm{\Phi }}^0-{\mathrm{\Phi }}_h^0\Vert }^2+{\Vert {\mathrm{\Psi }}^0-{\mathrm{\Psi }}_h^0\Vert }^2+{\Vert {\theta }^0-{\theta }_h^0\Vert }^2\right),$

where $\delta f^i=\left(f^i-f^{i-1}\right)/\mathrm{\Delta }t$.

Proof.

First, we subtract the first variational equation in $(\text{<a href="\#S2.E9" title="Equation 9 ‣ 2. Numerical approximation ‣ Numerical Analysis and Stability of the Moore-Gibson-Thompson-Fourier Model" class="ltx\_ref ltx\_font\_italic"><span class="ltx\_text ltx\_ref\_tag">9</span></a>})$ at time $t=t_n$ for a test function $\chi ={\chi }_h\in S_0^h\subset H_0^1(0,1)$ and the first discrete variational equation in $(\text{<a href="\#S2.E11" title="Equation 11 ‣ 2. Numerical approximation ‣ Numerical Analysis and Stability of the Moore-Gibson-Thompson-Fourier Model" class="ltx\_ref ltx\_font\_italic"><span class="ltx\_text ltx\_ref\_tag">11</span></a>})$ to obtain

(25) $(y_t^n-\delta y_h^n,{\chi }_h)+\frac{\gamma }{\alpha }(\nabla {\mathrm{\Phi }}^n-\nabla {\mathrm{\Phi }}_h^n,\nabla {\chi }_h)+\frac{\nu }{\alpha }(\nabla {\mathrm{\Psi }}^n-\nabla {\mathrm{\Psi }}_h^n,\nabla {\chi }_h)-$

$-\eta (\nabla {\theta }^n-\nabla {\theta }_h^n,\nabla {\chi }_h)=0,\text{for all }{\chi }_h\in S_0^h$

and so, we have

(26) $(y_t^n-\delta y_h^n,y^n-y_h^n)+\frac{\gamma }{\alpha }(\nabla {\mathrm{\Phi }}^n-\nabla {\mathrm{\Phi }}_h^n,\nabla (y^n-y_h^n))+\frac{\nu }{\alpha }(\nabla {\mathrm{\Psi }}^n-\nabla {\mathrm{\Psi }}_h^n,\nabla (y^n-y_h^n))$

$-\eta (\nabla {\theta }^n-\nabla {\theta }_h^n,\nabla (y^n-y_h^n))$

$=(y_t^n-\delta y_h^n,y^n-{\chi }_h)+\frac{\gamma }{\alpha }(\nabla {\mathrm{\Phi }}^n-\nabla {\mathrm{\Phi }}_h^n,\nabla (y^n-{\chi }_h))+\frac{\nu }{\alpha }(\nabla {\mathrm{\Psi }}^n-\nabla {\mathrm{\Psi }}_h^n,\nabla (y^n-{\chi }_h))$

$-\eta (\nabla {\theta }^n-\nabla {\theta }_h^n,\nabla (y^n-{\chi }_h)),\text{for all }{\chi }_h\in S_0^h$

Taking into account that

$(y_t^n-\delta y_h^n,y^n-y_h^n)$ $=(y_t^n-\delta y^n,y^n-y_h^n)+(\delta y^n-\delta y_h^n,y^n-y_h^n)=(y_t^n-\delta y^n,y^n-y_h^n)$

$+\frac{1}{2\mathrm{\Delta }t}\left({\Vert y^n-y_h^n-(y^{n-1}-y_h^{n-1})\Vert }^2+{\Vert y^n-y_h^n\Vert }^2-{\Vert y^{n-1}-y_h^{n-1}\Vert }^2\right)$

By the positivity of the terms ${\Vert y^n-y_h^n-(y^{n-1}-y_h^{n-1})\Vert }^2$, we get the following inequality

(27) $(y_t^n-\delta y_h^n,y^n-y_h^n)$ $\ge (y_t^n-\delta y^n,y^n-y_h^n)+\frac{1}{2\mathrm{\Delta }t}\left({\Vert y^n-y_h^n\Vert }^2-{\Vert y^{n-1}-y_h^{n-1}\Vert }^2\right),$

$(\nabla {\mathrm{\Phi }}^n-\nabla {\mathrm{\Phi }}_h^n,\nabla (y^n-y_h^n))$ $=({\mathrm{\Phi }}_t^n-\delta {\mathrm{\Phi }}_h^n,{\mathrm{\Phi }}^n-{\mathrm{\Phi }}_h^n)$

$=({\mathrm{\Phi }}_t^n-\delta {\mathrm{\Phi }}^n,{\mathrm{\Phi }}^n-{\mathrm{\Phi }}_h^n)+(\delta {\mathrm{\Phi }}^n-\delta {\mathrm{\Phi }}_h^n,{\mathrm{\Phi }}^n-{\mathrm{\Phi }}_h^n)$

$\ge ({\mathrm{\Phi }}_t^n-\delta {\mathrm{\Phi }}^n,{\mathrm{\Phi }}^n-{\mathrm{\Phi }}_h^n)$

$+\frac{1}{2\mathrm{\Delta }t}(\parallel {\mathrm{\Phi }}^n-{\mathrm{\Phi }}_h^n){\parallel }^2-\parallel {\mathrm{\Phi }}^{n-1}-{\mathrm{\Phi }}_h^{n-1}{\parallel }^2),$

$(\nabla {\mathrm{\Psi }}^n-\nabla {\mathrm{\Psi }}_h^n,\nabla (y^n-y_h^n))$ $=({\mathrm{\Psi }}^n-{\mathrm{\Psi }}_h^n,{\mathrm{\Phi }}_t^n-\delta {\mathrm{\Phi }}_h^n)$

$=({\mathrm{\Psi }}^n-{\mathrm{\Psi }}_h^n,({\mathrm{\Psi }}_t^n+\alpha {\mathrm{\Psi }}^n)-(\delta {\mathrm{\Psi }}_h^n+\alpha {\mathrm{\Psi }}_h^n))$

$=({\mathrm{\Psi }}^n-{\mathrm{\Psi }}_h^n,{\mathrm{\Psi }}_t^n-\delta {\mathrm{\Psi }}_h^n)+\alpha ({\mathrm{\Psi }}^n-{\mathrm{\Psi }}_h^n,{\mathrm{\Psi }}^n-{\mathrm{\Psi }}_h^n)$

$=({\mathrm{\Psi }}_t^n-\delta {\mathrm{\Psi }}_h^n,{\mathrm{\Psi }}^n-{\mathrm{\Psi }}_h^n)+\alpha {\Vert {\mathrm{\Psi }}^n-{\mathrm{\Psi }}_h^n\Vert }^2$

$\ge ({\mathrm{\Psi }}_t^n-\delta {\mathrm{\Psi }}^n,{\mathrm{\Psi }}^n-{\mathrm{\Psi }}_h^n)+\alpha \parallel {\mathrm{\Psi }}^n-{\mathrm{\Psi }}_h^n)\parallel {}^2$

$+\frac{1}{2\mathrm{\Delta }t}(\parallel {\mathrm{\Psi }}^n-{\mathrm{\Psi }}_h^n){\parallel }^2-\parallel {\mathrm{\Psi }}^{n-1}-{\mathrm{\Psi }}_h^{n-1}{\parallel }^2).$

Second, we subtract the second variational equation in $(\text{<a href="\#S2.E9" title="Equation 9 ‣ 2. Numerical approximation ‣ Numerical Analysis and Stability of the Moore-Gibson-Thompson-Fourier Model" class="ltx\_ref ltx\_font\_italic"><span class="ltx\_text ltx\_ref\_tag">9</span></a>})$ at time $t=t_n$ for a test function $\xi ={\xi }_h\in S_0^h\subset H_0^1(0,1)$ and the second discrete variational equation in $(\text{<a href="\#S2.E11" title="Equation 11 ‣ 2. Numerical approximation ‣ Numerical Analysis and Stability of the Moore-Gibson-Thompson-Fourier Model" class="ltx\_ref ltx\_font\_italic"><span class="ltx\_text ltx\_ref\_tag">11</span></a>})$ to obtain

(28) $({\theta }_t^n-\delta {\theta }_h^n,{\xi }_h)+\kappa (\nabla {\theta }^n-\nabla {\theta }_h^n,\nabla {\xi }_h)+\eta (\nabla y^n-\nabla y_h^n,\nabla {\xi }_h)=0,$

and so, we have

(29) $({\theta }_t^n-\delta {\theta }_h^n,{\theta }^n-{\theta }_h^n)+\kappa (\nabla {\theta }^n-\nabla {\theta }_h^n,\nabla ({\theta }^n-{\theta }_h^n))+\eta (\nabla y^n-\nabla y_h^n,\nabla ({\theta }^n-{\theta }_h^n))$

$=({\theta }_t^n-\delta {\theta }_h^n,{\theta }^n-{\xi }_h)+\kappa (\nabla {\theta }^n-\nabla {\theta }_h^n,\nabla ({\theta }^n-{\xi }_h))+\eta (\nabla y^n-\nabla y_h^n,\nabla ({\theta }^n-{\xi }_h)).$

Taking into account that

(30) $({\theta }_t^n-\delta {\theta }_h^n,{\theta }^n-y_h^n)$ $=({\theta }_t^n-\delta {\theta }^n,{\theta }^n-{\theta }_h^n)+(\delta {\theta }^n-\delta {\theta }_h^n,{\theta }^n-{\theta }_h^n)$

$\ge ({\theta }_t^n-\delta {\theta }^n,{\theta }^n-{\theta }_h^n)+\frac{1}{2\mathrm{\Delta }t}\left({\Vert {\theta }^n-{\theta }_h^n\Vert }^2-{\Vert {\theta }^{n-1}-{\theta }_h^{n-1}\Vert }^2\right).$

From (26)–(27) and using several times Young’s inequality $(\text{<a href="\#S2.E31" title="Equation 31 ‣ 2. Numerical approximation ‣ Numerical Analysis and Stability of the Moore-Gibson-Thompson-Fourier Model" class="ltx\_ref ltx\_font\_italic"><span class="ltx\_text ltx\_ref\_tag">31</span></a>})$

(31) $$ab\le \epsilon a^2+\frac{1}{4\epsilon }{b}^2,a,b\in ℝ,\epsilon \in {ℝ}^{+\ast }.$$

$(y_t^n-\delta y^n,y^n-y_h^n)+\frac{1}{2\mathrm{\Delta }t}\left({\Vert y^n-y_h^n\Vert }^2-{\Vert y^{n-1}-y_h^{n-1}\Vert }^2\right)+({\mathrm{\Phi }}_t^n-\delta {\mathrm{\Phi }}^n,{\mathrm{\Phi }}^n-{\mathrm{\Phi }}_h^n)+$

$+\frac{1}{2\mathrm{\Delta }t}(\parallel {\mathrm{\Phi }}^n-{\mathrm{\Phi }}_h^n){\parallel }^2-\parallel {\mathrm{\Phi }}^{n-1}-{\mathrm{\Phi }}_h^{n-1}{\parallel }^2)+({\mathrm{\Psi }}_t^n-\delta {\mathrm{\Psi }}^n,{\mathrm{\Psi }}^n-{\mathrm{\Psi }}_h^n)+\alpha \parallel {\mathrm{\Psi }}^n-{\mathrm{\Psi }}_h^n)\parallel {}^2$

$+\frac{1}{2\mathrm{\Delta }t}(\parallel {\mathrm{\Psi }}^n-{\mathrm{\Psi }}_h^n){\parallel }^2-\parallel {\mathrm{\Psi }}^{n-1}-{\mathrm{\Psi }}_h^{n-1}{\parallel }^2)-\eta (\nabla {\theta }^n-\nabla {\theta }_h^n,\nabla (y^n-y_h^n))\le$

$\le (y_t^n-\delta y_h^n,y^n-y_h^n)+\frac{\gamma }{\alpha }(\nabla {\mathrm{\Phi }}^n-\nabla {\mathrm{\Phi }}_h^n,\nabla (y^n-y_h^n))+\frac{\nu }{\alpha }(\nabla {\mathrm{\Psi }}^n-\nabla {\mathrm{\Psi }}_h^n,\nabla (y^n-y_h^n))$

$-\eta (\nabla {\theta }^n-\nabla {\theta }_h^n,\nabla (y^n-y_h^n))$

$=(y_t^n-\delta y_h^n,y^n-{\chi }_h)+\frac{\gamma }{\alpha }(\nabla {\mathrm{\Phi }}^n-\nabla {\mathrm{\Phi }}_h^n,\nabla (y^n-{\chi }_h))+\frac{\nu }{\alpha }(\nabla {\mathrm{\Psi }}^n-\nabla {\mathrm{\Psi }}_h^n,\nabla (y^n-{\chi }_h))$

$-\eta (\nabla {\theta }^n-\nabla {\theta }_h^n,\nabla (y^n-{\chi }_h)),\text{for all }{\chi }_h\in S_0^h.$

Next,

$\frac{1}{2\mathrm{\Delta }t}(\parallel y^n-y_h^n{\parallel }^2-\parallel y^{n-1}-y_h^{n-1}{\parallel }^2)+\frac{1}{2\mathrm{\Delta }t}(\parallel {\mathrm{\Phi }}^n-{\mathrm{\Phi }}_h^n){\parallel }^2-\parallel {\mathrm{\Phi }}^{n-1}-{\mathrm{\Phi }}_h^{n-1}{\parallel }^2)$

$+\alpha \parallel {\mathrm{\Psi }}^n{\mathrm{\Psi }}_h^n){\parallel }^2+\frac{1}{2\mathrm{\Delta }t}(\parallel {\mathrm{\Psi }}^n-{\mathrm{\Psi }}_h^n){\parallel }^2-\parallel {\mathrm{\Psi }}^{n-1}-{\mathrm{\Psi }}_h^{n-1}{\parallel }^2)-\eta (\nabla {\theta }^n-\nabla {\theta }_h^n,\nabla (y^n-y_h^n))$

$\le (y_t^n-\delta y_h^n,y^n-{\chi }_h)+\frac{\gamma }{\alpha }(\nabla {\mathrm{\Phi }}^n-\nabla {\mathrm{\Phi }}_h^n,\nabla (y^n-{\chi }_h))+\frac{\nu }{\alpha }(\nabla {\mathrm{\Psi }}^n-\nabla {\mathrm{\Psi }}_h^n,\nabla (y^n-{\chi }_h))$

$-\eta (\nabla {\theta }^n-\nabla {\theta }_h^n,\nabla (y^n-{\chi }_h))-(y_t^n-\delta y^n,y^n-y_h^n)-({\mathrm{\Phi }}_t^n-\delta {\mathrm{\Phi }}^n,{\mathrm{\Phi }}^n-{\mathrm{\Phi }}_h^n)$

$-({\mathrm{\Psi }}_t^n-\delta {\mathrm{\Psi }}^n,{\mathrm{\Psi }}^n-{\mathrm{\Psi }}_h^n),\text{for all }{\chi }_h\in S_0^h.$

It follows that

(32) $\frac{1}{2\mathrm{\Delta }t}(\parallel y^n-y_h^n{\parallel }^2-\parallel y^{n-1}-y_h^{n-1}{\parallel }^2)+\frac{\gamma }{2\alpha \mathrm{\Delta }t}(\parallel {\mathrm{\Phi }}^n-{\mathrm{\Phi }}_h^n){\parallel }^2-\parallel {\mathrm{\Phi }}^{n-1}-{\mathrm{\Phi }}_h^{n-1}{\parallel }^2)$

$+\frac{\nu }{2\alpha \mathrm{\Delta }t}({\Vert {\mathrm{\Psi }}^n-{\mathrm{\Psi }}_h^n\Vert }^2-{\Vert {\mathrm{\Psi }}^{n-1}-{\mathrm{\Psi }}_h^{n-1}\Vert }^2)+\nu {\Vert {\mathrm{\Psi }}^n-{\mathrm{\Psi }}_h^n\Vert }^2-\eta (\nabla {\theta }^n-\nabla {\theta }_h^n,\nabla y^n-\nabla y_h^n)$

$\le C(\parallel y_t^n-\delta y^n{\parallel }^2+\parallel y^n-y_h^n{\parallel }^2+\parallel {\mathrm{\Phi }}_t^n-\delta {\mathrm{\Phi }}^n{\parallel }^2+\parallel {\mathrm{\Phi }}^n-{\mathrm{\Phi }}_h^n){\parallel }^2+\parallel {\mathrm{\Psi }}_t^n-\delta {\mathrm{\Psi }}^n{\parallel }^2$

$+\parallel {\mathrm{\Psi }}^n-{\mathrm{\Psi }}_h^n){\parallel }^2+\parallel \nabla {\theta }^n-\nabla {\theta }_h^n{\parallel }^2+\parallel y^n-{\chi }_h{\parallel }^2+\parallel \nabla y^n-\nabla {\chi }_h{\parallel }^2+\parallel \mathrm{\Delta }{y}^n-\mathrm{\Delta }{\chi }_h{\parallel }^2)$

$+(\delta y^n-\delta y_h^n,y^n-{\chi }_h),\text{for all }{\chi }_h\in S_0^h.$

Proceeding with a similar approach for equations (29)-(30), we obtain the following estimates, for all ${\xi }_h\in S_0^h$,

(33) $\frac{1}{2\mathrm{\Delta }t}\left({\Vert {\theta }^n-{\theta }_h^n\Vert }^2-{\Vert {\theta }^{n-1}-{\theta }_h^{n-1}\Vert }^2\right)+\kappa {\Vert \nabla {\theta }^n-\nabla {\theta }_h^n\Vert }^2+\eta (\nabla y^n-\nabla y_h^n,\nabla {\theta }^n-\nabla {\theta }_h^n)$

$\le C(\parallel {\theta }_t^n-\delta {\theta }^n{\parallel }^2+\parallel {\theta }^n-{\theta }_h^n{\parallel }^2+\parallel \nabla {\theta }^n-\nabla {\theta }_h^n{\parallel }^2+\parallel \nabla y^n-\nabla y_h^n{\parallel }^2+\parallel {\theta }^n-{\xi }_h{\parallel }^2$

$+\parallel \nabla {\theta }^n-\nabla {\xi }_h{\parallel }^2)+(\delta {\theta }^n-\delta {\theta }_h^n,{\theta }^n-{\xi }_h).$

Combining estimates $(\text{<a href="\#S2.E32" title="Equation 32 ‣ 2. Numerical approximation ‣ Numerical Analysis and Stability of the Moore-Gibson-Thompson-Fourier Model" class="ltx\_ref ltx\_font\_italic"><span class="ltx\_text ltx\_ref\_tag">32</span></a>})$ and $(\text{<a href="\#S2.E33" title="Equation 33 ‣ 2. Numerical approximation ‣ Numerical Analysis and Stability of the Moore-Gibson-Thompson-Fourier Model" class="ltx\_ref ltx\_font\_italic"><span class="ltx\_text ltx\_ref\_tag">33</span></a>})$ it follows that, for all ${\chi }_h,{\xi }_h\in S_0^h$,

(34) $\frac{1}{2\mathrm{\Delta }t}(\parallel y^n-y_h^n{\parallel }^2-\parallel y^{n-1}-y_h^{n-1}{\parallel }^2)+\frac{\gamma }{2\alpha \mathrm{\Delta }t}(\parallel {\mathrm{\Phi }}^n-{\mathrm{\Phi }}_h^n){\parallel }^2-\parallel {\mathrm{\Phi }}^{n-1}-{\mathrm{\Phi }}_h^{n-1}{\parallel }^2)$

$+\frac{\nu }{2\alpha \mathrm{\Delta }t}({\Vert {\mathrm{\Psi }}^n-{\mathrm{\Psi }}_h^n\Vert }^2-{\Vert {\mathrm{\Psi }}^{n-1}-{\mathrm{\Psi }}_h^{n-1}\Vert }^2)+\nu {\Vert {\mathrm{\Psi }}^n-{\mathrm{\Psi }}_h^n\Vert }^2-\eta (\nabla {\theta }^n-\nabla {\theta }_h^n,\nabla y^n-\nabla y_h^n)$

$\frac{1}{2\mathrm{\Delta }t}\left({\Vert {\theta }^n-{\theta }_h^n\Vert }^2-{\Vert {\theta }^{n-1}-{\theta }_h^{n-1}\Vert }^2\right)+\kappa {\Vert \nabla {\theta }^n-\nabla {\theta }_h^n\Vert }^2+\eta (\nabla y^n-\nabla y_h^n,\nabla {\theta }^n-\nabla {\theta }_h^n)$

$\le C(\parallel y_t^n-\delta y^n{\parallel }^2+\parallel y^n-y_h^n{\parallel }^2+\parallel {\mathrm{\Phi }}_t^n-\delta {\mathrm{\Phi }}^n{\parallel }^2+\parallel {\mathrm{\Phi }}^n-{\mathrm{\Phi }}_h^n{\parallel }^2+\parallel {\mathrm{\Psi }}_t^n-\delta {\mathrm{\Psi }}^n{\parallel }^2$

$+{\Vert {\mathrm{\Psi }}^n-{\mathrm{\Psi }}_h^n\Vert }^2+{\Vert \nabla {\theta }^n-\nabla {\theta }_h^n\Vert }^2+{\Vert y^n-{\chi }_h\Vert }^2+{\Vert \nabla y^n-\nabla {\chi }_h\Vert }^2+{\Vert \mathrm{\Delta }{y}^n-\mathrm{\Delta }{\chi }_h\Vert }^2+{\Vert {\theta }_t^n-\delta {\theta }^n\Vert }^2$

$+\parallel {\theta }^n-{\theta }_h^n{\parallel }^2+\parallel \nabla {\theta }^n-\nabla {\theta }_h^n{\parallel }^2+\parallel \nabla y^n-\nabla y_h^n{\parallel }^2+\parallel {\theta }^n-{\xi }_h{\parallel }^2+\parallel \nabla {\theta }^n-\nabla {\xi }_h{\parallel }^2)$

$+(\delta y^n-\delta y_h^n,y^n-{\chi }_h)+(\delta {\theta }^n-\delta {\theta }_h^n,{\theta }^n-{\xi }_h).$

Multiplying the above estimates by $\mathrm{\Delta }t$ and summing up to $n$ we find that, for all ${\chi }_h,{\xi }_h\in S_0^h$,

(35) $\parallel y^n-y_h^n{\parallel }^2+\parallel {\mathrm{\Phi }}^n-{\mathrm{\Phi }}_h^n)\parallel {}^2+\parallel \mathrm{\Psi }{}^n-\mathrm{\Psi }{}_h{}^n\parallel {}^2+\parallel \theta {}^n-\theta {}_h{}^n\parallel {}^2\le$

$\le C\mathrm{\Delta }t{\displaystyle \sum _{i=0}^n}(\parallel y_t^i-\delta y^i{\parallel }^2+\parallel y^i-y_h^i{\parallel }^2+\parallel {\mathrm{\Phi }}_t^i-\delta {\mathrm{\Phi }}^i{\parallel }^2+\parallel {\mathrm{\Phi }}^i-{\mathrm{\Phi }}_h^i{\parallel }^2+\parallel {\mathrm{\Psi }}_t^i-\delta {\mathrm{\Psi }}^i{\parallel }^2$

$+\parallel {\mathrm{\Psi }}^i-{\mathrm{\Psi }}_h^i)\parallel {}^2+\parallel \nabla \theta {}^i-\nabla \theta {}_h{}^i\parallel {}^2+\parallel \nabla y{}^i-\nabla \chi {}_h{}^i\parallel {}^2+\parallel \theta {}_t{}^i-\delta \theta {}^i\parallel {}^2+\parallel \theta {}^i-\theta {}_h{}^i\parallel {}^2$

$+\parallel \nabla {\theta }^i-\nabla {\theta }_h^i{\parallel }^2+\parallel y^i-{\chi }_h^i{\parallel }^2+\parallel \delta y^i-\delta y_h^i{\parallel }^2+\parallel \nabla y^i-\nabla y_h^i{\parallel }^2+\parallel {\theta }^i-{\xi }_h^i{\parallel }^2+\parallel \nabla {\theta }^i-\nabla {\xi }_h^i{\parallel }^2)$

$+\mathrm{\Delta }t{\displaystyle \sum _{i=0}^n}\left((\delta y^i-\delta y_h^i,y^i-{\chi }_h^i)+(\delta {\theta }^i-\delta {\theta }_h^i,{\theta }^i-{\xi }_h^i)\right)$

$+C(\parallel y^0-y_h^0{\parallel }^2+\parallel {\mathrm{\Phi }}^0-{\mathrm{\Phi }}_h^0){\parallel }^2\parallel {\mathrm{\Psi }}^0-{\mathrm{\Psi }}_h^0{\parallel }^2+\parallel {\theta }^0-{\theta }_h^0{\parallel }^2).$

Finally, taking into account that

(36) $\mathrm{\Delta }t{\displaystyle \sum _{i=1}^n}(\delta y^i-\delta y_h^i,y^i-{\chi }_h^i)=$ $(y^n-y_h^n,y^n-{\chi }_h^n)+(y_h^0-z^1,y^1-{\chi }_h^1)+$

$+{\displaystyle \sum _{i=1}^{n-1}}(y^i-y_h^i,y^i-{\chi }_h^i-\left(y^{i+1}-{\chi }_h^{i+1}\right)),$

$\mathrm{\Delta }t{\displaystyle \sum _{i=1}^n}(\delta {\theta }^i-\delta {\theta }_h^i,{\theta }^i-{\xi }_h^i)=$ $({\theta }^n-{\theta }_h^n,{\theta }^n-{\xi }_h^n)+({\theta }_h^0-{\theta }^0,{\theta }^1-{\xi }_h^1)+$

$+{\displaystyle \sum _{i=1}^{n-1}}({\theta }^i-{\theta }_h^i,{\theta }^i-{\xi }_h^i-\left({\theta }^{i+1}-{\xi }_h^{i+1}\right))$

using again a discrete version of Gronwall’s inequality (see [6]) we obtain the desired a priori error estimates. ∎

The estimates provided in the above theorem can be used to obtain the convergence order of the approximations given by discrete problem $(\text{<a href="\#S2.E11" title="Equation 11 ‣ 2. Numerical approximation ‣ Numerical Analysis and Stability of the Moore-Gibson-Thompson-Fourier Model" class="ltx\_ref ltx\_font\_italic"><span class="ltx\_text ltx\_ref\_tag">11</span></a>})$. Hence, as an example, if we assume the additional regularity:

(37) $u\in H^4(0,T;L^2(0,1))\cap H^3(0,T;H^1(0,1))\cap C^2([0,T];H^3(0,1))$

$\theta \in H^2(0,T;L^2(0,1))\cap H^1(0,T;H^1(0,1))\cap C^0([0,T];H^1(0,1))$

we obtain the quadratic convergence of the algorithm applying some results on the approximation by finite elements (see [7]) and previous estimates already derived in [6]. We have the following.

Corollary 4.

Let $(y,\mathrm{\Phi },\mathrm{\Psi },\theta )$ be the solution of $(\text{<a href="\#S2.E8" title="Equation 8 ‣ 2. Numerical approximation ‣ Numerical Analysis and Stability of the Moore-Gibson-Thompson-Fourier Model" class="ltx\_ref"><span class="ltx\_text ltx\_ref\_tag">8</span></a>})$ and $(y_h,{\mathrm{\Phi }}_h,{\mathrm{\Psi }}_h,{\theta }_h)$ be that of the discrete system $(\text{<a href="\#S2.E11" title="Equation 11 ‣ 2. Numerical approximation ‣ Numerical Analysis and Stability of the Moore-Gibson-Thompson-Fourier Model" class="ltx\_ref"><span class="ltx\_text ltx\_ref\_tag">11</span></a>})$. Under the assumptions of Theorem 3, it follows that there exists a positive constant $C>0$, independent of the discretization parameters $h$ and $\mathrm{\Delta }t$, such that

(38) $\underset{0\le n\le N}{\max }\left\{{\Vert y^n-y_h^n\Vert }^2+{\Vert {\mathrm{\Phi }}^n-{\mathrm{\Phi }}_h^n\Vert }^2+{\Vert {\mathrm{\Psi }}^n-{\mathrm{\Psi }}_h^n\Vert }^2+{\Vert {\theta }^n-{\theta }_h^n\Vert }^2\right\}\le C(h^2+\mathrm{\Delta }{t}^2).$

3. Numerical Simulation

3.1. Numerical Convergence: error estimate with an exact solution

In a first example, our aim is to show the accuracy and efficiency of the proposed fully discrete example. Therefore, we will solve the problem:

(39)

with the following data:

(40) $T=1,\alpha =2\cdot {10}^{-2},\beta ={3.10}^{-3},\gamma ={10}^{-5},\eta ={10}^{-4},\kappa ={10}^{-5}.$

If we use the following initial conditions, for all $x\in (0,1)$,

(41) $u_0(x)=u_1(x)=u_2(x)=x^3{(1-x)}^3,{\theta }_0(x)=x^3{(1-x)}^3.$

considering homogeneous Dirichlet boundary conditions.

In the previous system of equations, the source terms $f_i,i=1,2$, can be easily calculate from the exact solution to the above problem and it has the form, for $(x,t)\in$ $[0,1]\times [0,1]:$

$$u(x,t)=e^t{x}^3{(1-x)}^3,\theta (x,t)=e^t{x}^3{(1-x)}^3.$$

Hence, for some values of the spatial and time discretization parameters, the approximated numerical errors given by (11) are shown in Table 1.

Fig. 1 illustrates how the error depends on the parameters $h^2$ and $\mathrm{\Delta }{t}^2$, demonstrating quadratic convergence. This confirms the theoretical results, assuming the solution meets certain regularity conditions. Moreover, Fig. 2 further validates these findings for different cases.

$h\downarrow \mathrm{\Delta }t\to$	0.02	0.01	0.005	0.0025	0.00125
0.02	0.143200720	0.052391942	0.024164654	0.017698156	0.017492619
0.01	0.100412690	0.028129534	0.008809871	0.003280523	0.001545261
0.005	0.090920463	0.023251807	0.006174570	0.001746539	0.000551093
0.0025	0.088623052	0.022107428	0.005592580	0.001442328	0.000384684
0.00125	0.088053433	0.021826038	0.005451913	0.001371276	0.000348277

Table 1. Computed numerical errors $\times {10}^{-4}$ for a final time $T=1$ and for some values of $h$ and $\mathrm{\Delta }t$.

Figure 1. Error behavior on the logarithmic scale.

Case	$h$	$\mathrm{\Delta }t$
Case 1: $\mathrm{\Delta }t\ne h$	$\frac{1}{40}$, $\frac{1}{80}$, $\frac{1}{160}$, $\frac{1}{320}$	$\frac{1}{50}$, $\frac{1}{100}$, $\frac{1}{200}$, $\frac{1}{400}$
Case 2: $\mathrm{\Delta }t=4h$	$\frac{1}{120}$, $\frac{1}{240}$, $\frac{1}{480}$, $\frac{1}{960}$	$\frac{1}{30}$, $\frac{1}{60}$, $\frac{1}{120}$, $\frac{1}{240}$
Case 3: $h$ fixed, $\mathrm{\Delta }t$ decreasing	$\frac{1}{300}$	$\frac{1}{54}$, $\frac{1}{108}$, $\frac{1}{216}$, $\frac{1}{432}$
Case 4: $\mathrm{\Delta }t$ fixed, $h$ decreasing	$\frac{1}{80}$, $\frac{1}{100}$, $\frac{1}{120}$, $\frac{1}{140}$	$\frac{1}{700}$

Table 2. Values of $h$ and $\mathrm{\Delta }t$ for different cases.

(a) Error behavior on the logarithmic scale for case 1.

(b) Error behavior on the logarithmic scale for case 2.

(c) Error behavior on the logarithmic scale for case 3.

(d) Error behavior on the logarithmic scale for case 4.

Figure 2. Error Behavior.

3.2. Discrete Energy: exponential decay

Now, we consider the system (8) with the following data :

(42) $$T=12,\alpha =10,\beta =2,\gamma =1,\eta =1,\kappa =1.$$

and the following initial conditions, for all $x\in (0,1)$,

(43) $$u_0(x)=u_1(x)=u_2(x)=x^3{(1-x)}^3,{\theta }_0(x)=x^2{(1-x)}^2\sin (x).$$

If we take the parameters $3h=\mathrm{\Delta }t=0.03$ and we use the following definition for the discrete energy $(\text{<a href="\#S2.E14" title="Equation 14 ‣ 2. Numerical approximation ‣ Numerical Analysis and Stability of the Moore-Gibson-Thompson-Fourier Model" class="ltx\_ref ltx\_font\_italic"><span class="ltx\_text ltx\_ref\_tag">14</span></a>})$ in Fig. 3(a) and Fig. 3(b) we represent discrete energy and discrete logarithm energy evolution of system $(\text{<a href="\#S1.E2" title="Equation 2 ‣ 1. Introduction ‣ Numerical Analysis and Stability of the Moore-Gibson-Thompson-Fourier Model" class="ltx\_ref ltx\_font\_italic"><span class="ltx\_text ltx\_ref\_tag">2</span></a>})$. We can clearly conclude that the discrete energy tends to zero and that an exponential energy decay is achieved.

(a) Natural scale behavior of ${ℰ}_h^n$

(b) Semi-log scale behavior of ${ℰ}_h^n$

Figure 3. Energy Behavior.

The numerical schemes were implemented using MATLAB.