"Unstable" solution of a system of two coupled first order ordinary differential equation

Question

Sorry for the long post. Maybe some of you will read the first few lines and already figure out what's going on, but in case it's not the case, I wanted to give as many details as possible. I came across a strange result while solving a system of two coupled first order ODE:

eqns = {
  vl pl'[t] == -sl (pl[t]-pc[t])
  vc pc'[t] == -sc pc[t] + vl pl'[t] + lback
  pl[0]==atm,
  pc[0]==atm
 }

For your curiosity, pc and pl represent the pressures in two volumes vc and vl; vc is pumped at a constant "volumetric" speed (liters/seconds) sc, while vl leaks into vc at a constant volumentric speed sl. There is also a constant contribution lback to the leak rate into vc.

The strange issue is that if I just ask Mathematica to solve the system (including initial conditions):

DSolve[eqns, {pc[t], pl[t]}, t]

and then plot the solution, it exhibits unexpected (and clearly wrong) behavior after a given time, for some (many) values of the parameters. Both pc and pl should just decrease almost exponentially until a constant value. Instead, after a certain time they seem to grow back and start to oscillate wildly.

After many many trials, I found out that solving the system through other methods does not produce the issue. For example: - not providing initial conditions to DSolve, and later using Solve to calculate the value of the constants. - "manually" obtaining a 2nd order ODE for one of the functions, solving (using DSolve, even with initial conditions) and then substituting back to find the other function

I also noticed the following. Suppose pcA is the solution for pc[t] obtained with the very first method (the bad one), and pcB is the solution obtained with one of the other two methods. pcA-pcB //Simplify gives 0. However, LogPlot[Evaluate[{pcA, pcB} /. vals], {t, 0, 1000}] (where vals is a list of rules to replace symbols with their numerical values) produces a plot in which the firs solution is weird (as described before) and the second looks fine.

Finally, LogPlot[Evaluate[{pcA, pcB} /.lback->0 /. vals], {t, 0, 1000}] also produces a plot in which the two solutions coincide and look fine, although obviously not what I want (because in my system lback is not 0).

Any idea what's going on?

----UPDATE----

After finding some inconsistency while trying to replicate my own examples, I finally figure out that the real difference between the weird solutions and the good one is that the good ones have been, directly or indirectly, "simplified" (using Simplify). In short: - if solving with initial conditions in DSolve, the solution is correct (i.e. not weird) is I issue the Simplify command on it before plotting - same applies for the case in which I solve without initial conditions, and calculate them later. If I don't Simplify (it turns out I was when I produced my example), I still get the same weird solution. - when solving going through the 2nd order ODE, Simplify doesn't seem to be necessary; probably I'm already forcing Mathematica to express the solution in a different form wrt the previous cases

I suppose this sheds quite a bit of light on the issue, but I didn't post this as an answer since it still doesn't clarify (to me) why the same solution should e plotted in different ways before and after being simplified... I suppose it will all boild down to a numerical issue, but it's a bit worrying. If I didn't know pretty well what to expect from the system, I could have concluded that the physical system would, indeed, become unstable after some time...

----UPDATE 2---- As requested, here is a complete sample code. Note that I had a sign wrong in my original post, but that doesn't make a difference.

vals = {
   atm -> 1000,
   sc -> 20,
   lback -> 10^-5,
   vc -> 200,
   sl -> 10^-8,
   vl -> 10^-6
   };
eqns = {
   vl pl'[t] == -sl (pl[t] - pc[t]),
   vc pc'[t] == -sc pc[t] - vl pl'[t] + lback,
   pl[0] == atm,
   pc[0] == atm
   };

solAuto = DSolve[eqns, {pc[t], pl[t]}, t];

LogPlot[Evaluate[{pc[t], pl[t]} /. solAuto /. vals], {t, 0, 1000}]

and here is the output

but if I issue the simplify command on the solution

