Open access peer-reviewed chapter - ONLINE FIRST

A Numerical Approach to Solving an Inverse Heat Conduction Problem Using the Levenberg-Marquardt Algorithm

By Tao Min, Xing Chen, Yao Sun and Qiang Huang

Submitted: May 9th 2019Reviewed: August 9th 2019Published: October 7th 2019

DOI: 10.5772/intechopen.89096

Downloaded: 42

Abstract

This chapter is intended to provide a numerical algorithm involving the combined use of the Levenberg-Marquardt algorithm and the Galerkin finite element method for estimating the diffusion coefficient in an inverse heat conduction problem (IHCP). In the present study, the functional form of the diffusion coefficient is an unknown priori. The unknown diffusion coefficient is approximated by the polynomial form and the present numerical algorithm is employed to find the solution. Numerical experiments are presented to show the efficiency of the proposed method.

Keywords

  • parabolic equation
  • inverse problem
  • Levenberg-Marquardt

1. Introduction

The numerical solution of the inverse heat conduction problem (IHCP) requires to determine diffusion coefficient from an additional information. Inverse heat conduction problems have many applications in various branches of science and engineering, mechanical and chemical engineers, mathematicians and specialists in many other sciences branches are interested in inverse problems, each with different application in mind [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15].

In this work, we propose an algorithm for numerical solving an inverse heat conduction problem. The algorithm is based on the Galerkin finite element method and Levenberg-Marquardt algorithm [16, 17] in conjunction with the least-squares scheme. It is assumed that no prior information is available on the functional form of the unknown diffusion coefficient in the present study, thus, it is classified as the function estimation in inverse calculation. Run the numerical algorithm to solve the unknown diffusion coefficient which is approximated by the polynomial form. The Levenberg-Marquardt optimization is adopted to modify the estimated values.

The plan of this paper is as follows: in Section 2, we formulate a one-dimensional IHCP. In Section 3, the numerical algorithm is derived. Calculation of sensitivity coefficients will be discussed in Section 4. In order to discuss on some numerical aspects, two examples are given in Section 5. Section 6 ends this paper with a brief discussion on some numerical aspects.

2. Description of the problem

The mathematical formulation of a one-dimensional heat conduction problem is given as follows:

ut=xqxux+fxt,xt0L×0T,E1

with the initial condition

ux0=u0x,0xL,E2

and Dirichlet boundary conditions

u0t=g1t,0tT,E3
u1t=g2t,0tT,E4

where fxt, u0x, g1t, g2tand qxare continuous known functions. We consider the problem (1)(4) as a direct problem. As we all know, if u0x, g1t, g2tare continuous functions and qxis known, the problem (1)–(4) has a unique solution.

For the inverse problem, the diffusion coefficient qxis regarded as being unknown. In addition, an overspecified condition is also considered available. To estimate the unknown coefficient qx, the additional information on the boundary x=x0, 0<x0<Lis required. Let the uxttaken at x=x0over the time period 0Tbe denoted by

ux0t=gt0tT.E5

It is evident that for an unknown function qx, the problem (1)(4) is under-determined and we are forced to impose additional information (5) to provide a unique solution pair uxtqxto the inverse problem (1)(5).

We note that the measured overspecified condition ux0t=gtshould contain measurement errors. Therefore the inverse problem can be stated as follows: by utilizing the above-mentioned measured data, estimate the unknown function qx.

In this work the polynomial form is proposed for the unknown function qxbefore performing the inverse calculation. Therefore qxapproximated as

qxqx^=p1+p2x+p3x2++pm+1xm,E6

where p1,p2,,pm+1are constants which remain to be determined simultaneously. The unknown coefficients p1,p2,,pm+1can be determined by using least squares method. The error in the estimate

Fp1p2pm+1=i=1nux0tip1p2pm+1gti2,E7

is to be minimized. Here, ux0tip1p2pm+1are the calculated results. These quantities are determined from the solution of the direct problem which is given previously by using an approximated qx^for the exact qx. The estimated values of pj, j=1,2,,m+1are determined until the value of Fp1p2pm+1is minimum. Such a norm can be written as

