Numerical Solution of a Nonlinear Inverse Heat Conduction Problem

(1)

Dissertation

Link¨oping Studies in Science and Technology

Numerical Solution of a Nonlinear Inverse

Heat Conduction Problem

Muhammad Anwar Hussain

LiTH - MAT - EX - 2010 / 10 - SE

Department of Mathematics

Linköping University, SE-581 83 Linköping, Sweden Linköping 2010

(2)

(3)

Numerical Solution of a Nonlinear Inverse

Heat Conduction Problem

Scientific Computing, Link¨opings Universitet

Muhammad Anwar Hussain

LiTH - MAT - EX - 2010 / 10 - SE

581 83 LINK ¨OPING

SWEDEN

Examensarbete: 30 hp

Level: D

Supervisor: Prof. Lars Eld´en,

Examiner: Prof. Lars Eld´en,

(4)

(5)

(6)

(7)

Abstract

The inverse heat conduction problem also frequently referred as the side-ways heat equation, in short SHE, is considered as a mathematical model for a real application, where it is desirable for someone to determine the temperature on the surface of a body. Since the surface itself is inaccessible for measurements, one is restricted to use temperature data from the inte-rior measurements. From a mathematical point of view, the entire situation leads to a non-characteristic Cauchy problem, where by using recorded tem-perature one can solve a well-posed nonlinear problem in the finite region for computing heat flux, and consequently obtain the Cauchy data [u, ux].

Fur-ther by using these data and by performing an appropriate method, e.g. a space marching method, one can eventually achieve the desired temperature at x = 0.

The problem is severely ill-posed in the sense that the solution does not depend continuously on the data. The problem solved by two different methods, and for both cases we stabilize the computations by replacing the derivative _∂t∂ in the heat equation by a bounded operator. The first one, a spectral method based on finite Fourier space is illustrated to supply an analytical approach for approximating the time derivative. In order to get a better accuracy in the numerical computation, we use cubic spline function for approximating _∂t∂ in the least squares sense.

The inverse problem we want to solve, by using Cauchy data, is a nonlin-ear heat conduction problem in one space dimension. Since the temperature data u = g(t) is recorded, e.g. by a thermocouple, it usually contains some perturbation in the data. Thus the solution can be severely ill-posed if the Cauchy data become very noisy. Two experiments are presented to test the proposed approach.

Keywords: inverse problem, ill-posed, Cauchy problem, heat conduction,

(8)

(9)

Acknowledgments

The entire work of this thesis was improved by the essential and gracious support of many individuals.

At first, I would like to express my deep respect and sincere gratitude to my supervisor Prof. Lars Eld´en for his invaluable guidance and generous attitude. His extensive programing knowledge and weighty discussion always influenced me to develop my creativity in the area of scientific computing. The reflection of my computing ability all over this work and the deliberation into the complete masters program would not have been possible without his willingly support and personal kindly interest.

I would like to thank Dr. Fredrik Berntsson whose keen computing skill affected me, and also some handy discussions provided me a good knowl-edge during the thesis work. I am really thankful to very generous PhD student Zohreh Ranjbar for her supporting mind, and for supplying me all the necessary ideas and techniques for the use of LaTeX.

I would like to thank my beloved teachers Prof. Irina Yakymenko, and Prof. Sven Stafstr¨om for their huge supporting role throughout my masters degree in Link¨oping University.

Finally, I would like to thank my honorable teacher Prof. Abdullahel Baqui who always helped me in my personal issues, and influenced me for developing a good career. I am very grateful to my family and friends for their unconditional support during the work.

Muhammad Anwar Hussain Link¨oping, Sweden

(10)

(11)

Introduction

For solving heat conduction problems numerically, it is important to iden-tify the problems as either well-posed or ill-posed. In the sense of Jacques Hadamard [18], for all admissible data, a well-posed problem should have the following properties:

A solution exists. The solution is unique.

The solution depends continuously on the data.

Mathematical models which violate these qualities are termed as ill-posed problems. In real applications, it is significant to specify which type of data that are admissible, which type of norm that should be applied for measuring continuity and exactly what it is meant by a particular solution [4].

Inverse Heat Conduction Problems, are often ill-posed models in science

and engineering. Such problems have been discussed in the last couple of decades by many authors [2, 5, 6, 7, 10]. It is sometimes necessary to determine, the surface temperature of a body from a measured temperature history at a given fixed location inside the body, when the surface itself is inaccessible for measurements. To clarify, inaccessible means locating a measurement device, e.g. a thermocouple, on the surface would disturb the measurements and the recorded temperature might be incorrect, or it could be difficult to set up a device because of high temperature. For all of these facts, it is always appropriate to compute the temperature from the internal measurements. A number of solution methods have been proposed in [4, 7, 10, 11, 12, 13] to solve these particular problems. Our main area of interest throughout this thesis paper is to solve nonlinear inverse heat conduction problems in one space dimension. We consider the nonlinear inverse heat conduction problem, where the temperature u(x, t) is determined for the region 0 < x < L from temperature measurement g(t) = u(L, t) and the

(14)

Hot region (Hr) gas or liquid Solid region (Sr) steel Hr - Sr interface u(0,t) = ? ux(0,t) = ? measured temperature and heat flux u(L,t) = g(t) and ux(L,t) = h(t) respectively measurement device e.g. thermocouple x=0 x=L x

Figure 1.1: Schematic diagram presenting a description of the inverse heat conduction problem for both linear and nonlinear case.

heat flux h(t) = ux(L, t), then u(x, t) satisfies

κ(u)ux_x= ρ(u)cp(u) ut, 0 < x < L, 0 < t < T,

ut(x, 0) = 0, 0 < x < L,

u(L, t) = g(t), 0 < t < T, ux(L, t) = h(t), 0 < t < T,

(1.1)

where the density ρ(u), the specific heat cp(u) and the thermal conductivity

κ(u) of the material vary with temperature, so that the problem (1.1) is nonlinear. In practice, since g is measured and h is computed from this data, there will be measurement errors, and we have only approximate data gδ

m and hδm satisfying,

||g − gδm|| < ǫg, and, ||h − hδm|| < ǫh,

where the positive constants ǫg and ǫh represents bounds on the

measure-ment errors.

The problem (1.1) can be considered as a Cauchy problem with appro-priate Cauchy data [u, ux], given on the line x = 1. Thus this problem is

often referred to as the sideways heat equation, so that we can distinguish it from other inverse heat conduction problems. Indeed, in most of the cases, Cauchy problems are well-known to be severely ill-posed [4, 18]. The prop-erties and investigation of ill-posedness is discussed further in section 2, and

(15)

3

a theoretical analysis of the problem, for the case of constant coefficients, is performed in section 5.

The problem, we are interest to solve in the differential equation form, (κ(u)ux)x = ρ(u)cp(u) ut, is a nonlinear problem. But the numerical

meth-ods, we mostly discuss in this paper, are based on the form κuxx = ut, in

order to get rid of the complexity of the nonlinearity in the theoretical de-scription. Indeed, the numerical implementation of these methods provide all necessary information that one can solve a typical linear problem, i.e. κuxx = ut and a problem with non-constant coefficients (κ(x, t)ux)x = ut,

as well as the nonlinear problem (1.1) dependent on temperature. Note that, for a more general case, e.g. differential equation (1.1), we can not use an integral equation1 of the form,

Z 1

0 k(t − τ)Ψ(τ)dτ = g δ m(t),

since the kernel k(t) is not explicitly known. Instead our purpose is to preserve the problem in the differential equation form and to solve it as an initial value problem by using a space marching method, i.e. a method of lines, is discussed in section 5. For this case we discretize the space variable with an equidistant grid point, which performed as hx in the numerical

experiment.

It was proved in [11, 12] that it is possible to implement space-marching methods efficiently, if the time derivative is approximated by a bounded operator. This is discussed in section 6.2. Several different methods for approximating the time derivative has been considered, and two of such methods are studied here. These are, a spectral approximation using finite Fourier space in the interval 0 ≤ x ≤ L and the approximation of time derivative of discrete functions by using cubic spline functions.

