вход по аккаунту


Rapid Solution of Stiff Differential Equations and Accurate Numerical Laplace Inversion of Steep and Oscillatory Functions using IMN Approximants.

код для вставкиСкачать
Dev. Chem Eng Mineral Process., l0(1/2), pp.143-164, 2002
Rapid Solution of Stiff Differential
Equations and Accurate Numerical Laplace
Inversion of Steep and Oscillatory Functions
using IMNApproximants
0. Taiwo* and R. King
Measurement and Control Group, Institute of Process and Plant
Technology, Technical University of Berlin, Hardenbergstr. 36a,
D-10623 Berlin, Germany.
lm approximants are a fast and convenient method of solving initial value problems
in linear stir differential algebraic equations as well as obtaining the numerical
inversion of Laplace transforms In the past it was impossible to use them to obtain
sufficiently accurate inversions of certain steep and highly oscillatory responses as
useable values of N had to be relatively small not only to ensure reliable evaluation of
I,, constants but also in order to avoid undue rounding errors in the computed
results However, the development of computer algebra systems such as Mathematica
which permit injinite precision computation has provided greater latitude for the
application of the method This work is an exposition of the potency of I m
approximants in accurately and cheaply inverting functions in the Laplace domain
whose time functions are steep, oscillatory or s t i r Previous applications of the
method to some examples in the literature led to wrong conclusions as the capabilities
of the method were not fully explored We show how to obtain very accurate results in
these circumstances using both the global and step-by-step methods The results of
using lMN
step-by-step technique to rapidly solve stifl diflerential equations of a large
staged process are also presented As Matlab gives inaccurate results, Mathematica
has been used to compute both the transferfunction and the analytical expressionfor
the time response of the plant Extended values of M and N as well as the ranges of
the corresponding Im constants are tabulatedfor IMNapproximants of full grade.
Im approximants have wide applicability in tackling practical probIems in
engineering and science. Their applications range fiom repeated computation of
closed-loop responses when designing control systems using the method of
* Authorfor correspondence (e-mail.femtaiwo@yahoo com), on leave@om Chemical
Engineering Dept, Obafemi Awolowo University, lle-ge, Nigeria
0 Taiwo and R. King
inequalities [l-51, the simulation of diffusion and reversible reaction in a sphere [ 6 ] ,
the identification of kinetic parameters in adsorption [7], the solution of differential
equations of a periodically cycled plate column [S], the design of circuits [9] and the
analysis of boundary conditions in both filled and empty tubes [lo]. The IMN
approximant numerically inverts the Laplace transform of the function describing a
phenomenon, either recursively or in a global manner rapidly yielding sufficiently
accurate results in typical applications.
The present work is to fiuther explore the untapped potentials of the technique for
solving process systems engineering problems. The first problem solved here is a
large-scale staged tridiagonal process with stiff differential equations. Such plants are
commonplace in the chemical plant. The particular example solved here is a benchmark and has been identified by Kim and Friedly [ 1 13 as an optimally designed staged
system. Im step-by-step technique is known to be immune to stiffness and has been
used to solve the resulting differential equations, while Mathematica has been used to
compute both the transfer function and the analytical expression for the time response
because Matlab gives inaccurate results. It has been decided to include the numerical
inversion of the sine wave as an example since this is a common test signal in
practical work. Moreover, it serves to highlight the way to obtain accurate results for
open- or closed-loop systems possessing oscillatory responses. Both the step-by-step
technique and higher lm approximants, applied in the global form, can solve this
problem. The other two examples have very steep responses. The first one is an openloop unstable process that can easily arise in the optimal steady-state design of
chemical reactors [12, 131. Either the step-by-step technique or higher Im
approximants, applied in global form, can provide a solution to this problem. The
other process with a steep response involves concentration profiles in an empty tube
at high Peclet Numbers. The recommended procedure for solving this problem is the
use of higher Im global method because the mathematical basis for the extension of
the step-by-step technique to partial differential equations is still being developed.
Using the notions of acuity and equivalent-width which have been explained in the
sequel, and using the derivation of the Im formula in the Appendix, it becomes clear
why higher N for the global method or small-step h values for the step-by-step
technique, is needed to obtain Im results of sufficient accuracy for functions which
vary rapidly at large t values.
Since Im approximants of full grade have certain properties (such as being
impervious to stiffness) which make them attractive for practical use, we tabulate both
the values of M and N and a sample of Im constants for Im approximants of full
grade up to N = 540. By using these higher approximants, it is not only possible to
accurately invert functions for which it had earlier been erroneously concluded that
the Im approximant did not produce accurate results, but also to check the accuracy
of the computed results for cases where no analytical solution exists. In view of the
large values of the Im constants involved for some of the examples solved here, and
in order to minimize rounding errors, infinite precision arithmetic available in
Mathematica has been exploited in the computations involving the global method.
Where applicable, the 1~ results have been compared with those obtained using
the method of Honig and Hirdes [I41 which is based on Fourier series expansion.
Higher Im results have been used to accurately estimate the error in the results of
both methods.
Rapid Soluiion of StiffDifferential Equations using IMNApproximanis
Zakian's Im approximants 115,161
Let f(t) have a Laplace transform:
where Re(s) > CT and f(t) is continuous fort 2 0 and of exponential order CT . The Im
approximant of f(t) (see the Appendix for the derivation) is an expression of the form:
IMN(f,t) = t " C K , L ( f , n l )
1 4
Many different sets Ki , a,may be defined, but for the purpose of this work, they
are assumed to satisfy the relation:
where eLfir is the (W)Pad6 approximant of e". Let
known that under certain conditions on f(t):
= min,{Re(a,)}.It is
f(t) Im(f,t) = O(PN+') ,t + 0'
whenever, ,&
> 0. Numerical tests have been performed [15] to determine all
(M, N) in 0 IM < N I 2 0 for which, ,&
> 0. Taiwo et al. [lo] carried out hrther
tests and produced a table of all (M, N) in 21 IN I60 (M < N) for which, ,&
In this work, because of the need to use still higher approximants in order to obtain
accurate results for our particular examples, new tests were done to determine all
(M, N) in 61 I NS 540 (M < N) for which &MN > 0 (see Table l), where M h in
Table 1 denotes the minimum value of M for the relationship to hold.
Table I . Minimum values of M (Mmi,Jfor 57 IN 5 540 such that
arefull grade, that is, satisfjr (3) and di,
57-71 72-88 89-107 108-129 130-155 156-183 184-214 215-249
288-329 330-376 3 77-426 427-48 1 482-540
0 Taiwo and R. King
A sample of the associated values of Ki , ai for some of Im approximants are also
tabulated (see Table 2). This table indicates that while Ki increases very rapidly with
N, the increase in the ai is more gradual. It is this rapid increase in Ki that has
necessitated the use of infinite precision arithmetic not only in the determination of
the constants but also in the actual inversion of the Laplace transfonns involving
higher lm approximants, so that the rounding errors involved in the computations
will be minimized.
Furthermore, if f(t) is of exponential order (r < 0 as t + 03, then:
f(t) - Im(f,t) = 0(tM”) ,t 3 00
Thus Im approximant is capable of good approximation at both large and small t
values; for a fixed N, a trade-off between accuracy at large t and accuracy at small t is
effected by changing M.
Equation (2) is useful in the numerical inversion of Laplace transforms.
Computationally, for N even, (2)may be reduced to an expression of the form.
Im(f,t) = 2t-’
since the constants Ki , cq occur in complex conjugate pairs. The constants Ki , ai are
obtained by evaluating the residues and poles in the partial fraction expansion of
e i i , see for example [ 171.
The inversion formulae give rise to a number of techniques for the solution of
initial value problems of the form
- = A x ( t ) , x(0) = x,
where A is an n x n constant matrix, x(t) is an n-vector and ~0 is the given initial
Firstly, the Laplace transform of the solution x(t) may be obtained explicitly from
(7) and the inversion formula (2) applied in a global manner. Secondly, (7) may be
Laplace transformed to give
where L(x,s) is the Laplace transform of the solution x(t). The vectors L(x, ai/t) may
be obtained by solution of the algebraic systems:
((ai/t)I-A) L(x, q/t) = x(O), i = 1,2,. ..,N
Rapid Solution of Stiff Differential Equations using lMNApproximants
for particular values o f t , and thus the inversion formula (2) may then be applied,
again in a global manner. When N is even, the formula (6) may be applied, whence
only N/2 systems (9) need to be solved.
However, the accuracy of Im at small t may best be exploited in a step-by-step
-f =
technique. Let f, r = 0, I, 2,.. be points on the half-line t 2 0 such that to = 0,
h and h > 0, where h depends, in general, on r. Let x, denote a numerical approximant
to x(G) Let X &) denote the Laplace transform:
Then replacing x(t) by x(t,+k) in (7) and taking Laplace transforms we obtain:
Clearly, the inverse of X,(s) evaluated at h is x(f+,). The Im approximant of x(h+h)
evaluated at h is:
On combining ( 9 ) and (lo), the following recursion formula is obtained:
((cq/t)I-A) Xi,, = x,, i = 1,2,...,N
x,~= h - ' x K , X , , r
r=0,1,2 ,...
The main work at each recursion lies in the solution of N uncoupled n x n algebraic
systems (13a). In fact, if N is even, the vectors KJ,, occur in complex conjugate pairs
and (13b) reduces to an expression of the form:
x,l= 2 h-'
Re(K,X,,,) r = 0, 1,2,...
hence only N/2 systems (13a) need to be solved. In practise, considerable economy is
achieved by keeping the step-length constant and calculating once, and for all, upper
and lower triangular factors of the matrix ((ai/t)I-A) for each i. Moreover, any special
structure of the A matrix, such as sparseness, is preserved in the matrices ((ai/t)1-A)
and can therefore be exploited. Further economy can be achieved by the use of small
N, while accuracy is controlled by the step-length.
Additionally, the use of small N reduces round-off error for the software
possessing only finite precision arithmetic.
(Note that all a, , K,tabulated occur in complex conjugatepairs. Only a member of thejirst and last pair for each N is listed and the
other is obtained as its conjugate. The i appearing after each imaginarypart denotes the imaginarypart, while yex denotes y x l fl )
Table 2. The constantsfor higher Im approximants
329 I
1 1
41 5.924816479 + 1.687446608i
0.9962499868 + 717.71782631
476.898505 1544592 + 1.6933944549i
0.02966058509 + 819.11877380754i
543.811 189651 + 1.69516286057i
10.6806122144463 + 900.71220844593
601.59671838312112+ 1.713923505287i
0.90853850863 + 914.61384771148i
607.435491 14196 + 1.6992723896593
0.979050379372 + 922.573287924i
612.7375075658+ 1.6995773882i
362.2670128305+ 1.685424511946i
, 0.0018196715
~ + 626.2887162478i
4.675746361 + 450.6290834331
30 1.025979203 + 1.6901349846i
~ 1 ~
Table 2 continued.
-1 -35056412684e156+ 1.14052022164e157i
1.935911333458156- 2.888 10794125753
-2.755357877585e179+ 2.295238846719e180i
-7.78264821296+ 6.0213286409283
-8.750206017174e205 + 6.959746977908e206i
-3.825317394664- 0.82421326817i
- 1.01698944621266e235+ 7.99420884155991e235i
-169477.1707896735 - 19784.0162945843
-1.4694175519657e260 + 1.005516542883e261i
-1.09362631293 - 9.71809142127i
4499014353491e262 + 3.42998 1760742e263i
10.52426252518- 0.070102682625i
-9.0528610644407e264 + 6.886360625796e265i
25 1.2 1622341391 - 2 19.0508366143i
-3.584050677799e129 + 2.8932789598119e130i
0 Taiwo and R. King
It has been shown [ 15, 161 that for Im approximants of full grade, the recursion
(13) has the useful property of A stability, i.e. when applied to the scalar problem:
- = m(t) x(0) = x,
for any Re(z) < 0 and using a constant step-length of any size, the recursion yields a
sequence {x,} which tends to zero as r -+ 0 0 , thus preserving the stability of (14).
Moreover, it is known (see, for example, [15]) that for any matrix A, all whose
eigenvalues have negative real part, the recursion yields a similarly converging
sequence when applied with constant step-length of any size to the problem (7). This
property is known as Z stability. It will be readily appreciated that a
recursion has considerable advantages over other methods, particularly when the
system (7) is stiff.
From the foregoing, the difference between the global and step-by-step methods is
explained as follows. For the global method the origin is fixed and the Im
approximant is calculated at a number of values o f t , whereby the distance from the
origin at which the approximant is calculated increases with t. However, for the stepby-step method, the origin is effectively shifted at the (r+l)th step to the point f, and
the Im approximant is calculated at a distance h from the new origin so that the
distance h from the new origin at which the approximant is calculated remains small.
It may be noted that both global and step-by-step techniques can be easily applied
to equations of the form:
- = A x @ )+ ~ ( t ) ,~ ( 0 =) X,
where u(t) is Laplace transformable and represents the input to the system. This is
dealt with in detail in [ 151.
As indicated in [lo], comparison of computed results with higher Im results
facilitates an estimation of the error for both the step-by-step and global methods.
Illustrative Examples
Example 1
Often separation equipment in the chemical plant are multi-stage systems having
linearized models with tridiagonal matrices. Such systems have been of interest to
many investigators [ 11, 18-23]. A number of such systems have characteristics which
render them difficult to analyse and control. Such characteristics include their large
scale, stifhess and ill-conditioning. Their large-scale nature has prompted
investigators [ l l , 18, 19, 221 to devise methods for reducing their dimensionality
while several robust control methods [20,21,231 have been proposed to minimize the
sensitivity caused by ill-conditioning. In this work, we propose the exploitation of Im
step-by-step method for their rapid and accurate simulation. For our example we take
the following model [ 1 1,221:
Rapid Solution of StiffDifferential Equations using IMNApproximanis
- = Ax(t) + bu(t), y(r) = cx(t), x(0)
= x0
where b = [ l , 0, 0,.... 0IT,c = [O,O, 0,...0,I],
... 0
... 0
... 0
. .
. .
0 ... - 2
0 ... 1
0 ... 0
u(t) is the plant input and y(t) is the output. This matrix was used by Evans et al. [24]
to test their eigenvalue algorithm, while Carmola and Chimowitz [22] used it to test
their model reduction technique. Kim and Friedly [l 11 noted that this A matrix
represents an optimally designed staged system.
We first consider the 50* order model. It was impossible to compute an accurate
transfer function model of this plant using Matlab, hence Mathematica can be used to
directly determine the transfer function from c(sI-Ay’b . However, as the A matrix is
tridiagonal, the transfer function can also be accurately and more rapidly computed
using Mathematica by separately determining the characteristic polynomial and the
numerator polynomial of the model. The numerator polynomial is determined from
fust principles and formulae for these are available [22]. The analytical expression for
the step response of the plant can be accurately determined using Mathematica by
decomposing the product of the transfer function and the step input into partial
fractions and Laplace inverting, since direct use of the inverse Laplace transform
function available in Mathematica did not produce a result because the personal
computer (with 128MB RAM) ran out of memory. For the determination of this
analytical expression, 40 digits of precision are needed in the computation in order to
obtain accurate results. From the analytical time response expression, it is noted that
the fractional response of the plant at t = 100 is 0.00893, signifying that the plant
exhibits a time delay of almost 100. Furthermore the response of the system is
99.99% complete at t = 2600. This is a stiff system with stiffness ratio 1 .O535x1O3.
We now investigate the results of numerically computing the step response of this
plant using 13,4 and
step-by-step approximants. Since the response of the plant is
practically zero for t < 100, we shall first investigate the performance of these
approximants for t = 50(50)2600 and consider other step-lengths later. Although we
know that the response at t = 50 is virtually zero, we include this point in the
experiment in order to test the goodness of the approximant for computing the
response of high order systems, with large relative degree (and hence large delay),
near the origin. It was found that with 13,4 approximant, although the absolute error
(2.0265~10-’)at t = 50 is small, the percentage error is 43.42 because the analytical
value of the response (4.667~10-’)at this point is small. It was therefore decided to
0 TaiwoandR King
use the Is,6 approximant to compute the single sample at t = 50, while the 13,4
approximant was used to compute the other samples for this step-length. The result of
doing this is given in the second row of Table 3. The row 13,4 (a) (see Table 3) is the
result obtained using 13,., approximant alone for t = 100(50)2600, since the response
for t < 100 is virtually zero, The row 13,4 (b) is the result obtained using 13,4
approximant alone for the interval t = 100(100)2600 since y(t < 100) =: 0. It is seen
that Im requires the smallest CPU time to do the computation. For this case and other
cases tabulated, the largest percentage error occurred during the first sample. This is
because the actual value of the response is small for small times. It was discovered
that the 13,4approximant gave sufficiently accurate results for other step-lengths,
larger or smaller than 50. It is worth pointing out that in performing the LU
factorization for this problem, the tridiagonal structure of the matrix was exploited to
reduce storage requirements and computing time [25].
Table 3. Maximum percentage errors and CPU times f o r the 5dhorder model of
Example I
and 15.6
13.4 (a)
Maximum percentage Error
Average CPU Time@)
0.43 1
The other methods with which Im approximant has been compared are the
various functions in Matlab specially written to handle stiff systems (these have
names beginning with ode) and the classical fourth-order Runge-Kutta method
(RK4TH). The Matlab hctions automatically determine the desired step-length,
although the user can observe the computed response at any desired value of time, so
in uniformity with the frrst case above the computed results were observed at
t = 50(50)2600. The best result fiom these functions was obtained using odel5s
followed by ode23T, ode23TB and ode23s. It is noted that the performance of ode23s
is much worse than those of the other three Matlab functions. On the other hand,
odel5s is more than 1.5 times faster than the other two. RK4TH required the largest
CPU time and the largest useable stable step-length is 0.685.
It was noted earlier that the Im method is immune to stiffness, hence one can use
large step-lengths in the computation. For t = 200(200)2600 and t = 300(300)2700,
13,4 approximant gave respectively the following maximum percentage errors: 1.064
and 0.416 with the corresponding CPU times of 0.1 1s and 0.055s.Still larger steps
were used with accurate results in all situations. Where one is required to compute
process responses for both short and long times for stable systems, it may be
expedient to use varying steps, with shorter ones near the origin and longer ones at
long times.
Rapid Solution of Stiff Differential Equations using IMNApproximnts
Further computations were also done for the 100" order model of (16). The
stifmess ratio in this case is 4.1 13x103,the response being 99.4% complete for t =
6000, and the fiactional response at t = 300 is 0.0014, signifying that the system
effectively has a time delay of 300. Here, 80 digits of computation are needed in order
to accurately determine the analytical expression for the time response. The computed
results were compared with the analytical ones at t = 200( 100)6000 and the relevant
results given in Table 4. As with the foregoing example, the 13,4 step-by-step
approximant requires the minimum CPU time, and the performance of Matlab
functions follows the trends described earlier. Further computations were carried out
using the I3,4 approximant with the following results. For t = 200(200)6000, although
the error (1.188~10-~)
of the computed result at t = 200 is small for this larger steplength, the percentage error (48.07) is large as the value of the response at this time
is small. Hence it was decided to compute this one sample using the Is,a
step-by-step approximant, while the remaining samples were computed using the
13,4 approximant with the maximum percentage error and CPU times given
respectively as 1.439 and 0.74s. For t = 300(300)6000,400(400)6000,500(500)6000
and 600(600)6000, the respective CPU times are OSOs, 0.39s, 0.33s and 0.27s with
the corresponding maximum percentage errors of 3.572,2.303,0.021 and 1.038.
Tabte 4. Maximum percentage errors and CPU times for the IOdhorder model of
Example 1
Ode 15s
Maximum percentage Error
1.42 1
Average CPU Time($
The following should be noted fiom the above results.
Unless otherwise stated, the largest percentage error occurred at the first sample
observed for all the methods. This is because the value of the response at this time
is small.
When it is required to compute process responses for both short and long times for
smooth stable hnctions using Im, it may be expedient to use varying steps, with
shorter ones near the origin and longer ones at long times.
The more widely separated the points are where responses are desired, the greater
the advantage of I ~ ~ over
R J other methods for solving initial value problems. It
gives sufficiently accurate results for typical problems in chemical engineering
even when large steps are used.
The step-by-step method is ideal when facilities are not available for infinite
precision computation. This is because low approximants give sufficiently
accurate results by utilizing the accuracy of the Im at small t, through the
appropriate control of h. Notice that an explicit expression is not needed for the
transfer function in order to numerically invert the transform in this case.
0 Taiwo and R. King
It is noted that a special procedure is not necessary in order to obtain accurate
results for this problem. This is because the response is smooth and monotonic, so
that the existence of a finite equivalent-width (27) (see a later section for the
explanation of this term) at long times for the small N used does not lead to
significant error as (23) gives sufficient accuracy in this case.
The Im techniques are easily programmed and programming expertise is not
required in order to implement it in different software packages.
Comments on Some IMNResults in the Literature
This section comprising Examples 2 and 3 has been included in order to properly
address some issues in the literature concerning what had been perceived as the
limitations of 1~ approximants. Additionally, it is felt that the functions involved are
readily encountered in analysis. The sine wave is a common test function and also
simulates the behaviour of a lightly damped open- or closed-loop system. Unstable
transfer functions are sometimes encountered in process engineering, principally in
the optimal steady-state design of chemical reactors [12, 131.
In the very useful book [26], three examples are given to illustrate numerical
Laplace inversion using I* approximants. Although all the results can be improved,
we deal only with the improvement of the results of the examples 9.17 and 9.18 [26,
pp.386-3871, since the results of the example 9.16 [26, p.3861 are sufficiently
acceptable. The necessary procedure required to improve the results of this latter
example even further is illustrated below,
Example 2
In this example, F(s) = l/(s2 + 1) which has the inverse f(t) = sin(t). On employing
the Ig,lo approximant, Rice and Do [26] obtain a result that is more than 25% in error
at t = 16.5. However, this result can be easily improved by using either, in a global
fashion, higher approximants or Im step-by-step method. We have reworked the
example for t = 0 31.0)20.5 and obtained the results given in Table 5 using Im global
method, ftom where it is clear that the 114,22 approximant gives a result accurate to
five decimal places and therefore of sufficient accuracy for practical purposes. If a
more accurate result is desired, then the 119,28 approximant could be used. The result
in this case (with accuracy of 11 decimal places) is more accurate than that obtained
using the method of Honig and Hirdes [14]. We have not manually set the ftee
parameters for this latter method here as it is felt that this result is already sufficiently
accurate for practical purposes. Notice that the time range used here is more than that
in the original problem. Using the 13,4 step-by-step method, the maximum error in the
computed results is 1.1 141x10” at t = 20.5 for a step h = 0.5s.
Table 5. Maximum errors in the computed results using the various approximants as
well as maximum errors by the method of Honig & Hirdes for Example 2.
Maximum Error, E,,
Location of E m , X E , , ~
3 . 9 8 106
Honig & Hirdes [ 141, automatic
setting of ftee parameters
Rapid Solution of Stiff Differential Equations using IMNApproximanis
Example 3
This is an unstable system (with unbounded response as t
transform: F ( s ) =
2s2 +3s-4
(s - 2)(s2 + 2s + 2)
with the Laplace
which has the inverse f(t) = e2t + e-%os(t)+ 2e-'sin(t). This is example 9.18 in [26,
p.3871, where it is indicated that the ratio of the analytical to the Im result is 2 . 4 1 ~ 1 0 ~
at their final time of 9.5. In demonstrating that higher Im approximants can be
applied to get a better result, we have performed the computations in the same time
interval as used in the ftrst example above. Note Erom Table 6 that the global 146.58
approximant results are accurate to six decimal places. This is probably sufficiently
accurate for practical purposes. However, if results of still higher accuracy are
, ~
needed, then higher approximants have to be used. For example, the 1 ~ 2approximant
results are accurate to sixteen decimal places. The method of Honig and Hirdes [ 141 is
inapplicable to systems having poles in the open right half plane.
Table 6. Maximum errors in the computed results using various approximantsfor
Erample 3
Maximum Error, E.,
7.7 1x 1 Q 8
4 . 9 7 1~0-'7
Location of Em20.5
step-by-step technique.
This problem has also been tackled using the 13,4and
For h = 0.125 s, they respectively give absolute errors of 0.1516 and 6 . 9 7 3 ~ 1 at
0 ~the
original final time of 9.5 s. Note that the value of the theoretical response at 9.5s is
lo*. These approximants were incapable of giving acceptable results at the
new final time. A higher approximant is needed and a smaller h in order to solve the
problem. But this was not pursued further as the problem has been solved above using
the global technique, Note that unlike other methods of numerical Laplace inversion,
Im approximant is applicable to systems having poles in the open right half plane.
For the foregoing two examples, the variable is changing relatively rapidly, so that
high focusing power or acuity is needed in order to obtain accurate results. As
explained later, this is achieved by choosing N large and h small (for step-by-step
method), thus ensuring that A ~ (a-- 1)
Example 4
This example concerns the dispersion model for an empty tube involving Danckwerts'
[27] correct boundary conditions, see also [lo, 281. The differential equation model of
this process is available in [27,28] and Seidel-Morgenstern [28] studied the effect of
the assumed boundary conditions on the computed exit response. In the present work,
we are interested in accurately determining the exit response when the correct
boundary condition is used, as well as providing novel results on required values of M
and N in order to get accurate results for different Peclet Numbers when Im
0 Taiwo and R. King
approximants are used for this problem. Here we employ only the global method as
the step-by-step technique has not been sufficiently developed for partial differential
equations. The Laplace transform of the nomalised concentration is given by:
F(x,s) =
A, exp(A,x)-3Ll exp(A, +a,x-A,)
~ ( ( 1-AI / P 4 A 2 exp(A.,)- (1 - A 2 / Pe)Al exp(A,,))
and Pe denotes the Peclet number for axial mass transport. This example is
particularly instructive in demonstrating the need to use higher 1~ approximants as
Pe (and correspondingly the steepness of the response) increases.
The smallest Peclet number considered here is 20 and the Laplace transform
inversion is done at T = 0.4(0.2)1.8, where T is dimensionless time [28]. From Table 7
we note that the results of the 19,1a approximant are accurate to 5 decimal places while
those of the 125,34 approximant are accurate to 12 decimal places. The results of using
the method of Honig and Hirdes [14] with automatic setting of f?ee parameters are
accurate to 11 decimal places. The next situation concerns Pe = 100 where the
Laplace inversion is done at T = 0.7(0.1)1.4. From Table 7, I14,22 approximant gives
results accurate to 4 decimal places while the 130,40 approximant results are accurate to
11 decimal places, and the results obtained using the method of Honig and Hirdes are
accurate to 12 decimal places. For Pe = 500, Laplace inversion is done at T =
0.9(0.02)1.08 with respective accuracies of 6 and 12 decimal places for 135, 46 and
approximants. Eleven decimal places of accuracy was obtained using the method of
Honig and Hirdes. With Pe = 10,000, the Laplace transform inversion is done at T =
0.97(0.01)1.03 giving results that are accurate to 6 and 14 decimal places for 114s,IM)
and 1268,288 respectively. When using the method of Honig and Hirdes, one has to
manually set the fiee parameters in order to obtain results accurate to 12 decimal
places. The difficulty encountered in this case is due to the step-like response (see
Figure 1). This is expected because for this very high Pe, the flow through the tube is
very close to ideal plug flow so that the exit composition response is virtually ~(T-I),
which is a delayed step response stepping at z = 1,
This is another example with steep response. By varying Pe values, we have
demonstrated that a steeper response demands higher acuity (see next section), which
means larger N, thus ensuring
(- - 1)
A(& t) ,so that as N increases the sifting
ability of the method is near perfect.
Rapid Solution of Stiff Differential Equations using IMN
Table 7. Maximum errors in the computed results using the various approximants as
well as the maximum errors when the method of Honig & Hirdes [I41 is
usedfor the dispersion problem at various Peclet numbersfor Example 4
Pe = 20
Honig & Hirdes, automatic
setting of fiee parameters
Pe =lo0
Honig & Hirdes, automatic
setting of fiee parameters
Pe =500
Maximum Error
1 O4
I. 14x lo4
I O-I2
Location of Em,
3.63x 1O5
1 O4
9.31 x lo4
1 . O ~ 1X0-13
5.1 8 x 10”
I .04
Pe= 10000
0 TaiwoandR King
Figure 1. Graph of composition against time in an empty tube with dispersion
Qualitative Analysis of Error in the Computed Results I291
The quantitative expressions for the error in Im results are given earlier. Those
expressions do not tell us why functions which have large rates of change at long
times must be solved using large N. In order to qualitatively analyse errors in Im
approximants, we note that the theory underlying the method is based on the sifting
property of the integral (A 6) in the Appendix. As N -, 00, perfect sifting occurs and
whenever t is a continuity point off, fN(t)-+ f(t). Equation (A.6) represents a linear
transformation TN as depicted in (A.8). The properties of this transformation are
entirely determined by the properties of its kernel A ~ (a.- - 1). Therefore the
properties of the truncation error depend entirely on the properties o f f and on the
properties of the kernel A N . As N + OD, certain properties of A,, approach those of
the delta function A(& t ) (defined by (A.7) with N = a)and TNbecomes the identity
transformation. One such property of the kernel is its area B defined by:
From (A.7) and (19) we obtain:
Rapid Solution of Stiff Diperential Equations using IMN
13= Z K , la,
which shows that the area of the kernel is a constant. Furthermore, for any value of N,
K , /a,= 1, thus ensuring that
the constants may be chosen to satisfy the condition
the area of the kernel is unity.
Another important property of the kernel is its focusing power or acuity. This
is the ability of the kernel to allow only values of f(h) in the immediate vicinity of
h = t to make a significant contribution to the integral in (A.6). To ensure a high
acuity, the constants ai, Ki are chosen such that the kernel is an optimal
approximation to the delta function A ( A , t ) . Thus for a given t, the kernel is a pulselike function of the variable A, attaining its pick value A,,,(t) at h = t. Thus the
equivalent-width of the kernel may be denoted by wN(t) and may be defined by:
Wdt) =
wiN( t )
The reciprocal of the equivalent width is adopted as a convenient measure of the
acuity of the kernel. The narrower the width, the higher the acuity. Thus TNbecomes
the identity transformation when the kernel has a pulse-like shape with 1-0 = 0 and
wN(t) -0 for all t. The more the quantities 1-13 and wN(t) deviate fiom zero, the
further TN departs fiom the identity transformation. However, when wN(t) is non-zero
and provided f(h) remains approximately constant during the interval:
= {a I t - w , ( t ) / 2 c a < t + W , 12
Then (A.6) may be replaced by the approximate expression:
which with (19) yields:
and therefore the corresponding truncation error is approximately (1-D)f(t). Thus if
f(h) is constant during the interval A N ( f ) the
, truncation error is almost entirely
proportional to (1-13), and any imperfection in the acuity of the kernel will make
negligible contribution to the truncation error. The foregoing discussion leads us to
consider the idea of relative acuity, a notion which involves both the acuity of the
kernel at h = t and the constancy o f f in the interval A,,, (t) . For a given f, the relative
acuity of the kernel at A = t will depend on the constancy of f in the interval A,,, (t) .
The greater the constancy, the higher the relative acuity.
0 Taiwo and R. King
When f undergoes substantial changes in the interval h N ( f the
) , nature of the
transformation TN is more clearly revealed if (A.6) is replaced by the following
approximate expression:
which involves the assumption that the kernel is approximately constant at the value
A N ( t ) during the interval A N ( t ) ,and is approximately zero outside this interval.
Equation (25)indicates that TNperforms a smoothing operation on f(h), that is fN(t) is
to a good approximation the average value of f(h) over the interval A N ( t ) .The
larger the kernel's width wN(t) ,the greater the smoothing action of TN,that is the less
will fN contain the rapid changes off. Thus TNmay be interpreted as the averager with
weighting function h,(t) . The pulse-like shape of the weighting function ensures
that the average of f is effectively taken over the interval A N( t ),while values of f(h)
outside this interval make negligible contribution to the average fN(t).
Note that TN represents a non-stationary system, that is a system with changing
properties such that if fN(t) = TN[f(h)] then the non-stationary characteristics imply
The non-stationary properties of the kernel of TNhave a great effect on the properties
of the truncation error. The exact nature of the non-stationary properties of the kernel
can be deduced from (A.7) where it is observed that the parameter t determines the
scale of the independent variable A. Hence the width of the kernel is directly
proportional to t, namely:
Thus for any finite N > 1, the acuity of the kernel is perfect at t = O+ and
progressively deteriorates with increasing t. Hence i f f undergoes changes of the same
magnitude in any interval of a given size, the relative acuity of the kernel will
decrease as t increases and therefore the truncation error will increase as t increases.
Zakian [29] estimated EN for particular situations.
Discussion and Conclusions
It has been shown in this work that Im techniques can be used to rapidly solve the
differential equations of large-scale systems which may be stiff as well as the solution
of oscillatory and steep responses. For the solution of the large-scale staged system, it
has been found that the step-by-step technique is faster than currently available
Matlab functions for solving such problems. Mathematica had to be used to compute
Rapid Solution of Stiff Diferential Equations using IMNApproxitnants
the transfer function of the large-scale system because Matlab gave inaccurate results.
This former software was also adapted to get the analytical expression of the time
response. For ordinary differential equations, we have shown that it is possible to
obtain accurate inversion of oscillatory functions by using either higher Im
approximants in a global manner or Im step-by-step techniques. Im global method
gives accurate results for steep responses either in ordinary or partial differential
equations. Through the concept of acuity and equivalent-width, explained above, and
the derivation of the Im formula given in the Appendix, we explain why the use of
high N or small h gives accurate IM results. In general, the use of higher
approximants facilitates the estimation of the error in the results. Another powerful
method of Laplace transform inversion is the method of Honig and Hirdes [14]. This
method, which unlike Im approximants is not applicable to unstable systems (having
poles in the open right half plane), has been applied to theexamples where it is
applicable yielding agreeable results. Although this method automatically gives a
rough estimate of the error, higher Im approximants have been exploited to
accurately estimate the errors for situations where there are no analytical Laplace
Table 1 indicates that the difference N-M- after being constant for a range of N
values increases in steps of one as N increases. When N =57, M- is 45, the
difference N-M,h being 12 in the range 57 5 N I 71. However, when N is increased
to 72, Mminis 59, with the difference N-Mh remaining at 13 in the range 72 SNI 88,
and so on. Another important observation is that the succeeding ranges of N values
for which N-Mm is constant increases with increasing N. For example, the range of
N values for which N-M- is 12 is 15 whereas the range of N values for which NMh is 25 is 59. These new observations are very useful in studying the properties of
Pad6 approximants, some of which are available (see Saff and Varga, 1975).
Table 2 provides the components of ai giving max(lct$ and the corresponding
components of Ki (which give min(lK, I)) as well as the components of a, giving
rnin(lai1) and the corresponding components of Ki (which give max(lKi I)) for various
values of M and N spanning the possible range of practical demand. This table shows
that both max(lail) and max(lKi I) increase with increasing N. For a given N both
max(lail) and max()Ki1) increase as M is increased, but the knowledge of max(lKi I),
for a given M and N gives a good idea of the order of magnitude involved as M is
varied with N fixed. In any case, based on long tern experience and fiom theoretical
consideration of ensuring good accuracy for short and long times for the global
method, M= Mh or close to M- is usually recommended. The importance of Table
2 is that it gives us an idea of how large the Im constants are so that one would know
beforehand how many digits to use in the computations to have a result of certain
number of digits of accuracy. For example, while MATLAB or double precision
arithmetic offered by FORTRAN would in general be good enough for computations
involving up to N = 16, it would be preferable to utilize extended or infinite precision
arithmetic for N > 16 in order to minimize rounding errors. For our case where
Mathematica has been used, knowing the range of values of the Im constants
facilitates the specification of the number of digits for an accurate determination of
the Im constants, and the subsequent numerical Laplace inversion. Consider for
example the need to use 1145,160 approximant, then at least 100 digits should be
0 Taiwo and R. King
specified in the computation of the Im constants and the subsequent numerical
Laplace inversion. The final number of digits specified will be determined by the
desired accuracy, the smallest value of time (t), and the largest value of F(ai/t), see
equation (2). Note that in using the Im method and for economy of computation, the
Im constants are computed once and stored for subsequent use.
0. Taiwo acknowledges the grant of the Alexander von Humboldt Fellowship,
1 Zakian, V and Al-Naib, U 1973 Design of dynamical and control systems by the method of
inequalities Proc IEE, 120, 1421-1427
2 Zakian, V 1979 New formulation for the method of inequalities Proc IEE, 126,579-584
3 Taiwo, 0 1980 Application of the method of inequalities to the multivariable control of binary
distillation columns Chem Eng Sci ,35,847-858
4 Zakian, V 1996 Perspectives on the principle of matching and the method of inequalities Int J
Contr~l,65, 147-175
5 Whidbome, J F and Liu, G P 1993 Critical Control Systems Theory, design and applications
Taunton Research Studies Press
6 Wu,Y T ,Zakian, V and Graves, D J 1976 Diffusion and reversible reaction in a sphere a numerical
study using Im approximants Chm Eng Sci, 31,153-162
7 Taiwo, 0 and King, R ,2000 Determination of kinetic parameters for the adsorption of a protein on
porous beads using symbolic computation and numerical Laplace inversion Submitted for
8 Edwards, M J 1977 Applications of M i a n ’ s 1~ and Jm approxitnants to the unsteady state solution
of the differential equations of aperiodically cycled plate column Chem Eng J , 13, 119-125
9 Atkinson, M P and Lang, S R 1972 A comparison of some inverse Laplace transform techniques for
use in circuit design Comp J , 15, 138-139
10 Taiwo, 0.Schultz, J and Krebs, V 1995 A comparison of two methods for the numerical inversion
of Laplace transforms Comput Chem Eng, 19,303-308
11 Kim, C and Friedly, 1 C 1974 Approximate dynamic modeling of large staged systems Ind Eng
Chem Proc Des Dev, 13, 177-181
12 Douglas, J M 1972 Process Dynamics and Control, Vol 1 & 2, Englewood Cliffs, New Jersey:
Prentice Hall
13 Fogler, H S 1999 Elements of Chemical Reaction Engineering Englewood Cliffs, New Jersey
Prentice Hall
14 Honig, G , and Hirdes, U 1984 A method for the numerical inversion of Laplace transforms, J
Comp Appl Math, 10,113-132
15 Zakian, V 1975 Properties of Im and J ~ ( N
approximants and applications to numerical inversion of
Laplace transforms and initial value problems J Math Anal Appl , 50,191-222
16 Zalrian, V 1975 Applications of Jm approximants to numerical initial value problems in linear
differential-algebraic systems J lnst Math Appl , 15,267-272
17 Zakian, V , and Edwards, M J 1978 Tabulation of constants for full grade IMN approximants Math
Comp, 32,519-531
18 Stoever, M A , and Georgakis, C 1982 Time domain order reduction of tridiagonal dynamics of staged
processes Chem Eng Sci ,37,687-698
19 Celebi, C , and Chimowitz, E H 1985 Analytic reducedader dynamic models for large equilibrium
staged cascades AIChE J , 3 1,2039-2051
20 Shimizu, K , Holt, B R , Morari, M , and Mah, R S H 1985 Assessment of control structllres for
binary distillation columns with secondary reflux and vaporization Ind Eng Chem Proc Des Dev ,
Rapid Solution of Stiff Differential Equations using IMNApproximants
Skogestad, S , Morari, M and Doyle, J C 1988 Robust control of ill-conditioned plants high-purity
distillation IEEE Trun Auro Control, 33, 1092-1105
Carmola, R E , and Chimowitz, E H 1990 Analysis of modal reduction techniques for the dynamics of
general tridmgonal systems Compur Chem Eng ,14 227-239
Prabhu, E S ,and Chidambaram, M 1991 Robust control of a distillation column by the method of
inequalities J Proc Control, 1, 171-176
Evans, D J , Shanehchi, J ,and Rick, C C 1985 A modified bisection algorithm for the determination
of the eigenvalues of a symmetric tridiagonal matrix Numer Math ,38,417-419
Cheney, W ,and Kincaid, D 1999 Numerical Mathematics and Computing Pacific Grove Brooks &
Cole Publishing, Co
Rice, R G I and Do, D D 1995 Applied Mathematics and Modeling for Chemical Engineers New
York John Wiley & Sons
Danckwerts, P V 1953 Continuous flow systems, Chern Eng Sci ,2. 1-9
Seidel-Morgenstem, A 1991 Analysis of boundary conditions in the axial dispersion model by
application of numerical Laplace inversion Chem Eng Sci , 46,2567-2571
Zakian, V 1969 Stiff equations solved by numerical inversion of Laplace transform,Report No 51,
Control Systems Centre, UMIST, Manchester, England
Saff, E B and Varga, R S 1975 On the zeros and poles of Pad6 approximants to e' Numer Murh ,
25, 1-14
Wolfram, S 1999 The Mathematica Book 4' ed Champaign & Cambridge, Wolfram Media and
Cambridge University Press
Received 24 November 2000; Accepted after revision 19 March 2001.
0 Taiwo and R. King
Appendix Derivation of the Iw Formula (2)
Using the property of the Dirac delta function, we have:
f ( t ) = 1 gj-(a)qA- ipa
Following Zakian [ 151, the function 6(- - 1) can be expressed in terms of the series:
6(- - 1) =
where the q and K, are defined by (3) for the case N =
If we denote the sequence 6, as the partial sum of 6 ,then:
which on substituting for 6, &om (A.3) into (A.4) yields:
which is the same as (2). For the purposes of the discussion in the main section of this
paper, we also define (AS) as:
Equation (A.6) represents a linear transformation TN,which converts f into fN,and
may therefore be written in the form
Без категории
Размер файла
1 129 Кб
using, rapid, imn, equations, numerical, approximants, stiff, oscillators, steel, solutions, differential, function, inversion, laplace, accurate
Пожаловаться на содержимое документа