Problem using WhenEvent to constrain solution

Note: This question has also been posted at the Wolfram Community

Problem: Simulate pressure in volume 1 and 2 for 1 second.

The circuit is as follows:

From this I set up the governing DE for both volumes and plot the solution:

Ad1 = 10/1000^2;

Cd1 = 0.67;

Cd2 = 0.67;

Cd3 = 0.67;

V1 = 10/1000;

V2 = 10/1000;

Rho = 875;

beta = 1000*10^6;

ps = 100*10^5;

Q3 = Q1 - Q5;

Q4 = Q2 + Q5;

s = NDSolve[{p1'[t] == (beta*Q3)/V1, p2'[t] == (beta*Q4)/V2,

p1[0] == p2[0] == 0, WhenEvent[p1[t] >= ps, p1[t] -> ps],

WhenEvent[p2[t] >= ps, p2[t] -> ps]}, {p1, p2}, {t, 0, 1}];

Plot[Evaluate[{p1[t]/10^5, p2[t]/10^5} /. s], {t, 0, 1}]


Mathematica aborts integration after t = 0.69 since it encounters complex solutions. This is due to p1 and p2 at some point getting larger than ps, which makes no sense, particularly since I have added 2 WhenEvents which shouldn't 'allow' p1 and p2 to be greater than ps (see code). Complex solutions can be avoided by adding Abs, however, then the solution seem to completely diverge:

Question: Why doesn't WhenEvent seem to 'work'?

PS.; I have obtained a more credible solution using Matlab and the same constrictions:

Matlab code:

close all;
[b][/b]
clear;

%Basic data

EndTime=1;

StepTime=1e-5;

ps=100*1e5;

Cd1=0.67;

Cd2=0.67;

Cd3=0.67;

V1=10/1000;

V2=10/1000;

rho=875;

beta=1000*1e6;

%Initialize

p1_initial=0;

p2_initial=0;

%Initially old values are simply set to current values

Time=0.0;

p1=p1_initial;

p2=p2_initial;

%Initialize counters so that plot data is only saved once pr. a number of

%time steps corresponding to ReportInterval

ReportCounter=0;

ReportInterval=10;

Counter=ReportInterval;

%Start time integration

while Time<EndTime

if p1>=ps

p1=ps;

end

if p2>=ps

p2=ps;

end

Q3=Q1-Q5;%Q3 to volume 1

Q4=Q2+Q5;%Q4 to volume 2

p1Dot=beta*Q3/V1;

p2Dot=beta*Q4/V2;

%report

if Counter==ReportInterval

Counter=0;

ReportCounter=ReportCounter+1;

Time_Plot(ReportCounter)=Time;

p1_Plot(ReportCounter)=p1*1e-5;

p2_Plot(ReportCounter)=p2*1e-5;

end;

%Time integrate

p1=p1+p1Dot*StepTime;

p2=p2+p2Dot*StepTime;

Time=Time+StepTime;

Counter=Counter+1;

end;

plot(Time_Plot,p1_Plot);

hold on;

plot(Time_Plot,p2_Plot,'r');

grid;


Matlab solution seems to be an OK fit with simulation results (from SimulationX):

-
Please first get rid of those nasty Subscript occurrences and do not use := to set values to simple constants use =. By the way C is reserved for something special in Mathematica, check the doc. – PlatoManiac Nov 27 '13 at 22:48
's' is being used a parameter and as the result of NDSolve, so that's an issue as well – chuy Nov 27 '13 at 22:54
Hi thanks for replying. I've cleaned up the code now, sorry I didn't do it before, but I encountered several new errors when I tried. Should be better now. C is a reserved character yes, but you can use it if you use a subscript to distinguish it. Changing 's' to some other name does not change the output errors. – Julian Nov 28 '13 at 9:30

