*3.2. Procedure for Solving Inverse Problem of TVIDE*

For the inverse problem in this paper, we specifically define the forcing term *F*(*x*, *t*) := *β*(*t*)*f*(*x*, *t*), where *β*(*t*) is a missing source function to be retrieved and *f*(*x*, *t*) is the given function. Thus, our considered time-dependent inverse TVIDE problem (1) becomes

$$u\_t(\mathbf{x},t) + \mathcal{L}u(\mathbf{x},t) = \int\_0^t \mathbf{x}\_1(\mathbf{x},\eta)u(\mathbf{x},\eta)d\eta + \int\_0^\mathbf{x} \mathbf{x}\_2(\mathbf{\tilde{x}},t)u(\mathbf{\tilde{x}},t)d\mathbf{\tilde{s}} + \beta(t)f(\mathbf{x},t),\tag{19}$$

where *u* is an approximate solution of *v* and the other parameters are defined as in (8). The initial and boundary conditions of (19) are (9) and (10), which satisfy the compatibility conditions. Now, we remove all spatial derivatives from (19) and use the shifted Chebyshev integration matrix (as explained in Section 2.1). Then, we obtain the following matrix equation, based on the same process as in (16), as

$$\begin{cases} \mathbf{A}^{n} + \tau \sum\_{j=0}^{n} \sum\_{i=0}^{n-j} (-1)^{i} \binom{n-j}{i} \mathbf{A}^{i+j} \mathbf{P}\_{n-j}^{(i)} - \frac{\tau^{2}}{2} \mathbf{A}^{n} \mathbf{K}\_{1}^{\langle m \rangle} - \tau \mathbf{A}^{n+1} \mathbf{K}\_{2}^{\langle m \rangle} \Big) \mathbf{u}^{\langle m \rangle} \\ + \mathbf{X}\_{n} \mathbf{d} - \tau \mathbf{A}^{n} \mathbf{f}^{\langle m \rangle} \boldsymbol{\beta}^{\langle m \rangle} = \frac{\tau^{2}}{2} \mathbf{A}^{n} \mathbf{K}\_{1}^{\langle 0 \rangle} \mathbf{u}^{\langle 0 \rangle} + \tau^{2} \sum\_{i=1}^{m-1} \mathbf{A}^{n} \mathbf{K}\_{1}^{\langle i \rangle} \mathbf{u}^{\langle i \rangle} + \mathbf{A}^{n} \mathbf{u}^{\langle m-1 \rangle}, \end{cases} \tag{20}$$

where *βm* = *β*(*tm*), **f***m* = [ *f*(*x*1, *tm*), *f*(*x*2, *tm*), *f*(*x*3, *tm*), ..., *f*(*xN*, *tm*)] and the other parameters in (20) are as defined in Section 3.1. However, the occurrence of missing data is caused by the given conditions being insufficient to ensure a unique solution to our inverse problem. Hence, an additional condition or observed data needs to be involved. Thus, we use an additional condition, regarding the aggregated solution of the system, in the following form:

$$\int\_{0}^{L} u(\xi, t) \, d\xi = \mathcal{g}(t), \; t \in [0, T], \tag{21}$$

where *g*(*t*) is the measured data at time *t*, which probably contains measurement errors. In order to illustrate the realistic phenomena of this problem, we assume that the measurement data of the aggregated solution *<sup>g</sup>*(*t*) involves some noise , which is denoted by *<sup>g</sup>*(*t*) (where *<sup>g</sup>*(*t*) <sup>−</sup> *<sup>g</sup>*(*t*)<sup>≤</sup> ) and define the noisy value by a random variable generated by the Gaussian normal distribution with mean *μ* = 0 and standard deviation *σ* = *p*|*g*(*t*)|, where *p* is the percentage of the noise to be input. Then, the additional condition (21) becomes

$$\int\_{0}^{L} u(\xi, t) \, d\xi = \mathcal{g}^{\mathfrak{c}}(t), \; t \in [0, T]. \tag{22}$$

Using the concept of FIM-SCP, the additional condition (22) at time *tm* can be written, in vector form, as

$$\int\_{0}^{L} u^{\langle m \rangle} \left( \xi \right) d\xi = \sum\_{n=0}^{N-1} c\_{n}^{\langle m \rangle} \int\_{0}^{L} \mathbf{S}\_{n} \left( \xi \right) d\xi = \sum\_{n=0}^{N-1} c\_{n}^{\langle m \rangle} \mathbf{S}\_{n} \left( L \right) d\xi := \mathbf{z} \mathbf{c}^{\langle m \rangle} = \mathbf{z} \mathbf{S}^{-1} \mathbf{u}^{\langle m \rangle} = \mathbf{g}^{\mathbf{c}} (t\_{m}), \tag{23}$$

where **z** = [*S*¯ <sup>0</sup>(*L*), *S*¯ <sup>1</sup>(*L*), *S*¯ <sup>2</sup>(*L*), ..., *S*¯ *<sup>N</sup>*−1(*L*)] and each *<sup>S</sup>*¯ *<sup>n</sup>*(*L*) is as defined in Lemma 1(iii). Finally, we can establish the following system of *m*th iterative linear equations for the inverse TVIDE problem (19) by utilizing (20) and (23), which has *N* + *n* + 1 unknown variables including **u***m*, **d**, and *βm*, as

$$
\begin{bmatrix}
\mathbf{H}^{\langle m\rangle} & \mathbf{X}\_{\text{ll}} & -\tau \mathbf{A}^{m} \mathbf{f}^{\langle m\rangle} \\
\mathbf{B} \mathbf{S}^{-1} & \mathbf{0} & \mathbf{0} \\
\mathbf{z} \mathbf{S}^{-1} & \mathbf{0} & \mathbf{0}
\end{bmatrix}
\begin{bmatrix}
\mathbf{u}^{\langle m\rangle} \\
\mathbf{d} \\
\boldsymbol{\mathcal{J}}^{\langle m\rangle}
\end{bmatrix} = \begin{bmatrix}
\mathbf{E}\_{2}^{\langle m\rangle} \\
\mathbf{F}^{\langle m\rangle} \\
\mathbf{g}^{c}(t\_{m})
\end{bmatrix},\tag{24}
$$

where **H***m* is the coefficient matrix of **u***m* defined in (20) and **E***m* <sup>2</sup> is the right-hand side column vector of (20). Before seeking an approximate solution **u***m* and source term *βm*, as we have mentioned, we must address that our inverse problem is ill-posed. When a noisy value is input into the system, it may cause a significant error. Hence, we need to stabilize the solution of (24) by employing the Tikhonov regularization method. We denote the linear system (24) by the simplified matrix equation as

$$\mathbf{R}\mathbf{y} = \mathbf{b}^c.\tag{25}$$

Applying the Tikhonov regularization method (6) in order to filter out the noise in the corresponding perturbed data, we can stabilize the numerical solution (25) by using (7). Thus, we have

$$\mathbf{y}\_{\lambda} = (\mathbf{R}^{\top}\mathbf{R} + \lambda\mathbf{I})^{-1}\mathbf{R}^{\top}\mathbf{b}^{c}.\tag{26}$$

Finally, we can receive the optimal regularization parameter *λ* by using Morozov's discrepancy principle combined with Newton's method, as described in Section 2.2. Thus, we can directly obtain the corresponding regularized solution by (26).