In order to solve the Cauchy problem for the heat equation (1.1) we essentially use industrial data T C1, T C2 and T C3, which are the functions of t, were recorded by a measurement device from a real application. For this case it is needed to arrange a particular settings, such that, we can first solve a test problem, e.g. model-1, to ensure that our approach for solving a sideways heat equation works well. As we accomplish the numerical procedures and experiments for model-1 successfully, we have all necessary information to determine the surface temperature at x = 0 by solving our main problem, e.g. model-2, see chapter 7.1. The problem is more difficult for model-2 than that of model-1, since we have no information about the temperature at x = 0. But, in the case of model-1, the setting is like that we only need to reconstruct the solution temperature at x = 110 mm. The entire situation can be seen graphically by the figure 7.1 in section 7.1.

(16)

The ill-posedness of the inverse problem can easily be seen by using the industrial data, since these data vectors are not so smooth, usually contain a large number of perturbation in the data, and also the thermal conductivity of the steel material vary with time, see chapter 7.

(17)

Chapter 2

Investigation of Ill-posedness

For investigating the ill-posedness we consider the linear form of the inverse problem (1.1), where it is assumed that the body is very large and the setting is one dimensional. To study the phenomenon, we consider the following heat equation problem in the quarter plane: Determine the temperature u(x, t) for 0 ≤ x ≤ 1 from temperature measurements g(t) = u(1, t), when u(x, t) satisfies

κuxx= ρcp ut, x > 0, t ≥ 0,

u(x, 0) = 0, _{x ≥ 0,}

u(1, t) = g(t), _{t ≥ 0, u|}_x→∞ bounded.

(2.1)

Since g(t) is measured, it is obvious that there will be some measurement errors. We now distinguish between the exact data g(t) and the measurement data gm(t) which is merely in the L2(R) norm. In order to illustrate the

distinction we have the following discussion: Let us assume that D is a closed interval, i.e., we define L2(D) for a ≤ D ≤ b to be the set of all functions f such that

Z

D |f(t)| 2

dt < ∞,

and we define the inner product and norm on L2(D) by hf1, f2i = Z D f1(t)f2(t)dt , kfk2 = hf, fi = Z D |f(t)| 2_{dt ,}

then there is a sequence {fn} that converses to f in norm, such that each

fn is continuous on D and vanishes outside some bounded set [14, chapter

3].

Thus we would actually have as data some function gm∈ L2(R) practically

contains some measurement errors, i.e., gm≡ gmδ, such that

kgm− gk = kgmδ(·) − u(1, ·)k ≤ ǫ, (2.2)

(18)

where ǫ > 0 is a constant. It is identified as a upper bound for the L2 norm on the measurement errors.

Indeed, one can use the measurement data gm to solve a well-posed

problem when x ≥ 1, also ux(1, ·) is determined. But for 0 ≤ x ≤ 1, we

would have the sideways heat equation2

κuxx = ut, x > 0, t ≥ 0,

u(x, 0) = 0, _{x ≥ 0,}

u(1, t) = gm(t), t ≥ 0, u|x→∞ bounded,

(2.3)

where ρ and cpare taken as unit value to make the investigation more simple.

The sideways heat equation (2.3) is ill-posed in the sense that the so-lution, if it exists, does not depend continuously on the data. We will investigate the situation using different approaches.

2.1 Analyze the Ill-posedness with Sine Function

We investigate the case where the data function g(t) is taken as a sine function with a fixed frequency ω such that for x = 1, we have

g(t) = bg(ω) sin(ωt) ∈ L∞(R). (2.4) The easiest way to solve the problem is by using the Ansatz

u(x, t) = C(x)eiωt+ D(x)e−iωt

in the sideways heat equation. Then we find the solution can be written

u(x, t) = bg(ω)e√2κω(1−x)sin ωt +

r ω

2κ(1 − x)

.

At x = 0, we obtain the surface temperature

u(0, t) = bg(ω)e√2κω sin(ωt +

r ω

2κ). (2.5)

We assume bg(ω) = 1/ω2, which provides, for ω → ∞ both g(t) and g′(t) tends to zero. But the amplitude of the solution u(0, t) tends to infinity as ω → ∞. Consequently we get a clear idea from the equation (2.5) is that the amplitude of oscillations are magnified by a factor e√2κω. This means that

the solution does not continuously depend on the data. The direct problem, that of computing u(1, ·) from u(0, 1), has a strong smoothing effect on the data. Thus, to obtain a solution that depends continuously on the data, the high frequency parts have to be removed.

The above example is analogous to that used by Hadamard [18] when he demonstrated the ill-posedness of the Cauchy problem for the Laplace equation.

(19)

2.2. Ill-posedness Analysis in the Fourier Space 7

2.2 Ill-posedness Analysis in the Fourier Space

In order to investigate the ill-posedness in the Fourier domain, we are only interested to find a solution in the finite time interval 0 ≤ t ≤ T . Since we have no idea about t ≥ T , we assume that the surface temperature drops dramatically to zero for all t ≥ T . Moreover, to simplify the analysis we define all functions to be zero for t < 0. Let bg be the Fourier transform of exact data function g,

b g(ξ) = _√1 2π Z _∞ −∞ g(t)e−iξtdt, _{−∞ < ξ < ∞.} (2.6) Now, suppose that the Fourier transform of F onR = (−∞, ∞) is

F[F (η)] = bF (ζ),

then differentiation takes bF (ζ) to iζ bF (ζ). Therefore, by taking Fourier trans-form of (2.1) we get the corresponding equations in the frequency domain [14, 28, chapter 7, 4]:

κbuxx(x, ξ) = iξbu(x, ξ), x > 0 and ξ ∈ R,

b

u(1, ξ) = bg(ξ), ξ ∈ R, b

u(·, ξ)|x→∞ bounded.

(2.7)

This is a family of ordinary differential equations parameterized by ξ. By using the Ansatz u(x) = cemx and plugging into (2.7) the general solution can be obtained in the form

b u(x, ξ) = C1(ξ)ex q iξ κ _{+ C} 2(ξ)e−x q iξ κ_, _(2.8)

where√iξ denotes the principal value of the square root, p iξ =    (1 + i) q |ξ| 2 , ξ ≥ 0, (1 − i) q |ξ| 2 , ξ < 0.

Since the real part of√iξ is positive, we can see that the exponential func-tion, the first term ex

q

iξ

κ _{in right side of (2.8), grows rapidly as x tends to}

infinity. Therefore, to find a bounded solution at infinity, we conclude that C1(ξ) = 0. The equation (2.8) now be simplified by using the boundary

condition x = 1, b u(1, ξ) = C2(ξ)e− q iξ κ, yields b g(ξ) = C2(ξ)e− q iξ κ_,

(20)

which leads to

C2(ξ) = e q

iξ κg(ξ)._b

Hence the equation (2.8) provides the solution

b

u(x, ξ) = e

q

iξ

κ(1−x)_bg(ξ). (2.9)

In order to obtain the solution (2.9) we have bound on the solution at infinity. Since the real part of√_{iξ is positive, and the solution b}u(x, ξ) is assumed to be in L2(R), we see that the exact data function, bg(ξ), must decays rapidly as ξ → ∞.

Since we measured the temperature by measurement device, there would have some measurement errors. Hence, we assume that

b

g_mδ(t) = g(t) + δ(t), (2.10)

where δ ∈ L2₍_{R) is a small measurement error. The Fourier domain if we}

use gm as data function, we get the following solution

b w(x, ξ) = ex q iξ κ(1−x)_b_{g(ξ) + b}_δ(ξ) _, This leads to b w(x, ξ) = bu(x, ξ) + ex q iξ κ(1−x)bδ(ξ). (2.11)

In the solution bw(x, ξ), since the measurement error bδ(ξ) do not have the same decay in frequency as much as the actual data bg(ξ), the solution bw(x, ξ) will not, in general, be in L2(R). Thus the ill-posedness can easily be ob-served in the blow-up of high frequency perturbations in the data when the equation (2.3) is solved numerically. Since the factor e√|ξ|/2κ magnify the frequency components of (2.9), thus it might be said that the degree of ill-posedness depends on the coefficient κ.

2.3 Standardization of the Region

In solving mathematical modelling problems using partial differential equa-tions it is almost always wise to rewrite the equaequa-tions into the dimensionless forms. By taking the proper choice of dependent and independent variables one can modify the original problem to the dimensionless form. Here we will scale both the time and space variables. We consider the heat equation in the linear form

ρcp

∂u ∂t = κ

∂2u

∂x2, x > 0, t ≥ 0, (2.12)

