Introduction to Numerical Approaches for Forward and Inverse Heat Transfer Problems

Abstract

We describe the formulation of forward and inverse problems in heat transfer and illustrate several techniques applicable to their solution.

1. Introduction

Heat flow is described by the heat (or diffusion) partial differential equation, which takes the form $\frac{\partial u}{\partial t}=k\frac{{\partial }^{2}u}{\partial x^2}$ with $u=u\left(x,t\right)$ representing the temperature at time $t$ and at spatial location $x$ (which may have one or more components). The constant $k$ is the ratio of the thermal conductivity in the material over the specific heat times the density of the material. Applying the conservation law principle to the heat flow across the rod yields the parabolic partial differential equation. If a heat source is present or the material is exothermic, a non-homogeneous version of the problem: $\frac{\partial u}{\partial t}-k\frac{{\partial }^{2}u}{\partial x^2}=f\left(x,t\right)$ with $f\left(x,t\right)\ne 0$ applies. The initial boundary value problem comes from assigning relevant initial and boundary conditions. For example, for a perfectly insulated rod as in Figure 1, the conditions could be $u\left(0,t\right)=u\left(L,t\right)=0$ . Given some initial profile of the temperature at time $t=0$ , $u\left(x,0\right)=h\left(x\right)$ , the analytical solution in 1D or 2D with a regular geometrical domain is given by separation of variables. Assuming $u\left(x,t\right)=F\left(x\right)G\left(t\right)$ , one obtains a system of two ODEs and from that, the series solution $u\left(x,t\right)={\sum }_{n=1}^{\infty }{u}_n\left(x,t\right)={\sum }_{n=1}^{\infty }{A}_{n}exp\left(-{\lambda }_n^{2}t\right)sin\left(\frac{\mathit{n\pi }}{L}x\right)$ with the coefficients $A_n$ satisfying the initial condition $u\left(x,0\right)={\sum }_{n=1}^{\infty }{A}_{n}sin\left(\frac{\mathit{n\pi }}{L}x\right)=f\left(x\right)$ and the eigenvalues ${\lambda }_n=\frac{\mathit{nk\pi }}{L}$ . The coefficients are evaluated via the sin series integration: $A_n=\frac{2}{L}{\int }_0^{L}f\left(x\right)sin\left(\frac{\mathit{k\pi x}}{L}\right)\mathit{dx}$ . For non-intregrable cases or more complicated geometries, numerical techniques for approximation are used [1].

Figure 1.

Perfectly insulated rod used a model of heat transfer, an exothermic slab, and sample modeled heat profile in 2D.

1.1 Derivation of homogeneous and non-homogeneous equations

We now discuss two simple scenarios which can be used to derive homogeneous and non-homogeneous forms of the heat Equation [2]. Heat is the flow of energy from a warmer to a cooler location. There are several types of heat transfer including conduction (flow of heat through stationary material), convection (flow of heat through fluids), and radiation (flow of heat through electromagnetic waves) [2]. In this chapter we are concerned mainly with heat conduction for which the differential heat model applies. We take $Q$ to represent the heat transfer rate in units of Watts. It is the rate at which heat is flowing through one of the slabs in Figure 2. When heat flows into and out of a mass, the temperature of the mass changes in accordance with the model: $Q_{\text{in}}-Q_{\text{out}}=\mathit{mc}\frac{\mathit{dT}}{\mathit{dt}}$ , where $T$ represents the temperature, $m$ the mass, and $c$ is the specific heat value of the material. It should be noted that in heat problems, $T$ and $u$ are ofter used interchangeably. We define $q=Q/A$ (heat over area) as the heat flux, which by Fourier’s law is proportional to the temperature gradient: $q=-\overline{k}\frac{\mathit{dT}}{\mathit{dx}}$ , where $\overline{k}$ is the thermal conductivity of the material [2]. Using $Q_{\text{in}}=Q\left(x\right)$ and $Q_{\text{out}}=Q\left(x+\mathit{\delta x}\right)$ along with heat flux $q=Q/A$ , we have $\mathit{Aq}\left(x\right)-\mathit{Aq}\left(x+\mathit{\delta x}\right)=\mathit{mc}\partial T/\partial t$ . Applying Fourier’s law yields:

Figure 2.

Heat conducting slabs without and with internal heat generation.

$$\overline{k}A\left[-{\left.\frac{\partial T}{\partial x}\right|}_x+{\left.\frac{\partial T}{\partial x}\right|}_{x+\mathit{\delta x}}\right]=\mathit{mc}\partial T/\partial t$$

Next, using Taylor’s approximation ${\left.\frac{\partial T}{\partial x}\right|}_{x+\mathit{\delta x}}\approx {\left.\frac{\partial T}{\partial x}\right|}_x+{\left.\frac{{\partial }^{2}T}{\partial x^2}\right|}_x\mathit{\delta x}$ and expressing $m=\mathit{\text{ρAδx}}$ in terms of slab density, area, and width, we obtain:

$$\overline{k}A\frac{{\partial }^{2}T}{\partial x^2}=\mathit{\rho c}\frac{\partial T}{\partial t}\Rightarrow k\frac{{\partial }^{2}T}{\partial x^2}=\frac{\partial T}{\partial t}$$

For the exothermic slab case exhibited in Figure 2 on the right, we obtain an inhomogeneous equation. We assume the slab is filled with a material which emits heat at a rate $\dot{q}$ with units $W/m^3$ (energy per unit volume). In this case, the in and out energy balance gives $Q\left(x\right)+\dot{q}V-Q\left(x+\mathit{\delta x}\right)=\mathit{mc}\partial T/\partial t$ , which with $q=Q/A,V=\mathit{A\delta x}$ , Fourier’s law and Taylor’s expansion, gives $\overline{k}A\frac{{\partial }^{2}T}{\partial x^2}\mathit{\delta x}+\dot{q}\mathit{A\delta x}=\mathit{mc}\frac{\partial T}{\partial t}\Rightarrow \overline{k}\frac{{\partial }^{2}T}{\partial x^2}+\dot{q}=\mathit{\rho c}\frac{\partial T}{\partial t}\Rightarrow k\frac{{\partial }^{2}T}{\partial x^2}+\frac{\dot{q}}{\mathit{\rho c}}=\frac{\partial T}{\partial t}$ . In the case of steady state, $\partial T/\partial t=0$ and a second order differential equation model results.

1.2 Numerical schemes

Discretization of the derivative operators yields the finite difference formulation. A number of different schemes are possible based on the difference schemes used to replace the derivatives in the equation resulting in either an explicit method (where the update for $u$ at the next time step is given by a direct calculation from the previous one) or an implicit method (where the update for $u$ is accomplished via a linear solve). For 2D heat transfer, the model becomes ${\alpha }^2\left(\frac{{\partial }^{2}u}{\partial x^2}+\frac{{\partial }^{2}u}{\partial y^2}\right)=\frac{\partial u}{\partial t}$ and again, the finite difference method can be used for discretization. The initial value problems are either forward or inverse, in the sense that they either predict the temperature profile at later time from the knowledge of the initial temperature distribution at time $t=0$ or in the inverse case, predict the temperature distribution at time $t=0$ from knowledge of the temperature at a future time. Sometimes, the initial boundary value problem may be prescribed extra conditions, such as Neumann or other derivate condition type and the inverse problem could instead seek the initial values for such expressions from future data [3]. In general, the inverse problem is more challening than the forward problem due to being ill-posed. The discretized operator matrix which can be assembled from the difference formulation is not well conditioned and as a result the problem is overly sensitive to slight changes in input values. Regularization techniques need to be used for the solution [4].

2. Forward problem

The forward problem is concerned by the prediction of heat profile in future time from the initial condition. We first consider the heat equation initial boundary value problem in 1D:

