Classroom Tips and Techniques: Diffusion with a Generalized Robin Condition Robert J. Lopez Emeritus Professor of Mathematics and Maple Fellow Maplesoft  
| Introduction | | | 
Recently, while perusing the backlog of unanswered questions on MaplePrimes, I came across an interesting boundary value problem in one-dimensional diffusion. In particular, the problem could be interpreted as describing heat transfer in a rod, with the left end insulated, and the right end subject to a generalized Robin condition. In the typical nonhomogeneous Robin condition, the "driving term" is known. In this problem, the nonhomogeneous term was the solution of an initial value problem for an ordinary differential equation that in turn, was driven by the temperature on the right end of the rod. Table 1 (where subscripts denote partial derivatives; and the prime, differentiation with respect to  ) states the problem in detail. The diffusion equation | 
| The initial temperature distribution | 
| The insulated left end | 
| Robin condition at right end | 
| Initial value problem for  | 
| Table 1 One-dimensional diffusion with a Robin condition whose nonhomogeneous term is determined by an initial value problem for an ordinary differential equation | 
Maple`s pdsolve command will solve one-dimensional evolution equations numerically (by finite difference techniques), and Maple`s dsolve command will solve (numerically) initial value problems for ordinary differential equations, no provision exists in Maple for the kind of coupling that occurs in the BVP given in Table 1. But the problem does yield to a simple finite difference scheme, as we will show below. In addition, since the BVP in Table 1 is linear with constant coefficients, it yields to a Laplace transform solution. While  and  , the Laplace transforms of  and  , respectively, are easily found, the expressions are sufficiently complex that inverting them analytically is apparently impossible. But at least two techniques for inverting the Laplace transform numerically succeed. The first is just a brute-force numerical evaluation of the Bromwich inversion integral; and the second, a more sophisticated numerical evaluation of this integral along a deformed contour of finite length. 
|
| Solution by a Finite Difference Scheme | | | 
Nearly every numerical analysis text that discusses the numeric solution of partial differential equations will exposit the simple explicit method for solving the one-dimensional heat equation. For the BVP in Table 1, we elected this technique on the grounds of its utter simplicity, even though there is a stepsize restriction for stability. 
For the sake of completeness, we summarize this technique here. Discretize with stepsizes  , and  , where  is required for stability. (We chose  , and with  , we had  .) Taking the discrete grid as  , and  , and using simple forward differences to approximate the derivatives in the one-dimensional heat equation, we have 
where  approximates  . Defining  , and isolating  , we have 
Notice that it takes three values of  at time  to determine one new value of  at time  . With this formula, the nodal values  , can be determined at time  , provided the nodal values are known at time  . When  , we have 
To obtain the value  , which falls outside the interval  , write the homogeneous Neumann condition  in terms of the forward difference 
so that  . At the right end of the interval, when  , we have 
To find a value for  , we discretize the derivative appearing in the Robin condition, obtaining 
where  approximates  . From this, we get 
so that we finally have 
Since  , values of  are computed in conjunction with a corresponding finite-difference solution of the ODE. For this, we write 
so that 
Of course,  . After a few experiments, it became apparent that the solution need not be extended beyond  , so  . Table 2 summarizes the Maple code that implements the finite-difference calculations. | • | Initialize problem parameters. |
| 
| | • | Initialize  . |
| 
| | • | Initialize 
|
| • | Advance  at 
|
| • | Advance  at "interior nodes" |
| • | Advance  at  |
| • | Advance 
|
| 
| Table 2 Code for finite-difference solution | There are some 8811 nodal values. Figure 1 is a graph of the solution surface  , drawn using only 891 such values. From the graph, it appears that the solution rapidly approaches a limiting value of 0.1. | Figure 1 Graph of solution surface  as determined by finite-difference solution | Figure 2 is a graph of  . It appears that  is rapidly asymptotic to a limiting value in the vicinity of 0.1. | Figure 2 Graph of  from finite-difference solution | 
(The code for Figures 1 and 2, and all other figures, is in the respective table holding the graphs. Likewise, cells tinted in yellow contain hidden input. To see such code, use the Table Properties dialog to select the Show Input option.) 
|
| Solution by Laplace Transform | | | | Obtaining the Laplace Transforms | | | 
To effect a solution of the BVP in Table 1 by the Laplace transform, transform the heat equation with respect to  , obtaining 
where subscripts indicate differentiation. Since  , we can write this ordinary differential equation as 
the general solution of which is 
The homogeneous Neumann boundary condition on the left leads to  , so that we can write the solution as 
The Robin condition on the right involves  so the IVP  , must be solved first. Taking the Laplace transform of the ODE, we get 
Since  , we have 
The Laplace transform of the Robin condition  leads to 
and hence to 
and  
Recalling that  , we can then determine  from the equation 
| (1) |
| 
(Maple code generating (1), and output in any other colored cell, is hidden in the cells. To see such code, use the Table Properties dialog to select the Show Input option.) 
The expression for  is surprisingly cumbersome, so we obtain it in Maple as 
| (2) |
| 
Hence, we now know  and  . For  we obtain 
| (3) |
| 
Fortunately, this can be written as 
| (4) |
| 
For  , we obtain 
| (5) |
|  
|
| Extracting Information with Tauberian Theorems | | | 
Under suitable conditions on the functions  and  , it is possible to extract "initial" and "terminal" values for these function from their Laplace transforms. Theorems justifying these results are typically called Tauberian, in contrast to Abelian theorems that provide information in the opposite direction. First, if  is of exponential order, and  is its Laplace transform, then  . If we apply this to  and  , we will have some evidence that these transforms are not a priori incorrect. Thus, we find and 
Both initial values are recovered. Of greater importance is the "final value theorem" that gives the limiting value of a function  in terms of the behavior of its transform as  . In particular, if  is bounded and its limit at  exists, then that limit is given by  . From the finite-difference numerical solution, we opined that both  and  tended to finite limits fairly rapidly. Since we don`t know analytically that these functions are well behaved at  , the following two limits are informative, but not conclusive. 
These limiting values are consistent with the finite-difference solution. 
|
| Numerical Inversion of the Laplace Transforms | | | | Justifying the Bromwich Integral | | | 
Sufficient conditions for the Bromwich (line) integral  
along  , yielding the function  from its Laplace transform  , include  being analytic in some right-half plane  and being  with  in that half plane. By the calculations at the end of the last section, both  and  are  . However, the Bromwich inversion integral will invert the transform  to 1, even though this transform does not satisfy all the sufficient conditions just stated. We implement this inversion integral with  because the zeros of 
the denominator of both  and  , appear to fall on the nonpositive real axis. We determine this not only from Figure 3, a graph of the denominator of  along the real line, but also by applying the Principle of the Argument to  on the circles  . The circles  are centered at  on the real line, and have a radius of  . Thus, as  increases, more and more of the half-plane  is enclosed. | Figure 3 Zeros of the denominator of  along the real line | Numerically evaluating a succession of integrals 
for increasing values of  shows that each such integral is zero. Since  is the number of zeros inside  and  is the number of poles, and  has no poles inside  , it then has no zeros inside either. 
|
| Evaluating the Bromwich Integral Directly | | | 
Direct numerical evaluation of the Bromwich integral generally has unsatisfactory results. Along a vertical line, the integrand oscillates, making it difficult to integrate to infinity. This is true for  , which, along the line  has the form. 
| (6) |
|  
The real part of (6) is an even function of  , whereas the imaginary part is odd in  . Thus, the Bromwich integral for  on any interval  can be reduced to 
As might be imagined, even a numeric evaluation of this integral is a challenge. If, instead, we integrate to a "large" finite upper limit (we used 10,000), and use one of the NAG library methods (_d01akc, designed for oscillating integrands on a finite interval), we can obtain the graph shown in Figure 4. | Figure 4 Numeric inversion of  by direct numeric evaluation of the Bromwich integral | 
Figure 4 uses ten computed points for the graph of  , and this takes several minutes in Maple. But the apparent asymptotic behavior is in agreement with the prediction of the Tauberian theorem in the Laplace transform section. Similar results for  are obtained in Figure 5 by numeric evaluation of the Bromwich integral on a uniform grid  , where  , and  . Since  has the same symmetry properties along  , the computations for Figure 5 are similar to those for Figure 4, except that the upper limit of integration is taken as  . It took more than two minutes to obtain Figure 4, and more than ten to obtain Figure 5. | Figure 5 Numeric inversion of  by direct numeric evaluation of the Bromwich integral | 
Like Figure 4, Figure 5 shows a limiting behavior for  consistent with the Tauberian prediction in the Laplace transform section. 
|
| Deforming the Bromwich Contour | | | 
The literature contains a number of methods for the numerical inversion of the Laplace transform. In [1], Talbot`s approach of deforming the Bromwich contour is implemented in variable-precision software. The inversion is given by 
where  
and  is the precision of the computing device. The relationship between  , and  is empirical, and allows for an increase in device precision as  increases. Our preliminary investigations suggest inverting  and  on  , so with  , the default ten digits in Maple are more than enough for our purposes. This allows us to use Maple`s built-in integrator rather than the trapezoid rule suggested in [1]. A graph of  obtained by this method appears in Figure 6. Drawn in about 15 seconds with 20 computed points, it is consistent with our Tauberian results and Figure 4. | Figure 6 Graph of  obtained by numeric evaluation of the Bromwich integral using a deformed contour | Figure 7 is obtained in about 80 seconds by applying the same algorithm to  . Again, the graph (drawn from 110 computed points) is consistent with our Tauberian results and with Figure 5. | Figure 7 Graph of  obtained by numeric evaluation of the Bromwich integral using a deformed contour | 
|
|
|
| References | | | 
[1] Abate J. and Valkó P. P., Multi-precision Laplace transform inversion, Int. J. Numer. Meth. Engng 2004, 60:979-993. 
|
Legal Notice: © Maplesoft, a division of Waterloo Maple Inc. 2010. Maplesoft and Maple are trademarks of Waterloo Maple Inc. This application may contain errors and Maplesoft is not liable for any damages resulting from the use of this material. This application is intended for non-commercial, non-profit use only. Contact Maplesoft for permission if you wish to use this application in for-profit activities. 
|
