# 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. Introduction 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 143 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. 144 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: N IMN(f,t) = t " C K , L ( f , n l ) 1 4 t 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): eMN = min,{Re(a,)}.It is - f(t) Im(f,t) = O(PN+') ,t + 0' whenever, ,& (4) > 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 >0 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 0 arefull grade, that is, satisfjr (3) and di, (N N-MN N-M- approximants 57-71 72-88 89-107 108-129 130-155 156-183 184-214 215-249 12 13 14 15 16 17 18 19 250-287 288-329 330-376 3 77-426 427-48 1 482-540 20 21 22 23 24 25 145 1 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 (5) 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. NI2 Im(f,t) = 2t-’ CRe{K,L(f,S)) i=l t 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 dx - = A x ( t ) , x(0) = x, dt (7) where A is an n x n constant matrix, x(t) is an n-vector and ~0 is the given initial condition. 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 146 (9) 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 (134 N x,~= h - ' x K , X , , r r=0,1,2 ,... ,=I 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: NI2 x,l= 2 h-' I Re(K,X,,,) r = 0, 1,2,... 4 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. 147 (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 3Q B 2 0 ~ 1456,480 1452,476 1450,464 1404,428 1354,376 I 1268,288 1223,238 479 1 475 I 329 I 1 1 375 I 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 187 , 0.0018196715 ~ ~ ~ + 626.2887162478i 287 4.675746361 + 450.6290834331 30 1.025979203 + 1.6901349846i ~ 1 ~ I I I 237 Table 2 continued. I I I -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 1 1 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: dr - = m(t) x(0) = x, dt 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 stable 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: c dr - = A x @ )+ ~ ( t ) ,~ ( 0 =) X, dt 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: 150 Rapid Solution of StiffDifferential Equations using IMNApproximanis lik - = Ax(t) + bu(t), y(r) = cx(t), x(0) dt = x0 where b = [ l , 0, 0,.... 0IT,c = [O,O, 0,...0,I], -2 1 1 -2 0 1 0 0 0 0 0 0 A= ... 0 ... 0 0 0 0 0 ... 0 . . . . 0 ... - 2 0 ... 1 0 ... 0 0 0 1 0 -2 1 -2 0 1 -2 1 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 151 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 Method and 15.6 13.4 (a) I3,4 I3,4 03) Odel5s Ode23T Ode23TB Ode23s RK4TH Maximum percentage Error 1.212 1.196 2.386 0.813 4.679 3.723 5.218 0.927 Average CPU Time@) 0.43 1 0.388 0.184 0.474 0.714 0.830 7.052 33.372 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. 152 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 (2.471~10-~) 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 Method 13.4 Ode 15s Ode23T Ode23TB Ode23s RK4TH Maximum percentage Error 1.398 0.944 4.723 3.799 5.718 1.42 1 Average CPU Time($ 1.45 1.48 2.15 2.70 30.38 158.24 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. 153 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. Method Maximum Error, E,, Location of E m , X E , , ~ 114.22 3 . 9 8 106 ~ 20.5 I19,28 7.58~10-‘~ 20.5 Honig & Hirdes [ 141, automatic 1.05~10-~’ 5.5 setting of ftee parameters IS4 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) - w) 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 Method I46,58 152.64 Maximum Error, E., 7.7 1x 1 Q 8 4 . 9 7 1~0-'7 Location of Em20.5 20.5 XE,, 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 1.784723~ 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) + A(A,t). t 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,,)) where: 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 A. means larger N, thus ensuring (- - 1) A(& t) ,so that as N increases the sifting t ability of the method is near perfect. - 156 Rapid Solution of Stiff Differential Equations using IMN Approximants 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 ~~~~ Method Pe = 20 14.10 19916 114~22 125~34 Honig & Hirdes, automatic setting of fiee parameters 410 Pe =lo0 I9,16 114~22 130~40 Honig & Hirdes, automatic setting of fiee parameters Pe =500 114.22 Maximum Error 2.79~ 1 O4 I. 14x lo4 3.90~ 2.36~10”~ 1.90~ I O-I2 ~ Location of Em, 1.4 1.8 1.8 1.8 1.8 3.63x 1O5 4.04~ 1 O4 9.31 x lo4 2.17~ lo-” 1 . O ~ 1X0-13 1.4 1.4 1.4 1.4 1.3 5.1 8 x 10” I .04 Pe= 10000 157 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 t 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: 158 Rapid Solution of Stiff Diperential Equations using IMN Approximants N 13= Z K , la, r=l which shows that the area of the kernel is a constant. Furthermore, for any value of N, N K , /a,= 1, thus ensuring that the constants may be chosen to satisfy the condition r=l 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 ) (21) 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: n,,,(t) = {a I t - w , ( t ) / 2 c a < t + W , 12 1 (22) 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. 159 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 that: 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 160 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 inversions. 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 161 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. Acknowledgement 0. Taiwo acknowledges the grant of the Alexander von Humboldt Fellowship, Germany. References 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 publication 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 , 24.852-858 162 Rapid Solution of Stiff Differential Equations using IMNApproximants 30 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 , 31 25, 1-14 Wolfram, S 1999 The Mathematica Book 4' ed Champaign & Cambridge, Wolfram Media and 21 22 23 24 25 26 27 28 29 Cambridge University Press Received 24 November 2000; Accepted after revision 19 March 2001. 163 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 t (A.1) I Following Zakian [ 151, the function 6(- - 1) can be expressed in terms of the series: t A 6(- - 1) = t m CK~exp(-ai r=l I -) t 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: Q and 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 164

1/--страниц