LogPlot[Evaluate[{pc[t], pl[t]} /. solAuto /. vals //Simplify], {t, 0, 1000}]

everything works fine...

Answer 1

Perhaps the simplest way to demonstrate Dr Hintze's observation about numerical instability is to examine the results of machine precision and arbitrary precision computations of the value of one of the functions. Machine precision corresponds to just under 16 digits of precision on my machine (or $MachinePrecision). (See Precision and Accuracy, as well as this tutorial.)

The following shows that at t -> 700, the computation loses more than 14 digits of precision. Since the computed machine-precision value is much less than the error, machine calculations in this neighborhood are useless.

pc[t] /. solAuto /. vals /. t -> 700.`
pc[t] /. solAuto /. vals /. t -> 700.`16
(*
  {-0.0008138084479616342`} 
  {0``1.902949081049885} 
*)

However, the simplified expression loses just over 2 digits in the arbitrary-precision calculation; the machine calculation loses a bit more, about 6, because arbitrary-precision calculations carry more guard digits. (Improving numerics is not the purpose of Simplify, only a coincidental side-effect in this case.)

expr = pc[t] /. solAuto /. vals // Simplify;
expr /. t -> 700.`
expr /. t -> 700.`16
(*
  {5.005628900673123`*^-7} 
  {5.005628901036847621902904972908550439687571275356483`13.852169564967438*^-7} 
*)

In the neighborhood of t -> 1000, the loss of precision is about 26 digits (16 - (-9.9)).

pc[t] /. solAuto /. vals /. t -> 1000.`
pc[t] /. solAuto /. vals /. t -> 1000.`16
(*
  {-2.5587709152932596`*^8} 
  {0``-9.918381068297375} 
*)

With such large errors, a plot should be expected to show large oscillations. A log plot might contain blank areas due to computed negative values of the functions, even though theoretically the exact values of the functions are positive.

Answer 2

(Please note that because of an error in my previous versions the answer was completely revised.)

Here is an analytical treatment of the problem which shows that the exact solutions can be stable or unstable, depending on the parameters.

The ODEs can be written in matrix form.

Let

p[t_] = {pl[t], pc[t]};
p0 = atm {1, 1};
q = {0, lback/vc};
m = {{-sl/vl, +sl/vl}, { -sl/vc, (sl - sc)/vc}}; (* corrected *)

Then the system of ODEs reads

Thread[p'[t] == m.p[t] + q && p[0] == p0]

(*
Out[458]= {Derivative[1][pl][t], 
   Derivative[1][pc][t]} == {(sl pc[t])/vl - (sl pl[t])/vl, 
   lback/vc + ((-sc + sl) pc[t])/vc - (sl pl[t])/vc} && {pl[0], 
   pc[0]} == {atm, atm}
*)

The solution is given by

p[t] -> - Inverse[m].q + MatrixExp[t m].(p0 + Inverse[m].q ); 

The character of the solutions is determined by the eigenvalues of the matrix m. These are given by

ev = Eigenvalues[m]

(*
Out[459]= {
(-sl vc - sc vl + sl vl 
    - Sqrt[-4 sc sl vc vl + (sl vc + sc vl - sl vl)^2])/(2 vc vl), 
(-sl vc - sc vl + sl vl + 
    + Sqrt[-4 sc sl vc vl + (sl vc + sc vl - sl vl)^2])/(2 vc vl)}
*)

It is not difficult to show that the eigenvalues can be classified into regions of a two dimensional diagram with vc/vl as the x-axis and sc/sl as the y-axis.

The region of unstable solutions is given by

cndun = vc/vl + sc/sl < 1;

The region of complex eigenvalues is given by

cndc = (1 - Sqrt[vc/vl])^2 < sc /sl < (1 + Sqrt[vc/vl])^2;

For each pair of the parameter quotients vc/vl and sc/sl we can deduct the type of solution in the property space (stable/ustable) x (monotonous/oscillating).

An example of an unstable oscillating solution is obtained for the parameters

{vc -> 9/2, vl -> 7, sc -> 9/5, sl -> 8}

If the real parts of the eivenvalues are negative the solutions are stable, and because they behave as Exp[ t Eigenvalue], they decay to the stable stationary value given by

ps = -Inverse[m].q

(*
Out[192]= {lback/sc, lback/sc}
*)

Observations

1) There are stable and unstable "physical" solutions which can be monotonous or oscillating depending on the paramters vc, vl, sc, sl

2) The "physical" time scales of the problem are given by the inverse real part of the eigenvalues, the frequencies by their imaginary parts. Times t far beyond the physical time scales are "physically" meaningless and should be considered artificial.

EDIT #1 (as of 23.09.15)

This following code let's us experiment with different parameters. It is a new improved version comprising the recent developments.

The code consists of the following parts

1) make a random choice of the parameters 2) tell the resulting type of solution 3) if the user wishes, solves the ODEs and plots the functions 4) plotting can be repeated with the same solution for differents time scales

