Ind. Eng. Chem. Res. 1993,32, 1798-1799
1798
W h y Standard Integration Packages Fail on a Near Index Problem Arthur W. Westerberg,'*+ Karl M. Westerberg3 Kenneth W. Westerberg,$ Yonsoo Chung,l and Boyd Safritt Department of Chemical Engineering and the Engineering Design Research Center, Carnegie Mellon University, Pittsburgh, Pennsylvania 15213, Department of Physics, Princeton University, Princeton, New Jersey, Lawrence Livermore National Laboratory, Livermore, California, and Korea Institute of Science and Technology, Seoul, Korea
We examine the analytical behavior of a "near index" problem which Chung and Westerberg used as an example to illustrate their newly proposed algorithm for solving such problems. The example problem has both a stiff and an unstable component, the latter causing both implicit (in principle) and explicit (in practice because of round-off error) numerical methods to diverge when attempting to track its known finite solution. We explain why the numerical method of C h u g and Westerberg, which corresponds to a perturbation expansion of the problem around the index problem to which it is close, has no difficulty in finding this solution. The Example Problem
differential equation:
In a recent paper two of the authors (Chung and Westerberg, 1992)presented a method to solve near index problems which can be discovered by the existence of small pivots during a Newton step in the solution algorithm. The method is the numerical equivalent of a perturbation expansion for solving such a problem. The paper illustrated the ideas with the following deceptively simple looking system of three differential-algebraic equations, which we shall call problem P1: x,' = x2
d2 1
which, because it is linear, we can solve by summing any particular solution with the homogeneous solution. 1'
= ~1,partieulm+ '1,homcgeneous
We find the homogeneous part by solving eq 1with a zero right-hand side. xl,homogeneoua
x i =y
$) = k , c o s h ( t ) + k2sinh(3) 1'
exp(
$)+
exp( -
'2
e1/2
x1 = ey + sin@)
When the value of the small parameter e is exactly zero, this problem cannot be solved by conventional numerical methods as it is an index 3 problem. They defined this problem to be a near index problem when t is small. This earlier paper argued that, when e became small, the problem becomes very stiff. In fact, the situation is more complex. An analysis of the problem shows it to have eigenvalues of fl/e1I2. For a positive e, one component in this problem does become stiff as e becomes very small; however, the other is always unstable, and it is the unstable part of the solution that caused LSODE (Hindmarsh, 1980) to fail to solve this problem without issuing a warning message. A second example in that paper involved tracking the levels of two tanks joined by a horizontal pipe with a valve in it. As the resistance of the valve approached zero, the tank levels became more and more coupled and began to act as a single state for the system. Here the problem simply became more and more stiff. LSODE failed when the resistance was extremely small, declaring that the problem was singular. In this case it does issue a warning message. The method in Chung and Westerberg stably computed a solution for both problems. As noted in the earlier paper, problem P1 has an analytical solution. One approach to find it is as follows. First eliminate y and x2 by using the second and third equations, getting the following second-order ordinary + Carnegie Mellon University. t Princeton University. 5 Lawrence Livermore National
(1)
Requiring XI to be real requires that c1 and c2 are real for e > 0 and that c1 is the complex conjugate of c2 for t < 0. Alternatively, we require that k l and k2/e1j2are real for e # 0. For e = 0,x1,homOgen~u. 0. The homogeneous solution has both a stable and an unstable component if e is positive. Both become very "fast" as e approacheszero. Thus, for small positive e, the homogeneous solution has a very stiff component and a very unstable component. If e is negative, the solution is oscillatory with the frequency increasing to infinity as e approaches zero. A particular solution can be found by solvingthe original model using a perturbation expansion as done in Chung and Westerberg. Rearranging eq 1,we write x , = sin(t)
+ exl"
which we can solve recursively to get x1 = (1- e
=-
+ e2 - e3 + ... + (-e)*')
1- (-e)" sin(t) l+€
sin(t) + enx1(2n)
+ cnxl(2n)
where ~ 1 (=~d2nxl/dt2n. ~ ) Assuming le1 < 1 and that the last term, e " z ~ ( ~tends ~), to zero as n gives us the following particular solution (which, when checked, does satisfy the original equations defining our model):
-
x l , p a r t i= ~
1 l+e sinW
The general solution for x1 is then the sum of the particular and homogeneous solutions.
Laboratory.
Korea Institute of Science and Technology. ,
0 1993 American Chemical Society
Ind. Eng. Chem. Res., Vol. 32, No. 8,1993 1799
+ Using the first and second equations for the originalmodel (Pl),we can also write the following as the corresponding solutions for x2 and y.
Upon checking,we fiid that these equationsare the general solution for our original problem provided that e # 0, -1, with the first term in each case being the particular solution. For e = 0, the solution is the particular solution only. For e = -1, the particular solution is x1 = -t COS(t)/2. Our interest is for \el