where cp is the specific heat, ρ is the density and κ is the thermal

(21)

2.3. Standardization of the Region 9

is the temperature inside the body and g(t) is the recorded temperature measured at a distance x = L from the surface of the body. In practice, g(t) can only be measured for a finite time interval 0 ≤ t ≤ T . We introduce dimensionless space coordinate ξ and time co-ordinate τ as follows:

ξ = x

L, and τ = t T.

Now, through the scaling, derivatives of original model can be transformed to derivatives of the dimensionless variables. By applying Euler’s chain rule we then have the required form

∂u ∂τ = K ∂2u ∂ξ2, (2.13) where K = _ρcκT pL2 is a constant.

For the measurement data, the function g is considered to be known on the interval 0 ≤ t ≤ 1. It can be observed that, when the value K decreases the degree of ill-posedness increases. Therefore for a smaller value of L we would have better stability properties of the problem. Indeed, the value of K also depends on T . Thus in measuring the data on a finite time interval one might deduce, that the length of interval can influence the stability properties of the problem. In fact, this is actually not the case. we can introduce a standard property of the Fourier transform

\

g(T t)(ξ) = 1 Tbg(

ξ T).

Based on this property we can make a conclusion, that any attempt to increase the stability by using a larger T will cause the scaled data function to be spread out in the frequency domain, consequently it would be difficult for someone to solve the problem.

It is noticeable, although we seek to recover u only for 0 ≤ x < 1, the problem specification includes the heat equation for x > 1, together with the boundness of the solution u(x, t) at infinity. Since we can obtain u for x > 1 by solving a well-posed quarter plane problem, also ux(1, ·) is determined.

Thus the problem (2.3) usually have the form on the line x = 1,

κuxx= ρcp ut, 0 < x < 1, 0 < t < 1,

u(1, t) = g(t), ux(1, t) = h(t),

u(x, 0) = 0,

(2.14)

(22)

where w(x, t) be the solution of well-posed problem κwxx = ρcp wt, 0 < x < 1, 0 < t < 1, w(0, t) = 0, wx(1, t) = h(t), w(x, 0) = 0, (2.15)

and v(x, t) be the solution of ill-posed problem

κvxx = ρcp vt, 0 < x < 1, 0 < t < 1,

v(1, t) = g(t) − w(1, t), vx(1, t) = 0,

v(x, 0) = 0.

(2.16)

Hence, from the above presentation it might be concluded that we can solve the inverse heat conduction problem (SHE) from a quarter plane problem by considering an appropriate length of the region, where we would have the measurement data g(t) and the heat flux h(t) by solving the well-posed problem, in the x-direction.

(23)

Chapter 3

Stabilizing Assumptions

The idea of stabilization is an important concept for any type of ill-posed problem. A stable approximation means small perturbations in the mea-surement data cause only small perturbations in the solutions. It means that the solution always depends continuously on the data. In the whole section we will discuss how to reconstruct the problem by imposing a priori bound on the solution data, and also how to get a solution that must con-tinuously exist on the data. We already found the solution operator defined by (2.9) is unbounded, that is, the problem is ill-posed. In order to make the solution stable we have to impose a priori bound M on the solution at x = 0, such that, ku(0, ·)k ≤ M. Thus, by considering some imprecisions in the matching of the data we can approximate the problem (2.3) by,

κuxx = ut, x > 0, t ≥ 0,

u(x, 0) = 0, _{x ≥ 0,} ku(1, ·) − gδm(·)k ≤ ǫ,

ku(0, ·)k ≤ M.

(3.1)

We will show that any two solutions of the equation (3.1), say u1 and u2,

from the given data g(t) at x = 1 with kui(1, ·) − gδm(·)k ≤ ǫ, and prescribed

bound ku(0, ·)k ≤ M at x = 0, satisfy the following inequality

ku1(x, ·) − u2(x, ·)k ≤ 2M1−xǫx, (3.2)

where ǫ, M > 0 are both known and ǫ << M.

The convexity inequality is considered for stabilizing the ill-posed con-tinuation problem when the data function g contains some perturbations. The inequality defined by (3.2) was first proved by [21]. Also a useful pre-sentation on convexity inequality can be found in [8]. Our prepre-sentation is based on that in [4]. Indeed, the degree of ill-posedness is stabilized by the use of such bound together with the analysis of resulting continuity with respect to the data. The estimate (3.2) is an approximation of Logarithmic

(24)

Convexity type. In the rest of this section, we will first define the concept of logarithmical convexity, and the proof of inequality (3.2) will be given later on.

3.1 The Proof

Definition 3.1.1. A function φ defined on a convex subset of a real vector space and taking positive values is said to be logarithmically convex if log φ(x) is convex function of x.

If φ > 0 and φ is log convex on an interval [Ω1, Ω2], then the inequality

log φ(x) ≤ (1 − γ(x)) log φ(Ω1) + γ(x) log φ(Ω2), x ∈ [Ω1, Ω2],

holds for γ(x) = Ω2−x

Ω2−Ω1 .

The inequality can be rewritten as

φ(x) ≤ φ(Ω1)(1−γ(x))φ(Ω2)γ(x), x ∈ [Ω1, Ω2]. (3.3)

Moreover, a non-negative C2 _{function φ is log convex if and only if}

φφ′′_{− (φ}′)2 _{≥ 0. This is a consequence of the fact that a function φ is convex} if its second order derivative is positive

(log(φ))′′ = φφ′′− (φ′)

2

φ2 ≥ 0, ⇐⇒ φφ′′− (φ′) 2_{≥ 0.}

In this paper, we use the log convex function for error estimates. Thus it is necessary to allow such functions to have zeros.

Lemma 3.1.2. Let φ(x) be continuous, nonnegative and log convex on the interval Ω. Then φ(x) > 0 on Ω or φ(x) ≡ 0 on Ω.

Proof. Suppose that φ(x) > 0, for some x ∈ Ω. We shall prove that φ > 0

in the whole interval. If φ has zeros in Ω, then due to the continuity there exists an s1 ∈ Ω such that

φ(s1) = 0, and φ(s2) > 0, s2 ∈ (s1, x].

Let ǫ > 0 be given and choose a ξ ∈ (s1, x) such that φ(ξ) < ǫ.

Since φ is a log convex and φ > 0 on [ξ, s1] we find that for s2 ∈ [ξ, s1], we

have

φ(s2) ≤ φ(ξ)ψ(s2)φ(x)1−ψ(s2)

≤ ǫψ(s2)_φ(x)1−ψ(s2)_,

where φ(s2) = x−s_x−ξ2 6= 0.

Since ǫ can be arbitrarily small we conclude that φ(s2) = 0, s2 ∈ (s1, x)

and since φ is continuous also φ(x) = 0, this contradicts the assumption that φ(x) > 0, hence φ has no zeros in Ω. This means that either φ = 0 or φ > 0 .

(25)

3.1. The Proof 13

Lemma 3.1.3. Let φ1 and φ2 are continuous, non-negative and log convex

on an interval Ω. Then φ1+ φ2 also has the same properties on the interval

Ω.

Proof. We shall prove the sum by applying H¨older’s inequality. We have to

show that the inequality (3.3) holds for every subinterval [Ω1, Ω2] ⊂ Ω. For

simplicity we may assume that [Ω1, Ω2] := [0, 1]. Let x ∈ (0, 1), then

φ1(x) + φ2(x) ≤ φ1(0)xφ1(1)1−x+ φ2(0)xφ2(1)1−x ≤nφ1(0)xp+ φ2(0)xp o1 pn φ1(1)(1−x)q+ φ2(1)(1−x)q o1 q ,

where 1_p +1_q = 1. If we substitute p = 1_x, we then get

φ1(x) + φ2(x) ≤ n φ1(0) + φ2(0) oxn φ1(1) + φ2(1) o_1−x .

Hence φ1+ φ2 is log convex.

We can now prove the inequality (3.2), and the proof presented here is a simplified version of a more general result by Levine [21, sec.3].

Theorem 3.1.4. Let u1 and u2 be two solutions of the problem (3.1). Then

ku1(x, ·) − u2(x, ·)k ≤ 2M1−xǫx, 0 ≤ x ≤ 1. (3.4)

Proof. For simplicity we use κ = 1, for equation (3.1), in the proof. Let

w = u1− u2, and φ(x) = kw(x, ·)k2. So we can define bw(x, ξ) at x = 1, using

