# TRANSIENT HEAT CONDUCTION IN ADJACENT QUADRANTS

Document Sample

```					                   TRANSIENT HEAT CONDUCTION IN ADJACENT QUADRANTS
SEPARATED BY A THERMAL RESISTANCE

Donald E. Amos a , James V. Beck b , Filippo de Monte c
a
Albuquerque, NM 87110, USA
b
Department of Mechanical Engineering , Michigan State University , East Lan sin g , MI 48824, USA
c
Dipartimento di Ingegneria Meccanica , Energetica e Gestionale, University of L ' Aquila
Località Monteluco, 67040 Roio Poggio, L ' Aquila , Italy
___________________________________________________________________________________

Abstract Two linear, transient heat conduction problems set in quadrants 1 and 2 of the (x,y) plane are
solved. In each problem, the quadrants have distinct, constant physical properties and are separated by
an infinitely thin thermal resistance along the y-axis. Each region is initially at zero temperature. In
Problem I, constant fluxes F1 and F2 are specified along the x-axis boundaries to complete the problem
definition; while in Problem II, constant temperatures B1 and B2 are specified.

An attempt at a solution to Problem I by classical application of the Laplace transform results in an
integral representation for the temperature in each quadrant. Unfortunately the integrals only converge in
the wedges y  x which are close to the x-axis. However, a modification in the path of integration into
the complex plane leads to a complete solution in terms of integrals of co-error functions of a complex
variable. Details on high accuracy numerical evaluation of error functions and quadratures are provided.

Problem II is solved by manipulating the solution of Problem I.

Explicit solutions for over 12 special cases, including some quarter plane problems, evolve in terms of
other functions by taking limits or specializing parameters. Numerical experiments compare the general
solution with a number of these cases.

Keywords Two regions           Heat Conduction     Unsteady State   Laplace Transform

1. Introduction

We solve two linear, transient heat conduction problems for Regions 1 and 2 corresponding to quadrants
1 and 2 of the (x,y) plane. In each problem, an infinitely thin resistance with heat transfer coefficient h
separates Regions 1 and 2 along the y-axis. The two problems differ only in the boundary conditions
along the x-axis. In Problem I, the flux ki Ti ( x, 0, t ) / y  Fi , i  1, 2 is specified along the x-axis
(y=0) for both regions; while in Problem II, the temperature Ti ( x, 0, t )  Bi , i  1, 2 is specified. Zero
initial temperatures and constant physical parameters are assumed for each problem. The common
elements for both problems are shown in Figure 1. In the past, there was high motivation to only
explore solutions which lead to computation in real arithmetic because only real algorithms were
available to compute the special functions and paths in the complex plane were used only for theoretical
purposes. The emphasis here is on a much neglected mathematical approach to a class of problems
where we do not avoid paths in the complex plane, but take advantage of their desirable properties.
More recently numerical evaluation of a number of special functions of a complex variable has shown
up in the literature and computational libraries, which will mitigate the need for only real arithmetic. As
always, these classical problems serve as benchmarks in code verification work to assess accuracies of

1
more complex numeric codes and may have direct physical applications. For code verification work it is
extremely helpful to have accurate tabulated results.

In Sections 2 and 3 boundary conditions along y=0 are supplied to complete the problem definition and
Problems I and II are solved. Sections 4 and 5 are devoted to special cases of Problems I and II by taking
limits or specializing parameters in equations leading to the solutions of Sections 2 and 3. Section 6
explains how the mathematical functions are computed. Section 7 compares some computational results
of the general formulae from Sections 2 and 3 with the special case formulae from Sections 4 and 5.
These results are presented in tabular form to help in code verification work. Since several of the
manipulations in Sections 2 and 3 could be considered as heuristic, a direct verification of the solutions
is shown the APPENDIX.

2. Problem I - Flux Conditions Specified on y=0

To complete the problem definition started in (1.1) and (1.2), we need to specify the flux on the
boundary y=0 for both regions,
T
(2.1)                          ki i ( x, 0, t )  Fi , i  1, 2 .
y
We start the solution by applying the Laplace transform to all equations,

L{Ti ( x, y, t )}  Ti ( x, y, p)
 2Ti  2Ti  p
 2  Ti , i  1, 2
x 2
y    i
T2                   T1
 k2                  k1                 h[T2  T1 ] x  0
x    x 0
x    x 0

T                F        T                        F1
 k2 2                2 ,  k1 1                    
y        y 0
p       y             y 0
p

Since one surmises that the heat flow is one-dimensional at large distances from the y-axis, the plan of
attack is to write solutions for x   and add regular solutions to satisfy boundary conditions at
x = 0. We continue by separating variables and drop the subscripts for this step,

T  X ( x)Y ( y ) .
Then,

2
1  2 X 1  2Y                  p
         q 2  0, q 
X x  2
Y y 2

and the separation step gives
1 2 X          1  2Y
 q2           2 .
X x 2          Y y 2
This leads to the pair of equations
2 X                                        2Y
 ( 2  q 2 ) X  0,                        2Y  0
x 2                                       y 2
with solutions
 2  q2
X  e x              , Y  sin( y), cos( y) .
 2 Gi   1 Gi
The solution for x   (transform for the solution of                                                   with boundary conditions (2.1)) is
y 2
 i t
Gi Fi
Gi  Ci e  qi y ,  ki              on y  0, qi                           p / i
y  p
Then, for x   ,
Fi
(2.2)                                                  Gi ( y, p)              e qi y , i  1, 2
pqi ki
Because the problem is linear, we can form the general solution from the separation of variables
solutions by summing linear combinations of solutions which (formally) satisfy the flux conditions (2.1)
along y = 0. Then, with coefficients Ai (  ) ,

Ti ( x, y, p)  Gi ( y, p)   Ai ( )e
 x qi2   2
cos( y ) d  , i  1, 2 .
0

Now, apply these solutions to the boundary conditions along x = 0 to determine A1 and A2 . The results
are
                                                             

 A1 ( )[k1 q1    h]cos( y )d    A2 ( )h cos( y )d   h G2  G1 
2    2
        
0                                                             0
                                

 A ( )h cos( y)d   A ( )[k                       q2   2  h]cos( y )d   h G2  G1  .
2
1                                  2           2                                          
0                                 0

These integrals are Fourier transforms which can be inverted according to the formulae
                                  

g ( y )   f ( ) cos( y )d  ,   f ( )   g ( y ) cos( y )dy
0
2                   0

and we have a pair of equations for A1 and A2 ,

 k q 2   2  h  A ( )  hA ( )  2h G ( y )  G ( y )  cos( y )dy
 1 1              1          2
  2
0
1      
(2.3)                                                                                      
2h
hA1 ( )   k2 q2   2  h  A2 ( ) 
 
2
G2 ( y )  G1 ( y )  cos( y )dy.

                               0

Using the first of the Fourier transform pair
                                                      
q        cos( y )              
e          cos( y )dy  2     , q 2         d   e  qy
 qy

0                        q  2
0 q 
2
2
we have

 F /k    F /k  1
R ( , p)   G2 ( y)  G1 ( y)  cos( y)dy   22 22  21 12 
                  
0                                    q2   q1    p
and the solution for A1 ( ) and A2 ( ) from (2.3) is
3
2h
A1 ( )         R ( , p)  k2 q2   2  / D
2

                           
2h
A2 ( )      R ( , p)  k1 q12   2  / D
                          
D  k1k2 q12   2 q2   2  hk1 q12   2  hk2 q2   2 .
2                              2

Then,
 x q1   2

 F /k    F /k 
2
2hk2       e
T1 ( x, y, p)  G1 ( y, p)                                      q2   2  22 22  21 1 2  cos( y )d  , x  0
2

p     0
D                          q2   q1   
(2.4)
2hk1  e x q2                                F /k    F /k 
2       2

p  D
T2 ( x, y, p)  G2 ( y, p)                                      q12   2  22 22  21 1 2  cos( y )d  ,         x0
0                                         q2   q1   
where
D  k1k2 q12   2 q2   2  hk1 q12   2  hk2 q2   2
2                              2

F
Gi ( y, p)  i e  qi y ,     qi  p /  i , i  1, 2 .
pqi ki
We have the Laplace transform of the solution. At this point, the classical approach would be to apply
the complex inversion integral
  i
1

L T j ( x, y , p ) 
1

 e T j ( x, y, p)dp
2 i   i
pt

to obtain a solution. This means identifying poles and branch points, selecting branch cuts and applying
Cauchy’s theorem to close the contour around the singularities. These cuts are generally taken in the left
half p-plane, Re(p)<0, so that pt gives exponential convergence for t>0. If the cuts are chosen properly,
the arcs connecting the inversion path   ip, -<p< to the extended contour vanish and only the
residues at the poles and the contributions from the two sides of the branch cuts remain. In the classical
cases, the contours on each side of a cut are conjugates whose difference yields a purely imaginary
result. Division by 2 i makes the result real. Examples of this procedure are contained throughout [9],
especially in Chapter XI and Appendix I. In the case at hand, the solutions for Problems I and II would
be represented as double integrals, making numerical evaluation difficult. In the developments to follow,
we get a direct inversion from tables, which limits the solution to one integral, but requires integration in
the complex plane with complex arguments to obtain the solution. We do not avoid complex contours
because the functions involved can be computed in the complex plane from accurate, numerically stable
algorithms. Further, the integrals converge exponentially over a large portion of the (x,y) plane.

To proceed, we simplify the inversion with the substitutions
1 iz
   p ,  j  1/  j ,                 j  1, 2                                cos( z )     e  eiz  .
2         
Then

1                  e  1 (  )   p
T1 ( x, y, p)  G1 ( y, p) 
     
      p3 / 2 ( p   )
H1 (  ) d  ,   x0


1                 e  2 (  )       p
T2 ( x, y, p)  G2 ( y, p) 
    
       p3 / 2 ( p   )
H 2 (  )d  ,    x0

where

4
Fj  j e y p /  j
G j ( y, p)                                , j  1, 2
kj       p3 / 2
 j ( )  x         2   2  iy 
j

h  k1 12   2  k2  2   2 
2
                                                                  
                                 h     1

1


k1 1  
2     2
k2  2  
2
  2       k1 1  

2
 2
k2  2   2
2


h        F2 / k2     F /k 
H j ( )                       2          21 1 2  .
 2   2  2       1   
2
kj       j

We have exposed each transform in terms of simple functions of p, and we can now invert. Tables of
inverse transforms give
 e p                
(2.5)                                                      L1  3 / 2   2 t ierfc     
 p 
                    2 t 
      e p  2 t e /(4t )  1            e   t    
2                                                      2

1
L  3/ 2                       2   erfc        2 erfc      t.
 p ( p  ) 
                                    2 t             2 t     
For computational purposes we manipulate this last equation into the form
     e p        1  2                                                                      
1
(2.6) L  3 / 2
 p ( p  )  
  e 2 t e ierfc    e

2              1 (   t )2
   
erfc    t  e erfc    
2

              
                 
using    /(2 t ) and ierfc( )  exp( 2 ) 1/ 2   erfc( ) . This form is suitable for large  (large h). For
small  (small h or large  ), we resolve the indeterminant form in (2.6) by expanding the expression
into a power series in  t :
(  t )n d n  z
             

d
(2.7) e(   t ) erfc    t  e erfc    ( t ) e z erfc  z                     e erfc  z   .
2                           2                                 2                                                   2

                 dz            n      z       n2       n!       dz           z 

Now, [1, 7.2.9 ] gives
d n  z2
e erfc  z   (2)n n ! i n erfc  z  e z , n  0 ,
2
(2.8)                           n              
dz
where i erfc(z ) is the iterated co-error function. With    /(2 t )
n

      e p             
  4te  (2 t ) [e i erfc( )]
 2       n2  2 n
L1  3/ 2
 p ( p  ) 
                       n2

for small  t . The ratio test using (6.3) shows an infinite radius of convergence. The explicit solution
can be summarized by

 U 1 ,  t V (  ) d  ,

2 t
T1 ( x, y, t )  G1 ( y, t )                                                          x0

2

(2.9)
 U  2 ,                      

T2 ( x, y, t )  G2 ( y, t ) 
2 t
    

t V1 (  ) d  , x  0

where

5
2 Fj  j t            y 
G j ( y, t )                   ierfc         ,            j  1, 2
kj              2  jt 
        
 j   j (  ) /(2 t ) ,  j (  )  x  2   2  i  y,
j                                                  2  1/  j ,
j

1                                1
  h ,          ( )                                                               ,
k2  2   2
2
k1 12   2

1                                kj 2   2                    F2 / k2     F /k 
V j ( )        H 3 j (  ) 
j
 2          21 1 2  ,
                    k1     k2
2
1
2
    2  
2
2
2 21   

 2
 e ierfc   

1  2
2 t e erfc    e

(   t )2
erfc    t  for 2 t 


1
4
               
 2
U ( ,  t )  e                                         
                        2 t  (2 t )n 2 [e i n erfc( )] for 2 t  .
2                    1

                              n2                                           4
The break point of ¼ was chosen so that the series would converge rapidly and, at the same time, not
allow too much loss of significance in the original expression when 2 t is small. In order to assess
the convergence of the integrals in (2.9), we examine the behavior of  j   j (  ) /(2 t ) in U and see
whether the relation [1, 7.2.14]
    2          
e z i n erfc( z )  O                 for z  , n  0 , arg z  3 / 4
2
(2.10)
  (2 z ) n 1 
can be applied to not only ensure convergence, but also to enhance convergence from the factor (2 z )n 1 .
 2                Re( 2 )
Now, with  j   j (  ) /(2 t ) in U ( ,  t ) , we note that | e                             j
| e               j
and convergence of (2.9)
depends on the sign of Re{ 2 }
j

 j 2  [ x 2  2  ( x 2  y 2 )  2 ]/(4t )  i [2 x y   2   2 ]/(4t )
j                                            j
(2.11)
j
Re{ 2 }  [ x 2  2  ( x 2  y 2 )  2 ]/(4t ) [( x 2  y 2 )  2 ]/(4t )  O(  2 )
j

j
If x  y , we have Re{ 2 }  0 and we get exponential convergence; while if x  y , Re{ 2 }  0 , and
j
we get exponential divergence. Notice also that Re{ 2 }  0 implies arg( 2 )   / 2 and
j

arg( j )   / 4 , making (2.10) applicable with  j  O (  ) . Then, with the first term of the series form
of U and n=2 in (2.10),
(2.12)            t  O(1/  ) , U ( j ,  t )  O e                          Re( 2 )
j
/
4
 , V ( )  O(1/ 
j
2
),
and it follows that each integrand with x  y is
 e Re( 2j ) 
(2.13)                     U ( j ,  t )V j (  )  O                for    .
 6 
              

For the case where x  y , Re{ j }   j x is still positive, but independent of  . The result is that the
2       2 2

convergence is only algebraic like O 1/                      6
 . In summary, the solution represented by (2.9) converges
only for the wedges y  x , leaving the solution for y  x as yet undetermined.

6
This division of the (x,y) plane, where one obtains only half of the solution, is rather unusual in heat
conduction problems though not unknown [5, 6.8.3 ]. The analysis to follow is devoted to obtaining the
complete solution. We employ a technique often used in functions of a complex variable to extend the
region of validity by regarding an integral as a contour integral in the complex plane and suitably
defining a contour to obtain the analytic continuation of the parameters to nearby regions. We carry out
the details as follows. If we re-define the contour from the real  -axis to the upper half plane, thus
adding a positive imaginary part to  , we can get an additional positive term added to the real part of
 2 (and  2 ) which will enhance the convergence for positive values of y. We move off the real line for
j        j

the contour C1  C2 shown in Figure 2. In this analysis, we avoid the branch points at   i  j , j  1,2
1       2
and always take contours within regions of analyticity. If we show that the integrals on C R and C R
vanish as R   , then we have, by Cauchy’s theorem, that the contour C1  C2 is equivalent to the
Fourier integral on the real line,

(2.14)                         

f j (  )d  =     
C1 C2
f j (  )d    f j (  )d    f j (  )d  ,
C1               C2
j  1,2

where
f j (  )  integrand in each of the integrals                 j  1, 2 .

In addition to showing that the integrals along C R and C R vanish, we also show that, while  values on
1       2

C1 are not the conjugate of  values on C 2 , the function values f j (  ) along C1 are the conjugates of
those along C2 , making the integral on C1  C2 real. To these ends, we parameterize C1 and C2 to get
C1 :   rei (  )   re i   r
C2 :   rei  r ,   cos  i sin   c  is

Notice that in the formulas for f j (  ) we need to consider only  and  2 :
2
C1 :     2  (r )2  (r )  conjugate of  2 for  on C2
(2.15)
 i  y  ry(i )  ryi  conjugate of  i  y for  on C2
The crucial result is
(2.16)                                j (r )   j (r ) or  j (r )   j (r ) .

7
We are now in a position to show that the integrals on C R and C R vanish as R   . For the immediate
1       2

 2           Re( 2 )
analysis we restrict the variables to x  y  0 and defer the case x  y . For C R , with | e
1                                                                              j
| e          j

and   Rei in (2.10) we get

(2.17)
 j 2   x 2 (  2   2 )  y 2  2  2i x y   2   2  (4t )
         j                                j
                   2
j   x  ( x 2  y 2 )  2i x y   2
2
                                             (4t )
Re j 2           2
j   x  ( x 2  y 2 )cos2  2 x y sin 2  R2 (4t ), , R  
2
                                                                
This latter relation shows that Re j   0 for x  y and 0       / 4 . Therefore the integral on C R
                      2                                               1

vanishes for R   . For C R we replace  by   Rei (  )   R e  i and note the results of (2.15). This
2

leads to Re j 2   0 also. Thus (2.14) is valid.


With (2.16) and the fact that all other functions in (2.9) are analytic and real on the positive real axis,
we use the reflection principle to move the conjugation symbol from the function argument to the
function itself to conclude that f ( r )  f (r ) . For the integrals, we have
0                                                                                                                            


C1C2
f j (  )d     f j (r ) dr   f j ( r ) dr   f j (r ) dr   f j ( r ) dr   f j ( r ) dr   f j ( r ) dr.
                         0                     0                        0                                 0                     0

Since the integral on C1  C2 is represented as the sum of complex conjugates of one another, we get
twice the real part along C2 and we can write the final solution for 0     / 4 as

 Re U 1 ,  t V (  ) 

4 t
T1 ( x, y, t )  G1 ( y, t )                                                                          dr ,       x0
        0
2
 r
(2.18a)
 Re U  2 ,                                

T2 ( x, y, t )  G2 ( y, t ) 
4 t
       0

t V1 (  )
  r
dr , x  0

where
2 Fj  j t            y 
G j ( y, t )                   ierfc        ,             j  1, 2
kj                2 jt 
       
x  2   2  iy 
j                                        ,  2  1/  j
j
j
2 t
1                    1
  h ,  (  )                                                        ,
k2        2
2
2
k1 12   2

kj 2   2              F2 / k 2     F /k 
Vj ( )                                                        21 1 2 
j
 2
k1 12   2  k 2          22   2   2       1   
2

 2
 e ierfc   
2 
1  2
2 t 
e erfc    e(   t ) erfc    t  for 2  t 
2



1
4
                   
U ( ,  t )  e                                    
                                                                               1
2  t  ( 2  t )n2 [e i n erfc( )] for 2  t 
2


                                 n2                                           4
  cos   i sin  .

8
Computationally, (2.18a) is best for small t, 0  t  1 . We re-write (2.18a) for large t by changing the
variable of integration from r to r t and compute for 1  t   :

 Re U 1 ,  t V (  ) 

4
T1 ( x, y, t )  G1 ( y, t )                                                      dr ,   x0
   0
2
 r
(2.18b)
 Re U  2 ,                  

T2 ( x, y, t )  G2 ( y, t ) 
4
    0

t V1 (  )
  r
dr , x  0

where

2 Fj  j t          y 
G j ( y, t )                 ierfc         ,      j  1, 2
kj              2  jt 
        

                                
 j  x  2   2  iy  /2 ,  2  1/( j t ) ,  t  h (  )
j                     j

  cos  i sin .
Functions V j (  ), U ( ,  t ) and  do not change and are still defined by (2.18a), but notice that  j ,
 2 , and  have new definitions for this form of the solution.
j

The analysis above applies to the case where | x | y  0 . To complete the analysis, we consider the case
for y>|x|. If y>|x|, the factor ( x 2  y 2 ) in (2.17) is negative and selecting  / 4       / 2 puts 2
j
and 2 in the second quadrant and again Re{ 2 }  0 . But in order to apply Cauchy’s theorem directly to
the case where    / 4 , we would have to make the path   Rei connecting the real axis to C 2 go
through the region 0     / 4 where Re j 2  could be negative for y>|x| and there could be

exponential divergence for R   . There is a way out of this dilemma by noticing that if    / 4 ,
cos2  0 , and the term involving ( x 2  y 2 ) in (2.17) is zero and there are no restrictions on either |x| or
y, other than they both be positive. By taking    / 4 , we could stop here and say that we have the
complete solution for y>0, but we want to take advantage of the higher rate of convergence offered by
the path where    / 4 and the term ( x 2  y 2 )cos 2  0 contributes. The way is to apply Cauchy’s
theorem to the sector composed of the 45o line   Rei / 4 , C 2 , and the arc CR :   Rei ,  / 4    
1

in the first quadrant and the mirror image in the second quadrant as shown in Figure 3. But the
j
observation made in the third sentence of this paragraph establishes what we need, namely Re{ 2 }  0 .
Then with R   , the integrals on the contour CR  CR of Figure 3 vanish. This establishes that the
1    2

contours   rei / 4 and C 2 are equivalent with the restriction y  x  0 and similarly for C1 in the
second quadrant. Thus, (2.18) also applies for y>|x| with  / 4     / 2 . The line    / 4 is a
crossover boundary between the cases | x | y  0 and y>|x|. Numerical experiments indicate that (2.18)
works quite well with

 / 6 for x  y

(2.19)                                                           .
 / 3 for x  y


Since the forms of  , U , and V j are the same here as they were for (2.9), the estimates (2.12) apply.
Therefore with the relation Re( z)  z and  2 defined by (2.18), the analog of (2.13) is
j

9
 e Re( 2j ) 
(2.20)                                          
Re U ( j ,  t )V j (  )  U ( j ,  t )V j (  )  O 
 6 
 for    .
              
One has to remember that these are ultimate convergence rates. For x and y close to zero,  must be
very large to see the exponential decay; in the interim, the behavior is only O(1/  ) which is the
3

ultimate rate for x=y=0. This influences greatly how one approaches a numerical integration.

Since  is arbitrary in each of the sets 0     / 4 or  / 4     / 2 , the solutions for different values
of  should be the same within each set. We can show this by taking rays C2 ( ' ) and C2 ( '' ) with
angles  ' and  " in the same set and form a closure with an arc of radius R and let R   . Then on
the arc
  Rei  ,  '     " with  ' and  " in [0,  / 4] or [ / 4,  / 2) but not both.

It is apparent from the product UV j in (2.18), with the replacements r by R and  by  in (2.20), that
the line integral along the arc vanishes rapidly as R   . With Cauchy’s theorem we are left with equal
integrals on C2 ( ' ) and C2 ( '' ) .

3. Problem II - Temperatures Specified on y=0

Problem II is the same as Problem I except for the boundary conditions along the x-axis, y=0. In this
problem, we complete (1.1) and (1.2) by specifying a temperature B j , j  1, 2 on each boundary in
place of the flux condition. Then, the boundary conditions on y=0 are
(3.1)                                T j ( x, 0, t )  B j , j  1, 2 .
While we can re-solve Problem II in the manner of Problem I, there is an elegant way to obtain the
solution in a simple manner. If we differentiate, with respect to y, the differential equations (1.1) and
boundary conditions (1.2) which contain y, we observe that the substitutions

T j ( x, y, t )                         Fj
(3.2)                   ˆ
T j ( x, y, t )                           ,        Bj                  ,   j  1, 2
y                                 kj
give
ˆ
1 T2                                              ˆ
1 T1
ˆ
 2T2                                               ˆ
 2T1            ,         x  0 , y>0
 2 t                                          1 t
ˆ
T2 ( x, y, 0)  0                                 ˆ
T1 ( x, y, 0)  0
T ˆ     B                                 ˆ
T       B
2                2                       1               1
y 0                                     y 0

ˆ
T2                       ˆ
T1
 k2                    k1                     ˆ    ˆ
 h[T2  T1 ] x 0
x     x 0
x    x 0

and we have the specifications for Problem II. Correspondingly, differentiation of the solution (2.18)
with respect to y along with the corresponding replacements yields the solution for Problem II. We start
by focusing on the solution for Problem I in (2.18a).

We need to differentiate (2.18a) with respect to y and apply (3.2). The formula for the derivative of the
integrand with respect to y can be computed from

10

Re U ( j ,  t )V3 j (  )  1
2 

U ( j ,  t )V3 j (  )  U ( j ,  t )V3 j (  )  , j  1, 2


                                   1  U  j                     U  j                 
y

Re U ( j ,  t )V3 j (  )    
2   j y
V3 j (  ) 
 j y
V3 j (  ) 

                                                   
 j     i      U
                U B ( ,  , t ) , U B ( ,  t )  erfc( )  e 2  t   t erfc(   t ).
2
,
y     2 t 
Then,

y
                             i
Re U ( j ,  t )V3 j (  )  U B ( j ,  t )V3 j (  )   U B ( j ,  t )V3 j (  )   /(2 t )
2                                                              

 U B ( j ,  t )V3 j (  )   U B ( j ,  t )V3 j (  )   /(2 t )
i
2                                                             

Now we have the difference of a complex number and its conjugate, which yields (2i) times the
imaginary part,

y
                                             
Re U ( j ,  t )V3 j (  )   Im U B ( j ,  t )V3 j (  )  /(2 t ) , j  1, 2 .     
Further, differentiation of G j from (2.18a) gives
G j                     Fj  y 
( y, t )   erfc         .
y               kj       2  jt 
        
Now (3.2) with these derivatives takes (2.18a) over to the solution for constant temperature boundary
conditions

2(1)3 j
(3.3a)       ˆ                 ˆ
T j ( x, y, t )  G j ( y, t ) 
   Im U ( j ,  t )V3 j (  )  dr, j  1, 2
0
B
                                          
  r

 y 
where                                       ˆ
G j ( y, t )  B j erfc          , j  1, 2
 2  jt 
        
 j ( )        x        2   2  i y                             1
j                                                   , 2 
j

2 t                       2 t
j
j
              1                         1       
  h (  ),            ( )=                                                    
 k1 12   2
                          k2           2   2 
2

kj 2   2                    B2           B   
Vj ( )                                                     2 1 2
j
 2
k1 12   2  k2           2   2  2      1   
2                2




e 2 erfc    e(  

t )2

erfc    t 
                 for 2  t 
1
4
U B ( ,  t )  e
2

                    
1
        2 t  (2 t ) n 1[e i n erfc( )]                             for 2  t 
2


                    n 1                                                              4
  cos   i sin  , 0     / 2 .
Here we have used (2.7) with (2.8) for the series form of U B and, as before, we take
 / 6 for x  y

                   .
 / 3 for x  y


11
This solution agrees with the Laplace transform solution. The regions of uniform convergence are the
same as those for Problem I, x  0, y  0 , with the integrand O(1/  ) for    .
4

As with Problem I, this expression is computationally best for small t, 0  t  1 . We re-write (3.3a) for
large t by replacing r with r t and compute for 1  t   :

2(1)3 j
(3.3b)       ˆ                ˆ
T j ( x, y, t )  G j ( y, t ) 
        Im U ( j ,  t )V3 j (  )  dr, j  1, 2
0
B

  r

where
 y 
ˆ
G j ( y, t )  B j erfc          , j  1, 2
 2  jt 
        

                    
 j  x  2   2  i  y / 2 ,  2  1/( j t ) ,  t  h (  )
j                       j

  cos   i sin  , 0     / 2 .
Functions V j (  ) , U ( ,  t ) and  do not change and are still defined by (3.3a), but notice that  j ,
B

 2 and  have new definitions for this form of the solution.
j

In this case there is a steady state solution formed by taking t  or letting  j  0 in (3.3b). In the
steady state case, the asymptotics for   r and r  0 are
  3 j
                     
2              2
e       j
e   j

V3 j (  )    3 j            and U ( j ,  t )  1  erf ( j ) 
B
 O             
2   r                                           ( j   t )     (   t )3 
 j           
kj                                                 C     1 1
j                B2  B1  ,    j  (  ) / 2 ,  t , C  h    ,   x  iy.
k1  k2                                                      k1 k2 
At first glance, it looks like there is a singularity at r  0 . However, the product written out explicitly in
powers of r,
 3 j                          3 j  3 j
U B ( j ,  t )V3 j (  )          1          O   3                  [  1/ C ]  O  r 2  ,
r         C                     r        
shows that the real part has the singularity, but the imaginary part, which is used in the integrand, is
analytic at r=0.

4. Special Cases for Problem I

Problems I and II are quite general and one expects simpler formulae, at least in some cases, when the
parameters are specialized. Not only are these formulas valuable for computation, but they serve as a
numerical check against the general solutions represented by (2.18) and (3.3). In some of the solutions
written below, reference is made to [2, Chapters 2 & 3, Folders 5, 10b, 21, and 22] where analytic and
computational forms of appropriate integrals are derived.

Thermal Barrier (h=0) on x=0

In this case there is no heat transfer between Regions 1 and 2 and, with h=0,   0 and U=0, the solution
(2.18) is one dimensional,

12
2 Fj  j t                y 
(4.1)                                  T j ( x, y , t )  G j ( y , t )                                 ierfc          , j  1, 2 .
kj                 2  jt 
        

Perfect Contact ( h   ) on x=0

If we take h   in the general solution (2.18a), then    and only the leading term of U remains.
Thus for perfect contact between Regions 1 and 2, we have a simplification for small t



 Re ierfc  V                            
4 t
(4.2)               T j ( x, y, t )  G j ( y, t )  ( 1)3 j                                                       3 j   (  )                dr ,   j  1, 2

j
  r
0

2 Fj  j t                y 
where                                               G j ( y, t )                           ierfc         ,               j  1, 2
kj                2  jt 
        
 j ( )
j                        ,  j ( )  x                    2   2  iy  ,  2  1/  j
j                  j
2 t

kj 2   2                    F2 / k2     F /k 
V j ( ) 
j
 2          21 1 2 
    2       1   
2
k1     k2
2
1
2                  2
2
2

  cos   i sin  , 0     / 2

The analysis leading to (2.18) is modified some by the limit h   . For at least one of x or y not zero,
 Re ( 2 )
/  ) with an interim rate of O(1/  ) . For x=y=0,we
4                                                            2
the ultimate rate of convergence is O(e                                     j

have an ultimate rate of only O(1/  2 ) .

The choices for  remain the same as that for (2.18a). A similar form exists for large t analogous to
(2.18b).

Equal Diffusivities

In the case where 1   2   , the inversion of the Laplace transform (2.4) becomes easy since q1  q2
and the radicals         q12   2 , q2   2 are equal. Then we have
2
q12   2  q2   2  q2   2 , q  p / 
2

and (2.4) yields
Fj e qy                          2hk3 j       
e x    q2   2
cos( y)                     F   F
 k k                                                    d   2  1  , j=1,2 .
3 j
(4.3) T j ( x, y, p)                  (1)
pqk j                               p        0                  q 2   2  h(k1  k2 )   q 2   2   k2 k1 
   1 2

Then
                  x q2   2
              
x
w                 
1                e              h           1          e      h            
L                                            L                                  
  k1 k2 q    h(k1  k2 )   q    
2   2                   2   2
  k1 k2 w  h  (k1  k2 )  w 
                          
                                                                           w p   2
 a p   
             2 t      e     b                     a             a 
1
          e   t U  0, , t   U  b, , t  
2
=             e               L 
k1  k2
    
 p  b p

k1  k2           2             2 

13
where
1   1 
a / 2  x /(2  ), b  h     ,                      U ( ,  , t )  e 2  t erfc( t   / t ) .
2

 k1 k2 
Now the solution, taking into account the coefficient 1/p, can be written as

2 Fj  t          y  2(1) k3 j
3 j
 F2 F1  t   a           a 

T j ( x, y, t )                                                               U  0, ,    U  b, ,     e   cos( y)d  d
2
ierfc                                   
kj             2 t     (k1  k2 )              k2 k1  0   2           2  0
But,
                                              b2
 1         
e
 2a
(4.4)                                                                cos(b)d                e       4a

0
2   a
gives
2 Fj  t          y  (1) k3 j 
3 j
 F2 F1  t   a           a  e
c /   2

T j ( x, y, t )               ierfc                                        U  0, ,    U  b, ,           d , j=1,2 .
kj              2 t     (k1  k2 )              k2 k1  0   2           2  

Now these integrals have been evaluated in [2, Chapters 2 & 3, Folder 22] with the designation
I 22 ( ,  ,  , t ) . The final result is

2 Fj  t          y     (1) k3 j   F2 F1    a
3 j
         a         
(4.5) Tj ( x, y, t )                   ierfc                            I 22  0, , c, t   I 22  b, , c, t   ,             j  1, 2
kj            2 t       (k1  k2 )  k2 k1    2                   2         
1      1 
a / 2  x /(2  ), b  h     , c  y /(2  ) .
 k1 k2 
From [2, Chapters 2 & 3, Folders 22 and 10b] we also have,

e  c w erfc(aw / 2)
2   2
 a                  a                                           1
(4.6)                        I 22  0, , c, t   2I1c  c, , T   2                       dw , T     .
 2                  2 
2
T           w                    t

This case can also be used for regions with common physical properties by taking k1  k2 . The formula
does not result in any great simplification. However, there is a simplification if we also take h   .

Common Physical Properties in and Perfect Contact Between Regions 1 and 2

This is an unsteady state , upper half plane problem with a jump in boundary values on y=0. In this case
we take 1   2   and k1  k2  k in (4.5) and let h   . Then b   , I 22  b, a / 2, c, t   0 , and (4.5)
with (4.6) goes over to

2 Fj  t          y  (1)3 j 
 F2  F1  I1c  c, , T  , T  1/ t ,
a
(4.7)           T j ( x, y, t )               ierfc                                                                    j  1, 2
k              2 t         k                    2 
a / 2  x /(2  ) , c  y /(2  )

where computational forms and properties of I1c ( ,  , T ) and I 22  ,  ,  , t  are found in [2, Chapters 2 &
3, Folders 10b and 22](see (4.10) for the continuity on x=0). Note also that this case can be computed
from (4.2).

14
Quarter Plane with Newtonian Cooling on x=0

If we take F2  0 and k2   in (4.5), then T2 ( x, y, t )  0 and we have the boundary condition
T1
k1             hT1 x 0
x x 0
on Region 1 which is the formula for Newtonian cooling on x=0. The formula for the equal diffusivities
case (4.5) goes over to
2 F1 1t        y  F1 1   a                        a         
(4.8)             T1 ( x, y, t )           ierfc                  I 0, , c, t   I 22  b, , c, t   ,
k1          2  t  k   22  2
                      2       
     1     1

a / 2  x /(2 1 ), b  h 1 / k1 , c  y /(2 1 ) .

As before we can make use of (4.6) and refer to [2, Chapters 2 & 3, Folders 22 and 10b] for properties
of these functions.

Quarter Plane with Zero Temperature on x=0

We consider (4.8) in the previous case further. If h   then b   and I 22  b, a / 2, c, t   0 . Using
(4.6), this gives a formula for the case where T(x,y,t)=0 on x=0:
2 F1 1t        y  2 F1 1 c  a 
(4.9)               T1 ( x, y, t )           ierfc             I1  c, , T  , T=1/ t
k1         2  t   k1      2 
   1 

a / 2  x /(2 1 ),         c  y /(2 1 ) .
T(0,y,t)=0 follows from

(4.10)                                     I1c (c, 0, T ) 
ierfc(cT ) .
T
Computation and properties of I1c are found in [2, Chapters 2 & 3, Folder 10b]

5. Special Cases for Problem II

As noted in Section 4., we expect some simpler results by specializing parameters. This section presents
these results for Problem II.

Thermal Barrier (h=0) on x=0

In this case there is no heat transfer between Regions 1 and 2 and, with h=0,   0 ,and U=0, the solution
(3.3) is one dimensional,
 y 
(5.1)                           ˆ                ˆ
T j ( x, y, t )  G j ( y , t )  B j erfc          , j  1, 2 .
 2  jt 
        

Perfect Contact ( h   ) on x=0

If we take h   in the general solution (3.3a), then    and only the leading term of U B remains.
Thus, for perfect contact between regions, we have a simplification for small t

15
3 j 
2(1)
(5.2)          ˆ                 ˆ
T j ( x, y, t )  G j ( y, t ) 
                Imerfc  V
0
j       3 j        (  )    
  r
dr ,   j  1, 2

 y 
where                                                     ˆ
G j ( y, t )  B j erfc          , j  1, 2
 2  jt 
        
 j ( )             x           2   2  i y                                 1
j                                                                  , 2 
j

2
j
2 t                              2 t                                        j

kj 2   2                        B2           B   
Vj ( )                                                                       2 1 2
j
 2
k1 12   2  k2                      2   2  2      1   
2                2

  cos   i sin  , 0     / 2 .

The choices for  remain the same as that in (3.3a). For large t, a reduction of (3.3b) exists also.

Equal Diffusivities

We can derive the Laplace transform for Problem II by differentiation of (4.3) with respect to y and
make replacements according to (3.2),
B j e qy                    2hk3 j      
 e x       q2   2
sin( y)
ˆ
T j ( x, y, p)                  (1)   3 j
 k k                                              d   B2  B1  , j=1,2 .
p                          p          0               q    h(k1  k2 )   q 2   2 
2        2
   1 2

The inverse Laplace transform on the second term is the same as in Problem I and one gets
2( 1)3 j k3 j                                              
ˆ ( x, y, t )  B erfc  y                                      a             a     
t

 B2  B1   U  0, ,    U  b, ,     e  sin( y )d  d .
2
(5.3) T j                           
 2 t    (k1  k2 )
j
0    2            2  0
Now differentiation of (4.4) with respect to b gives
                                                                   b2
 b                    
e
 2a
(5.4)                                                             sin(b)d                                  e       4a

0
4 a3 / 2
and (5.3) becomes
c2

 y                    (1)3 j k3 j y                            e 
 B2  B1   U  0, ,    U  b, ,   3 / 2 d
t
ˆ                                                                        a              a
T j ( x, y, t )  B j erfc                                                              
 2  t  2  (k1  k2 )             0    2            2  
1   1 
a / 2  x /(2  ), b  h     , c  y /(2  ), U ( ,  , t )  e 2  t erfc( t   / t ) .
2

 k1 k2 
These integrals, denoted as J 21 , have also been studied in [2, Chapters 2 & 3, Folder 21] and we can
express the solution as
 y             (1)3 j k3 j y
 B2  B1   J 21 (0, , c, t )  J 21 (b, , c, t )  , j  1, 2
ˆ                                                                                 a                   a
(5.5)     Tj ( x, y, t )  B j erfc                                                                                   
 2 t        2  (k1  k2 )                         2                   2         
1      1 
a / 2  x /(2  ), b  h     , c  y /(2  )
 k1 k2 
We also note that

 a                  a                                          1
J 21  0, , c, t   2I 5  c, , T   2  e  c w erfc(aw / 2)dw , T 
2 2
(5.6)
 2                  2           T                               t

16
where I 5 is computed in [2, Chapters 2 & 3, Folder 5]. Setting k1  k2  k does not produce much of a
reduction. However, it is interesting to note that if B1  B2 , we have identical one-dimensional
distributions on each side of the boundary x=0 (symmetry) and no heat flows across this boundary. It is
also interesting to note that in this case where B1  B2 , k and h play no role in the solution.

Common Physical Properties in and Perfect Contact Between Regions 1 and 2

This is an unsteady state, upper half plane problem with a jump in boundary values on y=0. In this case
we take 1   2   and k1  k2  k in (5.5) and let h   . Then b   , J 21  b, a / 2, c, t   0 , and (5.5),
with the help of (5.6), goes over to
 y       (1)3 j y
 B2  B1  I 5  c, , T  , T 
ˆ                                                                   a            1
(5.7)            T j ( x, y, t )  B j erfc                                                      ,                j  1, 2
 2 t       2                      2              t
a/2     x /(2  ) , c  y /(2  )
The computational forms and properties of I 5 ( ,  , T ) are found in [2, Chapters 2 & 3, Folder 5]. Note
also that this case can be computed from (5.2) with h very large. Using

(5.8)                                  I 5 (c, 0, T )     erfc(cT )
2c
we find that the temperature on x=0,
 y  (1)3 j                   y 
(5.9)                   ˆ
T j (0, y, t )  B j erfc              B2  B1  erfc        ,                    j  1, 2 ,
 2 t   2                      2 t 
is continuous with the value
ˆ                ˆ             1                y 
T1 (0, y, t )  T2 (0, y, t )  [ B1  B2 ]erfc       .
2                2 t 
The steady state solution can be written by taking t  and using the result
1          a
I 5 (a, b, 0)         tan 1   .
a            b
The steady state becomes
(1)      3 j
 y   
(5.10)                         ˆ
T j ( x, y, )  B j                    B2  B1  tan -1 
    ,

j  1, 2
                            x    
which can be written in a more symmetric form
 1          y                           1       y 
(5.11)           ˆ
T j ( x, y, )  B j 1  tan 1                    B3 j
1
 tan    , j  1, 2 .
 x                                     x
 
                                       
         

Quarter Plane with Newtonian Cooling on x=0

If we take k2   in (5.5), the formulas for equal diffusivities go over to
ˆ                        y        y[ B2  B1 ]           a                     a          
(5.12)               T1 ( x, y, t )  B1 erfc                        J 21 (0, 2 , c, t )  J 21 (b, 2 , c, t ) 
 2 t         2                                                 
ˆ                        y 
T2 ( x, y, t )  B2 erfc        
 2 t 
a / 2  x /(2  ), b  h  / k1 , c  y /(2  ) .
As before we can make use of (5.6) and refer to [2, Chapters 2 & 3, Folders 21 and 5] for properties of
these functions. Notice that if B2  0 , we have Newtonian cooling of Region 1 along x=0.
17
ˆ
The steady state case with t  gives T2 ( x, y, )  B2 and the boundary condition for Region 1
becomes
ˆ
T1
(5.13)                                     k1                            ˆ
 h[ B2  T1 ]          .
x                            x 0
x 0
The steady state formula for Region 1 is reduced by applying [2, Chapters 2 & 3, Folder 21]

2e2     e2 x  ( 2  x2 ) / t                               2          
J 21 ( ,  ,  , t ) 
 
  2  x2 e                 dx and J 21 (0,  ,  , ) 
 
tan 1  

using the relation tan 1 ( z )  tan 1 (1/ z )   / 2 . Partial fractions on the term 1/( 2  x 2 ) gives
2
J 21 ( ,  ,  , )          Im[e z E1 ( z )], z  2 (   i )
 
and the steady state form of (5.12) becomes
 x  2B                y  2[ B2  B1 ]
             
2B
ˆ
(5.14) T1 ( x, y, )  1 tan 1    2 tan 1                              Im eh / k1 E1  h / k1  ,   x  iy
            y                  x            
The computation of e En ( z ) , n  1 is described in [6] and can be computed from the algorithm in [7].
z

It is apparent that each term is the imaginary part of an analytic function in the cut plane and therefore
ˆ
each satisfies the Laplace equation. On y=0 we have T1 ( x ,0,  )  B1 since  is real and the imaginary
part of a real number is zero. Writing Im{e z E1 ( z )}  [e z E1 ( z )  e z E1 ( z )] /(2i), z  h / k1, differentiating with
respect to x, and setting x=0 shows that (5.13) is satisfied.

Also note that for h  0 , the case of perfect insulation on x=0, the asymptotics

eh / k1 E1(h / k1 )
   ln(h  / k1 )  i tan 1( y / x),  h 0
ˆ
in (5.14) give the steady state temperature for the quadrant, T ( x, y, )  B1 . In this context  is the
Euler constant.

Quarter Plane with Temperature B2 along x=0

We consider (5.12) with infinite conductivity in Region 2 further. If h   then b   and
J 21  b, a / 2, c, t   0 giving perfect contact between Regions 1 and 2. This gives the formulas:

ˆ                          y   y[ B2  B1 ]  a 
(5.15)              T1 ( x, y, t )  B1 erfc                    I 5  c, , T  , T=1/ t , j  1, 2
2   t                 2 
ˆ                           y 
T2 ( x, y, t )  B2 erfc        ,           a / 2  x /(2  ),   c  y /(2  )
2   t 
and with B2  0 we have the first quadrant temperature for zero boundary condition along x=0.
ˆ               ˆ
Computation and properties of I 5 are found in [2, Chapters 2 & 3, Folder 5]. T1 (0, y, t )  T2 (0, y, t ) on
x=0 follows from (5.8)

Taking t  gives the steady state temperature for the first quadrant from
1                       ˆ             2B         x  2B        y
(5.16) I5 ( ,  ,0)      tan 1   . Then      T1 ( x, y, )  1 tan 1    2 tan 1  
                                               y          x
with a temperature B2 on x=0. This agrees with (5.14) for h   since e E1 ( z ) 1/ z for z   .
z

18
6. Computational Notes

In order to numerically evaluate the temperature formulae, we need to compute integrands for complex
values of the arguments. Reference [2, Chapter 4] provides the critical FORTRAN codes. Subroutine
ZERF is used to compute the scaled function exp( z 2 )erfc( z ) directly. Subroutine DQUAD8 is used to
compute quadratures. At its most basic level, DQUAD8 uses an adaptive form of the 8-point Gauss
formula to compute (real) quadratures to high accuracy. But at its most advanced level, DQUAD8 [8],[2,
Chapters 1 & 3] uses a stepping procedure which will estimate truncation error on infinite intervals and
terminate the integration when a specified truncation error is reached. These routines are all double
precision routines computing with approximately 16 digits and returning accuracies O(10-13) on most
machines. Here the “Z” routines are double precision complex routines with complex numbers z carried
as double precision ordered pairs Z = (ZR,ZI) and manipulated with lower subroutines ZABS, ZMLT,
ZDIV, ZSQRT, etc. Reference [7] for (5.14) also uses these basic routines to do complex arithmetic.
2
For the sequence yn = e z inerfc(z) [3],[4] in the series part of (2.18) and (3.3), we coded the continued
fraction for the ratios rn  yn / yn1 (Re(z) > 0)
1/ 2
(6.1)                                  rN  0, rnN1 
N
                , n  N , N  1,...,1
z  nrnN
and normalize the ynN sequence with y1 :
2
yn  y1r0N r1N ...rnN  e z i n erfc( z ) , n  1, 2,..., N -1 ,   y1 
2
N
(6.2)

(Forward recurrence on the 3-term recurrence relation is not stable.) While the series in (2.18) and (3.3),
N                               N
1
2  t  ( 2  t )n2 yn , 2 t  ( 2  t )n1 yn
N                         N
,     2 t  ,      N  N
n2                             n1                                  4
require less than N  25 terms for errors O(425 )  O(1015 ) , N ( 120) needs to be considerably larger
in order for ynN , n<25, to achieve full accuracy over the right half z plane. Actually, as N increases, the
number of members ynN achieving full precision increases. But all members are not required to have full
precision because the factor (2  t ) n decreases rapidly in magnitude and shifts less accurate digits in
the product (2  t ) n ynN off of the floating point register when added to the accumulated sum. The
convergence of the continued fraction is slowest for z near the imaginary axis, but the convergence of
each series (see (6.5)) is helped by the decrease in magnitude of the approximation ynN to y n [3], [4],
e z 2 n
yn  e z i n erfc( z )                     , n  , z bounded .
2
(6.3)
2n (1  n / 2)
2
Since ez ierfc( z) is used in (2.18) and (4. 2), we developed a special subroutine called ZIERFC which
implemented the formulae
               e z 
2

ierfc( z )          1  z  e z erfc( z )           z 7
2


                                          
Re( z )  0     
(6.4)                                                   1 e z 
2
                                    (3 / 2) m
2 
                             (1)m                     z 7
              2  z m 0              z 2m
Re ( z )  0      ierfc( z )  2( z )  ierfc( z )

where (–z) is in the right half plane. The formula for Re z <0 is the analytic continuation into the left half
z plane with the right side of the equation computed from the formulas for Re( z)  0 . e z erfc( z ) can be
2

computed from subroutine ZERF [2, Chapter 4]. All sums are terminated on a relative error test,

19
(6.5)                                   term  accumulated sum 

where  is usually on the order of 0.5 x1014 .
2
Subroutine ZERF computes e z erfc( z ) by applying several formulas over regions in the right-half z plane
where rapid computation is possible. The computation uses the power series for small z , the continued
fraction for intermediate values of z , the asymptotic expansion for large values of z , and [1, 7.1.29] to
cover a region near the y-axis where the continued fraction is not efficient. An analytic continuation
formula is used to cover the left half plane.

7. Numerical Comparisons of (2.18) and (3.3) Against Special Cases

To test formulas (2.18) and (3.3) for the solutions of Problems I and II, we compared some numerical
evaluations over a significant range of variables against the special cases in Sections 2 and 3 where the
solutions are represented by different formulae. The results were in agreement with what we considered
to be the working accuracy with absolute (decimal place) errors O (10 13 ) or better. Some examples are
shown in TABLES 1-9 with TABLES 1-4 for Problem I and TABLES 5-9 for Problem II. TABLES
10 and 11 show full 10 digit solutions associated with TABLES 1 and 5.

In TABLES 1-8, we chose to display values of interest in columns for three values of y along a line
parallel to the x axis so that the discontinuity at x=0 can be seen when Regions 1 and 2 interact. Column
1 shows the x values in Regions 1 and 2 (values of x+0 and x-0 for the boundary x=0 are taken to be
0.10D  5 ). Columns 2, 4, and 6 show the values of the temperature at the point (x,y) as computed by
(2.18) or (3.3). Columns 3, 5,and 7 show the relative errors
(2.16)-Sc          (3.3)-Sc
or
Sc                Sc
between the formulas being compared with the special case value Sc at the point (x,y) taken as the
correct value; the absolute error (numerator) is displayed when the special case value Sc is zero. In
TABLE 9, we show how (3.3), as t gets large, approaches the steady state value Sc computed from
(5.14). Except for the exponential integral in (5.14), which is computed from [7], the special functions
needed for the value of Sc were computed from the FORTRAN codes in [2, Chapter 4].

Since we expect the results to be independent of  , numerical experiments comparing the results for
different values of  (    / 4 or    / 4 ) confirmed this expectation with relative errors O (10 14 )
consistently over wide ranges of parameters.

In some of the results, it looks like the accuracy decreases. For example, the entries in TABLE 3 for
t=1.0, x= 0.00D+0, and y= 2.0 are 0.92D-02 and 0.26D-12, with the relative error being O 1012  . In
this case the temperature is small, O 102  , and the relative error has dropped below our expected
precision of O 1013  . What we actually have is an absolute (decimal place) accuracy of O 1014  , two
leading zeros followed by approximately 12 significant digits.

20
PROBLEM I: COMPARISON OF (2.18) WITH SPECIAL CASES
TABLE 1:       COMPARISON OF (2.18) WITH (4.5)
F1 =1.00 F2 =2.00 α1 =0.50 α2 =0.50 k1 =0.75 k2 =0.50 h=1.30
t =   1.0                                                              t = 2.0
x          y =   0.5             y =   1.0              y =   2.0            x        y = 0.5           y = 1.0           y =                2.0
-0.10D+01   0.15D+01   0.23D-14   0.63D+00   0.58D-14    0.63D-01   0.34D-13 -0.10D+01 0.26D+01 0.14D-14 0.14D+01 0.25D-14 0.35D+00                0.16D-13
-0.50D+00   0.14D+01   0.21D-14   0.57D+00   0.58D-14    0.57D-01   0.72D-13 -0.50D+00 0.23D+01 0.76D-15 0.13D+01 0.21D-14 0.32D+00                0.16D-13
-0.10D-05   0.11D+01   0.10D-14   0.46D+00   0.39D-14    0.46D-01   0.77D-13 -0.10D-05 0.19D+01 0.46D-15 0.11D+01 0.10D-14 0.27D+00                0.14D-13
0.10D-05   0.84D+00   0.79D-15   0.36D+00   0.31D-14    0.37D-01   0.64D-13  0.10D-05 0.15D+01 0.59D-15 0.88D+00 0.18D-14 0.22D+00                0.12D-13
0.50D+00   0.65D+00   0.29D-14   0.29D+00   0.79D-14    0.30D-01   0.90D-13  0.50D+00 0.12D+01 0.14D-14 0.73D+00 0.36D-14 0.19D+00                0.18D-13
0.10D+01   0.57D+00   0.43D-14   0.25D+00   0.98D-14    0.26D-01   0.56D-13  0.10D+01 0.11D+01 0.29D-14 0.63D+00 0.49D-14 0.17D+00                0.24D-13
t =   3.0                                                              t = 4.0
x          y =   0.5             y =   1.0              y =   2.0            x        y = 0.5           y = 1.0           y =                2.0
-0.10D+01   0.34D+01   0.79D-15   0.21D+01   0.15D-14    0.72D+00   0.66D-14 -0.10D+01 0.40D+01 0.11D-14 0.27D+01 0.18D-14 0.11D+01                0.56D-14
-0.50D+00   0.31D+01   0.29D-15   0.19D+01   0.12D-14    0.65D+00   0.72D-14 -0.50D+00 0.36D+01 0.73D-15 0.24D+01 0.15D-14 0.10D+01                0.49D-14
-0.10D-05   0.25D+01   0.35D-15   0.16D+01   0.14D-15    0.56D+00   0.44D-14 -0.10D-05 0.31D+01 0.29D-15 0.21D+01 0.00D+00 0.86D+00                0.31D-14
0.10D-05   0.21D+01   0.43D-15   0.14D+01   0.99D-15    0.48D+00   0.55D-14  0.10D-05 0.25D+01 0.35D-15 0.18D+01 0.25D-15 0.76D+00                0.25D-14
0.50D+00   0.17D+01   0.13D-14   0.12D+01   0.25D-14    0.42D+00   0.97D-14  0.50D+00 0.21D+01 0.83D-15 0.15D+01 0.14D-14 0.67D+00                0.50D-14
0.10D+01   0.15D+01   0.24D-14   0.10D+01   0.37D-14    0.37D+00   0.11D-13  0.10D+01 0.19D+01 0.15D-14 0.14D+01 0.23D-14 0.60D+00                0.71D-14

TABLE 2:       COMPARISON OF (2.18) WITH (4.7)
F1 =1.00 F2 =2.00 α1 =0.50 α2 =0.50 k1 =0.75 k2 =0.75 h=1.0D+40

t =   1.0                                                              t = 2.0
x          y =   0.5             y =   1.0              y =   2.0            x        y = 0.5           y = 1.0           y =                2.0
-0.10D+01   0.10D+01   0.11D-14   0.42D+00   0.30D-14    0.42D-01   0.19D-13 -0.10D+01 0.17D+01 0.51D-15 0.98D+00 0.79D-15 0.24D+00                0.88D-14
-0.50D+00   0.94D+00   0.83D-15   0.39D+00   0.30D-14    0.39D-01   0.41D-13 -0.50D+00 0.16D+01 0.00D+00 0.90D+00 0.61D-15 0.22D+00                0.81D-14
-0.10D-05   0.79D+00   0.14D-15   0.33D+00   0.15D-14    0.34D-01   0.38D-13 -0.10D-05 0.14D+01 0.32D-15 0.80D+00 0.28D-15 0.20D+00                0.62D-14
0.10D-05   0.79D+00   0.28D-15   0.33D+00   0.15D-14    0.34D-01   0.38D-13  0.10D-05 0.14D+01 0.32D-15 0.80D+00 0.83D-15 0.20D+00                0.72D-14
0.50D+00   0.64D+00   0.14D-14   0.28D+00   0.44D-14    0.29D-01   0.54D-13  0.50D+00 0.12D+01 0.11D-14 0.69D+00 0.22D-14 0.18D+00                0.11D-13
0.10D+01   0.57D+00   0.22D-14   0.24D+00   0.52D-14    0.26D-01   0.32D-13  0.10D+01 0.11D+01 0.17D-14 0.62D+00 0.29D-14 0.16D+00                0.14D-13
t =   3.0                                                              t = 4.0
x          y =   0.5             y =   1.0              y =   2.0            x        y = 0.5           y = 1.0           y =                2.0
-0.10D+01   0.23D+01   0.19D-15   0.15D+01   0.46D-15    0.50D+00   0.31D-14 -0.10D+01 0.28D+01 0.48D-15 0.19D+01 0.83D-15 0.78D+00                0.30D-14
-0.50D+00   0.21D+01   0.21D-15   0.14D+01   0.00D+00    0.47D+00   0.30D-14 -0.50D+00 0.26D+01 0.35D-15 0.17D+01 0.64D-15 0.72D+00                0.25D-14
-0.10D-05   0.19D+01   0.71D-15   0.12D+01   0.73D-15    0.43D+00   0.14D-14 -0.10D-05 0.23D+01 0.00D+00 0.16D+01 0.14D-15 0.67D+00                0.13D-14
0.10D-05   0.19D+01   0.35D-15   0.12D+01   0.55D-15    0.43D+00   0.38D-14  0.10D-05 0.23D+01 0.00D+00 0.16D+01 0.00D+00 0.67D+00                0.12D-14
0.50D+00   0.16D+01   0.82D-15   0.11D+01   0.17D-14    0.39D+00   0.62D-14  0.50D+00 0.20D+01 0.44D-15 0.14D+01 0.78D-15 0.61D+00                0.31D-14
0.10D+01   0.15D+01   0.14D-14   0.97D+00   0.24D-14    0.35D+00   0.73D-14  0.10D+01 0.18D+01 0.98D-15 0.13D+01 0.14D-14 0.56D+00                0.44D-14

TABLE 3:    COMPARISON OF (2.18) WITH (4.8)
F1 =1.00 F2 =0.00 α1 =0.50 α2 =0.50 k1 =0.75 k2 =1.0D+40 h=1.30
t =    1.0                                                                      t =   2.0
x          y =   0.5             y =    1.0             y =   2.0               x          y =   0.5             y =   1.0             y =   2.0
0.00D+00   0.26D+00   0.14D-13   0.99D-01    0.52D-14   0.92D-02   0.26D-12     0.00D+00   0.39D+00   0.34D-14   0.20D+00   0.12D-13   0.45D-01   0.83D-13
0.50D+00   0.42D+00   0.63D-14   0.17D+00    0.22D-13   0.16D-01   0.24D-12     0.50D+00   0.66D+00   0.30D-14   0.35D+00   0.97D-14   0.80D-01   0.61D-13
0.10D+01   0.49D+00   0.75D-14   0.20D+00    0.18D-13   0.20D-01   0.10D-12     0.10D+01   0.80D+00   0.50D-14   0.44D+00   0.91D-14   0.10D+00   0.51D-13
0.15D+01   0.52D+00   0.62D-14   0.22D+00    0.13D-13   0.22D-01   0.50D-13     0.15D+01   0.88D+00   0.54D-14   0.49D+00   0.86D-14   0.12D+00   0.30D-13
0.20D+01   0.53D+00   0.23D-14   0.22D+00    0.40D-14   0.22D-01   0.21D-13     0.20D+01   0.91D+00   0.53D-14   0.52D+00   0.73D-14   0.13D+00   0.28D-13
0.30D+01   0.53D+00   0.00D+00   0.22D+00    0.00D+00   0.23D-01   0.00D+00     0.30D+01   0.93D+00   0.36D-15   0.53D+00   0.63D-15   0.13D+00   0.42D-15
t =    3.0                                                                      t =   4.0
x          y =   0.5             y =    1.0             y =   2.0               x          y =   0.5             y =   1.0             y =   2.0
0.00D+00   0.47D+00   0.19D-14   0.27D+00    0.73D-14   0.84D-01   0.40D-13     0.00D+00   0.54D+00   0.29D-14   0.33D+00   0.68D-14   0.12D+00   0.28D-13
0.50D+00   0.81D+00   0.22D-14   0.48D+00    0.68D-14   0.15D+00   0.34D-13     0.50D+00   0.93D+00   0.25D-14   0.58D+00   0.69D-14   0.22D+00   0.23D-13
0.10D+01   0.10D+01   0.35D-14   0.62D+00    0.68D-14   0.20D+00   0.26D-13     0.10D+01   0.12D+01   0.42D-14   0.76D+00   0.65D-14   0.30D+00   0.20D-13
0.15D+01   0.11D+01   0.43D-14   0.71D+00    0.67D-14   0.24D+00   0.28D-13     0.15D+01   0.13D+01   0.42D-14   0.89D+00   0.61D-14   0.35D+00   0.22D-13
0.20D+01   0.12D+01   0.43D-14   0.76D+00    0.60D-14   0.26D+00   0.19D-13     0.20D+01   0.14D+01   0.42D-14   0.96D+00   0.65D-14   0.39D+00   0.16D-13
0.30D+01   0.12D+01   0.30D-14   0.80D+00    0.43D-14   0.28D+00   0.10D-13     0.30D+01   0.15D+01   0.40D-14   0.10D+01   0.54D-14   0.43D+00   0.11D-13

TABLE 4:       COMPARISON OF (2.18) WITH (4.9)
F1 =1.00 F2 =0.00 α1 =0.50 α2 =0.50 k1 =0.75 k2 =1.0D+40 h=1.0D+40
t =    1.0                                                                      t =   2.0
x          y =   0.5             y =    1.0             y =   2.0               x        y =     0.5             y =   1.0             y =   2.0
0.00D+00   0.00D+00   0.00D+00   0.86D-15    0.86D-15   0.24D-14   0.24D-14     0.00D+00-0.33D-15     0.33D-15   0.33D-15   0.33D-15   0.25D-14   0.25D-14
0.50D+00   0.30D+00   0.58D-14   0.11D+00    0.22D-13   0.98D-02   0.32D-12     0.50D+00 0.43D+00     0.29D-14   0.21D+00   0.10D-13   0.45D-01   0.85D-13
0.10D+01   0.45D+00   0.57D-14   0.18D+00    0.14D-13   0.17D-01   0.97D-13     0.10D+01 0.69D+00     0.37D-14   0.36D+00   0.75D-14   0.82D-01   0.54D-13
0.15D+01   0.50D+00   0.46D-14   0.21D+00    0.99D-14   0.21D-01   0.42D-13     0.15D+01 0.82D+00     0.43D-14   0.45D+00   0.73D-14   0.11D+00   0.28D-13
0.20D+01   0.52D+00   0.19D-14   0.22D+00    0.30D-14   0.22D-01   0.18D-13     0.20D+01 0.88D+00     0.44D-14   0.50D+00   0.60D-14   0.12D+00   0.24D-13
0.30D+01   0.53D+00   0.00D+00   0.22D+00    0.00D+00   0.23D-01   0.00D+00     0.30D+01 0.92D+00     0.24D-15   0.53D+00   0.42D-15   0.13D+00   0.21D-15
t =    3.0                                                                      t =   4.0
x        y =     0.5             y =    1.0             y =   2.0               x          y =   0.5             y =   1.0             y =   2.0
0.00D+00-0.67D-15     0.67D-15   0.00D+00    0.00D+00   0.21D-14   0.21D-14     0.00D+00   0.00D+00   0.00D+00   0.00D+00   0.00D+00   0.17D-14   0.17D-14
0.50D+00 0.51D+00     0.18D-14   0.28D+00    0.70D-14   0.82D-01   0.47D-13     0.50D+00   0.57D+00   0.24D-14   0.33D+00   0.76D-14   0.12D+00   0.31D-13
0.10D+01 0.84D+00     0.30D-14   0.49D+00    0.60D-14   0.15D+00   0.27D-13     0.10D+01   0.95D+00   0.35D-14   0.59D+00   0.58D-14   0.22D+00   0.22D-13
0.15D+01 0.10D+01     0.37D-14   0.63D+00    0.56D-14   0.21D+00   0.28D-13     0.15D+01   0.12D+01   0.34D-14   0.78D+00   0.56D-14   0.30D+00   0.22D-13
0.20D+01 0.11D+01     0.33D-14   0.72D+00    0.53D-14   0.24D+00   0.17D-13     0.20D+01   0.13D+01   0.37D-14   0.89D+00   0.57D-14   0.36D+00   0.14D-13
0.30D+01 0.12D+01     0.25D-14   0.79D+00    0.35D-14   0.27D+00   0.83D-14     0.30D+01   0.15D+01   0.35D-14   0.10D+01   0.46D-14   0.42D+00   0.99D-14

21
PROBLEM II: COMPARISON OF (3.3) WITH SPECIAL CASES
TABLE 5: COMPARISON OF (3.3) WITH (5.5)
B1 =1.00 B2 =2.00 α1 =0.50 α2 =0.50 k1 =0.75 k2 =0.50 h=1.30
t =   1.0                                                              t = 2.0
x          y =   0.5             y =   1.0              y =   2.0            x        y = 0.5           y = 1.0           y =                2.0
-0.10D+01   0.12D+01   0.18D-15   0.62D+00   0.36D-15    0.87D-01   0.22D-14 -0.10D+01 0.14D+01 0.32D-15 0.91D+00 0.49D-15 0.29D+00                0.31D-14
-0.50D+00   0.12D+01   0.19D-15   0.58D+00   0.96D-15    0.80D-01   0.71D-14 -0.50D+00 0.13D+01 0.17D-15 0.85D+00 0.79D-15 0.27D+00                0.29D-14
-0.10D-05   0.99D+00   0.11D-15   0.49D+00   0.90D-15    0.69D-01   0.82D-14 -0.10D-05 0.11D+01 0.00D+00 0.73D+00 0.76D-15 0.24D+00                0.33D-14
0.10D-05   0.78D+00   0.14D-15   0.41D+00   0.81D-15    0.60D-01   0.64D-14 0.10D-05 0.92D+00 0.00D+00 0.63D+00 0.18D-15 0.21D+00                 0.12D-14
0.50D+00   0.66D+00   0.33D-15   0.36D+00   0.11D-14    0.53D-01   0.71D-14 0.50D+00 0.79D+00 0.14D-15 0.55D+00 0.40D-15 0.19D+00                 0.15D-14
0.10D+01   0.63D+00   0.35D-15   0.33D+00   0.50D-15    0.48D-01   0.27D-14 0.10D+01 0.75D+00 0.30D-15 0.51D+00 0.22D-15 0.17D+00                 0.17D-14
t =   3.0                                                              t = 4.0
x          y =   0.5             y =   1.0              y =   2.0            x        y = 0.5           y = 1.0           y =                2.0
-0.10D+01   0.15D+01   0.15D-15   0.11D+01   0.42D-15    0.45D+00   0.17D-14 -0.10D+01 0.15D+01 0.14D-15 0.11D+01 0.19D-15 0.56D+00                0.98D-15
-0.50D+00   0.14D+01   0.00D+00   0.98D+00   0.68D-15    0.41D+00   0.22D-14 -0.50D+00 0.15D+01 0.00D+00 0.11D+01 0.42D-15 0.52D+00                0.13D-14
-0.10D-05   0.12D+01   0.18D-15   0.86D+00   0.52D-15    0.37D+00   0.21D-14 -0.10D-05 0.13D+01 0.18D-15 0.94D+00 0.00D+00 0.47D+00                0.83D-15
0.10D-05   0.99D+00   0.11D-15   0.74D+00   0.00D+00    0.33D+00   0.83D-15 0.10D-05 0.10D+01 0.00D+00 0.82D+00 0.00D+00 0.43D+00                 0.65D-15
0.50D+00   0.86D+00   0.00D+00   0.66D+00   0.17D-15    0.30D+00   0.13D-14 0.50D+00 0.89D+00 0.00D+00 0.73D+00 0.30D-15 0.39D+00                 0.11D-14
0.10D+01   0.81D+00   0.14D-15   0.61D+00   0.00D+00    0.28D+00   0.79D-15 0.10D+01 0.84D+00 0.13D-15 0.68D+00 0.33D-15 0.36D+00                 0.11D-14

TABLE 6: COMPARISON OF (3.3) WITH (5.7)
B1 =1.00 B2 =2.00 α1 =0.50 α2 =0.50 k1 =0.75 k2 =0.75 h=1.0D+40
t =   1.0                                                              t = 2.0
x          y =   0.5             y =   1.0              y =   2.0            x        y = 0.5           y = 1.0           y =                2.0
-0.10D+01   0.12D+01   0.18D-15   0.61D+00   0.36D-15    0.86D-01   0.18D-14 -0.10D+01 0.14D+01 0.16D-15 0.90D+00 0.37D-15 0.29D+00                0.25D-14
-0.50D+00   0.11D+01   0.19D-15   0.56D+00   0.99D-15    0.78D-01   0.60D-14 -0.50D+00 0.13D+01 0.00D+00 0.83D+00 0.67D-15 0.26D+00                0.27D-14
-0.10D-05   0.93D+00   0.00D+00   0.48D+00   0.70D-15    0.68D-01   0.71D-14 -0.10D-05 0.11D+01 0.20D-15 0.72D+00 0.46D-15 0.24D+00                0.29D-14
0.10D-05   0.93D+00   0.00D+00   0.48D+00   0.70D-15    0.68D-01   0.71D-14 0.10D-05 0.11D+01 0.20D-15 0.72D+00 0.00D+00 0.24D+00                 0.15D-14
0.50D+00   0.71D+00   0.31D-15   0.39D+00   0.14D-14    0.58D-01   0.82D-14 0.50D+00 0.85D+00 0.13D-15 0.61D+00 0.37D-15 0.21D+00                 0.17D-14
0.10D+01   0.64D+00   0.35D-15   0.34D+00   0.65D-15    0.51D-01   0.31D-14 0.10D+01 0.77D+00 0.29D-15 0.54D+00 0.00D+00 0.19D+00                 0.22D-14
t =   3.0                                                              t = 4.0
x          y =   0.5             y =   1.0              y =   2.0            x        y = 0.5           y = 1.0           y =                2.0
-0.10D+01   0.15D+01   0.15D-15   0.10D+01   0.42D-15    0.44D+00   0.15D-14 -0.10D+01 0.15D+01 0.14D-15 0.11D+01 0.19D-15 0.56D+00                0.99D-15
-0.50D+00   0.14D+01   0.16D-15   0.97D+00   0.69D-15    0.41D+00   0.20D-14 -0.50D+00 0.14D+01 0.00D+00 0.11D+01 0.42D-15 0.52D+00                0.85D-15
-0.10D-05   0.12D+01   0.00D+00   0.85D+00   0.26D-15    0.37D+00   0.18D-14 -0.10D-05 0.12D+01 0.18D-15 0.93D+00 0.00D+00 0.48D+00                0.58D-15
0.10D-05   0.12D+01   0.00D+00   0.85D+00   0.00D+00    0.37D+00   0.89D-15 0.10D-05 0.12D+01 0.18D-15 0.93D+00 0.24D-15 0.48D+00                 0.58D-15
0.50D+00   0.92D+00   0.12D-15   0.72D+00   0.46D-15    0.33D+00   0.15D-14 0.50D+00 0.96D+00 0.00D+00 0.80D+00 0.42D-15 0.43D+00                 0.10D-14
0.10D+01   0.83D+00   0.27D-15   0.64D+00   0.17D-15    0.30D+00   0.11D-14 0.10D+01 0.87D+00 0.13D-15 0.71D+00 0.31D-15 0.39D+00                 0.14D-14

TABLE 7:       COMPARISON OF (3.3) WITH (5.12)
B1 =1.00 B2 =0.00 α1 =0.50 α2 =0.50 k1 =0.75 k2 =1.0D+40 h=1.30
t =    1.0                                                                      t =   2.0
x          y =   0.5             y =    1.0             y =   2.0               x          y =   0.5             y =   1.0             y =   2.0
0.00D+00   0.35D+00   0.51D-14   0.15D+00    0.55D-15   0.19D-01   0.44D-13     0.00D+00   0.39D+00   0.69D-14   0.20D+00   0.51D-14   0.55D-01   0.19D-13
0.50D+00   0.54D+00   0.23D-14   0.25D+00    0.42D-14   0.33D-01   0.29D-13     0.50D+00   0.60D+00   0.35D-14   0.35D+00   0.27D-14   0.97D-01   0.11D-13
0.10D+01   0.60D+00   0.93D-15   0.30D+00    0.15D-14   0.41D-01   0.79D-14     0.10D+01   0.68D+00   0.65D-15   0.42D+00   0.12D-14   0.13D+00   0.95D-14
0.15D+01   0.61D+00   0.54D-15   0.31D+00    0.71D-15   0.44D-01   0.28D-14     0.15D+01   0.71D+00   0.63D-15   0.46D+00   0.12D-14   0.14D+00   0.41D-14
0.20D+01   0.62D+00   0.18D-15   0.32D+00    0.18D-15   0.45D-01   0.15D-15     0.20D+01   0.72D+00   0.62D-15   0.47D+00   0.94D-15   0.15D+00   0.24D-14
0.30D+01   0.62D+00   0.00D+00   0.32D+00    0.00D+00   0.45D-01   0.00D+00     0.30D+01   0.72D+00   0.00D+00   0.48D+00   0.23D-15   0.16D+00   0.71D-15
t =    3.0                                                                      t =   4.0
x          y =   0.5             y =    1.0             y =   2.0               x          y =   0.5             y =   1.0             y =   2.0
0.00D+00   0.40D+00   0.53D-14   0.23D+00    0.39D-14   0.79D-01   0.14D-13     0.00D+00   0.41D+00   0.48D-14   0.24D+00   0.23D-14   0.95D-01   0.85D-14
0.50D+00   0.63D+00   0.41D-14   0.39D+00    0.24D-14   0.14D+00   0.94D-14     0.50D+00   0.64D+00   0.40D-14   0.41D+00   0.18D-14   0.17D+00   0.60D-14
0.10D+01   0.71D+00   0.62D-15   0.48D+00    0.12D-14   0.19D+00   0.55D-14     0.10D+01   0.73D+00   0.61D-15   0.51D+00   0.11D-14   0.23D+00   0.47D-14
0.15D+01   0.75D+00   0.74D-15   0.52D+00    0.13D-14   0.21D+00   0.61D-14     0.15D+01   0.77D+00   0.72D-15   0.56D+00   0.99D-15   0.26D+00   0.55D-14
0.20D+01   0.76D+00   0.73D-15   0.55D+00    0.10D-14   0.23D+00   0.20D-14     0.20D+01   0.79D+00   0.71D-15   0.59D+00   0.94D-15   0.29D+00   0.17D-14
0.30D+01   0.77D+00   0.43D-15   0.56D+00    0.59D-15   0.24D+00   0.10D-14     0.30D+01   0.80D+00   0.56D-15   0.61D+00   0.55D-15   0.31D+00   0.90D-15

TABLE 8:       COMPARISON OF (3.3) WITH (5.15)
B1 =1.00 B2 =0.00 α1 =0.50 α2 =0.50 k1 =0.75 k2 =1.0D+40 h=1.0D+40
t =      1.0                                                                    t =     2.0
x        y =     0.5           y =      1.0           y =     2.0               x          y =   0.5           y =     1.0           y =     2.0
0.00D+00-0.11D-15     0.11D-15-0.72D-15      0.72D-15-0.11D-14     0.11D-14     0.00D+00   0.00D+00   0.00D+00-0.50D-15     0.50D-15-0.10D-14     0.10D-14
0.50D+00 0.43D+00     0.91D-15 0.17D+00      0.60D-14 0.20D-01     0.47D-13     0.50D+00   0.46D+00   0.36D-15 0.23D+00     0.39D-14 0.56D-01     0.20D-13
0.10D+01 0.57D+00     0.98D-15 0.27D+00      0.17D-14 0.35D-01     0.93D-14     0.10D+01   0.63D+00   0.88D-15 0.36D+00     0.15D-14 0.10D+00     0.12D-13
0.15D+01 0.60D+00     0.55D-15 0.30D+00      0.92D-15 0.42D-01     0.28D-14     0.15D+01   0.69D+00   0.64D-15 0.43D+00     0.13D-14 0.13D+00     0.47D-14
0.20D+01 0.61D+00     0.00D+00 0.31D+00      0.18D-15 0.44D-01     0.16D-15     0.20D+01   0.71D+00   0.47D-15 0.46D+00     0.96D-15 0.14D+00     0.23D-14
0.30D+01 0.62D+00     0.00D+00 0.32D+00      0.00D+00 0.45D-01     0.00D+00     0.30D+01   0.72D+00   0.00D+00 0.48D+00     0.23D-15 0.16D+00     0.71D-15
t =      3.0                                                                    t =     4.0
x          y =   0.5           y =      1.0           y =     2.0               x          y =   0.5           y =     1.0           y =     2.0
0.00D+00   0.22D-15   0.22D-15-0.33D-15      0.33D-15-0.89D-15     0.89D-15     0.00D+00   0.33D-15   0.33D-15-0.11D-15     0.11D-15-0.72D-15     0.72D-15
0.50D+00   0.47D+00   0.23D-15 0.25D+00      0.32D-14 0.79D-01     0.16D-13     0.50D+00   0.48D+00   0.00D+00 0.26D+00     0.26D-14 0.94D-01     0.10D-13
0.10D+01   0.65D+00   0.68D-15 0.40D+00      0.12D-14 0.14D+00     0.69D-14     0.10D+01   0.67D+00   0.67D-15 0.43D+00     0.91D-15 0.17D+00     0.60D-14
0.15D+01   0.72D+00   0.77D-15 0.49D+00      0.13D-14 0.19D+00     0.67D-14     0.15D+01   0.74D+00   0.75D-15 0.52D+00     0.86D-15 0.23D+00     0.65D-14
0.20D+01   0.75D+00   0.59D-15 0.53D+00      0.11D-14 0.22D+00     0.22D-14     0.20D+01   0.77D+00   0.86D-15 0.57D+00     0.98D-15 0.27D+00     0.17D-14
0.30D+01   0.77D+00   0.29D-15 0.56D+00      0.60D-15 0.24D+00     0.92D-15     0.30D+01   0.80D+00   0.56D-15 0.60D+00     0.55D-15 0.30D+00     0.92D-15

22
TABLE 9:       COMPARISON OF (3.3) WITH (5.14)
B1 =1.00 B2 =0.00 α1 =0.25 α2 =0.25 k1 =0.50 k2 =1.0D+40 h=1.3
y =    0.5                                                                      y =   1.0
t          x =   0.0             x =    0.5             x =   1.0               t          x =   0.0             x =   0.5             x =   1.0
0.10D+01   0.24D+00   0.29D+00   0.43D+00    0.33D+00   0.47D+00   0.39D+00     0.10D+01   0.67D-01   0.68D+00   0.13D+00   0.70D+00   0.15D+00   0.74D+00
0.10D+02   0.33D+00   0.35D-01   0.62D+00    0.42D-01   0.73D+00   0.54D-01     0.10D+02   0.19D+00   0.11D+00   0.39D+00   0.12D+00   0.51D+00   0.14D+00
0.10D+03   0.34D+00   0.36D-02   0.64D+00    0.44D-02   0.76D+00   0.57D-02     0.10D+03   0.21D+00   0.12D-01   0.43D+00   0.13D-01   0.58D+00   0.15D-01
0.10D+04   0.34D+00   0.36D-03   0.64D+00    0.44D-03   0.77D+00   0.57D-03     0.10D+04   0.21D+00   0.12D-02   0.44D+00   0.13D-02   0.59D+00   0.15D-02
0.10D+05   0.34D+00   0.36D-04   0.64D+00    0.44D-04   0.77D+00   0.57D-04     0.10D+05   0.21D+00   0.12D-03   0.44D+00   0.13D-03   0.59D+00   0.15D-03
0.10D+11   0.34D+00   0.36D-10   0.64D+00    0.44D-10   0.77D+00   0.57D-10     0.10D+11   0.21D+00   0.12D-09   0.44D+00   0.13D-09   0.59D+00   0.15D-09
y =    2.0                                                                      y =   3.0
t          x =   0.0             x =    0.5             x =   1.0               t          x =   0.0             x =   0.5             x =   1.0
0.10D+01   0.18D-02   0.98D+00   0.36D-02    0.99D+00   0.44D-02   0.99D+00     0.10D+01   0.82D-05   0.10D+01   0.17D-04   0.10D+01   0.21D-04   0.10D+01
0.10D+02   0.76D-01   0.34D+00   0.17D+00    0.35D+00   0.24D+00   0.37D+00     0.10D+02   0.31D-01   0.60D+00   0.70D-01   0.61D+00   0.10D+00   0.62D+00
0.10D+03   0.11D+00   0.41D-01   0.25D+00    0.43D-01   0.36D+00   0.45D-01     0.10D+03   0.72D-01   0.88D-01   0.16D+00   0.90D-01   0.25D+00   0.93D-01
0.10D+04   0.12D+00   0.42D-02   0.26D+00    0.44D-02   0.38D+00   0.46D-02     0.10D+04   0.79D-01   0.92D-02   0.18D+00   0.94D-02   0.27D+00   0.97D-02
0.10D+05   0.12D+00   0.42D-03   0.26D+00    0.44D-03   0.38D+00   0.47D-03     0.10D+05   0.79D-01   0.93D-03   0.18D+00   0.94D-03   0.27D+00   0.97D-03
0.10D+11   0.12D+00   0.42D-09   0.26D+00    0.44D-09   0.38D+00   0.47D-09     0.10D+11   0.79D-01   0.93D-09   0.18D+00   0.94D-09   0.27D+00   0.97D-09

10 DIGIT PRESENTATIONS
OF THE SOLUTIONS (2.18) and (3.3)
TABLE 10:     EVALUATION OF (2.18) WITH THE PARAMETERS OF TABLE 1
F1 =1.00 F2 =2.00 α1 =0.50 α2 =0.50 k1 =0.75 k2 =0.50 h=1.30
t = 1.0                                                       t = 2.0
x         y = 0.5          y = 1.0          y = 2.0           x         y = 0.5          y = 1.0          y = 2.0
-0.100D+01 0.1519883451D+01 0.6311959660D+00 0.6316310428D-01 -0.100D+01 0.2577708151D+01 0.1443631170D+01 0.3536145164D+00
-0.500D+00 0.1393924227D+01 0.5703391290D+00 0.5652226754D-01 -0.500D+00 0.2337772550D+01 0.1295296929D+01 0.3157839640D+00
-0.100D-05 0.1115372864D+01 0.4601177026D+00 0.4615250782D-01 -0.100D-05 0.1913464028D+01 0.1073483363D+01 0.2659978191D+00
0.100D-05 0.8387905502D+00 0.3597786297D+00 0.3715728239D-01  0.100D-05 0.1517066613D+01 0.8814742515D+00 0.2247044539D+00
0.500D+00 0.6530896414D+00 0.2862976787D+00 0.3024410924D-01  0.500D+00 0.1234194264D+01 0.7335985411D+00 0.1915136906D+00
0.100D+01 0.5691168253D+00 0.2457264541D+00 0.2581688475D-01  0.100D+01 0.1074237197D+01 0.6347090471D+00 0.1662933224D+00
t = 3.0                                                       t = 4.0
x         y = 0.5          y = 1.0          y = 2.0           x         y = 0.5          y = 1.0          y = 2.0
-0.100D+01 0.3372716514D+01 0.2125815536D+01 0.7246461665D+00 -0.100D+01 0.4028515414D+01 0.2714323506D+01 0.1106616676D+01
-0.500D+00 0.3052853410D+01 0.1908924555D+01 0.6495885258D+00 -0.500D+00 0.3647594392D+01 0.2442650065D+01 0.9960293566D+00
-0.100D-05 0.2536165936D+01 0.1607897543D+01 0.5568388024D+00 -0.100D-05 0.3063213311D+01 0.2080891485D+01 0.8639034407D+00
0.100D-05 0.2065872254D+01 0.1352528596D+01 0.4811912736D+00  0.100D-05 0.2539372964D+01 0.1777483928D+01 0.7571119023D+00
0.500D+00 0.1721413939D+01 0.1151843921D+01 0.4193581247D+00  0.500D+00 0.2149785577D+01 0.1536311542D+01 0.6690279584D+00
0.100D+01 0.1508171869D+01 0.1007249934D+01 0.3693196976D+00  0.100D+01 0.1895838229D+01 0.1355195915D+01 0.5953030788D+00

TABLE 11:     EVALUATION OF (3.3) WITH THE PARAMETERS OF TABLE 5
B1 =1.00 B2 =2.00 α1 =0.50 α2 =0.50 k1 =0.75 k2 =0.50 h=1.30
t = 1.0                                                       t = 2.0
x         y = 0.5          y = 1.0          y = 2.0           x         y = 0.5          y = 1.0          y = 2.0
-0.100D+01 0.1215462816D+01 0.6152747943D+00 0.8651621503D-01 -0.100D+01 0.1409106343D+01 0.9096599864D+00 0.2887843611D+00
-0.500D+00 0.1163055654D+01 0.5753373822D+00 0.7990410391D-01 -0.500D+00 0.1341539137D+01 0.8466976579D+00 0.2664133143D+00
-0.100D-05 0.9874016041D+00 0.4919585167D+00 0.6927746003D-01 -0.100D-05 0.1146315499D+01 0.7336723826D+00 0.2356964050D+00
0.100D-05 0.7815741113D+00 0.4124188405D+00 0.5998230907D-01  0.100D-05 0.9243614237D+00 0.6297186967D+00 0.2099005464D+00
0.500D+00 0.6644714116D+00 0.3568329302D+00 0.5289787982D-01  0.500D+00 0.7942123313D+00 0.5543685131D+00 0.1894226069D+00
0.100D+01 0.6295333035D+00 0.3302079888D+00 0.4848980574D-01  0.100D+01 0.7491675277D+00 0.5123936275D+00 0.1745085757D+00
t = 3.0                                                       t = 4.0
x         y = 0.5          y = 1.0          y = 2.0           x         y = 0.5          y = 1.0          y = 2.0
-0.100D+01 0.1494587627D+01 0.1056061789D+01 0.4467518186D+00 -0.100D+01 0.1544979193D+01 0.1146438142D+01 0.5637302445D+00
-0.500D+00 0.1420881201D+01 0.9825899130D+00 0.4130623142D+00 -0.500D+00 0.1467960579D+01 0.1067027175D+01 0.5223600131D+00
-0.100D-05 0.1218548312D+01 0.8573981347D+00 0.3692579552D+00 -0.100D-05 0.1261952985D+01 0.9352467766D+00 0.4700374956D+00
0.100D-05 0.9909044413D+00 0.7437079207D+00 0.3329918809D+00  0.100D-05 0.1031401823D+01 0.8163439964D+00 0.4270328546D+00
0.500D+00 0.8560158487D+00 0.6602467352D+00 0.3037889748D+00  0.500D+00 0.8940634275D+00 0.7284903973D+00 0.3921511762D+00
0.100D+01 0.8068782317D+00 0.6112654842D+00 0.2813293053D+00  0.100D+01 0.8427176845D+00 0.6755497530D+00 0.3645710220D+00

References

[1] Abramowitz S, Stegun IA (1965) Handbook of Mathematical Functions, AMS 55, Dover
Publications Inc., New York, 1046 pp

[2] Amos, DE, Walker, G (2011), "Integrals Related to Heat Conduction and Diffusion,"
https://thermalhub.org/resources/459 , www.beckeng.com , and

[3] Amos DE, Burgmeier JW (1973) Computation with Three-Term Linear Recursion Relations,
SIAM Review, Vol. 15, No. 2, April: 335-351

[4] Gautschi W (1967) Computational Aspects of Three-Term Recurrence Relations, SIAM Review,
Vol. 9, No. 1, January: 24-82

23
[5] Beck JV, Cole KD, Haji-Sheikh A, Litkouhi B (1992) Heat Conduction Using Green’s Functions,
Hemisphere Press, Washington D.C., 523 pp

[6] Amos DE (1990) Computation of Exponential Integrals of Complex Argument, ACM Trans. Math.
Software, Vol. 16, No. 2, June: 169-177

[7] Amos DE (1990) Algorithm 683, A Portable Subroutine For Exponential Integrals of Complex
Argument, ACM Trans. Math. Software, Vol. 16, No. 2, June: 178-182

[8] Amos DE (1978) Evaluation of Some Cumulative Distribution Functions by Numerical Quadrature,
SIAM Review, Vol. 20, No. 4, October: 778-799

[9] Carslaw HS, Jaeger JC, (1948) Conduction of Heat in Solids, Oxford Univ Press, London, 386pp

APPENDIX

Verification of T1 (x,y,t) and T2 (x,y,t) as Solutions to Problem I

We can write the terms containing x, y and t in (2.18a) in the form W  t (U1  U2  U3 ) where
U m , m  1, 2,3 are the terms of U and the real part of any complex number  as Re( )  (   ) / 2 . It
is a straight forward computation to verify that not only t Um , m  1, 2,3 satisfy the heat equation but
also each of their complex conjugates and hence Re(W) satisfies the heat equation also.

However, some analysis is required for the boundary conditions along the boundaries x=0 and y=0. The
uniform convergence of the integrals in T1 ( x, y, t ) and T2 ( x, y, t ) along with the uniform convergence of
the derivative integrals makes differentiation under the integral sign valid. The derivations presented
below provide the model for computing derivatives.

Differentiation of T1 (x,y,t) and T2 (x,y,t)

In the sections to follow, we will need to differentiate T1 ( x, y, t ) and T2 ( x, y, t ) with respect to x and y.
This means taking derivatives of the integrand in (2.18a),
(A.1)                                    
Re U ( j ,  t )V3 j (  )  U ( j ,  t )V3 j (  )  U ( j ,  t )V3 j (  )  , j  1, 2 .
1
2                                                      

Now with the chain rule, we have
                                  1  U ( j ,  t )  j                 U ( j ,  t )  j              
x
                         
Re U ( j ,  t )V3 j (  )  
2        j        x
V3 j (  ) 
 j         x
V3 j (  ) 

(A.2)                                                                                                                 
 U  j
                        

 Re            V3 j (  ) 
  j x
                        

and direct differentiation of the appropriate quantities in (2.18a) gives

24
U
 U B ( ,  , t ) ,                     U B ( ,  t )  erfc( )  e 2                    t  2t
erfc(   t )

(A.3)
 j           2   2  x                      3 j
2   2                   j                  3 j
2   2
                              ( 1)                                                   ( 1)
j                                             j                                                           j
,
x                2 t         x                              2 t                       x                                2 t
Then

(A.4)
x
                                        
Re U ( j ,  t )V3 j (  )  Re U B ( j ,  t )(1)2 j  2   2 V3 j (  )  /(2 t ) .
                           j


And similarly for y,

                                 1  U  j                U  j              
y

Re U ( j ,  t )V3 j (  )  
2   j y

V3 j (  ) 
 j y
V3 j (  ) 

(A.5)                                                                                              
 j     i       j       i
        ,               .
y     2 t       y       2 t
Then,

y
                           i

Re U ( j ,  t )V3 j (  )  U B ( j ,  t )V3 j (  )   U B ( j ,  t )V3 j (  )   /(2 t )
2                                                               

 U B ( j ,  t )V3 j (  )   U B ( j ,  t )V3 j (  )   /(2 t )
i
2                                                               

Now we have the difference of a complex number and its conjugate, which yields (2i) times the
imaginary part,

(A.6)
y
                                                     
Re U ( j ,  t )V3 j (  )   Im U B ( j ,  t )V3 j (  )  /(2 t ) , j  1, 2 .                           
Boundary Conditions on x=0

In this part of the verification process, we show that the flux boundary conditions on x=0,
T2                     T1
(A.7)                                              k2                     k1                 h[T2  T1 ] x 0 ,
x     x 0             x    x 0

are satisfied by (2.18a). Here,
 j ( )       x      2   2  iy 
j                                                               2  1/  j
j
                                   ,           j
2 t                         2 t
and we denote
 j ( )              i  y                                                                 1                            1        
0                                , j  1, 2 ,                   =h ( )=h                                                           ,
2 t                   2 t                                                         k 2   2                      k1    12   2 
x 0                                                                     2 2                                            
Then,
h[T2  T1 ] x 0 

 Re U  0 ,  t V1 (  )  r dr                                                                         
                                                                       
h[G2  G1 ] 
4h t
         0
  0
4h t
Re U  0 ,  t V2 (  )                    r
dr

and note for x = 0, y  x  0 and  / 4     / 2 applies. Reduction with
 F /k     F /k 
V1 (  )  V2 (  )   22 2 2  21 1 2 
 2    1   
gives

25
                                                           
  F / k                F1 / k1  


 Re U  0 , 
4h t            
(A.8)            h[T2  T1 ] x 0  h[G2  G1 ]                                                         t         2       2
           2 
   dr .
   0                           
2
2
2
1      r
2

Restating U in (2.18a) and using (A.3), we have
1                                          U ( ,  t )
U ( ,  t )  ierfc                                        U B ( ,  t ) and                                U B .
2 t                                               
Then
1  1 U ( ,  t ) 
(A.9)                                U ( ,  t )  ierfc                                           
 2 t
               

and (A.8) becomes

2h
       
                         F2 / k2              F1 / k1   
h[T2  T1 ] x 0  h[G2  G1 ]                              Re 2         t ierfc  0                                                   dr
    0                                

2
2
2
12   2  
           r
(A.10)
2
        
   1 U                   F /k     F /k    

    Re      (  ) 
( 0 ,  t )  22 2 2  21 1 2  
 2  
         1     r
 
dr
0        
If we let
Fj / k j
(A.11)                                                  g j (  )  2 t ierfc  0 (  )                               , j  1, 2
2   2
j
in the first integral of (A.10), then this integral can be written as

1                    1                     1
(A.12)            Re [ g   2    g1 ]  r dr                     [ g2  g1 ]d   2 C [ g2  g1 ]d   2 C C [ g2  g1 ]d  .
2 C1
0                                                                            2                    1 2

Now, g j (  ) has poles at   i  j and ierfc( z ) is analytic in the whole z plane. This suggests applying
the residue theorem to a closed contour composed of C1 (   ) , C2 ( ) , and a large circular arc of
radius R connecting C2 and C1 , where   Rei ,        . We also note at this point that
y  x  0, implies  / 4     / 2 . In order for the analysis to be successful, we have to show that the
line integral along the arc vanishes as R   . To this end, we examine (2.10) with n=1 and
z   0  iy  /(2 t ) and verify that (2.10), written in terms of  0 and  ,

(A.13)

ierfc( 0 )  O e  0 /  0
2 2


arg  0  Rei       / 2 and  / 4          3 / 4
             
applies. Now, (A.13) shows that    / 2 is within the range arg( 0 )   / 4 . Therefore (2.10) provides
the estimate with e 0 / 0 being the dominant factor. The relevant formulas in terms of  are
2
2

y 2 R 2 ei (2  )                  y 2 R 2 cos 2
(A.14)                       0 
2
and Re( 0 )  
2
,
4t                                   4t
which, for  in the range given by (A.13), puts 2 in the left half plane where cos 2  0 and
Re( 0 )  0 . We use the worse case when 2  2   / 2 or 3 / 2 in (A.14) to obtain Re( 0 )  0 for
2                                                                                       2

the estimates
e 0               e Re( 0 )
2                  2
 1 
                       O 2                and        ierfc( 0 )  O 1/ R 2  .
       2
0                 2
0              R 
With the factor (  2   2 )1 in (A.11) g j (  )  O(1/ R4 ) on the circular arc   Rei ,        .
j

26
Integration of g j (Rei ) on  between  and    gives an extra factor of R from d   i Rei d
and the line integral goes to zero faster than O (1/ R 3 ) for R   .

With the vanishing of the line integral connecting C2 and C1 , we are left with the line integral on
C1 + C2 which, according to the residue theorem, is equal to 2 i times the sum of the residues at the
poles within the contour. Then with (A.12),

1                      2 i
 Re [ g2  g1 ]  r dr  2 C C [ g2  g1 ]d   2 sum of residues at   i 1 ,   i 2 
0                                  
                                     
1       2

            y  1 F2                        y  1 F1  
  i  2 t ierfc                     2 t ierfc                    G ( y , t )  G1 ( y , t )
            2  t  2  2 i k2              2  t  2 1i k1  2 2
               2                              1            
keeping in mind that  j  1/  j , j  1, 2 . This integral, when multiplied by its coefficient 2h /  ,
cancels the leading term in the expression for h[T2  T1 ] x 0 in (A.9) and we have
2
         
        1      U                           F /k     F /k    
    (  )  0
(A.15)        h[T2  T1 ] x 0                 Re         ( ,                                     t )  22 2 2  21 1 2       dr .
0                                                           2  
         1      r
 
Now we want to compare direct differentiation of each T j ,                                                 j  1, 2, with h[T2  T1 ] x 0 in (A.15). Then,
T j              2(1)3 j (2 t )                      U
                 (1)3 j (k j )                        
(A.16) k j                                                     
Re     ( 0 ,  t )                   2   2 V3 j (  ) 
        dr
x    x 0
                      0                         2 t
j
   r

where we have used (A.2) and (A.3). But

k j  2   2 V3 j (  )
j

k  1   12   2      k    2
2

2   2  F / k
2   2     F /k  1
 21 1 2 
 2
k1 12   2  k2               2   2  2       1    2 t
2                 2
2 t
1      F2 / k2     F /k 
                2          21 1 2 
2 t  (  )  2       1   
2

and (A.16) becomes
T j                          2
      
   1 U                   F /k     F /k    
k j
x

        Re    (  ) 
( 0 ,  t )  22 2 2  21 1 2  
 2  
         1     r
 
dr .
x 0                      0      
Now we recognize the right side of this equation as h[T2  T1 ] x 0 from (A.15). Therefore, (A.7) is
satisfied.

Boundary Conditions on y=0

In this part of the verification process, we show that the flux boundary conditions on y=0,
T j
(A.17)                                                                  k j                  Fj ,       j  1, 2 ,
y     y 0

are satisfied by (2.18a). We start the derivation by denoting

 j ( )                       x  2   2  i y                            x 2   2
0 
j                                            j
                                                     =                   .
2 t        y 0
2 t                                   2 t
y=0

Differentiating (2.18a) with respect to y and applying (A.6) yields
27
T j               G j           2(1)3 j 
y
( x, 0, t ) 
y
(0, t ) 
       Im U ( 0 ,  t )V3 j (  ) 
0
B
                                  r
dr , j  1, 2

for x  y, or x  0 . Since we have shown the equivalence of all paths C2 ( ) in the interval 0     / 4 ,
 can be taken to be zero with   1 . This makes the argument of the Im(*) function purely real. Since
the imaginary part of a real number is zero, we have
T j                        G j
k j          ( x, 0, t )  k j          (0, t ) .
y                          y
Computation of G j / y from G j in (2.18) gives
G j                     y 
k j      ( y , t )  F j erfc         
y                       2  jt 
        
and with y=0 we get (A.17).

Initial Conditions for t  0

The final part of the verification procedure examines the way in which the initial conditions are
satisfied. While direct evaluation at t=0 is not possible, the limiting form T j ( x, y, 0)  0 , j  1, 2 is
considered an acceptable replacement for T j ( x, y, 0)  0 , especially for semi-infinite regions (e.g.
T ( x, t )  T0 erfc[ x /(2 t )] , 0  x   for the semi-infinite slab initially at zero temperature). Thus, with
this understanding, (2.18a) satisfies the initial conditions of Problem I.

In summary, we have verified that the differential equations, boundary conditions, and initial conditions
as stated in (1.1), (1.2) and (2.1) are satisfied by (2.18a). In other words, (2.18a) solves Problem I and
the combination of (2.18a) and (2.18b) provides the numerical means for computation.

Uniqueness of the solution can be surmised from a standard energy argument. Assume that two distinct
solutions exist and consider the difference D  0 , t>0. Because of the linearity of the problem, D
satisfies the diffusion equations and the internal boundary conditions on x=0 with zero flux values on
y=0. Since there is no energy source, D cannot change from its initial value of zero. Consequently the
two solutions must be equal, contrary to the assumption that they were distinct. A general version of this
argument can be found in [9] with citations to the literature that discuss the assumptions necessary to
make solutions unique.

28

```
DOCUMENT INFO
Shared By:
Categories:
Tags:
Stats:
 views: 24 posted: 11/16/2011 language: English pages: 28