FP=UPGTUPG,E8

where PT=p1p2pm+1denotes the vector of unknown parameters and the superscript Tabove denotes transpose. The vector UPGTis given by

UPGT=ux0t1pgt1u(x0t2p)gt2u(x0tnp)gtn.E9

FPis real-valued bounded function defined on a closed bounded domain DRm+1. The function FPmay have many local minimum in D, but it has only one global minimum. When FPand Dhave some attractive properties, for instance, FPis a differentiable concave function and Dis a convex region, then a local maximum and problem can be solved explicitly by mathematical programming methods.

3. Overview of the Levenberg-Marquardt method

The Levenberg-Marquardt method, originally devised for application to nonlinear parameter estimation problems, has also been successfully applied to the solution of linear ill-conditioned problems. Such a method was first derived by Levenberg (1944) by modifying the ordinary least-squares norm. Later Marquardt (1963) derived basically the same technique by using a different approach. Marquardt’s intention was to obtain a method that would tend to the Gauss method in the neighborhood of the minimum of the ordinary least-squares norm, and would tend to the steepest descent method in the neighborhood of the initial guess used for the iterative procedure.

To minimize the least squares norm (8), we need to equate to zero the derivatives of FPwith respect to each of the unknown parameters p1p2pm+1, that is

FPp1=FPp2==Fppm+1=0.E10

Let us introduce the sensitivity or Jacobian matrix, as follows:

JP=UTPPT=up1x0t1pup2x0t1pupm+1x0t1pup1x0t2pup2x0t2pupm+1x0t2pup1x0tnpup2x0tnpupm+1x0tnp,E11
orJij=upjx0tip=ux0tippj,i=1,2,n,j=1,2,m+1.E12

The elements of the sensitivity matrix are called the sensitivity coefficients, the results of differentiation (10) can be written down as follows:

2JTPUPG=0.E13

For linear inverse problem the sensitivity matrix is not a function of the unknown parameters. The Eq. (13) can be solved then in explicit form:

P=JTJ1JTG.E14

In the case of a nonlinear inverse problem, the matrix Jhas some functional dependence on the vector p. The solution of Eq. (13) requires an iterative procedure, which is obtained by linearizing the vector UPwith a Taylor series expansion around the current solution at iteration k. Such a linearization is given by

UP=UPk+JkPPk,E15

where UPkand Jkare the estimated temperatures and the sensitivity matrix evaluated at iteration k, respectively. Eq. (15) is substituted into (14) and the resulting expression is rearranged to yield the following iterative procedure to obtain the vector of unknown parameters P:

Pk+1=Pk+JkTJk1JkTGUPk.E16

The iterative procedure given by Eq. (16) is called the Gauss method. Such method is actually an approximation for the Newton (or Newton-Raphson) method. We note that Eq. (14), as well as the implementation of the iterative procedure given by Eq. (16), require the matrix JTJto be nonsingular, or

JTJ0,E17

where is the determinant.

Formula (17) gives the so called identifiability condition, that is, if the determinant of JTJis zero, or even very small, the parameters pj, for j=1,2,,m+1, cannot be determined by using the iterative procedure of Eq. (16).

Problems satisfying JTJ0are denoted ill-conditioned. Inverse heat transfer problems are generally very ill-conditioned, especially near the initial guess used for the unknown parameters, creating difficulties in the application of Eqs. (14) or (16). The Levenberg-Marquardt method alleviates such difficulties by utilizing an iterative procedure in the form:

Pk+1=Pk+JkTJk+μkΩk1JkTGUPk,E18

where μkis a positive scalar named damping parameter and Ωkis a diagonal matrix.

The purpose of the matrix term μkΩkis to damp oscillations and instabilities due to the ill-conditioned character of the problem, by making its components large as compared to those of JTJif necessary. μkis made large in the beginning of the iterations, since the problem is generally ill-conditioned in the region around the initial guess used for iterative procedure, which can be quite far from the exact parameters. With such an approach, the matrix JTJis not required to be non-singular in the beginning of iterations and the Levenberg-Marquardt method tends to the steepest descent method, that is, a very small step is taken in the negative gradient direction. The parameter μkis then gradually reduced as the iteration procedure advances to the solution of the parameter estimation problem, and then the Levenberg-Marquardt method tends to the Gauss method given by (16). The following criteria were suggested in literature [13] to stop the iterative procedure of the Levenberg-Marquardt method given by Eq. (18):