the idea of equation (2.9), as

b

w(x, ξ) = e√iξ(1−x)w(1, ξ).b By applying Parseval relation we determine that

φ(x) = Z _∞ −∞| bw(x, ξ)| 2_{dξ =}Z ∞ −∞ e√2|ξ|(1−x)_{| b}_{w(1, ξ)|}2dξ. Using the concepts, described in definition (3.1.1), we have

φ(x)φ′′_{(x) − (φ}′(x))2 = Z _∞ −∞| bw(x, ξ)| 2 dξ Z ∞ −∞(2 |ξ|)| bw(x, ξ)| 2 dξ − Z _∞ −∞ p (2 |ξ|)| bw(x, ξ)|2dξ2 _{≥ 0,} where the last inequality follows from the Cauchy-Schwarz inequality. This means that kw(x, ·)k2 is log convex. Thus we obtain

ku1(x, ·) − u2(x, ·)k2≤ ku1(0, ·) − u2(0, ·)k2

_1−x

ku1(1, ·) − u2(1, ·)k2

x

(26)

If we impose the priori bound kui(0, ·)k ≤ M for i = {1, 2} on the solution,

and apply the condition kui(x, ·)−gδm(·)k ≤ ǫ, then by the concept of triangle

inequality we finally obtain

ku1(x, ·) − u2(x, ·)k ≤ 2M1−xǫx,

which completes the proof of (3.4).

The above theorem proves that we have a stable dependence on the data for the problem (3.1) in the sense (3.2), and so we call (3.1) a stabilized problem. The inequality (3.4) is sharp (cf. 4.1), thus we can not expect to find a numerical method which satisfy for a better error estimate for approx-imating solutions of (3.1). But one difficulty of (3.1) is that it does not have uniqueness. Therefore, in order to find a useful numerical procedure it is convenient to modify the problem by choosing some parameters dependence on both ǫ and M.

(27)

Chapter 4

Regularizations

In this section we will study how to stabilize the sideways heat equation by cutting off high frequencies from the solution. It means that the ill-posed problem will be regularized. A variety of numerical and analytical tech-niques have been developed, e.g. [7, 17], to find stable solutions for the inverse heat conduction problems. We usually in common have some meth-ods like Tikhonov regularization method, mollification method, Alifanov’s iterative regularization method and the method of Beck. Here we consider a regularization method in the Fourier space by cutting off high frequencies from the solution. Similar approach of the method can be found in [4, 25]. We introduce the problem in Fourier space parameterized by ξ,

κbuxx(x, ξ) = iξbu(x, ξ), x > 0 and ξ ∈ R,

b

u(1, ξ) = bg(ξ), ξ ∈ R,

b

u(·, ξ)|x→∞ bounded.

(4.1)

It has a solution in the Fourier domain, already illustrated in section 2, of the form u(x, t) = √1 2π Z _∞ −∞ e q iξ κ(1−x)_b_g(ξ)eiξt_dξ. _(4.2)

Since the principal value of √iξ in the term e

q

iξ

κ(1−x) has a positive real

part, small error in high frequency components can blow-up and completely destroy the solution for 0 ≤ x < 1. A natural way to stabilize the problem is to remove all high frequencies from the solution and instead consider (4.1) only for |ξ| < ξmax. Then we have the regularized solution defined by the

relation

b

w(x, ξ) = Xmaxu(x, ξ),b

which yields to the form

w(x, t) = _√1 2π Z _∞ −∞ e q iξ κ(1−x)_b_g(ξ)X max eiξtdξ, (4.3) Hussain, 2010. 15

(28)

where Xmax is the characteristic function of the interval [−ξmax, ξmax].

In real applications, instead of having exact data g we would actually have an access only for some approximation gm of g, where the perturbation

can be defined

δ = g − gmδ,

and let wδ_{(x, t) be the regularized solution defined by (4.3) with the}

mea-sured data g_mδ. The difference between exact solution u(x, t) and regularized solution can be divided into two parts,

ku(x, ·) − wδ(x, ·)k ≤ ku(x, ·) − w(x, ·)k + kw(x, ·) − wδ(x, ·)k = RT + RX,

(4.4)

where RT and RX are defined

In the regularization process these two terms behave differently. The trunca-tion error RT, is a measure of the effect of the introduced cut off frequency.

By increasing ξmax we include a larger part of the spectrum in the solution,

and thus we reduce the truncation error. We have also propagated data error RX, is generated by the inexact data. If we increase ξmax, then RX

will eventually dominate the solution. Thus the total error depends on the eliminating part of the frequency, ξmax.

In fact, we would have some optimal cut off frequency, but we can not determine it explicitly since it depends on unavailable information about the exact data g, e.g. the smoothness and on the errors in the data. However, it is possible to give a parameter choice rule which, given some a priori information about the solution, leads to a convergent scheme [17, section 3], that is,

ku(x, ·) − wδ(x, ·)k → 0, if kδk → 0.

For solving the ill-posed problems, such a scheme is referred to as a regu-larization method in [17]. It is necessary to investigate the quality of our method using the perturbed data gδ_m. In the following sections we will dis-cuss about ‘parameter choice rule’ by deriving an error estimate for the approximate solution (4.3).

4.1 Error Estimate

In general, it is impossible to determine the solution of an ill-posed problem without knowing any information of the error level [1, 30], to get the error

(29)

4.1. Error Estimate 17

estimate for an approximate solution [1, 29] and even to get a rate of conver-gence for approximate solutions on the whole infinite dimensional Banach space [1]. So it would be convenient to have all available a priori information about the exact solution for constructing a set of rules with special prop-erties. In this section we illustrate such a bound on the difference between unbounded operator u(x, ·) and bounded operator w(x, ·) of (4.2) and (4.3) respectively. We assume that we have an a priori bound ku(0, ·)k ≤ M on the exact solution of the problem (4.1), where δ ≥ 0 be an error level of u(1, ·), is supposed to be known.

Now we will investigate some error bound conditions using the regular-ized solutions, say, w1 and w2 of (4.3) by the following lemmas.

Lemma 4.1.1. Suppose that we have two regularized solutions w1 and w2

defined by (4.3) with data g1 and g2, satisfying kg1− g2k ≤ ǫ. If we select

ξmax= 2κ(log(M/ǫ))2, then we get the bound

kw1(x, ·) − w2(x, ·)k ≤ M1−xǫx. (4.5)

Proof. In order to obtain L2estimates for the difference {w1(x, ·)−w2(x, ·)},

we may introduce Parseval relation. Then we have

kw1(x, ·) − w2(x, ·)k2 = Z _∞ −∞|w 1(x, t) − w2(x, t)|2dt = Z _∞ −∞| b w1(x, ξ) − bw2(x, ξ)|2dξ = k bw1− bw2k2 = Z ξmax −ξmax e q iξ κ(1−x)_{( b}_g 1− bg2) 2 dξ ≤ e q iξmax κ (1−x) 2 k bg1− bg2k2 = e q 2_ξmax κ (1−x)_kg 1− g2k2.

Using ξmax= 2κ(log(M/ǫ))2, we finally obtain

kw1− w2k ≤ M1−xǫx.

From Lemma 4.1.1 we observe that the regularized solution (4.3) contin-uously depends on the input data. In the next lemma we will examine the difference between (4.2) and (4.3) using the exact data g(t).

Lemma 4.1.2. Suppose, u and w be the solutions defined by (4.2) and (4.3)

with the same exact data g, and let ξmax = 2κ(log(M/ǫ))2. If ku(0, ·)k ≤ M

then

(30)

Proof. In order to obtain L2 estimates for the difference {u(x, ·) − w(x, ·)}, we may introduce Parseval relation and consider the fact, that the solutions coincide for ξ ∈ [−ξmax, ξmax]. we have

ku(x, ·) − w(x, ·)k2 = Z _∞ −∞|u(x, t) − w(x, t)| 2_dt = Z _∞ −∞|bu(x, ξ) − bw(x, ξ)| 2 dξ = Z |ξ|>ξmax e q iξ κ(1−x)_bg(ξ) 2 dξ ≤ e−x q 2_ξmax κ Z |ξ|>ξmax e q iξ κ_b_g(ξ) 2 dξ ≤ e−x q 2_ξmax κ Z |ξ|>ξmax |bu(0, ξ)|2dξ ≤ e−x q 2ξmax κ kbu(0, ·)k2.