$$u_t={\mathit{ku}}_{\mathit{xx}},0<x<1,$$ E1

$$u\left(x,0\right)=f\left(x\right),0<x<1,$$ E2

$$u\left(0,t\right)=u\left(1,t\right)=0,t>0$$ E3

with constant $k\in \mathrm{ℝ}$ . Via a transformation, this problem can be extended to nonzero boundary conditions. For (3), separation of variables gives the solution as the superposition of solutions of the form:

$$\tilde{u}\left(x,t\right)=e^{-{\pi }^2\mathit{kt}}sin\left(\mathit{\pi x}\right)$$

A quick calculation shows that $\left[\frac{\partial }{\partial t}-k\frac{{\partial }^2}{\partial x^2}\right]\tilde{u}\left(x,t\right)=0$ . In order to form the well-known explicit scheme, we use the forward in time and central in space difference operators:

$${\frac{\partial u}{\partial t}}_{\mathrm{fd}}\left(x_i,t_j\right)=\frac{u\left(x_i,t_{j+1}\right)-u\left(x_i,t_j\right)}{l}$$ E4

$${\frac{{\partial }^{2}u}{\partial x^2}}_{\mathrm{fd}}\left(x_i,t_j\right)=\frac{u\left(x_{i+1},t_j\right)-2u\left(x_i,t_j\right)+u\left(x_{i-1},t_j\right)}{h^2}$$ E5

which results in the explicit finite difference scheme:

$$u_{i,j+1}=u_{i,j}+\frac{\mathit{kl}}{h^2}\left(u_{i+1,j}-2u_{i,j}+u_{i-1,j}\right)$$ E6

This can be re-written as $u_{i,j+1}=\left(1-2\alpha \right)u_{i,j}+\alpha \left(u_{i+1,j}+u_{i-1,j}\right)$ with $\alpha =\frac{\mathit{kl}}{h^2}$ , with $l$ the time step ( $\mathrm{\Delta }t$ ) and $h$ the spatial step size ( $\mathrm{\Delta }x$ ). Von Neumann stability analysis yields the condition $l\le \frac{h^2}{2k}$ . For non-homogeneous cases the general model is $\frac{\partial u}{\partial t}=k\frac{{\partial }^{2}u}{\partial x^2}+f\left(x,t\right)$ and the finite difference analysis proceeds in a similar way, with an explicit scheme of the form: $u_{i,j+1}=u_{i,j}+\frac{\mathit{kl}}{/}{h}^2\left(u_{i+1,j}-2u_{i,j}+u_{i-1,j}\right)+\mathit{lf}\left(x_i,t_j\right)$ .

If we use the backward time finite difference approximation for $\frac{\partial u}{\partial t}$ , we get: $\frac{u_{i,j+1}-u_{i,j}}{k}=\frac{u_{i+1,j+1}-2u_{i,j+1}+u_{i-1,j+1}}{h^2}$ , yielding the implicit scheme: $-\alpha u_{i-1,j+1}+\left(1+2\alpha \right)u_{i,j+1}-\alpha u_{i+1,j+1}=u_{i,j}$ which can be solved for time $j+1$ from the information at time $j$ using a linear system with a tridiagonal matrix with elements $\left(-\alpha ,\left(1+2\alpha \right),-\alpha \right)$ . The implicit scheme is numerically stable for any $\alpha$ but is costlier than the explicit scheme. Another example, is the so called Crank–Nicolson scheme based on the time derivative approximation at the midpoint of the time stamp. It takes the form:

$$\frac{u_{i,j+1}-u_{i,j}}{l}=k\frac{\left(u_{i+1,j+1}-2u_{i,j+1}+u_{i-1,j+1}\right)+\left(u_{i+1,j}-2u_{i,j}+u_{i-1,j}\right)}{2h^2}$$

Again using $\alpha =\frac{\mathit{kl}}{h^2}$ , we arrive at the scheme:

$$-\alpha u_{i+1,j+1}+2\left(1+\alpha \right)u_{i,j+1}-\alpha u_{i-1,j+1}=\alpha u_{i+1,j}+2\left(1-\alpha \right)u_{i,j}+\alpha u_{i-1,j}$$

Note that this takes the form of the linear system $\mathit{Au}=b$ , where $A$ is a tridiagonal matrix of terms $\left(-\alpha ,2\left(1+\alpha \right),-\alpha \right)$ , as before, $u$ containing the heat distribution at time $j+1$ and $b$ containing the heat distribution information at time $j$ . A linear system solves then yields the new time information from previous data.

In 2D, the standard explicit scheme arises from the discretization:

$$\frac{T_{i,j}^{n+1}-T_{i,j}^n}{\mathrm{\Delta }t}=k\left(\frac{T_{i,j+1}^n-2T_{i,j}^n+T_{i,j-1}^n}{{\left(\mathrm{\Delta }x\right)}^2}+\frac{T_{i+1,j}^n-2T_{i,j}^n+T_{i-1,j}^n}{{\left(\mathrm{\Delta }y\right)}^2}\right)$$

Like in the 1D scheme, suitable finite difference approximations yield the implicit method.

A practically useful direction to consider is the use of non-standard finite difference (nsfd) schemes [5]. Based on (5), we define the discretized operators:

$${\overline{d}}_{t}f\left(x,t\right)=f\left(x,t+l\right)-f\left(x,t\right)$$ E7

$${\overline{d}}_x^{2}f\left(x,t\right)=f\left(x+h,t\right)-2f\left(x,t\right)+f\left(x-h,t\right)$$ E8

and proceed to evaluate $\left({\overline{d}}_t-\mathit{kl}\frac{{{\overline{d}}^2}_x}{h^2}\right)\tilde{u}\left(x,t\right)$ for $\tilde{u}\left(x,t\right)=e^{-{\pi }^2\mathit{kt}}sin\left(\mathit{\pi x}\right)$ . The standard fd scheme does not set this expression to zero, and hence this analytically valid solution does not satisfy the original partial differential equation. However, based on the calculations below, we can make a replacement for $\frac{\mathit{kl}}{h^2}$ , which would force the equation to zero.

$${\overline{d}}_t\left[\tilde{u}\left(x,t\right)\right]=sin\left(\mathit{\pi x}\right)\left[e^{-{\pi }^{2}k\left(t+l\right)}-e^{-{\pi }^2\mathit{kt}}\right]=sin\left(\mathit{\pi x}\right)e^{-{\pi }^2\mathit{kt}}\left(e^{-{\pi }^2\mathit{kl}}-1\right)=\tilde{u}\left(x,t\right)\left(e^{-{\pi }^2\mathit{kl}}-1\right)$$ E9

$${\overline{d}}_x^2\left[\tilde{u}\left(x,t\right)\right]=e^{-{\pi }^2\mathit{kt}}\left[sin\left(\pi \left(x+h\right)\right)-2sin\left(\mathit{\pi x}\right)+sin\left(\pi \left(x-h\right)\right)\right]$$ E10

$$=e^{-{\pi }^2\mathit{kt}}\left[\frac{1}{2i}\left(e^{\mathit{i\pi }\left(x+h\right)}-e^{-\mathit{i\pi }\left(x+h\right)}\right)-2sin\left(\mathit{\pi x}\right)+\frac{1}{2i}\left(e^{\mathit{i\pi }\left(x-h\right)}-e^{-\mathit{i\pi }\left(x-h\right)}\right)\right]$$ E11

$$=e^{-{\pi }^2\mathit{kt}}\frac{1}{2i}\left[e^{\mathit{i\pi x}}{e}^{\mathit{i\pi h}}-e^{-\mathit{i\pi x}}{e}^{-\mathit{i\pi h}}-4isin\left(\mathit{\pi x}\right)+e^{\mathit{i\pi x}}{e}^{-\mathit{i\pi h}}-e^{-\mathit{i\pi x}}{e}^{\mathit{i\pi h}}\right]$$ E12