Fpk+1<ε1,E19
JkGUpk<ε2,E20
pk+1pk<ε3,E21

where ε1,ε2and ε3are user prescribed tolerances and denotes the Euclidean norm. The criterion given by Eq. (19) tests if the least squares norm is sufficiently small, which is expected in the neighborhood of the solution for the problem. Similarly, Eq. (20) checks if the norm of the gradient of Fpis sufficiently small, since it is expected to vanish at the point where Fpis minimum. The last criterion given by Eq. (21) results from the fact that changes in the vector of parameters are very small when the method has converged. Generally, these three stopping criteria need to be tested and the iterative procedure of the Levenberg-Marquardt method is stopped if any of them is satisfied.

Different versions of the Levenberg-Marquardt method can be found in the literature, depending on the choice of the diagonal matrix Ωkand on the form chosen for the variation of the damping parameter μk. In this paper, we choose the Ωkas

Ωk=diagJkTJk.E22

Suppose that the vector of temperature measurements G=gt1gt2gtnare given at times ti, i=1,2,,nand an initial guess P0is available for the vector of unknown parameters P. Choose a value for μ0, say, μ0=0.001and k=0. Then,

Step 1. Solve the direct problem (1)(4) with the available estimate Pkin order to obtain the vector UPk=ux0t1pku(x0t2pk)u(x0tnpk).

Step 2. Compute FPkfrom the Eq. (8).

Step 3. Compute the sensitivity matrix Jkfrom (12) and then the matrix Ωkfrom (22), by using the current value of Pk.

Step 4. Solve the following linear system of algebraic equations, obtained from (18): JkTJk+μkΩkΔPk=JkTGUPkin order to compute ΔPk=Pk+1Pk.

Step 5. Compute the new estimate Pk+1as Pk+1=Pk+ΔPk.

Step 6. Solve the exact problem (1)(4) with the new estimate Pk+1in order to find UPk+1. Then compute FPk+1.

Step 7. If FPk+1FPkreplace μkby 10μkand return to step 4.

Step 8. If FPk+1FPk, accept the new estimate Pk+1and replace μkby 0.1μk.

Step 9. Check the stopping criteria given by (19). Stop the iterative procedure if any of them is satisfied; otherwise, replace kby k+1and return to step 3.

4. Calculation of sensitivity coefficients

Generally, there have two approaches for determining the gradient; the first is a discretize-then-differentiate approach and the second is a differentiate-then-discretize approach.

The first approach is to approximate the gradient of the functional by a finite difference quotient approximation, but in general, we cannot determine the sensitivities exactly, so this method may led to larger error.

Here we intend to use differentiate-then-discretize approach which we refer to as the sensitivity equation method. This method can be determined more efficiently with the help of the sensitivities

uk=upk,k=1,2,m+1.E23

We first differentiate the flow system (1)(4) with respect to each of the design parameters p1p2pm+1, to obtain the m+1continuous sensitivity systems: for k=1,2,,m+1

ukt=xp1+p2x++pm+1xmukx+xk1uxukx0=0uk0t=0ukLt=0.E24

There have (m+2) equations, we can make them in one system equation and use the finite element methods to solve the system of equation. Here, we give the vector form of the equation as follow:

P1Ut+Γ=FUx0=U0xU0t=G1tULt=G2t,E25

where

U=uu1u2um+1,Γ=p1+p2x++pm+1xmuxp1+p2x++pm+1xmu1xuxp1+p2x++pm+1xmu2xxuxp1+p2x++pm+1xmum+1xxmux,F=fxt000.E26
U0x=u0x000T,G1t=g1t000T,G2t=g2t000T.E27

We use the Galerkin finite element method approximation for discretizing problem (25). For this, we multiply the Eq. (25) by a test function v:0LR, vV0H010Land integrate the obtained equation in space form 0to L. We obtain the following equation:

0LUxttvxdx0LΓvxdx=0LFxtvxdx,E28

integrating by parts gives

0LΓvxdx=Γvx0L0LΓvxtxdx.E29

We can change the first derivative in time and the integral. We have v0=0=vL, because vV0. This leads to an equivalent problem to P1: t>0,find Uxtsatisfying

ddt0LUxtvxdx+0LΓvxtxdx=0LFxtvxdx,E30

for all vV0H010L. To simplify the notation we use the scalar product in L20L

fg=0Lfxgxdx.E31

We also can define the following bilinear form:

aUv=0LΓvxtxdx=0Lp1+p2x++pm+1xmuxvxdx0Lp1+p2x++pm+1xmu1xvxuxvxdx0Lp1+p2x++pm+1xmu2xvxxuxvxdx0Lp1+p2x++pm+1xmum+1xvxxm+1uxvxdx.E32

Finally, we obtain with this notations the weak problem of P1:

P2ddtUvL2+aUv=FvL2Ux0=U0xU0t=G1tULt=G2t.E33

4.1 Space-discretization with the Galerkin method

In this section, we search a semi-discrete approximation of the weak problem P2, using the Galerkin finite element method. This leads to a first order Cauchy-problem in time.

Let Vhbe a Nx+1dimensional subspace of Vand V0,h=VhV0. Then the following problem is an approximation of the weak problem, find uh,u1,h,u2,hum+1.hVhthat satisfies:

ddtUhvh+aUhvh=FvhUhx0=U0,hxUh0t=G1tUhLt=G2t,E34

for all vhV0,h. where Uh=uhu1,hu2,hum+1.hT.

The choice of Vhis completely arbitrary. So we can choose it the way that for later treatment, it will be as easy as possible. For example, we subdivide the interval 0Linto partitions of equal distances h:

0=a1<a2<<aNx<aNx+1=L,,ai=i1hE35
Vh=vhC00L:vhai,ai+1P1i=1Nx,E36
V0,h=vhVh:vh0=vhL=0.E37

Note, that the finite dimension allows us to build a finite base for the corresponding space. In the case of V0,h, we have: φii=2Nxwhere i=2Nx.

φix=0x[a0,ai1]xhi2xai1aiixhxaiai+10xai+1aNx+1,E38

while we add for Vhthe two functions φ1and φNx+1defined as:

φ1x=1xh,ifxa1a20,ifxa2aNx+1,E39
φNx+1x=0,ifxa1aNxxhNx+1,ifxaNxaNx+1,E40

so that we can write Uhas a linear combination of the basic elements:

Uhxt=j=1Nx+1U˜jtφjx,E41
U0,hxt=j=1Nx+1U0xjφjx,E42

where U˜1t=G0tand U˜Nx+1t=G1t. Using that ais bilinear form and that Eq. (34) is valid for each element of the base φii=2Nx, we obtain

j=1Nx+1ddtU˜jtφjφi+j=1Nx+1U˜jtaφjφi=Fφi,E43

i=2Nx. This equation can be written in a vector form. For this we define the vectors u,u0and Fwith components

Fit=FφiL2,ujtu˜jt,u0,j=u0xj,E44

and matrices Mand Aas

mijφiφjL2,aijaφiφj,E45

Note that M,ARNx1×Nx+1,uRNx+1,andFRNx1. So that (43) is equal to the Cauchy problem

Mddxut+Aut=Ftut0=u0,E46

the Crank-Nicolson method can be applied to (46) at time tk, resulting in

Muk+1ukΔt+12Auk+1+12Auk=12Fk+Fk+1,E47

where uk=utk, Fk=Ftk,k=0,1,..

The Eq. (47) can be written in simple form as

M+Δt2Auk+1=MΔt2Auk+Δt2Fk+Fk+1,E48

the algebraic system (48) is solved by Gauss elimination method.

5. Numerical experiment

In this section, we are going to demonstrate some numerical results for uxtqxin the inverse problem (1)(5). Therefore the following examples are considered and the solution is obtained.