Finally, we use the a priori error bound ku(0, ·)k ≤ M together with the frequency ξmax = 2κ(log(M/ǫ))2 , the unknown solution leads to the error

bound

ku − wk ≤ M1−xǫx.

Now we have all information to formulate the major part of this section. Throughout the following theorem, we will observe that the regularized so-lution wδ(x, t) is actually an approximation to the actual solution u(x, t). From the relation (3.2) we can also conclude that the approximation er-ror depends continuously on the measurement erer-ror. Similar survey can be viewed in [4].

Theorem 4.1.3. Let u be the solution (4.2) with exact data g and wδ the

solution (4.3) with measured data gδ

m. If we have a bound ku(0, ·)k ≤ M

and the measured function g_mδ _{satisfies kg − g}_mδ_{k ≤ ǫ and also if we choose}

ξmax= 2κ(log(M/ǫ))2, then we get the error bound

ku − wδk ≤ 2M1−xǫx. (4.7)

Proof. Let wδ be the solution described by (4.3) with exact data g. Using

the inequality

|p + q|2 ≤ |p|2+ |q|2+ 2|p||q| ≤ 2(|p|2+ |q|2), together with the previous two Lemmas, we obtain

(31)

4.1. Error Estimate 19

It is noticeable, that, the sharpness of the inequality (4.6) depends only on the spectral contents of the exact solution. Now

ku(x, ·)k2 ≈ e−x q 2ξmax κ Z |ξ|≥ξmax |bu(0, ξ)|2dξ,

if the support of bu(0, ξ) is localized close to the cut-off frequency, conse-quently the estimate, in Lemma 4.1.2, can be very tight. On the other hand, Z ξmax −ξmax e q iξ κ(1−x)_{(bg − bg}δ m) 2 dξ ≈ e q 2_ξmax κ (1−x)_{kbg − bg}δ mk,

if the support of the measurement error (bg − bgδ

m) is localized close to ξmax.

Hence, we make a statement that, the Lemma 4.1.2 necessarily gives a sharp inequality, and for any case, we can not expect a better error estimate than that in 4.1.3. We can also find a handy discussion on sharpness of the estimate (4.7) in [21, 7].

Next, we will discuss some error bounds on the heat flux ux for the

domain 0 < x < 1. From (4.2), we also obtain a formula for the gradient ux(x, t) by differentiating with respect to x inside the integral. We may

have, ux(x, t) = − 1 √ 2π Z _∞ −∞ r iξ κ e q iξ κ(1−x)_g(ξ)e_b iξt_dξ. _(4.8)

From (4.2) and (4.8) we find a relation between h and g,

bh(ξ) = − r

iξ κ bg(ξ)

We obtain error bounds in the L2-norm for temperature gradient h := ux(1, ·) to the approximated measured data hm at x = 1.

Corollary 4.1.4. Let 0 < x < 1. Let u and wδ are defined as in Theorem

4.1.3. Suppose that we have a bound kux(0, t)k ≤ Md, and that the measured

temperature gradient hm satisfies kh − hδmk < ǫd, and that we choose

ξmax= 2κ(log(Md/ǫd))2, then we get the error bound

kux− wδxk ≤ 2M1−xd ǫxd. (4.9)

In the sense of H¨_{older continuity on L}2 space, we conclude that a com-puted temperature gradient depends continuously on the data.

(32)

(33)

Chapter 5

Numerical Implementation

for the SHE

The problem of solving the sideways heat equation (SHE) the solution is severely ill-posed in the sense that the solution does not depend continuously on the data. The situation is actually caused by the time derivative which is an unbounded operator. Therefore we can obtain well-posed behavior of an ill-posed problem by substituting _∂t∂ with a bounded approximation. Thus it is convenient to establish a numerical implementation of the problem by discretizing space and time separately.

In this section we consider the time variable with equidistant grid points, and search for an approximation of _∂t∂ by using discretized Fourier Trans-form of the problem. Numerical methods, based on approximating the time derivative, have been considered by Eldén [12, 11], by Eldén and Regińska and Berntsson [13], see also [4].

5.1 A Method of Lines

The space discretization scheme should be selected independently of the time discretization. Therefore it is assumed that a time grid and a matrix representing the _∂t∂ have been selected. We illustrate in this section, that the resulting problem can be solved numerically, essentially as an ordinary differential equation in the space variable. In our implementation, the time derivative is computed by a spectral approximation or a spline function, and will be discussed more details later in this section. Now, we consider the “initial-value” problem in the linear form

u κux x = 0 κ−1I ∂ ∂t 0 u κux , 0 < x < 1, _{0 ≤ t ≤ 1,} (5.1a) with initial-boundary values

u(1, t) = gδ_m(t), ux(1, t) = hδm(t), 0 ≤ t ≤ 1, (5.1b)

(34)

u(x, 0) = 0. 0 < x < 1. (5.1c) The initial values for ux can be computed numerically by solving the

well-posed quarter plane problem with data u(1, t) = gδ_m_{(t) for x ≥ 1. In a formal} manner the solution of the system (5.1) can be written

u κux = e−B(1−x) g_mδ(t) κhδ m(t) ,

where B is an unbounded operator defined as B = 0 κ−1I ∂ ∂t 0 .

We already analyzed that the sideways heat equation (2.3) is severely ill-posed. The same phenomenon can be observed in the frequency domain by replacing (5.1a) with a family of ordinary differential equations parameter-ized by ξ. It is easily verified that the matrix corresponding to B is

b Bξ= 0 κ−1 iξ 0 .

The eigenvalues of b_Bξ are ±

q

iξ

κ. Therefore the spectrum of the operator B

is unbounded in the left-half plane. It implies that the operator e−B(1−x) is also unbounded and consequently solving (5.1) is an ill-posed problem, high frequency noise can blow up and completely destroy the numerical solution. Even if the data are filtered [7, 23], so that, the high frequency perturbations are removed, the problem is still ill-posed: rounding errors introduced in the numerical solution will be magnified and make the accuracy of the numerical solution deteriorate as the initial value problem is integrated. The methods have been investigated in different papers [12, 11, 13, 4] for solving the sideways heat equation numerically, where to obtain a well-posed problem the unbounded operator _∂t∂ have to be bounded. Thus we essentially replace

∂

∂t with a bounded approximation.

Although the differential equation is valid for all t ∈ R but, in real applications, data are only available for a finite period of time. In our case we normalized the data t ∈ [0, 1]. Thus by discretizing (5.1a) on an equidistant grid {tk}n−1_k=0 ∈ [0, 1], we obtain an initial-value problem for an

ordinary differential equation, U κUx x = 0 κ−1_I D 0 U κUx , _{0 ≤ x ≤ 1,} (5.2)

with initial data,

(35)

5.2. Time Derivative Discretization 23

where the matrix D is discretization of the time derivative; U = U (x) and Ux = Ux(x) are semi-discrete representations of the solution and its

deriva-tive. That is,

U (x) ≈ u(x, t1), u(x, t2), . . . , u(x, tn)T,

and,

Ux(x) ≈ ux(x, t1), ux(x, t2), . . . , ux(x, tn)T,

and, Gδ

m and Hmδ are vectors containing samples from gmδ and hδm

respec-tively on the grid. Thus (5.2) can be considered as a method of lines [16, 20]. The discretized system (5.2) can be solved using some standard method. We use an explicit fourth order Runge-Kutta method for solving the system of equations (5.2) which works well. It can also found in [12] that it is sufficient to use an explicit method, e.g. ode45 in Matlab, based on RK method.

It is noticeable that, the approach, a method of lines, is easy to generalize to more general equations. We can solve the initial value problem for any κ, where κ can only be a constant, it can only depend on x, it can also depend on both x and U . For a non linear system of ordinary differential equations, i.e., κ := κ(x, U ), and taking V = κ(x, U )Ux, we have

U V x = 0 κ(x, U )−1I D 0 U V , _{0 ≤ x ≤ 1,} (5.3)

with initial data,

U (1) = Gδ_m, and V (1) = κ(1, Gδ_m)H_mδ,

where κ(x, U ) = diag κ(x, U1), . . . , κ(x, Un) is a diagonal matrix. For

sim-plicity, we assume in this discussion for the nonlinear problem (1.1) that ρ(U ) = cp(U ) = 1.

