@digerdiga

Executive summary: if you have a fixed CO>1/(2*k_OHCO), the next term in the expansion in NO2 has order 3/2. So looking for a Taylor series is futile.

Also, the code that you posted breaks for me (Maple 2017.3) at op_OH_0 := op(2, eq_OH_0). Extracting parts of an expression in this way is rather fragile. numer(eq_OH_0) works though.

(The link is below, but I can't embed the worksheet in the post; I'm getting "Failed to load the worksheet /maplenet/convert/rootof.mw".)

rootof.mw

@Carl Love

Ah, of course, a nested procedure does the trick nicely. I was thinking, Can this be done with seq((x,t)->td(x,t)[i], i=1..3)? Nope, seq doesn't work that way with scoping constructs. I ended up with seq(unapply('td'(x,t)[i],[x,t]), i=1..3), which wasn't very satisfactory because I was back to uneval.

P.S. In fact, SignalProcessing:-FFT has custom code only for one-dimensional arrays, for higher dimensions it just goes to DiscreteTransforms:-FourierTransform.

@acer

I think the documentation could state more clearly that the functional form of plot/fsolve/Maximize can be used to delay the evaluation of the first argument, and then the function will be evaluated only for numerical values of the arguments (plot of x->rand(x..x+1)()).

Here [(x, t) -> td(x, t)[1], (x, t) -> td(x, t)[2], (x, t) -> td(x, t)[3]] will do the job, but it's more tricky if one wants to generate the list programmatically. And it's not very efficient to call td(x, t) three times. If plot/fsolve accepted vector-valued procedures, that would be a cleaner solution.

@Carl Love

But simplify receives diff(..., t), not diff(..., x). The integral does not depend on t, so zero is actually correct :).

@vv

I don't think it's been fixed. I'm on 2017.3, I believe 2017.3 already includes the fix from August 11. At least the examples from pdsolve_adjustment.mw all run fine.

@Markiyan Hirnyk

This one diverges though, so what should evalf/Int return?

@acer

Both Split() and epsilon hints are very useful, thanks. In fact, I don't even need epsilon close to 10^(-digits), a result correct to just a few digits is already enough to analyze the behavior of the integral when the limit of integration approaches 2 (I do need a higher setting of digits for that, but epsilon can stay fixed).

It's a bit unfortunate that Maple actually does most of the job but doesn't return an answer and doesn't give any indication that it happened because the answer hadn't met the required tolerance. If, roughly speaking, the default setting of epsilon allows the loss of only one significant digit, that seems too stringent.

@Christopher2222

Right, but point #3 is about numerical integration. int evaluates the integral, so no numerical integration is taking place.

Symbolic int is the last example in point #1. The value 1.587-1.262*I is wrong.

@Ronan

I tried to do a variation of the integral discussed in your post numerically, the result probably belongs here too:

evalf(Int(1/(1-4*JacobiSN(s, 1/2)^2), s = 0 .. InverseJacobiSN(1/2, 1/2)-10^(-10))-
  InverseJacobiSN(1/2, 1/2));
                                     -1.*10^(-10)

The integrand has a pole at InverseJacobiSN(1/2, 1/2), but the numerical result approaches InverseJacobiSN(1/2, 1/2) instead of increasing. Might indicate that there is a wrong symbolic transformation somewhere.

@digerdiga

It's the other way around: without uneval (?' opens the help page for it), the first argument of eval is being evaluated before a and x receive numeric values, hence the unexpected result.

If you know in advance that you have two positive roots, you can just do this:

plot([x -> r(x, 1)[1], x-> r(x, 1)[2]], -.4 .. .2);

This way r(x, 1) is evaluated only for numeric values of x. Note that you can't plot a procedure that returns a list, you have to write out the components.

@digerdiga

This is the issue:

r(x, a);
                               []

For symbolic arguments, the test in select never gives true. This would work:

eval('r(x, a)', [x = .1, a = 1]);
  [0.55884077383873353962589443719771779995909073393621 + 0. I,
   0.18621290461186062353425226698348130683442111573059 + 0. I]

And also there is a standard way of making a procedure stay unevaluated for non-numeric inputs:

r := proc(x, a)
 if not [x,a]::list(numeric) then return 'r(x,a)' end if;
 # other stuff
end proc

But without knowing more about your actual problem, I don't see why you don't want to use RootOf, which is happy to stay unevaluated for symbolic parameters.

@vv

As an aside, the integrals are doable by changing x^2 to t:

ee:=(value@SI)(1, [x^2+y^2+z^2 - 1, 5*x^2-z^2 - 1], [x,y,z]);

applyrule(
  'int'(e::anything, r::anything) =
    '(value@simplify@PDETools:-dchange)(
      `if`(is(op([2, 1], r) > 0), {x = sqrt(t)}, {x = -sqrt(t)}),
      Int(e, r))',
  ee);
        (16/3*I)*sqrt(2)*sqrt(3)*EllipticPi(1/3, (1/3)*sqrt(3)*sqrt(5))-(16/3*I)*sqrt(2)*sqrt(3)*
        EllipticK((1/3)*sqrt(3)*sqrt(5))-(16/3*I)*sqrt(2)*sqrt(3)*EllipticPi((1/5)*sqrt(3)*
        sqrt(5), 1/3, (1/3)*sqrt(3)*sqrt(5))+(16/3*I)*sqrt(2)*sqrt(3)*EllipticF((1/5)*sqrt(3)*
        sqrt(5), (1/3)*sqrt(3)*sqrt(5))+2*sqrt(6)*EllipticK((1/9)*sqrt(6))+(10/9)*sqrt(6)*
        EllipticPi(4/9, (1/9)*sqrt(6))-(12/5)*sqrt(6)*EllipticE((1/9)*sqrt(6))

evalf(%, 20);
        11.146863418040527043-2.781662233*10^(-9)*I

That's for the surface shown by RegionPlot3D, which I think would be defined as

ImplicitRegion[
  x^2 + y^2 + z^2 == 1 && 5 x^2 - z^2 <= 1 ||
  x^2 + y^2 + z^2 <= 1 && 5 x^2 - z^2 == 1, {x, y, z}]

(doesn't matter if the inequalities are strict or non-strict).

@digerdiga

This would work:

r := c -> max(select(t -> is(Im(t) = 0), [sol(c)])); # or is(abs(Im(t)) < eps)
plot('r(c)', c = -10 .. 10);

But instead you can use RootOf to select the first real root of p(-x):

plot(-RootOf(subs(x = -_Z, (x-1)*(x-2)*(x-3)+c), index = real[1]), c = -10 .. 10);

Then you don't have to worry about the spurious imaginary part or about which radical expression corresponds to the real root for a given c.

If you say you get a wrong solution for your problem of degree 8, you'll need to post the actual code. Most likely a wrong solution branch is being selected somewhere.

@Vortex

Sorry for the delay, you probably have found the error already -- if you look at the error message, it complains about ".99 = .001 .. .999". That should have been "t = .001 .. .999". t got assigned the value .99 somewhere in your session.