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.
The Robin Boundary Condition
The general one-dimensional Robin boundary condition is of the form , where is a known function, is a bounding edge of the interval over which the accompanying BVP has been set, and and are constants. In the BVP of Table 1, , but is not a "known" function. The normal derivative, denoted by , becomes for the BVP in Table 1.
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
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
(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
Hence, we now know and . For we obtain
Fortunately, this can be written as
For , we obtain
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
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.
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.
?