Note also that, if D represents a discrete approximation of the time derivative, then D is a bounded operator [12, 11, 13]. Thus the problem of solving (5.3) is, in principle, a well-posed initial value problem. However, the original problem is ill-posed, and for a fine discretization (5.3) can extremely be ill-conditioned [19].

In the following sections we will discuss the approximation of the time derivative D by using some methods.

5.2 Time Derivative Discretization

Methods for approximating the time derivative, the matrix D in (5.3), have been considered in different papers. For example, discretization using finite differences was considered in [11, 12], Galerkin approximation using wavelets in [4, 13], a spectral approximation in [25, 13] and a method based on the approximation by cubic splines in [4]. In this section we have a discussion on the approximated derivative based on a spectral method and cubic splines.

(36)

5.2.1 A Spectral Method in the Finite Fourier Space

A spectral approximation is the approach, to write the numerical solution as a finite expansion in terms of some known type of functions, where the basis functions are associated with the spectrum of some operator [15]. Now, to deal with periodic solutions, a natural choice of basis functions can be trigonometric polynomials. This leads to Fourier methods, such as, Discrete Fourier Transform (DFT) that provides the decomposition of a vector in Rn _{in terms of the frequency components. Thus we can use the DFT to}

eliminate all high frequencies, i.e., |ξ| > ξmax, where ξmaxis the cut-off level.

A use of this method was proposed in different papers, e.g. [4, 13, 25]. We explain the DFT for approximating the derivative in details in the following paragraphs.

The DFT can be used to obtain very accurate approximations of the derivative of smooth functions. It is a linear mapping that operates on com-plex n-dimensional vectors in much the same way that the Fourier transform operates on functions onR [14, chapter 7]. To motivate it, we consider the problem of numerical approximation of Fourier transforms.

In order to use the Fourier transform

b g(ξ) = √1 2π Z _∞ −∞ g(t)e−iξtdt,

in numerical computations, we must approximate it by something that in-volves only a finite number of algebraic calculations performed on a finite set of data. This work is done in several steps.

First, we replace the integral over [−∞, ∞] by the integral over a finite interval. In other words, we assume that g vanishes outside some bounded interval [r, r + δ]. It will be convenient to assume that r = 0, which can always be achieved by replacing g(t) by g(t − r) [14, chapter 7].

Next, we try to calculate bg(ξ), not at every ξ ∈ R but only at a finite sequence of points contained in some bounded interval [−C, C]. We intro-duce the trigonometric interpolation of 2π-periodic functions u(x, ·) that are known in the finite interval ∈ [0, 1] at grid points tk = kht, such that,

uk = u(x, tk) for k = 0, 1, . . . , n − 1, where ht is the step size in the time

scale. Assuming that n is an even number with nht = 2π, we write the

approximation e u(x, t) = _√1 2π n/2−1_X k=−n/2 b ukeiξkt, (5.4)

where ξk= 2πk. Thus under the interpolation condition

e

(37)

the coefficients are given by [15, chapter 12],

b uk= 1 √ 2π n−1 X k=0 uke−iξktkht, k = −n/2, −n/2 + 1, . . . , n/2 − 1. (5.5)

This is the Discrete Fourier Transform. If we differentiate the interpolating function (5.4), and denote eut(x, t) = ev(x, t), we then have

e v(x, t) = _√1 2π n/2−1_X k=−n/2 b vkeiξkt= 1 √ 2π n/2−1_X k=−n/2 iξkbukeiξkt,

The trigonometric interpolation is very accurate for smooth functions, both for the function itself and for its derivatives [15, chapter 12]. If u(x, ·) is known only at the grid points, we want to compute approximations of ∂u_∂t at the grid points as well. Hence the algorithm is

1. Compute buk, k = −n/2, −n/2 + 1, . . . , n/2 − 1 by the DFT

2. Compute bvk= ikbuk, k = −n/2, −n/2 + 1, . . . , n/2 − 1

3. Compute vk, k = 0, 1, . . . , n−1 by the inverse discrete Fourier transform

In general, when the discrete Fourier transform is carried out by using the fast Fourier transform, the operation count of the numerical differentiation is O(n log n). However, we compute the derivative by executing the above steps. Hence the regularized solution (4.3) can be computed by solving (5.2), with the matrix D ≡ DF defined by

DF = FHΛFF,

where F is the Fourier matrix [22], and ΛF is a diagonal matrix defined

as, ΛF = diag(iξk), corresponding to differentiation of the trigonometric

interpolant, where the frequency components with |ξ| > ξmax are explicitly

set to zero. Thus the diagonal elements of ΛF are

(ΛF)k,k =

