From: spellucci@mathematik.tu-darmstadt.de (Peter Spellucci) Subject: Re: System of 2 differential equations and 1 algebraic equation Date: 25 Feb 1999 15:22:26 GMT Newsgroups: sci.math.num-analysis In article <36D52AAF.F3EAF8F5@elis.rug.ac.be>, Bart Desoete writes: |> Hello, |> |> I am trying to solve the following system of 2 differential equations |> and 1 algebraic equation |> for x(t), y(t) and u(t) in the interval [0,T]: |> |> dx/dt = I(u) - I(x) |> dy/dt = (y-u) * I'(x) |> (y-u) * I'(u) = I(u) - I(x) |> |> where I(x) = (ln(1+exp(x/2)))^2 and I'(x) is the derivative of I(x) with |> respect to x. Both I(x) |> and I'(x) are monotonously increasing functions of x. There are two |> boundary conditions: |> x(0)=x_0 and x(T)=x_f, where T, x_0 and x_f have some specific values. |> |> I have tried to solve this problem numerically using Maple. It doesn't |> seem to be trivial for |> Maple to solve this system of equations simultaneously, even if I use |> only initial value |> conditions. A separation of the 2 differential equations and the |> algebraic equation, such that |> during the integration of the differential equations for each value of x |> and y the corresponding |> value of u is calculated, also doesn't seem to work. whether or not such a DAE is difficult to solve or not depends on the so called "index" of the system. somewhat crude this may be defined as the number of diffentiations of the given system needed to transform it into an explicit DE. If the derivative with respect to u of the third equation, i.e. (y-u)*I"(u)-2I'(u) is always nonzero , then you would have a system of index 1 and for such systems your imagined solution (transformation into an initial value problem, solving for u during integration) should work. There is however the trouble that the shooting method itself may fail if the interval [0,T] is not very small. for the solution of an initial value problem for a DAE of index 1 DASSL from netlib/slatec (http://www.netlib.org) is an excellent solution. however DASSL solves not the boundary value problem. you could try exolicit differentiation of the third equation with respect to t and solving the system by multiple shooting e.g. using codes like bvpsol of Deuflhard (http://elib.zib.de ) or colnew of ascher and bader (http://www.netlib.org in the library ode). then you need to add the algebraic equation for the left endpoint as an additional boundary condition (hence you have 3 such). there is the problem that you must provide a reasonable initial estimate for the solution. not an easy task good luck peter ============================================================================== From: trhoffend@mmm.com (Tom Hoffend) Subject: Re: Differntial Algebraic Equations DAEs Date: 10 Mar 1999 22:18:53 GMT Newsgroups: sci.math.num-analysis Toshiro K. Ohsumi (ohsumit@coffeepot.cs.rpi.edu) wrote: : Perhaps the best introductory source on solving DAEs (and ODEs) is the recent : book by Uri M. Ascher and Linda R. Petzold, "Computer Methods for Ordinary : Differential Equations and Differential-Algebraic Equations", published by : SIAM (1998). I hope this helps. Good luck, To answer the specific question about Runge-Kutta methods, there is also `The Numerical Solution of Differential-Algebraic Systems by Runge-Kutta Methods' by Ernst Hairer, Christian Lubich, and Michel Roche, Springer-Verlag Lecture Notes in Mathematics number 1409 (1989). I haven't looked at the newer Siam book yet, so I don't know how much (if any) overlap there is. There are some neat examples in the back of the book by Hairer et. al. TRH -- Thomas R. Hoffend Jr., Ph.D. EMAIL: trhoffend@mmm.com 3M Company 3M Center Bldg. 236-2A-06 My opinions are my own and not St. Paul, MN 55144-1000 those of 3M Company. Opinions expressed herein are my own and may not represent those of my employer. ============================================================================== From: rusin@vesuvius.math.niu.edu (Dave Rusin) Subject: Re: System of 2 differential equations and 1 algebraic equation Date: 11 Mar 1999 08:51:57 GMT Newsgroups: sci.math.num-analysis I got around to reading this article just before it rolled off my system. I hope this is not repeating what has been said already. Bart Desoete wrote: > I am trying to solve the following system of 2 differential equations > and 1 algebraic equation > for x(t), y(t) and u(t) in the interval [0,T]: > > dx/dt = I(u) - I(x) > dy/dt = (y-u) * I'(x) > (y-u) * I'(u) = I(u) - I(x) > > where I(x) = (ln(1+exp(x/2)))^2 and I'(x) is the derivative of I(x) with > respect to x. Both I(x) > and I'(x) are monotonously increasing functions of x. There are two > boundary conditions: > x(0)=x_0 and x(T)=x_f, where T, x_0 and x_f have some specific values. You probably want to solve for x, y, and u as functions of t, but you can get better results by thinking globally. Let's see, you have three equations in 4 unknowns. The solution set ought to be a curve, or more precisely (since you have two derivatives) a 2-parameter family of curves. The two boundary conditions will then pin down the particular curve winding through xyut-space. Near most points on the curve, any of these variables can be thought of as a function of any of the others (e.g., you'd like to think of x as a function of t .) Your differential equations are autonomous, that is, t is only implicitly mentioned. It follows that you can avoid references to t, giving just two equations in the other three unknowns, which describe the projection of the curve to xyu-space. Indeed, dividing the first equation by the second and substituting in the third gives the equations dy/dx = I'(x)/I'(u) (y-u)*I'(u) = I(u)-I(x). Here the presence of the "dy/dx" suggests viewing this as a system to be solved for two functions y(x) and u(x) . You can now solve the last equation for y, of course, expressing it in terms of x and u(x) . Differentiating this expression for y gives dy/dx = (\partial y/\partial x) + (\partial y/\partial u)(du/dx) so we can equate to the other expression for dy/dx and get a single differential equation relating x and u. That is, we have described the projection of the curve from xyu-space to the xu-plane, viewing the result as the graph of a function u(x) of x, say. Surprisingly, this last differential equation admits a symbolic solution. To obtain it, I first used the transformations x=2*ln(1/X-1), u=2*ln(1/U-1) (in which we want X, U in (0,1). ); then I(x)=(ln(X))^2 and I'(x)=(ln(X))(X-1). It turns out to work best to take this variable U as the parameter marking points on the curve. Then u is given in terms of U as above, and x and y are easily obtained once we solve the differential equation, which I compute to be dX/dU = 2 2 3 2 2 2 1/4 X (3 ln(U) U - 3 ln(U) - ln(U) U - ln(X) + ln(X) U + ln(U) U ln(X) )/ (ln(X) ln(U) U (-1 + U)). The substitution X=exp(z) leaves no transcendental terms in z, and if we multiply by 2z we can substitute z^2 = z2 to get a _linear_ differential equation in z2 . By some miracle, a symbolic solution can be obtained without much effort, although as it turns out Maple V5R4 makes a bad mistake in one of them (see a related thread in sci.math.symbolic). The general solution to the linear DE turns out to be 2 1/2 z2 = ln(U) + C (ln(U) (U - 1)) from which we recover 2 1/2 1/2 X = exp((ln(U) + C (ln(U) (U - 1)) ) ) and thus obtain explicit formulas for x, y, and u in terms of the curve parameter U and the undetermined constant C. [Note there is no branch-of-sqrt problem for U in (0,1) . ] So there is one parameter U from which we may obtain X, and then u,x,and y. This gives an explicit parameterization of each curve satisfying the original equations. If you want to follow along on your own, use this Maple expression for the parameterized point [x,y,u]: [2*ln(exp((ln(U)^2+C*(ln(U)*(U-1))^(1/2))^(1/2))-1), 2*ln(1/U-1)-C/(ln(U)*(U-1))^(1/2)), 2*ln(1/U-1)]; The curve quickly reaches the (positive) asymptote x=y=u as U -> 0. As U->1, the points become essentially [-2a, -C exp(a), -a] with increasing a. Now, the parameter U varies with t (and so then x, y, and u do, too). The initial conditions give 2 equations involving three unknowns U(0), U(T), and C. You can determine a third equation, and thus determine the value of C, from one of the original differential equations: for example, dx/dt = I(u) - I(x) may be rewritten dt = [ (dx/dU) / ( I(u) - I(x) ) ] dU, which we may integrate from t=0 to t = T, i.e., from U=U(0) to U=U(T). This leads to the equation T = (-1/C)* int( (dx/dU) / (ln(U)*(U - 1))^(1/2) ), U=U(0) .. U=U(T) ) involving the three unknowns U(0), U(T), and C. Unfortunately, this integrand (dx/dU)/(I(x)-I(u)) is a terrible function of U, and unlikely to have a closed-form antiderivative. So what seems easier is to simply guess C, compute U(0) and U(T) from the boundary conditions, then integrate and see how close we come to T. Guess a new C, recompute U(0) and U(1), and then you should be ready for something like the bisection method, say, to determine C. For example, I looked to see how to get a curve with x(0)=2, x(1)=5. Starting with C=1, just a few iterations eventually suggest C=6.5 is a pretty good fit; for this C, the parameterization [x,y,u], above, comes close to x=2 and x=5 when U=.76230 and U=.33988, respectively. We may then compute an integral from U=.76230 to any other U in (0,1) to compute t as a function of U. This lets us compute a table of points [x,y,u,t] on the original curve, and we find that t=0 when x=2 and t=1 (almost) when x=5. So the final result is neither particularly efficient nor elegant. But it has the virtue that the long-term trends in the curve can be described precisely, subject only to the initial determination of C -- there are not, for example, any propogation errors resulting from solving the initial differential system using Euler's method or whatever. dave