Example 1. Consider (1)(4) with

ux0=sinx,0x1,E49
u0t=0,0t1E50
u1t=sin1et,0t1,E51
fxt=sinxx24+x2+1etx+12cosxetsinxet,0x10t1E52

We obtain the unique exact solution

qx=1+0.5x+0.25x2E53

And

uxt=sinxet.E54

We take the observed data gas

gt=u0.5t=sin0.5et0t1E55

The unknown function qxdefined as the following form

qx^=p1+p2x+p3x2,E56

where p1,p2,p3are unknown coefficients.

Table 1 shows how the Levenberg-Marquardt algorithm can find the best parameters after 12 iterations when it is initialized in four different points.

Starting point0.5 0.5 0.51 1 110 10 1050 50 50
Iteration 120.999729028233135
0.499885876453067
0.252009862457275
0.999729028233183
0.499885876453056
0.252009862457315
0.999729028233194
0.499885876453057
0.252009862457325
0.999729028307261
0.499885876454169
0.25200986249336
Error F8.7564944405 × 10−148.7564944427 × 10−148.7564944420 × 10−148.7564944420 × 10−14

Table 1.

Performance of the algorithm when it is run to solve the model using three different parameters guesses.

Figures 14 show the fitness of the estimated parameters and the rate of convergence.

Figure 1.

All the initial values for the parameters are set 0.5.

Figure 2.

All the initial values for the parameters are set 1.

Figure 3.

All the initial values for the parameters are set 10.

Figure 4.

All the initial values for the parameters are set 50.

Figures 58 show the comparison between the inversion results q^xand the exact value qx:

Figure 5.

The comparison chart with all the initial values for the parameters is set 0.5.

Figure 6.

The comparison chart with all the initial values for the parameters is set 1.

Figure 7.

The comparison chart with all the initial values for the parameters is set 10.

Figure 8.

The comparison chart with all the initial values for the parameters is set 50.

Table 2 shows the values of qjΔxand ujΔx0.5in x=jΔxwith the all the initial values are set 1.

NumericalExactNumericalExact
jqjΔxqjΔxujΔx0.5ujΔx0.5
00.999729028233183100
11.052237714503061.05250.06055931902391730.0605520280601669
21.109786598022091.110.1205117976117860.120499040271796
31.172375678790261.17250.1792570590785210.179242065904716
41.240004956807581.240.2362074490803440.236194164064666
51.312674432074041.31250.2907939432508690.290786288212692
61.390384104589651.390.3424718283616250.342472971890064
71.473133974354411.47250.3907261148970890.390737778838824
81.560924041368311.560.4350766304105870.435098463062163
91.653754305631361.65250.4750827175305320.475111787267016
101.751624767143551.750.5103474067133680.510377951544573

Table 2.

The values of qjΔxand ujΔx0.5in x=jΔxwith the all the initial values being set to 1.

Example 2. Consider (1)(4) with

ux0=xex,0x1,E57
u0t=tet,0t1E58
u1t=1+te1t,0t1E59
fxt=et+xt+xet+xexet+xt+xet+x+ex2et+xt+xet+x,0x10t1E60

We obtain the unique exact solution

qx=ex,E61

And

uxt=x+text.E62

We take the observed data gas

gt=u0.5t=0.5+te0.5t0t1.E63

The unknown function qxdefined as the following form

qx^=p1+p2x+p3x2+p4x3+p5x4+p6x5+p7x6+p8x7, where p1,p2,,p7,p8are unknown coefficients.

Table 3 shows how the Levenberg-Marquardt algorithm can find the best parameters after 20 iterations when it is initialized in four different points.