(

iξmax, when |ξk| < ξmax,

0, _{when |ξ}k| ≥ ξmax,

where the ξk are defined as in (5.4). This will filter the data and remove

the high frequency parts of the measurement error from the solution. The product of F and a vector can be computed in O(n log n) operations using the Fast Fourier Transform(FFT). Thus multiplication by DF can be carried

out as a FFT, followed by a multiplication by a diagonal matrix, and finally an inverse FFT. This work requires in total O(n log n) operations.

5.2.2 Approximating the Derivative by a Spline Function In this section we discuss another method to approximate time derivatives of discrete functions using “cubic spline approximations”. In some cases, it

(38)

β1 β2 β3 β4 B 1 B 0 B2 B -1 B3 B4 B5 β5 β6 β7 β0 β−1 β−2 β−3

Figure 5.1: The basis functions B_j3_{(t), for j = −1, . . . , 5. Every cubic spline,} defined on the interval [β0, β4] can be written, uniquely, as a linear

combi-nation of the basis functions B_j3.

is difficult to approximate the operator _∂t∂ using Fourier series, e.g. when the truncation is not periodic, instead cubic spline function can be a bet-ter choice. It often provides betbet-ter approximate result compared to using Fourier series.

The accurate estimations of the derivative _∂t∂ require that, the cubic splines for approximating u(x, t) must be all piecewise cubic functions. The cubic splines work as smooth curves with continuous first and second deriva-tives, and are associated to go through a set of points. In this process the whole interval is divided into a set of subintervals. Moreover, each spline curve of every subinterval is a third-degree polynomial, which are normally different cubic polynomial curve. These pieces of curves are so well matched where they are glued for every subinterval. We now illustrate the numer-ical procedure for the approximations of _∂t∂ using cubic spline. A similar approach of the method can be found in [4].

To make the procedure more easier, we substitute the function u(x, t) to u(t) for each forward step in the sideways heat equation, where u(t) is known only for the discretized grid points tk = kht ∈ [0, 1] as uk = u(tk),

k = 0, 1, . . . , n − 1, and ht is the time step size. Now, in order to obtain

an accurate estimation of the derivative ∂u_∂t we first find a cubic spline s(t) that approximates u(t), and then setting ∂u_∂t _≈ ∂s_∂t. To explain the work we

(39)

assume that, the points {β−3, · · · , βN +2} be a set of coarse grid, and s(t) is

a cubic spline function defined on the interval [β₋₃, βN +2], such that it has

natural breaks at {βj} and also s, ∂s_∂t , and ∂

2_s

∂t2 are continuous and smooth

for all breaks, i.e., in each subinterval [βj, βj+1] for j = −3, . . . , N + 1, the

function s(t) is a cubic polynomial. The set of all cubic splines can be represented using basis functions, B-spline, on the interval [β0, βN −1] such

that s(t) = N X j=−1 cjBj3(t − βj), (5.6)

where B_j3(t) are the B-spline basis functions, which is nonzero on the interval [β_j−2, βj+2], and can be seen in figure 5.1.

To fitting the coarse grid {βj} with respect to the discretize grid {tk}, it

is convenient to select them in such a way, so that, t0 = β0 and tn−1= βN −1.

Thus u(t) can be approximated by cubic spline s(t) in applying least-squares sense, which lead to solve the problem

min

c_∈RN +2kAc − bk2, (5.7)

where the elements of the matrix A and the vector b can be defined as

(A)i,j = B3j(ti), and (b)i = u(ti).

Since the basis functions B3

j have compact support with the banded matrix

A has only four non-zero elements on each row, thus only O(n) operations are needed for computing the vector c.

In the previous section 5.1 for solving the initial value problem (5.3), we need to compute derivatives of the vector U (t), that contain the approximate temperatures {u(x, tk)}. The technique presented in the previous paragraph

is used for this purpose.

It is remarkable, that, when we compute the approximation of the deriva-tive _∂t∂ by using a discrete operator D and by executing differentiation of ap-proximating cubic splines, a number of B-splines basis functions are used to control the accuracy of the computed derivatives and also the norm of opera-tor D. In fact, for solving the discrete initial value problem (5.3), if we com-pute the derivatives very accurately by using too many basis functions with respect to errors in the initial data, the solution can extremely be unstable. Thus, to stabilize the computations the grid parameter ∆β = min(βj+1−βj)

in the coarse grid {βj} acts as a regularization parameter [17, 19]. The

(40)

(41)

Chapter 6

Numerical Analysis and

Solution Procedures for

Well-posed Nonlinear

Problems

In this section we first present a discussion on the properties of stiff problems, and then we will describe solving procedures of a well-posed nonlinear initial value problem numerically by using implicit methods.

The problem u′ _{= au+f (t) is called stiff if Re(ah∗) is large and negative,} where h∗ is some measure of the time scale of the problem [26, chapter 1].

Systems of first-order equations, e.g. time dependent partial differential equation, heat equation or diffusion equation, of the form u′ = f (t, u) nor-mally show this stiffness behavior when the jacobian matrix fu has some

eigenvalues whose real parts are large and negative. To motivate it, if A is diagonalizable, A = P−1DP , the linear system u′ = Au + f (t) reduces to a set of uncoupled equations of z′ = λizi+ hi(t), where the zi and hi are

the components of z = P u and h = P f , respectively, and the λi are the

eigenvalues of A. Similar discussion on stiff problems can be found in [26, chap 1].

Although Euler’s method and other explicit multi step methods can be used to solve a ordinary differential equation which is stiff, but it will require excessively small time steps for to stabilize the problem numerically. Instead we will search some other multi step methods which are well designed for stiff problems. These are known as implicit methods. By using such methods a stiff problem, linear or nonlinear algebraic equation, can be solved with sufficiently large time steps. For the one dimensional problem it can easily be handled, since the involved matrices are tridiagonal. It is suggested that for solving a stiff problem numerically some 2nd order implicit methods, such as, Adams-Moulton method also well known as Crank-Nicolson’s method,

(42)

and Matlab solver ode23s which is based on a modified Rosenbrock formula, can be very useful.

Next, to solve the well-posed problem we consider a nonlinear PDE in 1D (heat or diffusion equation) in the general form:

∂u ∂t = ∂ ∂x κ(u) ∂u ∂x , _{0 ≤ x ≤ 1,} _{t ≥ 0,} (6.1a) with the Dirichlet boundary conditions:

u(0, t) = α(t), and u(1, t) = β(t), (6.1b)

and an initial condition:

u(x, 0) = 0 (6.2)

It is noticeable, in real applications, we most of the time deal with the problems which are situated in some other interval for x-domain (for exam-ple, 0 ≤ x ≤ L or, L1 ≤ x ≤ L2). For all of such cases, we first scale the

problem into the dimensionless form, as already discussed in section 2.3, and later on we will present the implementation of scaling in solving the problem numerically, see chapter 7.

Indeed, the initial-boundary value problem become nonlinear due to the nonlinearity of the governing differential equations, or of the boundary con-ditions, or both. Most physical problems in science and engineering are actually nonlinear. In the case of linear problems one can always solve a well-posed problem (i.e., when κ does not depend on u) and reduce the problem to one with the initial condition (6.2) above but that does not work for nonlinear problems. More natural may be to assume that the problem is at steady state. Which means that

ut(x, 0) = 0. (6.3)

In the steady state case, there is no time variation in the solution. Therefore the steady state solution satisfies the boundary value problem (BVP):

d dx κ(u) du dx = 0, (6.4a)

with the boundary conditions:

u(0) = α0, and u(1) = β0, (6.4b)

In order to find a initial steady state temperature u0(x, t) for the well-posed

problem (6.1), we first solve the BVP (6.4) by using finite difference method and next, we will solve (6.1) by discretizing space and time separately.

(43)

6.1. Space Domain Discretization 31

6.1 Space Domain Discretization

In this section we will discretize the space domain to deal with the right hand side of (6.1a), the boundary conditions, and also the steady state initial condition. In order to execute these, we use the finite difference method by which the continuous variable x is discretized and the derivatives are replaced by considering the difference approximations. Similar approach for a nonlinear problem is suggested in [9]. In the following discussion we transform the ODE problem (6.4) into an algebraic system of equations, which is in fact, can also be used for the left side of (6.1a).

We first discretize the x-interval and introduce a numbering of grid points: x_{0} x_{1} x_{i−1} x_{i} x_{N+1}x_{N} x i xN+1 0 hx h1/2 h1/2 1 x x N x 1 xi−1 xi+1 x 0

where hxis constant step size over the whole discretization, and hx/2 = h1/2

in the figure. Therefore the number of inner points can be defined by the relation

hx(N + 1) = 1,

in the interval [0, 1]. By introducing a second-order formula and by using the central difference approximation for the first derivative, we have

du dx(xi) = u(xi+ hx) − u(xi− hx) 2hx + O(h 2 x). (6.5)

Applying this approximation in (6.4a) we can derive

d dx κ(u) du dx (xi) = 1 hx κ(u)du dx xi+hx/2− κ(u) du dx xi−hx/2 + O(h2x) = 1 hx κ(uxi+hx/2) uxi+hx− uxi hx − − κ(uxi−hx/2) uxi− uxi−hx hx + O(h2x).

At an arbitrary grid point xi, we define uxi = ui, which leads to

d dx κ(ui) dui dx = 1 hx κ(ui+1/2) ui+1− ui hx − − κ(ui−1/2) ui− ui−1 hx + O(h2x).

(44)

The function κ contains the values ui+1/2and ui−1/2for u which are actually

not given in the grid points. Therefore from algebraic point of view

λa+1/2=

1

2(λa+1+ λa),

we then have a simple symmetric approximation with second order accu-racy throughout the grid points by setting them to the mean values of the neighboring points, that is, by linear interpolation,

ui+1/2 = 1 2(ui+1+ ui), ui−1/2= 1 2(ui+ ui−1). Thus we obtain d dx κ(ui) dui dx = 1 hx κ ui+1+ ui 2 ui+1− ui hx − − κ ui+ u₂ i−1ui− u_h i−1 x + O(h2x).

This relation is valid at all inner points, i.e., i = 1, 2, . . . , N . In order to get a matrix form, denoting by u the vector with components u1, u2, . . . , uN

then the above relation obtain to the form d dx κ(u) du dx = 1 h2 x tridiag κ ui−1+ ui 2 , −κ ui−1+ ui 2 + + κ ui+1+ ui 2 , κ ui+1+ ui 2 u. (6.6)

Now, for simplicity, we describe (6.6) to the standard form d dx κ(u) du dx = A κ(u∗)u+ b, (6.7)

where the vector u∗ is different from the vector u for every forward step that can be formulated as

u∗ = α(·) + u(1) 2 , u(1) + u(2) 2 , . . . , u(N − 1) + u(N) 2 , u(N ) + β(·) 2 T ,

and it contains (N + 1) elements, i.e., u∗ = (u∗₁, u∗₂, . . . , u∗_N, u∗_{N +1})T. Nor-mally, it is reduced to be length N when it perform in the matrix operations for A, where A a is symmetric, positive definite, and N x N tridiagonal matrix defined as

A = 1 h2

x

tridiag κ(u∗_i+1) _{, −}κ(u∗_i) + κ(u∗_i+1) , κ(u∗_i) ,

and b is N -elements vector including the boundary values at x0 and xN +1

such that b= 1 h2 x κ(u∗ 1)α(·), 0, . . . , 0, κ(u∗N +1)β(·) T .

(45)

6.2. Time Domain Integration 33

To clarify with the notations, the vector elements κ(u∗₁) and κ(u∗_{N +1}) are actually taken in every forward step by computing the thermal conductivity function κ(u) using u∗_{. For computing the initial steady state temperature}

u(x, 0) from (6.4), we take from the vectors α(·) and β(·) the first two elements as α0 and β0 respectively. Since the problem is nonlinear we can

use Newton-raphson iterative method to solve (6.4), where in each iteration the jacobian, i.e., the matrix A is tridiagonal. There are some built-in solver in Matlab language, e.g. “ fsolve ”, using this solver we can solve the problem (6.4).

Further, we apply all the ideas of space discretization to compute the temperature u(x, t) ≡ u(x, 0), a steady state temperature, in the time do-main [0, 1].

6.2 Time Domain Integration

When dealing with initial-boundary value problems, it is useful to introduce the Method of Lines (MoL), is also called semi discretization, where only the space derivatives are discretized (the time derivatives are kept). In brief, it means that u(xi, t) ≈ ui(t) where ui(t) is a time-dependent solution function

is associated to the space grid point xi. Mathematically, we can define the

illustration in the matrix form du

dt = A(u)u + b(t), u(0) = u0, (6.8) where A is a tridiagonal matrix and b(t) is a vector containing the bound-ary values. The initial vector u(0) is to be determined using (6.4). An effective approach to handle the solving procedure of the ODE (6.8) is that, first, discretize the space domain to yield a nonlinear system of algebraic equations, and then these equations to be solved at each time step by using some higher order (more than first order) implicit methods often used for the stiff problems. In solving the well-posed problem we will apply two dif-ferent methods, which both have second order accuracy in the x and the t directions. For the first one, we implement a recursive formula based on the trapezoidal method, e.g. Crank-Nicolson implicit method for the nonlinear system of equations, and for the next one, we use “ode23s” Matlab solver based on a modified Rosenbrock formula. The discretization of space deriva-tives are already illustrated in the previous section. Now, by applying these two methods in (6.8), we compute the solution u(x, t) in the time interval [0, 1].

6.2.1 Crank-Nicolson Method for Nonlinear Problems The discretization of space domain for (6.8) has already implemented by applying finite difference approximation of second order, consequently we

(46)

have a nonlinear system of algebraic equations. Now, to linearize the equa-tions we introduce an approach, where we first compute the coefficients of the system at the previous k level using Newton-Raphson iterative method. Meanwhile, we would have a recursive formula by which we can compute the coefficients at the level k + 1 by lagging it one step behind, i.e., by eval-uating the coefficients at the previous level k. A presentation in details on Newton-Raphson iterative process and lagging properties by one time step can be found in [32, chapter 3, 9]. Now to construct the recursive formula by the aid of Newton-Raphson iteration scheme, we consider the system of algebraic equations in the vector form

F(x) = 0, (6.9a)

where we need to find x = (x1, . . . , xN)T such that the system of equations,

F1(x), . . . , FN(x)

T

= 0, is satisfied. The Taylor series expansion of the system of equations (6.9a) is considered as

F(xj+1) = F(xj) + (∂F ∂x)(x

j+1_{− x}j_{) + . . .} _(6.9b)

Now for solving F(xj+1) = 0 for xj+1, the truncated Taylor series gives

xj+1= xj_{− (}∂F ∂x) −1_F(xj_), _(6.9c) where J_≡ ∂F(x j₎ ∂x =    ∂F1 ∂x1 . . . ∂F1 ∂xN .. . ... ... ∂FN ∂x1 . . . ∂FN ∂xN    .

Equation (6.9c) defines the Newton-Raphson method of iterations. Here the superscript j denotes the values obtained on the j’th iteration and j + 1 indicates the values on the j + 1’th iteration.

In order to make the procedure more simpler to develop the desired recursion formula, we consider the problem (6.8) to the system of form

ut= A(u)u. (6.10)

We now introduce a Trapezoid method for the ordinary differential equation

ut= F(u), (6.11)

where the Trapezoid method can be derived from the trapezoid rule for the integration. It has a simple form

tk+1 = tk+ htt,

uk+1 = uk+ht 2 F(t

(47)

6.2. Time Domain Integration 35

where htis the step length for time variable. Equivalently,

uk+1₋ht 2 F(t

k+1_{, u}k+1₎_{= u}k₊ ht

2 F(t

k_{, u}k₎_. _(6.12)

From (6.12) we can see that it is an implicit formula. Plugging (6.10) into (6.12), the system of equations obtain the form

uk+1₋ht 2 A(u k+1_)uk+1_{= u}k₊ ht 2 A(u k_)uk_, equivalently uk+1₋ht 2 A(u k+1_)uk+1 _{= G,} _(6.13) where G = uk+ ht 2 A(uk)uk is known.

Thus, to find uk+1 _{we must solve the system of equations}

F_{(u) = u −}ht

2A(u)u − G = 0, (6.14)

where, for simplicity, we have omitted superscripts. Newton’s method now can be written uj+1= uj_{− F}′(uj)−1F (uj), (6.15) where F′_p,q= ∂Fp ∂uq ; p = 1, . . . , N, q = 1, . . . , N. Now, since u′ =      u1 u2 .. . uN      ′ =      1 0 . . . 0 0 1 . . . 0 .. . ... ... ... 0 0 . . . 1     = I, and G is independent of u, we have

F′_{(u) = I −}ht 2 A(u)u ′ = I − ht 2 (A(u)) ′ u+ A(u)I.

If we assume that A(u) is a slowly varying function of u, we can replace A(u)′ by 0, and then we have

F′_{(u) ≈ I −}ht 2 A(u).

Omitting A(u)′ amounts to using Newton’s method with an approximate jacobian. It also means that we ignore the nonlinearity in this integration step. So with this approximation, there is no need to perform more than one Newton iteration.

(48)

Let u0 be the starting value in the numerical solution of (6.14). Then u1 = u0_{− I −}ht 2A(u 0₎−1 (I − h₂tA(u0))u0_{− G}0 = u0_{− Iu}0_{+ I −}ht 2A(u 0₎−1_G0_,

which is the same as

u1 _{= I −}ht 2 A(u

0₎−1 _{I +}ht

2A(u

0₎_u0_.

Thus, the integration of the differential equation can be written

uk+1_{= I −}ht 2A(u

k₎−1 _{I +} ht

2A(u

k₎_uk_, _(6.16)

which is nothing but Crank-Nicolson finite difference method used for solving nonlinear heat equations. It is a second order method in both the x- and the t- direction, implicit in time, and is numerically stable specially for stiff problems. Even if the nonlinearity is ignored in each particular integration step, it is taken care of over the whole integration by “lagging it behind one step”.

We now introduce the boundary conditions to give the above recurrence equation in the more general form

uk+1 _{= I −}ht 2A(u

k₎−1 _{(I +} ht

2A(u

k_))uk_{+ b}k_, _(6.17)

where the first and last elements of b contains the terms, the boundary values at x = 0 and x = 1 multiplied with κ(u)’s first and last elements respectively, and the rest of the elements of b are zeros. During the computation, A and bare usually changed in every next step. We have a more details about A, b and κ in the previous section. Finally, by applying the recursion formula (6.17) we can solve the nonlinear system of equations (6.8).

It is noticeable, since the boundary values are time dependent, there-fore in each time step it is needed to take care of those values. The work is normally done by introducing one dimensional interpolation, e.g. using the Matlab built-in function “interp1”, where a function uses polynomial techniques and fitting the supplied data with polynomial functions between data points and then evaluating the appropriate function at the desired interpolation points.

For solving nonlinear system, such as (6.8), the computation time is longer than that resulting from any linear system, because the temperature dependent thermal conductivity κ(u), density ρ(u) and specific heat cp(u)

need to be evaluated at each time level k + 1 from the previous time level k. Moreover, if the matrix A is full instead sparse, it would have comparatively

Numerical Solution of a Nonlinear Inverse Heat Conduction Problem

Dissertation

Numerical Solution of a Nonlinear Inverse

Heat Conduction Problem

Numerical Solution of a Nonlinear Inverse

Heat Conduction Problem

Abstract

Acknowledgments

Contents

Chapter 1

Introduction

Chapter 2

Investigation of Ill-posedness

2.1

Analyze the Ill-posedness with Sine Function

2.2

Ill-posedness Analysis in the Fourier Space

2.3

Standardization of the Region

Chapter 3

Stabilizing Assumptions

3.1

The Proof

Chapter 4

Regularizations

4.1

Error Estimate

Chapter 5

Numerical Implementation

for the SHE

5.1

A Method of Lines

5.2

Time Derivative Discretization

Chapter 6

Numerical Analysis and

Solution Procedures for

Well-posed Nonlinear

Problems

6.1

Space Domain Discretization

6.2

Time Domain Integration