The equations are

eqns = {
   vl pl'[t] == -sl (pl[t] - pc[t]), 
   vc pc'[t] == -sc pc[t] + vl pl'[t] + lback,
   pl[0] == atm,
   pc[0] == atm};

Classification of solutions according to the eigenvalues

stable = vc/vl + sc/sl > 1;
complex = (1 - Sqrt[vc/vl])^2 < sc /sl < (1 + Sqrt[vc/vl])^2;

We shall set the parameters to some random values chosen from

r := Rationalize[1 + 10*Random[], 1/10] (* avoiding zero *)

rr6 := Array[r &, 6]  (* all 6 parameters *)

(* ----------- entry to next choice of parameters *)
rr6s = rr6; (* store parameter choice in rr6s *)
rep = Thread[{sc, vc, sl, vl, atm, lback} -> rr6s]
Print[If[stable /. rep, "stable", "unstable"], ", ", 
 If[complex /. rep, "oscillating", "monotonous"]]

(*
Out[236]= {sc -> 3/2, vc -> 5, sl -> 43/4, vl -> 64/7, atm -> 29/3, lback -> 7/4}
*)

unstable, oscillating

(* ----------- entry to next plot *)
sol = DSolve[eqns /. rep, {pc[t], pl[t]}, t] // FullSimplify;
pp = {pc[t], pl[t]} /. sol
Plot[pp, {t, 0, 20}, PlotRange -> All]

(*
Out[239]= {{7/6 + (17 E^(
    863 t/2560) (1566911 Cos[(Sqrt[1566911] t)/2560] - 
      1631 Sqrt[1566911] Sin[(Sqrt[1566911] t)/2560]))/3133822, 
  7/6 + (17 E^(
    863 t/2560) (1566911 Cos[(Sqrt[1566911] t)/2560] - 
      863 Sqrt[1566911] Sin[(Sqrt[1566911] t)/2560]))/3133822}}
*)

Answer 3

!! EDITED !!

Dr. Wolfgang Hintze wrote: What you observe is a numerical instability due to limited accuracy. Try using only rational numbers for all parameters. The numerical instability should vanish.

To fix it , you can use a Rationalize[ ] function.

  Rationalize[{vc = 200, vl = 10^-6, sc = 20, sl = 10^-8, lback = 10^-5,
  atm = 1000}, 0];

  eqns = {vl*pl'[t] == -sl*(pl[t] - pc[t]), 
  vc*pc'[t] == -sc*pc[t] + vl*pl'[t] + lback, pl[0] == atm, 
  pc[0] == atm};

  sol = DSolve[eqns, {pc[t], pl[t]}, t];

  LogPlot[Evaluate[{pc[t], pl[t]} /. sol], {t, 0, 1000}, 
  PlotLegends -> {"pc[t]", "pl[t]"}, AxesLabel -> {t, {pc[t], pl[t]}}, 
  PlotRange -> All]

Solution with Maple: