Method of Lines and Finite Difference Schemes with the Exact Spectrum for Solution the Hyperbolic Heat Conduction Equation

Abstract. This paper is concerning with the 1-D initial–boundary value problem for the hyperbolic heat conduction equation. Numerical solutions are obtained using two discretizations methods – the finite difference scheme (FDS) and the difference scheme with the exact spectrum (FDSES). Hyperbolic heat conduction problem with boundary conditions of the third kind is solved by the spectral method. Method of lines and the Fourier method are considered for the time discretization. Finite difference schemes with central difference and exact spectrum are analyzed. A novel method for solving the discrete spectral problem is used. Special matrix with orthonormal eigenvectors is formed. Numerical results are obtained for steel quenching problem in the plate and in the sphere with holes. The hyperbolic heat conduction problem in the sphere with holes is reduced to the problem in the plate. Some examples and numerical results for two typical problems related to hyperbolic heat conduction equation are presented.


Introduction
Problems involving the hyperbolic heat conduction equation arise in modelling intensive steel quenching, laser pulse duration and other processes [2,4]. This paper is concerning with the corresponding 1-D initial-boundary value problem. Analytical and numerical solutions are obtained using the finite difference scheme (FDS) and a difference scheme with the exact spectrum (FDSES) [5].
The second order differential operator in space is approximated by a secondorder FDS on the uniform grid. We employ the method of lines for the discretization in time and solve the corresponding equations by continuous and discrete Fourier methods. The discrete spectral method leads to a transcendent equation for the eigenvalues of the FDS. We derive a new transcendent equation to compute the two largest eigenvalues and the corresponding eigenvectors.
We find that the number of distinct eigenvalues depends on a special parameter Q = Lσ 1 σ 2 /(σ 1 + σ 2 ), where σ 1 , σ 2 are the heat transfer coefficients in the boundary conditions and L is the length of the interval.
We define the difference scheme with the exact spectrum (FDSES) as follows. Using the eigenvalues D and eigenvectors P of the finite difference matrix A, we write the Jordan form of the matrix A = P DP T and replace the discrete eigenvalues on the diagonal of D with the first eigenvalues of the differential operator.
Some applications are presented in Section 6. Two problems for the hyperbolic heat conduction equation are considered, featuring boundary conditions of the first kind and of the third kind. Spectral problems for the difference and differential operators are solved. We also present analytical and numerical results for four related problems: boundary value problems for ordinary differential equations (arising in the stationary case) and initial-boundary value problems for the standard heat transfer process, wave and hyperbolic heat conduction equations. The advantages of the FDSES method for the case of the first kind boundary conditions are demonstrated by comparison with FDS for several numerical examples. Some examples and numerical experiments are presented. Numerical solutions in the time for these problems are obtained by the MATLAB solver "ode15s".

Mathematical Models
We consider the following 1-D hyperbolic heat conduction problem in a plate: wherek is the heat conductivity, t f is the final time, τ is the relaxation time (τ < 1), T l , T r , T 0 (x), f (x, t) are given temperatures and a source function, α 0 , α 1 are the heat transfer coefficients (for boundary conditions of the first kind α 0 = α 1 = ∞).
In a sphere with a hole defined in the spherical coordinates by 0 < r 0 < r < R, assuming radial symmetry, with homogeneous boundary conditions (BC) we obtain the following 1-D hyperbolic heat conduction problem: 2) where r is the polar radius, R, r 0 are the radius of the sphere and hole. The transformation V = T r, x = r − r 0 reduces the problem (2.2) to the following hyperbolic heat conduction problem in the plate: ∂x .

Solution of the Problem (2.1) with Homogeneous BC of the First Kind
If the surrounding temperature T l , T r in (2.1) is constant and α 0 = α 1 = ∞ then using the transformation we obtain a problem with homogeneous BCs of the first kind and V (x, 0) We consider the uniform grid in the space x k = kh, k = 0, N , N h = L. Using the standard finite difference approximation for partial derivatives of the second order with respect to x, we obtain an initial value problem for the system of ordinary differential equations (ODEs) of the second order: where A is the standard 3-diagonal matrix of order M = N −1 with the elements The corresponding discrete spectral problem has the following solution [6]: eigenvalues µ k = 4 h 2 sin 2 kπ 2N and the orthonormal eigenvectors p k with elements p k i = 2 N sin πik N , i, k = 1, M . From the spectral problem AP = P D it follows that the matrix A is represented in the form A = P DP, where P = P −1 is the symmetric orthogonal matrix with the columns being the eigenvectors of A, elementwise p i,k = p k i , i, k = 1, M , and D is the diagonal matrix with the elements d k,k = µ k , k = 1, M .
The solution of the spectral problem for the corresponding continuous (differential) problem where δ k,m is the Kronecker delta. For the discrete problem the integral in the scalar product (p k , p m ) is approximated by the trapezoidal formula.
The difference scheme with the exact spectrum (FDSES) is obtained by replacing A withÃ = PDP , whereD is the diagonal matrix with the first M eigenvalues λ k of the differential operator (− ∂ 2 ∂x 2 ) as elements, i.e.d k,k = λ k . Note thatÃ is a full matrix.
The FDSES method is more stable than the method of finite difference approximation with central difference (FDS), because the eigenvalues are larger λ k > µ k . We can construct analytical solutions of (3.1) using the spectral representation of matrix A in the form A = P DP. The transformation W = P U decouples the system of ODEs The solution of this system is given by where κ k = 0.25/τ 2 −kd k /τ . If 4kd k τ > 1, then the hyperbolic functions are replaced with the trigonometrical and the parameter κ k with k d k /τ −0.25/τ 2 .
If κ k = 0, then we have We can also construct analytical solutions by using the Fourier method in the following form: where p k (x) are the or-thonormal eigenfunctions for the differential operator (− ∂ 2 ∂x 2 ) with homogeneous boundary conditions, w k (t) is the solution of (3.3), with For a given function q(x) ∈ C n [0, L] the Fourier coefficients b k can be estimated The 3-diagonal matrix A can be represented by the difference operator [6] Ay = Using the scalar product of two vectors [y 1 , ) it can be proved, that the operator A is symmetric and [Ay, y] ≥ 0 [6].
The solution of the corresponding spectral problem can be given in the following form: where (y n , y m ) = L 0 y n (x)y m (x) dx = δ n,m and λ n are the positive roots of the transcendental equation: , n = 1, 2, 3, . . . .
The discrete spectral problem Ay n =µ n y n , n=1, . . . , N +1 has the solution [6] where p n are the positive roots of the following transcendental equation The constants C 2 n = [y n , y n ] can be obtained in the following form h sin(2p n L) cos(p n h) 4 sin(p n h) + σ 1 sin 2 (p n L) cos(p n h).
In the limit case it follows that µ n → λ 2 n as h → 0. The experiments with MATLAB (the spectral problem is solved with MAT-LAB solver "eig") show that the first roots p n , n = N − 1 belong to the interval ((n − 1)π, nπ), µ n ≤ 4 h 2 and can be obtained from (4.2), but the two last roots p N , p N +1 can not be obtained from this transcendental equation (µ n ≥ 4 h 2 , n = N, N + 1). In this case special solutions with µ ≥ 4 h 2 can be obtained in the form µ = 4 h 2 cosh 2 ph 2 . Then y j = (−1) j [C 1 sinh(px j ) + C 2 cosh(px j )]. The constants C 1 , C 2 are determined from the difference equations by j = 0, j = N and the values p N , p N +1 are obtained from the new transcendental equation. Depending on the fixed parameter Q = Lσ 1 σ 2 /(σ 1 + σ 2 ) we can obtain one (Q < 1) or two (Q ≥ 1) roots from the following transcendental equation: coth(p n L) = sinh 2 (p n h) + h 2 σ 1 σ 2 h(σ 1 + σ 2 ) sinh(p n h) , n = N, N + 1.
The expression Q = 1 is obtained from this equation in the limit case when p n → 0.
The corresponding eigenvalues and orthonormal eigenvectors are The constants C 2 n in this case are given by Then we have the orthonormal eigenvectors y n , y m for all n, m = 1, N + 1.
For Q = 1 the single eigenvalue is µ N = 4/h 2 , (p N = 0) and the corresponding components y N j of the orthonormal eigenvector y N are obtained as the limit (p N → 0) of the expression y N j /C N in the following form: The spectral problem for σ 1 = 0, σ 2 = α has the solution Then we have orthonormal eigenvectors y n , y m , where the discrete scalar product is defines as [y n , y m ] = δ n,m for all n, m = 1, . . . , N + 1. The solution of the corresponding spectral problem for the differential equation −y (x) = λ 2 y(x), x ∈ (0, L), y (0) = 0, y (L) + αy(L) = 0, is given in following form where (y n , y m ) = L 0 y n (x)y m (x) dx = δ n,m and λ n are the positive roots of the transcendental equation tan(λ n L)λ n = α, n = 1, 2, 3, . . .. In the scalar product [y n , y m ] the integral (y n , y m ) is approximated using the trapezoidal formula. In the limit case, as h → 0, we have p n → λ n for n = 1, . . . , N .

The Solution of the Problem (2.1) with Homogeneous BC of the Third Kind
Similarly we obtain from (2.1) for T l = T r = 0 the initial value problem for the system of ODEs (3.1), where the order of column-vectors is M = N + 1. The matrix A can be represented in the form A = P DP T , where the columns of the matrix P and the diagonal matrix D contain M orthonormal eigenvectors y n and eigenvalues µ n , n = 1, M . From P T P = E it follows that P −1 = P T . We can also use the Fourier method in the form T ( For the FDSES the matrix A is represented in the form A = P DP T , where the diagonal matrix D contains the first M eigenvalues d n = λ 2 n , n = 1, . . . , M of the differential operator (− ∂ 2 ∂x 2 ) respectively. It is possible to replace the column of the matrix P with the first eigenfunctions y n j = y n (x j ), n, j = 1, M of the continuous problem. Then the eigenvectors (y n , y m ) or column-vectors of the matrix P is not orthogonal in the discrete scalar product [y n , y m ] (P T P = E) and the analytic form (3.3) is not valid.
If the heat conductivityk(x) is a functions of x then we can obtain the discrete eigenvectors and eigenvalues using the MATLAB operator [P, D] = eig(A), where A is the matrix of the difference operator with second order approximation. The orthogonal eigenvectors p n should be normalized with respect to the scalar product [y n , y n ]. The eigenvalues of the differential problem can be obtained with the MATLAB solver "bvp4". However, this algorithm is of high computational complexity.
The FDSES method performs well for nonlinear problems, when the diffusion term in (2.1) is given in the form ∂ 2 ∂x 2 φ(T ). Then for dφ(T ) dT =k(x, t): In this case from (3.1) we have the nonlinear system of ODEs with the term Aφ(U ), where φ(U ) = (φ(u 1 (t)), φ(u 2 (t)), . . . , φ(u M (t))) T , A = P DP T . The nonlinear heat transfer equation with power functions φ(T ) is considered in [3]. Similarly we can consider the analytical solutions of (3.1) using the spectral representation of matrix A in the form A = P DP T . From the transformation W = P T U we obtain the decoupled system of ODEs (3.2), where W (0) = P T U 0 ,Ẇ (0) = P T V 0 , G(t) = P T F (t) are the column-vectors of N + 1 order. The solution of this system is given in the form (3.3). In this case, in order to compute the scalar product the first and last components of vectors U 0 , V 0 , F (t) are divided by √ 2, but w 1 (t), w N +1 (t) need to be multiplied by √ 2. Similarly we can consider the analytical solutions of (2.3).

Some Examples and Numerical Results
Boundary value problems for ODE with the homogeneous BC of the first kind (stationary problem). In the stationary case from (2.1) (k = 1) we have the boundary value problem with BCs of the first kind Then from (3.1) we have a system of linear algebraic equations AU = F, where U, F are the column-vectors with the elements u k ≈ V (x k ), f k = f (x k ), k = 1, M , M = N − 1. Using the transformation W = P U or U = P W we can decouple system (3.2) into the form DW = G, (G = P F ) with the solution We obtain FDS method by setting d k = µ k and FDSES method by If The boundary value problem is solved by f (x) = 12x 2 C 0 + σ 1 sin(x), C 0 = (σ 1 cos(L)+σ 1 σ 2 sin L+σ 2 )/(4L 3 +σ 2 L 4 ). The exact solution of the differential problem is u(x) = −x 4 C 0 + 1 + σ 1 sin(x).
The solution of the problem is obtained with MATLAB in the form U = In this case the first and last components of vector F are divided with √ 2, but elements of the solution vector U u 1 , u N +1 are multiplied by √ 2. In the Table 1 we can see the last eigenvalues p N , p N +1 , µ N , µ N +1 obtained from (4.1), (4.2) and the maximum norm of the error of the finite difference scheme (δ(FDS )) and FDSES (δ(FDSES )).

Conclusions
• The hyperbolic heat conduction problems can be numerically modelled using FDS and FDSES methods.
• The advantages of the FDSES method for the problem with constant heat conductivity and with the first kind boundary conditions (BC) are demonstrated via several numerical examples in comparison with the FDS method.
• In the case of the BC of the third kind the eigenvectors of the differential and discrete problems are different and the FDSES is obtained by using the orthonormal eigenvectors of the discrete problem.
• For nonlinear diffusion term ∂ 2 ∂x 2 φ(T ) it is possible to obtain the algorithm of FDSES by numerical solutions of a system of nonlinear ODEs.
• For the discrete Fourier method numerical oscillations are present, when the initial and boundary conditions are discontinuous. For larger N these oscillations disappear. The method of FDSES is free of oscillations.
• For the discrete spectral problem new transcendental equation and eigenvectors are obtained. All eigenvectors for the finite difference operator are obtained analytically.