In article <ggjWj.3412$J16.1931 at (no spam) newssvr23.news.prodigy.net>,
Monty Hall <monty at (no spam) hall.com> writes:
Hello,
When I assemble a system that I want to simulate I get a matrix equation
of the following form. My input vector is a list of function
derivatives. I'm solving for pipeline flow rates using a linear model -
something very similar to a electrical circuit. Would the matrix
equation be classified as a differential algebraic equation?
A x b
--------------- ------ -
a11 a12 a13 a14 m''(t) c1
a21 a22 a23 a24 * m'(t) = c2
a31 a32 a33 a34 m(t) c3
a41 a42 a43 a44 p c4
!!!!! this looks quite strange: maybe you intended here to write this in the
original meaning but had the system already transformaed in the first order
system form. but even then the system would look like
A(y)y' + B(y) = F(t)
y(t)=(y1(t),y2(t),p(t)) with y1(t)=m(t) and y2(t)=m'(t)
The above was written abbreviated. If I had 3 pipelines, my "x" vector
would be
m1''(t)
m2''(t)
m3''(t)
m1'(t)
m2'(t)
m3'(t)
m1(t)
m2(t)
m3(t)
...
Currently, I place the initial conditions into the "b" input vector -
which includes the initial/current value to the diff eq's - and multiply
by a LU decomposed "A" matrix to solve to the "x" vector. At that point
I accumulate the solutions - an euler integration in appearance, place
accumulated values into "b" vector and solve again for each time step.
"euler" performed at (no spam) each time step:
----------------------------------
m(t1) = m(t0) + m'(t1) placed into "b" for next iteration
m'(t1) = m'(t0) + m''(t1) placed into "b" for next iteration
m''(t1) = solution from "x"
obviously you forgot the factor delta_t on the right hand side
and took Euler implicit?
What are my choices for a numerical solution? Any examples on the web
that solves a similar system? It looks like a DAE, from what I can tell
from the web, but DASPK solver is a fortran program. I would a
it _is_ a DAE, hence your Euler solution would fail for sure :
a DAE can be considered as a stiff ODE with infinite stiffness,
not a good point to start with Euler
with DASSL or DASPK you are on the rihgt track
c/c++/java version w/ source if possible. I would very much like to
keep the above form w/ a minimal of symbolic rearrangement (avoid having
to hand craft ODE's from the above equation for each new piping network
that is created..) Any help appreaciated, TIA.
Regards,
Monty
concerning your wishes for C++:
this is from Hairer' homepage
There is a folder, written by Blake Ashby "bmashby (at) stanford.edu",
which contains C++ versions of the nonstiff integrator DOPRI5
and of the stiff integrator RADAU5.
It can be unfolded with "tar xzf IntegratorT.tgz".
radau5 allows the implicit form with singular mass matrix (your case)
Hairer has also specail codes for constrained mechanical systems
look here:
http://www.unige.ch/~hairer/software.html
hth
peter