The problem with WhenEvent has to do with the OP's DE. For an event to be detected, there has to be a point at which the condition is crossed, that is, changes from False for t < t0 to True for t > t0. NDSolve then applies a root-finding algorithm to approximate the value of t0 at which the event occurs. In your DE, the solution p1[t] theoretically is always less than or equal to ps. If it is ever greater, it is because of error due to approximation; further, the equations become non-real when p1[t] > ps. So on the one hand p1[t] should not cross ps, and practically if in computing a step it does, the value for p1[t] at that step becomes complex and so p1[t] > ps becomes nonsense. It is neither True nor False. This makes finding the value of t where p1 crosses from being less than ps to being greater than ps impossible. The same goes for p2[t]. In short, you cannot use WhenEvent in this way.

How to use WhenEvent

The way around it is to construct a condition that will switch from False to True. In this case, checking that p1[t] (resp. p2[t]) has no imaginary part does the trick.

s = NDSolve[{p1'[t] == Identity[(beta*Q3)/V1],
p2'[t] == Identity[(beta*Q4)/V2],
p1[0] == 0, p2[0] == 0,
WhenEvent[Im@p1[t] != 0, p1[t] -> ps],
WhenEvent[Im@p2[t] != 0, p2[t] -> ps]}, {p1, p2}, {t, 0, 1},
Method -> "StiffnessSwitching"]

(* {{p1 -> InterpolatingFunction[{{0., 1.}}, <>],
p2 -> InterpolatingFunction[{{0., 1.}}, <>]}} *)

Plot[Evaluate[{p1[t]/10^5, p2[t]/10^5} /. s], {t, 0, 1}]


Original way (avoiding WhenEvent)

Another way is to clip the value of ps - p1[t] etc. to always be nonnegative. Since theoretically they always are, this is valid.

Clear[clip];
clip[x_?NonNegative] := x;
clip[x_?NumericQ] := 0;

d1 = 10/1000^2;
Cd1 = 0.67;
Cd2 = 0.67;
Cd3 = 0.67;
V1 = 10/1000;
V2 = 10/1000;
Rho = 875;
beta = 1000*10^6;
ps = 100*10^5;
Q3 = Q1 - Q5;
Q4 = Q2 + Q5;

s = NDSolve[{
p1'[t] ==(beta * Q3)/V1,
p2'[t] ==(beta * Q4)/V2,
p1[0] == 0, p2[0] == 0},
{p1, p2}, {t, 0, 1}, Method -> "StiffnessSwitching"]

(* {{p1 -> InterpolatingFunction[{{0.,1.}},"<>"],
p2 -> InterpolatingFunction[{{0.,1.}},"<>"]}}  *)


Further explanation

The important feature of the DE in this case is a branch point, where p1[t0] == p2[t0] == ps. If, using these initial conditions, one were to integrate backward from t0, one has to choose which branch to follow, the constant solution, p1[t] == p2[t] == ps, or the non-constant solution. When doing the forward integration, the non-constant solution runs into the constant one.