$$=e^{-{\pi }^2\mathit{kt}}\frac{1}{2i}\left[e^{\mathit{i\pi x}}\left(e^{\mathit{i\pi h}}+e^{-\mathit{i\pi h}}\right)-e^{-\mathit{i\pi x}}\left(e^{-\mathit{i\pi h}}+e^{\mathit{i\pi h}}\right)-4isin\left(\mathit{\pi x}\right)\right]$$ E13

$$=e^{-{\pi }^2\mathit{kt}}\frac{1}{2i}\left[2cos\left(\mathit{\pi h}\right)\left(e^{\mathit{i\pi x}}-e^{-\mathit{i\pi x}}\right)-4isin\left(\mathit{\pi x}\right)\right]$$ E14

$$=e^{-{\pi }^2\mathit{kt}}\frac{1}{2i}\left[4icos\left(\mathit{\pi h}\right)sin\left(\mathit{\pi x}\right)-4isin\left(\mathit{\pi x}\right)\right]$$ E15

$$=2e^{-{\pi }^2\mathit{kt}}sin\left(\mathit{\pi x}\right)\left[cos\left(\mathit{\pi h}\right)-1\right]=2\tilde{u}\left(x,t\right)\left[cos\left(\mathit{\pi h}\right)-1\right]$$ E16

Setting $\left({\overline{d}}_t-\frac{\mathit{kl}}{h^2}{{\overline{d}}^2}_x\right)\tilde{u}\left(x,t\right)=0$ , we obtain:

$$\tilde{u}\left(x,t\right)\left[e^{-{\pi }^2\mathit{kl}}-1\right]-\frac{2\mathit{kl}}{h^2}\tilde{u}\left(x,t\right)\left[cos\left(\mathit{\pi h}\right)-1\right]=0\Rightarrow \left[e^{-{\pi }^2\mathit{kl}}-1\right]=\frac{2\mathit{kl}}{h^2}\left[cos\left(\mathit{\pi h}\right)-1\right]$$

Hence, the nsfd scheme follows if we make the following replacement in (6):

$$\frac{\mathit{kl}}{h^2}=\frac{1}{2}\left[\frac{e^{-{\pi }^2\mathit{kl}}-1}{cos\left(\mathit{\pi h}\right)-1}\right]$$

When $h,l\ll 1$ , the nsfd substitution coincides with the standard finite difference scheme. We use the small $\mid q\mid$ Taylor series approximations $e^q\approx 1+q$ and $sin\left(q\right)\approx q$ to obtain:

$$\frac{1}{2}\left[\frac{e^{-{\pi }^2\mathit{kl}}-1}{cos\left(\mathit{\pi h}\right)-1}\right]=\left[\frac{e^{-{\pi }^2\mathit{kl}}-1}{-2{sin}^2\left(\frac{\mathit{\pi h}}{2}\right)}\right]\approx \frac{1}{2}\left[\frac{1-{\pi }^2\mathit{kl}-1}{-\frac{{\pi }^2{h}^2}{2}}\right]=\frac{\mathit{kl}}{h^2}$$

By this it follows that for small $h,l$ , the stability constraint $\frac{\mathit{kl}}{h^2}\le 1$ is satisfied in the nsfd scheme when $l\le \frac{h^2}{2}$ just like in the original scheme. This also holds for large $h$ , as long as $h$ is chosen so that $\mid cos\left(\mathit{\pi h}\right)-1\mid$ is not very small.

2.1 Numerical experiments

Below, we present some sample results for the case of constant $k\in \mathrm{ℝ}$ in (3), where we found that the described nsfd scheme often outperforms the fd scheme at the same time step. In Figure 3, we illustrate the solutions with fd and nsfd methods with the initial condition, $u\left(x,0\right)=sin\left(\mathit{\pi x}\right)+x$ for problem 1, with boundary conditions $u\left(0,t\right)=0$ and $u\left(1,t\right)=1$ and with the initial condition $u\left(x,0\right)=12sin\left(9\mathit{\pi x}\right)-7sin\left(4\mathit{\pi x}\right)$ and zero homogeneous boundary conditions for problem 2. For the second example, errors vs. step size are also compared. In general, we observe more accurate results with the use of the nsfd scheme.

Figure 3.

Comparison of fd and nsfd scheme solutions for the 1-D heat equation for problems 1 and 2 above.

An interesting case is that of heat transfer with an exothermic material, which can be modeled with a heat source distributed throughout [3]. Considering an infinitesimal segment as in Figure 1 yields the energy balance equation: $\dot{q}+\mathit{\alpha dx}-\left(\overline{q}+\frac{d\dot{q}}{\mathit{dx}}\mathit{dx}\right)=0\Rightarrow \frac{d\dot{q}}{\mathit{dx}}=\alpha$ . It follows that the heat model takes the form $\frac{1}{\alpha }\frac{\partial u}{\partial t}=\frac{{\partial }^{2}u}{\partial x^2}+\frac{\dot{q}}{k}$ . The steady state behavior does not depend on time so that $\frac{\partial u}{\partial t}=0$ and an ordinary differential equation model $\frac{d^{2}u}{{\mathit{dx}}^2}+\frac{\alpha }{k}=0$ results. Upon a double integration, the general solution is $u\left(x\right)=-\frac{\alpha }{2k}{x}^2+\mathit{Ax}+B$ , which with the boundary conditions $T\left(0\right)=T\left(L\right)=T_w$ becomes $u\left(x\right)=-\frac{\alpha }{2k}{x}^2+\mathit{Lx}+T_w$ , a parabola. We can simulate this steady state example with a finite difference scheme for the time derivative, yielding the temperature profile shown in Figure 4.

Figure 4.

Temperature distribution inside a thin slab with an exothermic material with homogeneous boundary conditions as simulated by an fd scheme.

3. Inverse problem

The inverse problem consists of predicting the initial condition from the observation of the heat profile at some later time $t$ . In general, the problem is ill-posed and regularization is necessary for its solution [4]. Let us look at a simple example in 1-D. Starting from the initial boundary value problem on $x\in \left[0,L\right]$ :

$$\frac{\partial u}{\partial t}=k\frac{{\partial }^{2}u}{\partial x^2}$$ E17

$$u\left(x,0\right)=f\left(x\right),u\left(0,t\right)=0=u\left(L,t\right)$$ E18

In this case, separation of variables can be used and after a well-known sequence of steps, the following series solution can be derived [3]:

$$u\left(x,t\right)=\sum _{n=1}^{\infty }{A}_{n}sin\left(\frac{\mathit{n\pi x}}{L}\right)exp\left(-k{\left(\frac{\mathit{n\pi }}{L}\right)}^{2}t\right)$$

It follows that $u\left(x,0\right)={\sum }_{n=1}^{\infty }{A}_{n}sin\left(\frac{\mathit{n\pi x}}{L}\right)$ which is a sin series; that is, $A_n=F\left[u\left(x_0\right)\right]=\frac{2}{L}{\int }_0^{L}u\left(x,0\right)sin\left(\frac{\mathit{n\pi x}}{L}\right)$ . It follows that the sin coefficients of the solution at $t=T$ are given by the sin coefficients at $t=0$ multiplied by $exp\left(-k{\left(\frac{\mathit{n\pi }}{L}\right)}^{2}T\right)$ . The values of $u$ at $t=T$ can then be computed via the inverse sin transform applied to the product. In a similar way, it follows that from the transform of the solution at $t=T$ , $F\left[u\left(x,T\right)\right]$ , we can in principle go backwards, to obtain $u\left(x,0\right)=F^{-1}\left[F\left[u\left(x,T\right)\right].\ast exp\left(k{\left(\frac{\mathit{n\pi }}{L}\right)}^{2}T\right)\right]$ with the $.\ast$ representing component-wise operation [3]. However, with the negative sign inside the exponential gone, there is potential for blow up of any large magnitude terms. The issue can be addressed by performing filtering on the large coefficients, setting the largest magnitude portion to zero or some maximum value. The magnitude cut off value can be selected based on quartile statistics of the nonzero absolute values.

