@FormatB

Here's the same Explore in a worksheet. See if that works better.

rat.mw

@Axel Vogt

The main idea is quite simple, abs(x)<y can be transformed into {x^2<y^2,0<y}, so then it's a polynomial system, for which Maple knows how to generate the cylindrical algebraic decomposition. For our problem, the CAD will consist of bounded sectors {a<x<b,f(x)<y<g(x)}, and an iterated integral can be computed.

So the bulk of the code just does the conversion from {a<x<b,f(x)<y<g(x)} to int(1, y=f(x)..g(x), x=a..b). You probably want to throw away lower-dimensional components right away (like vv's code does), and in general you also need to take care of unbounded regions.

@mah_mheasen

plot calls f (the argument will be the value for x5), so you can do the same inside a loop or inside seq(). The computed values are stored in the variable sols.

seq(f(x5), x5 = .1 .. .5, .1);
    [300.9500875, 516.4703171, 105.9894141, 1469.316944], 
    [283.1990384, 183.4187181, 211.6322728, 1155.942231], 
    [505.7638611, 123.3106999, 433.2315225, 1192.981049], 
    [undefined, undefined, undefined, undefined], 
    [undefined, undefined, undefined, undefined]

op(sols);
    table([.3 = [505.7638611, 123.3106999, 433.2315225, 1192.981049],
      .2 = [283.1990384, 183.4187181, 211.6322728, 1155.942231],
      .4 = [undefined, undefined, undefined, undefined],
      .5 = [undefined, undefined, undefined, undefined],
      .1 = [300.9500875, 516.4703171, 105.9894141, 1469.316944]])

@ecterrab

I only noticed an unrelated issue with evalf:

timelimit(60, MeijerG([[1], [0]], [[], []], .2));
Error, (in evalf/hypergeom) time expired

Also, this has a removable singularity at _k1=0, so it cannot be directly passed to evalf:

convert(MeijerG([[0, 1], []], [[], []], z), Sum);
   Sum(((-2*_k1^2-2*_k1)*Psi(_k1)+(-_k1^2-_k1)*ln(z)-3*_k1-2)/(z^(1+_k1)*_k1*GAMMA(2+_k1)^2),
     _k1 = 0 .. infinity)+1

(Symbolically it's fine, separating the term _k1=0 would be more like a quality of life improvement.)

@ecterrab

I disagree. Take MeijerG([[1/2], []], [[], [1]], z), abs(z)>1. The integrand is GAMMA(1/2+y)/GAMMA(y)*z^y, asymptotic to z^y*sqrt(y) for large y. There is only one contour over which the integral converges -- the left loop. The value of the integral is not zero. How do you define MeijerG to make it zero?

@Axel Vogt

Right, the curves in the plot coincide on -1..1, it's exactly the same issue again. MeijerG is discontinuous on the unit circle when p=q and m+n<=p. So you can construct rather a lot of examples like this.

@Axel Vogt

You're correct, it just means that the series representation has the same issue.

Let's trace what is being computed here.

term := (s -> subs(op([2, 1], s) = k, op(1, s)))(convert(MeijerG([[a], [b]], [[c], [d]], z), Sum));

res := -residue(GAMMA(1-a+y)*GAMMA(c-y)/(GAMMA(b-y)*GAMMA(1-d+y))*z^y, y = c+k) assuming k::nonnegint;

simplify(term-res);
                               0

So the series is the sum of the residues over the "right" (as opposed to "left") poles. So you can skip the intermediate step in your reasoning and just compute the sum over the right poles of GAMMA(1-a+y). Of which there are none, so the sum is zero.

The catch is that the contour integral does not converge (for this kernel) over the right loop if abs(z)>1. You have to choose a contour where the integral converges, that's part of the definition. That will be the left loop, and the integral will be equal to the sum of the left residues.

For the sake of having a pretty picture:

Download fsolve.mw

@mah_mheasen

I just transformed the system a bit, because in its original form it seems to be more sensitive to the choice of the starting point or the bounding interval in fsolve, and then chose a value of x5 for which there is a solution. Could have taken x5=1 or x5=2 just as well.

@Markiyan Hirnyk

It's a question of the assumptions. If x and y are real, {sqrt(x)=y} is equivalent to {x=y^2,y>=0}. x>=0 follows from x=y^2, just like Kitonum says. But {sqrt(x)=y} isn't equivalent to {x=y^2,x>=0}. So my reply above about what's equivalent to what wasn't quite accurate.

If x and y in your system are complex, that should be stated from the start. There are additional solutions then.

@Markiyan Hirnyk

I don't think so. The extra condition rhs(sys[2]) >= 0 already says that the value of the square root is nonnegative. That is exactly equivalent to your condition (which says that what's under the square root is nonnegative), regardless of whether we're solving over the reals or over the complexes.

Note though that we are indeed solving over the reals. For the transformed system Maple is using its semi-algebraic solver, which computes the solutions over the real numbers. If you want complex solutions (so y and 4*x^4+4*x^2*y+1/2 still have to be real, but other quantities don't have to be), I think you'll have to rewrite the system in terms of the real and imaginary parts.

@roubeur

I'm simply showing how you can get a result that looks like (Layer1 and Layer3 and ...) and is logically equivalent to the input. Without knowing exactly how a "layer" is defined, I cannot suggest anything more.

@roubeur

You're correct, I did mean a conjunction every term of which is either a variable or the negation of a variable (aka a literal). That's how I interpreted your words that no "or" operation should exist in the simplified equation.

If you want to generate a CNF (that is, if your layers are in fact disjunctions of literals), Normalize(ee, form = CNF) does that in no time. But it's clearly not a minimal length CNF in this case. I don't think Maple has the functionality to generate a minimal length CNF.

Some progress can be made with additional information. E.g., if there are some restrictions on how the variables can be grouped into the layers.

To find a disjunction I simply did this:

vv := indets(ee, symbol);
proc() local i,j,sat;
  for i from 1 to nops(vv) do
    for j from i+1 to nops(vv) do
      sat:=[Satisfy(subs({vv[i]=false,vv[j]=false},&not ee) &and subs({vv[i]=false,vv[j]=true},ee)
        &and subs({vv[i]=true,vv[j]=false},ee) &and subs({vv[i]=true,vv[j]=true},ee))];
      if sat <> [] then return sat[1] end if
end do end do end proc();

@roubeur

Unless I'm missing something, your original expression is not a single conjunction of literals:

ee:=((((((((((((((OD &and (NPi &or NEXTotp)) &and &not ((((((((POinoBlaze &and &not (OD &and
  ((POLYGi1 &and POBIASPdrawing2) &or  (POLYGi1 &and POBIASMdrawing2)))) &or (OD &and
  (spolyp2 &or spolym2))) &and &not (OD &and ((POLYGi2 &and POBIASPdrawing4) &or
  (POLYGi2 &and POBIASMdrawing4)))) &or (OD &and (spolyp4 &or spolym4))) &and &not
  (OD &and ((POLYGi3 &and POBIASPdrawing6) &or (POLYGi3 &and POBIASMdrawing6)))) &or
  (OD &and (spolyp6 &or spolym6))) &and &not BLACKBOXall) &and &not POLYACTsign)) &and
  (((BULK &and &not (SUBS1 &and &not MSUB0)) &and &not (NWELi &and &not MKRMSKREG)) &and &not PWbki))
  &and &not ((NACTALL0 &and &not  NWRNW) &and &not  PWELL1)) &and &not OD225dg) &and &not OD218dg)
  &and &not ((VNPNMSIPW &and &not EMPTY) &or (VNPNMSIPW25 &and &not EMPTY))) &and &not RPOi) &and &not
  ((NACTALL0 &and PWELL1)  &and  POi)) &and &not (((OD &and (PPi &and &not (SEALRINGi  &and &not
  PMDMY))) &and &not POi) &and PWELL1)) &and &not VTHNi) &and &not VTLNi) &and ((DCOng &or DCOmsk) &or
  DCOdg)) &and &not MKRHPA) &and &not EMPTY;

BooleanSimplify(subs({BULK = true, DCOng = false, EMPTY = false, MKRHPA = false, MSUB0 = false,
    NEXTotp = true, NPi = false, NWELi = false, NWRNW = true, OD = true, OD218dg = false,
    OD225dg = false, PMDMY = false, POLYGi1 = false, POLYGi2 = false, POLYGi3 = false, POi = true,
    PPi = false, PWELL1 = true, PWbki = false, RPOi = false, SUBS1 = false, VTHNi = false,
    VTLNi = false, spolym2 = false, spolym4 = false, spolym6 = false, spolyp2 = false, spolyp4 = false,
    spolyp6 = false, BLACKBOXall = true, MKRMSKREG = true, NACTALL0 = false, POBIASMdrawing2 = false,
    POBIASMdrawing4 = false, POBIASMdrawing6 = false, POBIASPdrawing2 = false, POBIASPdrawing4 = false,
    POBIASPdrawing6 = false, POLYACTsign = true, POinoBlaze = false, SEALRINGi = true,
    VNPNMSIPW = false, VNPNMSIPW25 = false}, ee));
                   Logic:-&or(DCOdg, DCOmsk)

Are you sure the original formula is correct.

@tomleslie

1a. The numerator is 1, the denominator has absolute value less than 1. So the result must have absolute value greater than 1. We don't want to unnecessarily extend the resulting interval by including -1..1.

1b. I don't think Maple allows variables as the interval bounds. And in any case, writing INTERVAL(-1..x, x..1) doesn't imply that x has to be between -1 and 1.

2. Documented at ?evalrC.

3. The "bounded variables" are a separate thing, I can't decypher what Maple is doing for those:

evalr(1/INTERVAL(x, -1 .. 1));
       INTERVAL(1/INTERVAL(x, 0 .. 1), 1/INTERVAL(x, -1 .. 0))

So apparently those are two different x's now, and the resulting interval is from something positive to something negative. So either it's empty or it crosses zero. Doesn't seem to make sense either way.