2
$\begingroup$

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

enter image description here

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...

enter image description here

Improve this question
$\endgroup$
6
  • 1
    $\begingroup$ Possibly related: mathematica.stackexchange.com/q/7109 -- But how can you expect accurate help? Problems with code usually require (complete) code that reproduces the problem. $\endgroup$ Commented Sep 17, 2015 at 1:10
  • $\begingroup$ Fair enough. Updated my question... $\endgroup$ Commented Sep 18, 2015 at 4:06
  • 3
    $\begingroup$ @ Giacomo Ciani: This problem is linear and admits an analytic solution in terms of trigonometric and exponential functions. For non negative paramters vl..sc the solutions are stable. What you observe is a numerical instability due to limited accuracy. Try using only rational numbers for all parameters. The numerical instability should vanish. $\endgroup$ Commented Sep 18, 2015 at 6:45
  • $\begingroup$ Note for closers: I voted to leave the question open since there are 3 good answers and it may be useful for reference later on. But I am not committed to this. $\endgroup$ Commented Sep 20, 2015 at 20:12
  • $\begingroup$ @ Dr Hintze: I knew the problem admitted analytic solutions. I was just wondering what was going on with the plot... I think you hinted at it (although my parameters were rational numbers!), and Micheal E2 gave the final answer. $\endgroup$ Commented Sep 20, 2015 at 20:51

3 Answers 3

Reset to default
3
$\begingroup$

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.

Improve this answer
$\endgroup$
2
  • $\begingroup$ Thanks for the thorough explanation. My complain at this point (not to you) is: if you can easily check when limited precision produces completely inaccurate results, why can't Mathematica do it itself, or at least throw a warning? $\endgroup$ Commented Sep 20, 2015 at 20:54
  • $\begingroup$ @GiacomoCiani MachinePrecision is for fast calculation without such checking. For a tested algorithm in which the error loss is understood, one wants to be able to calculate as fast as possible. Continually keeping track of error propagation costs time, which many people would find objectionable. Arbitrary precision does keep track, although I think you might have to explicitly check (with Precision) if the precision loss is significant. -- And thanks for the accept. $\endgroup$ Commented Sep 20, 2015 at 21:01
1
$\begingroup$

(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.

enter image description here

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}} *)

enter image description here

Improve this answer
$\endgroup$
0
$\begingroup$

!! 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]

enter image description here

Solution with Maple:

enter image description here

Improve this answer
$\endgroup$
6
  • $\begingroup$ With your values, my code (the one I just added to my original answer) give a weird behavior after t ~ 800. It does not oscillate, but pl starts growing exponentially while pc drops sharply towards 0... $\endgroup$ Commented Sep 18, 2015 at 4:10
  • 1
    $\begingroup$ ……Your Rationalize doesn't have any effect here, it's completely outside of those Set. Your plot is "better" only because this time the parameters are directly introduced into DSolve rather than plugged in by ReplaceAll afterwards so the expressions now are simpler and the numerical instability is somewhat relieved. However, the plot is still wrong, see my comment above. $\endgroup$ Commented Sep 18, 2015 at 8:32
  • $\begingroup$ Yes, Mathematica is packed with bugs. mathematica.stackexchange.com/questions/10820/… $\endgroup$ Commented Sep 18, 2015 at 9:35
  • 1
    $\begingroup$ @xzczd Rationalize wouldn't have any effect, even if it were applied properly. The parameters are already rational. $\endgroup$ Commented Sep 18, 2015 at 10:30
  • $\begingroup$ …Sorry, I should partly take back my words. Though the Rationalize[] is redundant, the plot is correct. It seems that I encountered a bug of LogPlot. Try this: #[Evaluate[{pc[t], pl[t]} /. solAuto /. vals], {t, 0, 1000}, PlotRange -> All, WorkingPrecision -> 16] & /@ {Plot, LogPlot} $\endgroup$ Commented Sep 18, 2015 at 11:04

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.