A more general approach for inversion is to utilize the inverse problem formulations of the matrix heat propagation formulation. For example, for the standard implicit or Crank-Nicholson scheme, we have the system $\mathit{Au}=b$ , where $A$ is a tridiagonal matrix of terms $\left(-\alpha ,2\left(1+\alpha \right),-\alpha \right)$ . This means that $u^1=A^{-1}b\left(u^0\right)$ , $u^2=A^{-1}b\left(u^1\right)=A^{-2}b\left(u^0\right)\Rightarrow b\left(u^0\right)=A^2{u}^2$ . This then implies that the initial condition at time $t=0$ , $u^0$ can be recovered from the application of the inverse matrix raised to a power $n_t$ corresponding to the time point. As another example, if we discretize $u_t={\mathit{ku}}_{\mathit{xx}}$ as:

$$\frac{\partial u\left(t\right)}{\partial t}-\frac{k}{h^2}\left[u_{i+1}^{j+1}\left(t\right)-2u_i^{j+1}\left(t\right)+u_{i-1}^{j+1}\left(t\right)\right]=0=\frac{u_i^{j+1}-u_i^j}{\mathrm{\Delta }t}-{\mathit{Mu}}^{j+1}\left(t\right),$$

with $M=\frac{k}{h^2}\text{tridiag}\left(-1,2,-1\right)$ , for $j=0,1,\dots ,n_t-1$ . Using $u\left(x,0\right)=U_0$ , we obtain the relations:

$$u^0=U_0;u^{j+1}={\left(I+\mathrm{\Delta }\mathit{tM}\right)}^{-1}{u}^j={\left(I+\mathrm{\Delta }\mathit{tM}\right)}^{-\left(j+1\right)}{U}_0\Rightarrow {\overline{L}}^{-\left(j+1\right)}{U}_0=u^{j+1}$$

This means that $U_0$ can in theory be obtained right from ${\left(I+\mathrm{\Delta }\mathit{tM}\right)}^{n_t}{u}^{n_t}$ . The matrix needs not be formed explicitly, as the finite differencing scheme can be repeatedly applied to $u^{n_t}=u\left(T\right)$ . However, as shown in Figure 5, raising the typical differencing matrix to a power results in non-linear singular value decay. Hence, the application of the inverse of matrix ${\left(I+\mathrm{\Delta }\mathit{tM}\right)}^{-n_t}$ to the approximant for $b=u\left(T\right)$ , the heat profile at time $T$ , would result in the blow-up of any numerical errors due to the ill-conditioning. We must instead consider the inverse formulation ${\mathit{LU}}_0=b$ with $L={\overline{L}}^{-n_t}$ and utilize a regularization scheme.

Figure 5.

Tridiagonal structure of differencing operator matrix and singular value decay of operator raised to a power.

This can be accomplished by replacing the naive solution $L^{-1}b$ by the regularized formulation ${min}_u\left(\parallel \mathit{Lu}-b{\parallel }_2^2+\lambda \parallel u{\parallel }_2^2\right)$ , which for large enough $\lambda >0$ , alleviates the problem due to the fast decay of singular values. This well-known Tikhonov regularization formulation essentially acts to remove the influence of small singular values on the solution, which would otherwise blow up in value upon inversion. The regularization parameter $\lambda$ should be chosen with care to provide regularization but not to move too far away from the naive inverse operator. For this reason, $\lambda$ is usually selected based on an L-curve criterion [6], whereby the parameter is initially set at some fraction of $\parallel L^{T}b{\parallel }_{\infty }$ and then iteratively lowered to a lower fraction (with the values either linearly or logarithmically spaced). At each $\lambda$ value, the solution to $\left(L^{T}L+\mathit{\lambda I}\right)u_{\lambda }=L^{t}b$ is obtained (e.g. with a Conjugate Gradient scheme) using the solution at the previously used $\lambda$ as the initial guess. Typically, the progression to lower $\lambda$ stops once the value of $\parallel {\mathit{Lu}}_{\lambda }-b{\parallel }_2/N$ reaches a plateau. It is possible to employ a predictor-corrector scheme to iteratively improve the solution. In this approach, the approximation to $U_0$ can be obtained as above, then the forward operator can be applied to propagate this approximated initial heat profile at $t=0$ to time $t=T$ and compare with observations [4]. We can compare the residual at time $T$ and attempt some corrections to $U_0$ based on this residual such that the propagated solution at $t=T$ best matches observations. In the case of only a finite subset of observations, interpolation techniques can be utilized to complete the profile at $t=T$ across all values of $x$ [7] prior to regularization to obtain an approximate smooth and continuous $u\left(T\right)$ heat profile.

3.1 Numerical experiments

As an example, consider the 1D cosine based initial heat profile with $k=1$ and zero boundary conditions. The profile and subsequent integration results are shown in Figure 6. It is clear that the temperature in the rod will approach zero at all locations with increasing time. The inverse problem consists of approximately recovering the initial profile $u\left(x,0\right)=0.5cos\left(2\pi /65x\right)+1)$ , with support on $\left[N/4,3N/4\right]$ . For this, we take the sin transform of the solution at $t=T$ and multiply by $exp\left(k{\left(\frac{\mathit{n\pi }}{L}\right)}^{2}T\right)$ . In Octave, this can be accomplished with the code:

Figure 6.

Homogeneous 1D profile and change of temperature with time.

fuT=fft (u (:, n)); % temperature profile at time T

R=exp (k.*k ′*T);

filter (R);

init_cond_approx=real (ifft (fuT.*R));

The term $R$ may contain large entries, which will blow up when the inverse fft is performed. For this reason, the filtering of small entries is necessary. This can be accomplished simply by replacing all large entries of $R$ with either zero or some max value; e.g. taking the cutoff to be all values above $1e6$ and replacing them either with 0 or $1e2$ . The cutoff can be estimated with standard statistical outlier techniques, by taking the quartiles $Q_1,Q_2,Q_3$ of the $\mid R\mid$ values, computing $\mathit{IQR}=Q_3-Q_1$ and filtering all entries of $R$ above $Q_3+\mathit{\text{αIQR}}$ with $\alpha \ge 2$ by setting them to the $Q_2$ value. The process is illustrated in Figure 7, for the same example as above with filtering for absolute value of entries above $1e1,1e3,1e5$ . It can be observed that filtering only the entries of $R$ above $1e5$ in absolute value leads to a good approximation of the initial condition. Tikhonov regularization can be utilized when it is possible to build up the finite difference propagation operator ( $L$ ), either as an explicit matrix, or as a function. An example is shown to the right of Figure 7, with an initial triangular temperature profile that is approximately recovered by a simple continuation scheme. In the same Figure, we also show the profiles of $u\left(T\right)$ both as observed and as propagated from the approximated initial values. The smoothing properties of the heat kernel [3] give rise to the observed transformation between $t=0$ and $t=T$ .

Figure 7.

Initial condition reconstruction with different filtering levels and illustration of inversion for initial condition using Tikhonov regularization with parameter tuned via the L-curve approach.

4. Conclusions

This article discusses numerical approaches for the heat equation, both for the forward and inverse problems with examples for homogeneous and non-homogeneous cases. For the forward heat propagation modeling problem, we have presented a simple non-standard finite difference formulation, which satisfies the differential model and appears to result in improved performance over the standard finite difference scheme with similar magnitude time and spatial step sizes. For the inverse problem, we have described approaches based on Fourier filtering and Tikhonov regularization.