Re: Best Numerical EDO method for this problem?


In article <763b6b9c-b170-43db-823f-686336fa7afc@xxxxxxxxxxxxxxxxxxxxxxxxxxx>,
Ketakopter@xxxxxxxxx writes:
On 24 mayo, 11:16, Ketakop...@xxxxxxxxx wrote:
Hi all,

I have to solve numerically the following system:

dm
-- =3D f1(m,T)
dt

dT
-- =3D f2(m,T)
dt

where f1 and f2 are:

f1(m,T) =3D M(in) - M(out)
=A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 T =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =
=A0mRT dV
f2(m,T) =3D A*{ - - ( M(in) - M(out) ) - --- -- + B*M(in) - C*T*M(out) }
=A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 m =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =
=A0 V =A0dt

and M(in) and M(out) are in and out fluxes, which are also functions
of m and T:

M(in) =A0=3D A1 * sqrt[D*(E - (mRT/V))] * Phi(EV/mRT)
M(out) =3D A2 * sqrt[2(m/V)*(mRT/V - F)] * Phi(mRT/FV)

The Phi function is a compressibility correction function and is quite
complicated, but behaves almost as linear. All A,A1,A2,B... are
constants, and V is a known function of time. Please forgive me if you
don't properly see the equations.

The solution is straightforward I believe, as everything is well
defined. Now, I'd like to know which method to use so as to optimize
complexity and accurateness. Which method would you recommend?

Thanks in advance,

keta :)

Uh, I meant ODE on the title :P Also, just wanted to note that A1 and
A2 are not constants, but rather functions of time, and not
differentiable in some points (they're valve areas, which open and
close instantly). Phi is also not differentiable.

I'm going for a Runge-Kutta method of order 4. Is the Butcher array
the same for a system than for a scalar problem? If not, what are the
best coefficients for this system?

yes , Runge-Kutta for systems works word by word as in the scalar case.
but I would discourage the use of the "classical" scheme. better use a scheme
with build in error control like dormand-prince or shampine-lawrence.
(order 4 and 5, 6 stages).
then it is possible to use variable stepsizes dependent on the
prescribed local error.
due to the discontinities you must integrate stepwise, nether integrate
through a jump discontinuity (this causes large errors at least locally)
if you cannot locate the jump discontinuities beforehand you must
use a numerical scheme for locating them during the integration process.
the integrators in the nag library and matlab's ode-suite have such
features: you must be able to characterize the occurence of the discontinuity
by a smooth function which changes sign there.
and whether the use of a nonstiff integrator works at all depends on the
properties of your system which cannot be judged from the information
above. best you try it with an explicit scheme first and if you detect
that the adaptive stepsize scheme chooses very small steps, then try an
implicit scheme.
you find adaptive integrators in http://www.netlib.org/ode
hth
peter

.