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