Starting point0.1 0.1 0.1 0.1
0.1 0.1 0.1 0.1
0.5 0.5 0.5 0.5
0.5 0.5 0.5 0.5
1 1 1 1
1 1 1 1
2 2 2 2
2 2 2 2
Iteration 201.01536263526644
0.896348846894057
0.954285303464511
–0.890298938193057
1.40131927153131
–0.871276408882294
0.183785623507722
0.0359103726343979
1.01536263500695
0.896348850692318
0.954285278486704
–0.890298849338373
1.40131909032315
–0.871276197318301
0.183785492186491
0.0359104061952515
1.01536263525763
0.896348847022403
0.954285302637587
–0.890298935334171
1.40131926588117
–0.871276402487896
0.18378561965322
0.0359103735933848
1.01536263525905
0.896348846999736
0.954285302790922
–0.890298935876618
1.40131926696099
–0.87127640370648
0.183785620380814
0.0359103734148875
Error F7.89749200363512×10117.89749200363504×10117.89749200353888×10117.8974920035389×1011

Table 3.

Performance of the algorithm when it is run to solve the model using four different parameters guesses.

Figures 9, 10, 11, 12 show the fitness of the estimated parameters and the rate of convergence.

Figure 9.

All the initial values for the parameters are set 0.1.

Figure 10.

All the initial values for the parameters are set 0.5.

Figure 11.

All the initial values for the parameters are set 1.

Figure 12.

All the initial values for the parameters are set 2.

Figures 13, 14, 15, 16 show the comparison between the inversion results q^xand the exact value qx:

Figure 13.

The comparison chart with all the initial values for the parameters is set 0.1.

Figure 14.

The comparison chart with all the initial values for the parameters is set 0.5.

Figure 15.

The comparison chart with all the initial values for the parameters is set 1.

Figure 16.

The comparison chart with all the initial values for the parameters is set 2.

Table 4 shows the values of qjΔxand ujΔx0.5in x=jΔxwith the all the initial values are set 1.

NumericalExactNumericalExact
jqjΔxqjΔxujΔx0.5ujΔx0.5
01.0153626352576310.3033429626440880.303265329856317
11.113781680590131.105170918075650.3292019646771260.329286981656416
21.227656949593991.221402758160170.3474928822248860.347609712653987
31.355490213058741.3498588075760.3593475687026780.359463171293777
41.49737221492651.491824697641270.3658087113211590.365912693766539
51.654328284152061.648721270700130.3677923782088570.367879441171442
61.827850568693931.822118800390510.3660912184541820.366158192067887
72.019634990464642.013752707470480.3613877614737960.361433054294643
82.231541020068792.225540928492470.3542678692730630.354291330944216
92.465792370156872.459603111156950.3452330236180590.345235749518249
102.725436706223332.718281828459050.3347126048031750.334695240222645

Table 4.

The values of qjΔxand ujΔx0.5in x=jΔxwith the all the initial values are set 1.

6. Conclusions

A numerical method to estimate the temperature uxtand the coefficient qxis proposed for an IHCP and the following results are obtained.

  1. The present study, successfully applies the numerical method involving the Levenberg-Marquardt algorithm in conjunction with the Galerkin finite element method to an IHCP.

  2. From the illustrated example it can be seen that the proposed numerical method is efficient and accurate to estimate the temperature uxtand the coefficient qx.

Acknowledgments

The work of the author is supported by the Special Funds of the National Natural Science Foundation of China (Nos. 51190093 and 51179151). The author would like to thank the referees for constructive suggestions and comments.

Conflict of interests

The authors declare that there is on conflict of interests regarding the publication of this article.

Download

chapter PDF

© 2019 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution 3.0 License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

How to cite and reference

Link to this chapter Copy to clipboard

Cite this chapter Copy to clipboard

Tao Min, Xing Chen, Yao Sun and Qiang Huang (October 7th 2019). A Numerical Approach to Solving an Inverse Heat Conduction Problem Using the Levenberg-Marquardt Algorithm [Online First], IntechOpen, DOI: 10.5772/intechopen.89096. Available from:

chapter statistics

42total chapter downloads

More statistics for editors and authors

Login to your personal dashboard for more detailed statistics on your publications.

Access personal reporting

We are IntechOpen, the world's leading publisher of Open Access books. Built by scientists, for scientists. Our readership spans scientists, professors, researchers, librarians, and students, as well as business professionals. We share our knowledge and peer-reveiwed research papers with libraries, scientific and engineering societies, and also work with corporate R&D departments and government entities.

More About Us