You don't necessarily have to use digits=25. But it might be interesing to find out if the computation becomes stable when epsilon is held "large" while digits is raised from 10 towards 15, and then upwards from there.

Other approaches might involve trying to find a better representation of the integral for that problematic range.

You asked how epsilon works. It is estimated, according to the nature of the particular quadrature rule. For example,if you have a quadrature rule of order n (exactly correct results for polynomials up to degree n) then the error at stepsize h can be estimated using the result calculated at stepsize h along with the result calculated at stepsize h/2. This is because you know that the error term is roughly equal to constant K times h^(n+1), where K depends on the (n+1)th derivative of the integrand (and under some assumptions about the derivatives higher than nth being bounded). The refinement of the grid allows comparison of two results, which allows for estimation of the error. Sometimes you can nicely show how it works by plugging the (n+1)th order truncated Taylor series into the quadrature formula.

If I understand you rightly, you are wondering whether it is the evalf/Int computation in which numerical stability might be occuring.

But when you raise Digits high above 15-16, then it becomes too slow.

Note that evalf(Int(...)) takes separate optional keyword parameters to specify either accuracy (epsilon) or working precision (lowecase digits). By default, they both depend on the Digits environment variable. When Digits is rasied then the tolerance epsilon is taken to be a very small number close to 10^(-Digits) or so. And that means the the solver must work very hard. But by supplying them both you can raise the working precision -- to attain greater stability and robustness of individual integrand evaluations -- while simultaneously keeping the desired accuracy tolerance in check -- to keep it from being too slow. It's not always the case, but sometimes it is mainly the striving to attain greater accuracy that makes the quadrature solvers churn.

So, as merely one example, you could try it like this.

  evalf(Int(...., digits=25, epsilon=1.0e-5 ) )

If I understand you rightly, you are wondering whether it is the evalf/Int computation in which numerical stability might be occuring.

But when you raise Digits high above 15-16, then it becomes too slow.

Note that evalf(Int(...)) takes separate optional keyword parameters to specify either accuracy (epsilon) or working precision (lowecase digits). By default, they both depend on the Digits environment variable. When Digits is rasied then the tolerance epsilon is taken to be a very small number close to 10^(-Digits) or so. And that means the the solver must work very hard. But by supplying them both you can raise the working precision -- to attain greater stability and robustness of individual integrand evaluations -- while simultaneously keeping the desired accuracy tolerance in check -- to keep it from being too slow. It's not always the case, but sometimes it is mainly the striving to attain greater accuracy that makes the quadrature solvers churn.

So, as merely one example, you could try it like this.

  evalf(Int(...., digits=25, epsilon=1.0e-5 ) )

That last alternative looks interesting, Joe. In LR's examples, things are simple. But this whole frontend as a way to get derivatives w.r.t non-names comes up a lot.

My first (unposted) idea was this.

  Matrix(3,(i,j)->frontend(diff,[rhs([rossler_sys][j]),[x,y,z](t)[i]]));