In the OP's straightforward implementation of Euler's method, this difficulty is ignored, except for constraining the values of p1, p2. (That implementation is easily translated into Mathematica, by the way.) Depending on the options used, Mathematica's rather sophisticated NDSolve stops the integration when the branch point is reached (approximately). This could be considered correct behavior in some contexts. One runs into similar issues in DEs that can be integrated in finite terms. Consider the analogous DE $y' = \sqrt{1-y}$, whose solutions are $y(t) = (4-(t-C)^2)/4$, which is valid only for $t \le C$ (for $t > C$, $y'(t) < 0$), and the constant solution $y(t) = 1$. There is something to be said for keeping the branches distinct. One can trace the trajectory backwards and forwards on the non-constant solutions. But on the constant solution, the history of the trajectory cannot be determined.

-
thanks alot. I understand all your arguments, although to me, it seems strange that there is no better way to tell mathematica p1<=ps. This is the only time I have ever felt that Matlab does a better job than Mathematica, simply because of the more convenient way of adding boundary conditions. Anyway, I just added the Abs@ to avoid complex solutions and clipped the plot at 100 bar. Still feels like this is a 'cheating' way of doing it, but it does the job I guess. Thank you very much. – Julian Dec 4 '13 at 16:02
@Julian I found a way to use WhenEvent (see edit). I feel rather that Mathematica does a better job here, but just not the job that you need for your application (at least not easily). There are mathematical reasons why Mathematica does what it does. I tried to hint at them above, but I can't really say a lot more without sheer guesswork at what M is actually doing. – Michael E2 Dec 4 '13 at 22:41
Yes I saw it, and I thought it was the most straightforward way of solving this problem. I used it in a solution myself (see below). Thanks a lot. – Julian Dec 4 '13 at 23:18

You could try specifying the event differently (this will trigger when the absolute difference between p1 and p2 is smaller than some value):

Ad1=10/1000^2;
Cd1=0.67;
Cd2=0.67;
Cd3=0.67;
V1=10/1000;
V2=10/1000;
Rho=875;
beta=1000*10^6;
ps=100*10^5;

Q3=Q1-Q5;
Q4=Q2+Q5;

s = NDSolve[{p1'[t] == switch[t] (beta*Q3)/V1,
p2'[t] == switch[t] (beta*Q4)/V2, p1[0] == p2[0] == 0, switch[0] == 1,
WhenEvent[Abs[p2[t] - p1[t]] <= 10, switch[t] -> 0]},
{p1, p2}, {t, 0, 1}, DiscreteVariables -> {switch}]

Plot[Evaluate[{p1[t]/10^5, p2[t]/10^5} /. s], {t, 0, 1}]


-
very nice, elegant solution. Thanks alot – Julian Dec 4 '13 at 18:47

2 problems: You were comparing, in the WhenEvent, solution, which had complex value at that t, to real numbers. I used Abs. If this does not work for you, you can use Re, but can't compare complex number to real number using >.

Second, your system is stiff, need to use StiffnessSwitching to help NDSolve.

d1 = 10/1000^2;
Cd1 = 0.67;
Cd2 = 0.67;
Cd3 = 0.67;
V1 = 10/1000;
V2 = 10/1000;
Rho = 875;
beta = 1000*10^6;
ps = 100*10^5;
Q3 = Q1 - Q5;
Q4 = Q2 + Q5;
s = NDSolve[{p1'[t] == (beta*Q3)/V1, p2'[t] == (beta*Q4)/V2,
p1[0] == 0, p2[0] == 0,
WhenEvent[Abs@p1[t] >= ps, p1[t] -> ps],
WhenEvent[Abs@p2[t] >= ps, p2[t] -> ps]
}, {p1, p2},
{t, 0, 1}, MaxSteps -> 10000, Method -> "StiffnessSwitching"
];
Plot[Evaluate[{p1[t]/10^5, p2[t]/10^5} /. s], {t, 0, 1}]


-
thanks for reply. It's possible to avoid the complex numbers all together and plot the solution up to t=1 (see edit). The problem then is, that the 'WhenEvent' doesn't seem to have any effect at all, as p2 goes way beyond the 'allowed' value. – Julian Nov 28 '13 at 13:58
WL result It looked like the Matlab image to me. Can you post the text of the full Matlab code which you say works so one can compare easier? Hard to read code from an image. – Nasser Nov 28 '13 at 14:57
Absolutely. It's up now. – Julian Nov 28 '13 at 18:24
@Julian I can't tell now why it failes for over 1 second. But looking at your matlab code, in there, you are not using ode45, but using simple Euler steps. So your Matlab solution itself might not be accurate. If all else fails, you can code it in Mathematica the same way as you did it in Matlab, i.e manual integration. What happens is that the ps-p1[t] becomes negative at one point, and sqrt returns complex then. But you really need to first write down a mathematical model of your equations, for your own benefit. – Nasser Nov 28 '13 at 19:18
I've tried using Method -> "ExplicitEuler". In that case, it yields the exact same curve as before. – Julian Nov 28 '13 at 19:47

Here is another way to solve this issue.

I know, it is written everywhere, but this is a common mistake to think that the orifice equation is

Cd * Ad * Sqrt[ 2/Rho * ( ps - p1[t] ) ]


and could give complex results ! There is no physics going imaginary in the real world. The equation is just wrong. When the pressure reverse, the flow reverse, at least.

So a better equation that stands for the reversibility of flow requires a signed square root function :

SignedSqrt[x_] := Sign[x] * Sqrt[ Abs[x]]


With that function you get rid of the complex issue, without having to make any hypothesis about the flow direction.

The other thing is that this equation is not realistic around null flow. Around null flow, in reality, there is still a resistive effect, that is not described by this equation. It is known to be laminar flow, while the previous equation describe only the turbulent part, which happen over the reynolds critical number. The equation is of the form

ModifiedSignedSqrt[x_] := If[ Abs[x]< 1, x, Sign[x] * Sqrt[ Abs[x]] ]


(it switches between flow types around 1, in this simplified case, but must be adapted to the real transition point).

Without this resistive effect the solver try to solve an undamped system dynamics around zero, which is numerically impossible without approximation.

So, if your goal is to simulate a physical system, a simple solution is to use an equation of flow more physical, this will solve all issues : no complex, no stiffness, no need of event handler...

-

Ok, so this is my attempt at a short discussion and solution (after getting a little bit wiser).

Due to these equations being stiff (Wikipedia) as pointed out by @Nasser, certain numerical methods have difficulty at tracking the solution, and thus, approximations of the orifice equations yield complex solutions at some point when p1[t]>ps and/or p2[t]>ps.

There are several workarounds to this (I also had some help from the Wolfram Community). The one I find to be the most straightforward is the one pointed out by @Michael E2 by letting WhenEventdetect imaginary parts of the orifice equations and constrain p1[t]and p2[t] to a max value of ps.

To be honest, I still don't quite understand why the original WhenEvents didn't 'work'. The first time I solved this I used Matlab, and all I had to do was constrain p1 and p2 like this:

if p1>=ps
p1=ps;
end
if p2>=ps
p2=ps;
end


That's it. My understanding was that WhenEventwould operate similarly, as sort of an override, or a boundary condition that was non-negotiable so to speak.

Anyway, here's my own solution which gives the correct curves.

s = NDSolve[{p1'[t] == (beta*Q3)/V1, p2'[t] == (beta*Q4)/V2, p1[0] == 0, p2[0] == 0, WhenEvent[Im@(p1[t] - p2[t]) != 0, p1[t] -> ps], WhenEvent[Im@(p1[t] - p2[t]) != 0, p2[t] -> ps], WhenEvent[Im@(ps - p2[t]) != 0, p2[t] -> ps], WhenEvent[Im@(ps - p1[t]) != 0, p1[t] -> ps]}, {p1, p2}, {t, 0, 1}, Method -> {"StiffnessSwitching"}]; Thanks a lot guys.

-
The reason why the original WhenEvent doesn't work: For WhenEvent to work, the condition has to evaluate to True or False. If p1[t] is complex, then p1[t] >= ps does not evaluate to True or False. You could force the issue with TrueQ: For example, WhenEvent[! TrueQ @ Quiet[p1[t] < ps], p1[t] -> ps] works. (The Quiet keeps it from complaining about comparing a complex number.) – Michael E2 Dec 5 '13 at 1:58
But p1[t] and p2[t] are not complex? It's the difference between them that becomes complex, due to errors in approximation. From the orifice equation Cd*Ad*Sqrt[(2/Rho)*(ps-p1[t])], at some point p1[t]gets larger than ps and voila, we have complex numbers. My understanding was, that adding the WhenEvent[p1[t]>=ps,p1[t]->ps] would tell Mathematica that "psis always larger or equal to p1[t], period", and therefore avoid this issue. So, when at some step when p1[t] does get greater than ps, I thought that WhenEventwould go from Falseto Trueand that's that. – Julian Dec 5 '13 at 10:28
So, I think it's strange that there is no more convenient way of adding these boundary conditions (as in Matlab). It all comes down to my understanding of how WhenEventactually works I guess. And thanks a lot @MichaelE2 for that last workaround also. Worked perfectly. – Julian Dec 5 '13 at 10:33