I used to think that subs-fore-and-back or `freeze` & `thaw` (so I wouldn't have to think up names for the subs) were not as nice as frontend for this job. But with frontend, there is always the nagging worry that it will freeze or not freeze the right things. Unless the example is "obviously" simple.

Consider the example, where someone says they want the derivative of y(y(t)^2) w.r.t y(t). Using your last alternative, I could get this.

> frontend(diff,[y(y(t)^2),y(t)],[{'Not(identical(x(t),y(t),z(t)))'},{}]);
                                         2
                              2 D(y)(y(t) ) y(t)

That gives something that looks useful. But if the frontend arguments have to be that specific, maybe it's easier after all to use double-subs (freeze and thaw by hand).

That last alternative looks interesting, Joe. In LR's examples, things are simple. But this whole frontend as a way to get derivatives w.r.t non-names comes up a lot.

My first (unposted) idea was this.

  Matrix(3,(i,j)->frontend(diff,[rhs([rossler_sys][j]),[x,y,z](t)[i]]));

I used to think that subs-fore-and-back or `freeze` & `thaw` (so I wouldn't have to think up names for the subs) were not as nice as frontend for this job. But with frontend, there is always the nagging worry that it will freeze or not freeze the right things. Unless the example is "obviously" simple.

Consider the example, where someone says they want the derivative of y(y(t)^2) w.r.t y(t). Using your last alternative, I could get this.

> frontend(diff,[y(y(t)^2),y(t)],[{'Not(identical(x(t),y(t),z(t)))'},{}]);
                                         2
                              2 D(y)(y(t) ) y(t)

That gives something that looks useful. But if the frontend arguments have to be that specific, maybe it's easier after all to use double-subs (freeze and thaw by hand).

Could you give more detail than "many occasions"? For example, do you have an example which replicates the problem (even intermittently, if not consistently)?

The term Software Change Request (SCR) is a euphemism for a bug report.

See here, which I see when logged in in the "navigation" section of the left-pane.

The term Software Change Request (SCR) is a euphemism for a bug report.

See here, which I see when logged in in the "navigation" section of the left-pane.

Did you try the 11.00 or the updated point release? (For Vista, there was the same problem that you describe with file save/open stuff, in the original 11.00 release.)

Would it work if you had just eval(phi,ans) instead, at that point in the loop?

Would it work if you had just eval(phi,ans) instead, at that point in the loop?

> data:=[13.,17.,19.]:

> eqs:={x-z,x+sin(Q*y)-3,y-z-2,y-w,v*x-y}:

> fcn := proc(a,b,c,d,e)
>     c-d-e/tan(a*b);
> end proc:

> expr:=x-y-z/tan(v*w):

> for i from 1 to nops(data) do
>     sol[i]:=fsolve(eval(eqs,Q=data[i]),{w,v,x,y,z});
>     expr_ans[i]:=eval( expr, sol[i] );
>     fcn_ans[i]:=fcn( op(eval([v,w,x,y,z],sol[i])) );
> end do:

> eval(expr_ans);
         table([1 = -1.300706592, 2 = 4.527178744, 3 = -1.699096815])

> eval(fcn_ans);
         table([1 = -1.300706592, 2 = 4.527178744, 3 = -1.699096815])

> eval(sol);
table([1 = {v = 1.782415881, w = 4.556185333, x = 2.556185333, y = 4.556185333,

    z = 2.556185333}, 2 = {v = 1.519017240, w = 5.853436543, x = 3.853436543,

    y = 5.853436543, z = 3.853436543},

    3 = {v = 1.981981500, w = 4.036698248, x = 2.036698248, y = 4.036698248,

    z = 2.036698248}

    ])

> expr_ans[2];
                                  4.527178744

> sol[2];
{v = 1.519017240, w = 5.853436543, x = 3.853436543, y = 5.853436543,

    z = 3.853436543}

> data:=[13.,17.,19.]:

> eqs:={x-z,x+sin(Q*y)-3,y-z-2,y-w,v*x-y}:

> fcn := proc(a,b,c,d,e)
>     c-d-e/tan(a*b);
> end proc:

> expr:=x-y-z/tan(v*w):

> for i from 1 to nops(data) do
>     sol[i]:=fsolve(eval(eqs,Q=data[i]),{w,v,x,y,z});
>     expr_ans[i]:=eval( expr, sol[i] );
>     fcn_ans[i]:=fcn( op(eval([v,w,x,y,z],sol[i])) );
> end do:

> eval(expr_ans);
         table([1 = -1.300706592, 2 = 4.527178744, 3 = -1.699096815])

> eval(fcn_ans);
         table([1 = -1.300706592, 2 = 4.527178744, 3 = -1.699096815])

> eval(sol);
table([1 = {v = 1.782415881, w = 4.556185333, x = 2.556185333, y = 4.556185333,

    z = 2.556185333}, 2 = {v = 1.519017240, w = 5.853436543, x = 3.853436543,

    y = 5.853436543, z = 3.853436543},

    3 = {v = 1.981981500, w = 4.036698248, x = 2.036698248, y = 4.036698248,

    z = 2.036698248}

    ])

> expr_ans[2];
                                  4.527178744

> sol[2];
{v = 1.519017240, w = 5.853436543, x = 3.853436543, y = 5.853436543,

    z = 3.853436543}

See here.

(Why doesn`t that reply show up in Mapleprimes Search results against the terms shake and evalr. I found it using Google itself. There a few users whose posts do not seem to be returned in Search result, based purely on words unique to their own replies. I wonder if that is deliberate.)

Context can be everything. If I recall correctly, the earlier queries were about obtaining equations/results for integrals with a mix of symbolics and floats, and not immediately about some end-goal of doing numerical quadrature to obtain purely floating-point results. So, IntregrationTools:-Expand might have been the answer to the question then asked, but the context in which it was asked may have not reflected the actual final goal.

For example, here your integrand contained symbolic M. Even though you said it represented a constant numeric value, your question did seem to be all about getting an integration result with M a symbol.