acer

33255 Reputation

29 Badges

20 years, 246 days
Ontario, Canada

Social Networks and Content at Maplesoft.com

MaplePrimes Activity


These are replies submitted by acer

@Ali Guzel  You didn't supply the details I described as being more important.

I changed your Post into a Question.

But your posting is missing important details. What version are you trying now? More importantly, what did the intermediate results look like in Maple 9.5 or whenever it last worked as desired?

@Hullzie16 For this plotting task it is much faster to use implicitplot rather than a procedure which calls fsolve.

The idea is to first compute the implicit plot where equation R holds (S vs C), and then to transform that plot result, using the T expression as a formula for the new vertical component.

Also, there appear to be two roots, for significant portions of your parameter space.

In this attachment I show Explore working on the implicitplot. I also include the two "curves" generated by taking either the first or second root that fsolve can find, just to show that together they correspond to the implicitplot result.

In the code below the procedures CC (and CC2 and CCboth) are now only used to first show that the new technique is ok. For the faster implicitplot approach (and its Explore) all the calls CC(S,mu,Z,Pr) in expression T are never executed or called(!). Instead, they are replaced by the vertical component y-values from the implicitplot of R, when transforming it. That's what's going on in the second line of procedure Timpl below.

The plot recomputation in the Explore assembly is so much faster with this approach that I can get rid of the continuous option on the Sliders, and get a nice smooth action.

restart;

T := 2^(2/3)*(16*CC(S,mu,Z,Pr)^2*Pi^2 + 2*Pi^2*Z^2 + 16*Pi*S*CC(S,mu,Z,Pr) + 3*S^2)*Pr^(1/3)/(4*(4*CC(S,mu,Z,Pr)*Pi + S)^(2/3)*(2*Pi^2*Z^2 + 4*Pi*S*CC(S,mu,Z,Pr) + S^2)^(2/3)*CC(S,mu,Z,Pr)^(1/3)*Pi^(2/3));

(1/4)*2^(2/3)*(16*CC(S, mu, Z, Pr)^2*Pi^2+2*Pi^2*Z^2+16*Pi*S*CC(S, mu, Z, Pr)+3*S^2)*Pr^(1/3)/((4*CC(S, mu, Z, Pr)*Pi+S)^(2/3)*(2*Pi^2*Z^2+4*Pi*S*CC(S, mu, Z, Pr)+S^2)^(2/3)*CC(S, mu, Z, Pr)^(1/3)*Pi^(2/3))

Tf := unapply(T,S,mu,Z,Pr):

Tf2 := unapply(subs(CC=CC2,T),S,mu,Z,Pr):

R := mu = 2^(2/3)*S*(16*C^2*Pi^2 - 2*Pi^2*Z^2 - S^2)*Pr^(1/3)/(4*(4*C*Pi + S)^(2/3)*(2*Pi^2*Z^2 + 4*C*Pi*S + S^2)^(2/3)*C^(4/3)*Pi^(2/3));

mu = (1/4)*2^(2/3)*S*(16*C^2*Pi^2-2*Pi^2*Z^2-S^2)*Pr^(1/3)/((4*C*Pi+S)^(2/3)*(2*Pi^2*Z^2+4*C*Pi*S+S^2)^(2/3)*C^(4/3)*Pi^(2/3))

Rf := unapply(R,S,mu,Z,Pr):

CCboth:=proc(S,mu,Z,Pr) option remember; local res;
  if not [S,mu,Z,Pr]::list(numeric) then return 'procname'(args); end if;
  res := fsolve(Rf(S,mu,Z,Pr),C=0.01..100,maxsols=2);
  if nops([res])=1 then res,undefined; else res; end if;
end proc:
CC := proc(S,mu,Z,Pr)
  if not [S,mu,Z,Pr]::list(numeric) then return 'procname'(args); end if;
  CCboth(S,mu,Z,Pr)[1];
end proc:
CC2 := proc(S,mu,Z,Pr)
  if not [S,mu,Z,Pr]::list(numeric) then return 'procname'(args); end if;
  CCboth(S,mu,Z,Pr)[2];
end proc:

 

CCboth(1,0,0.5,0.2);
CC(1,0,0.5,0.2);
CC2(1,0,0.5,0.2);

.1938622551, undefined

.1938622551

undefined

CCboth(4,0.5,0.25,0.5);
CC(4,0.5,0.25,0.5);
CC2(4,0.5,0.25,0.5);

.4611869603, 3.243089599

.4611869603

3.243089599

P:=plots:-implicitplot(eval(R,[mu=0.5,Z=0.25,Pr=0.5]),
                       S=0..16,C=0..100,signchange=false);

Timpl := proc(mu,Z,Pr) local P; uses plots, plottools:
  P := implicitplot(eval(R,[:-mu=mu,:-Z=Z,:-Pr=Pr]),
                    S=0..16,C=0..100,signchange=false);

  display(transform((x,y)->[x,eval(eval(T,[CC(:-S,:-mu,:-Z,:-Pr)=y]),
                                   [:-S=x,:-mu=mu,:-Z=Z,:-Pr=Pr])])(P),
          view=[0..15,0..1.5]);

end proc:

 

P1 := Timpl(0.5,0.25,0.5);

# Recall,.. this is not especially fast
P2 := plot(Tf(S,0.5,0.25,0.5),S=0..15,view=[0..15,0..1.5],
           adaptive=false,numpoints=49,color=blue);

# The second if faster than the first, due to remembered values.
P3 := plot(Tf2(S,0.5,0.25,0.5),S=0..15,view=[0..15,0..1.5],
           adaptive=false,numpoints=49,color=magenta);

 

Explore(Timpl(mu,Z,Pr),
        parameters=[[mu=0 .. 1.0, minorticks=0.25],
                    [Z=0..0.5, minorticks=0.1],
                    [Pr=0..1.0, minorticks=0.25]],
        initialvalues=[mu=0.5,Z=0.25,Pr=0.5], adaptview=false);

 

 

Download Explore_Plot_Problem_ac4.mw


Without the above comparison of approaches, this implicitplot approach alone is quite short, and utilizes just your original expression forms for T and R.

T := 2^(2/3)*(16*CC(S,mu,Z,Pr)^2*Pi^2 + 2*Pi^2*Z^2 + 16*Pi*S*CC(S,mu,Z,Pr)
     + 3*S^2)*Pr^(1/3)/(4*(4*CC(S,mu,Z,Pr)*Pi + S)^(2/3)*(2*Pi^2*Z^2
     + 4*Pi*S*CC(S,mu,Z,Pr) + S^2)^(2/3)*CC(S,mu,Z,Pr)^(1/3)*Pi^(2/3)):

R := mu = 2^(2/3)*S*(16*C^2*Pi^2 - 2*Pi^2*Z^2 - S^2)*Pr^(1/3)/(4*(4*C*Pi + S)^(2/3)
          *(2*Pi^2*Z^2 + 4*C*Pi*S + S^2)^(2/3)*C^(4/3)*Pi^(2/3)):

Timpl := proc(mu,Z,Pr) uses plots, plottools:
  display(transform((x,y)->[x,eval(eval(T,[CC(:-S,:-mu,:-Z,:-Pr)=y]),
                                   [:-S=x,:-mu=mu,:-Z=Z,:-Pr=Pr])])(
            implicitplot(eval(R,[:-mu=mu,:-Z=Z,:-Pr=Pr]),
                         :-S=0..16,:-C=0..100,'signchange'=false)),
          'view'=[0..15,0..1.5]);
end proc:

Explore(Timpl(mu,Z,Pr), 'initialvalues'=[mu=0.5,Z=0.25,Pr=0.5], 'adaptview'=false,
        parameters=[[ mu = 0 .. 1.0, 'minorticks'=0.25 ],
                    [ Z = 0..0.5, 'minorticks'=0.1 ],
                    [ Pr = 0..1.0, 'minorticks'=0.25 ]]);

 

 

Download Explore_Plot_Problem_ac4B.mw

@Hullzie16 I prefer assignments over Equation Labels, and I wanted a quick way to shoehorn the R expression into the procedure CC.

There are other ways, eg. using eval(R, [:-S=S,:-mu=mu, etc]) as the argument to fsolve, so that the arguments of CC got put into the R expression when CC were called.  Explore_Plot_Problem_ac2.mw

And here's another variant, using procedures instead of expressions&eval.
Explore_Plot_Problem_ac3.mw

@GFY I wasn't so confused by your use of Eval (vs eval) as I was that you were multiplying D[1] by something. The name D is special, and denotes the differential operator (for use on procs). If your using D[1] as a mere name then I suggest using another name.

On the subject of D, note this comparison of derivates in diff/Diff expression form and in D operator form. Note also that you can utilize Diff instead of diff, versus Eval instead of eval. I'm not saying that D would serve you better, as I don't know your full set of use-cases. But you might find it useful to learn about such, and conversions.

D_diff_ex.mw

(And the D operator can also handle multivariate procs. That's all related to why I found the product involving D[1] confusing. See also the Help-page.)

ps. Like sand15, I simply used lprint of what was actually being sent to algsubs (ie. the result of expand(eqn1)), to see what was going on. And now that's in your tool-set.

I have seen a report of it working in M2025 and M2026 on Windows 11.

@GFY I see now. Somehow I'd got the impression that there might be some other, more involved example.

Yes, this is a key purpose of IntegrationTools:-Expand, which is why I mentioned it in my Answer.

Your example-at-hand shows that I may have been too worried about unwanted extra expansion of subterms in the integrand. I was considering an unfortunate example such as,

    goo := Int( (1-cos(t)^2)*(1-sin(2*x)^2), x=a..b );
    IntegrationTools:-Expand(goo);

But if that kind of thing isn't causing you problems, then all the better.

ps. In your last Reply you seem to be using the word "summation" to refer to an explicit, finite sum of terms, ie. the result of explicitly adding a concrete, finite number of terms. That is not a good use of that word, especially since in your original Question you used around Sigma notation Sigma and a name n as upper end-point, for your query 1).

You wrote, "...using Expand does indeed make it possible to transform the inert integral of a summation into a summation of inert integrals."  But that's not true at all. Transforming an integral of a summation into a summation of integrals is the topic of your other query 1). The Expand command makes it possible to transform an integral of an explicit sum of terms into an explicit sum of integrals.

@GFY You haven't supplied your more complicated, explicit examples, so it's difficult to be sure about what difficulties you might have had with part 2.

It's not clear to me whether you are starting with the integrand in a state that is factored into f(t) and g(x) multiplicands, and whether the stock command attempts were just muddying up the integrand too much for easy recovery.

It's quite possible that a dedicated routine would be much better, once in f(t)*g(x) form. I mean, that could be a reusable, short proc -- something nice in an init file. Would that be useful?

@GFY He has used a Code Edit Region.

This is a link to its Help page in a more recent version than your Maple 2022.

In your version your can access its Help page by query for Topic CodeEditRegion,

@emendes 

I realize that there were already comments that described that proc as being less efficient. But I don't understand the point of showing something flawed in ways that can be remedied. (That's a bit like a straw man...)

Also, another variant given by janhardo in the same Reply contains,

   [seq(seq(subsop(i = model[i] + t, model), t = convert(missing[i], list)), i = 1..n)];

but that lacks a check on whether missing[i] is of type `+`. If it's not, and if it's a single product term, then conversion to a list will produce the wrong result. That's a bug in the code.

Also, not all of these variants seem to give the same results for your example problem.

Parts of the revision shown by janhardo is inefficient. It illustrates a classic programming problem.

I have seen AI generated code (reposted by janhardo), make this mistake before, in other postings.

result := [op(result),
           subsop(i = model[i] + t, model)];

One of the problems with repeatedly augmenting a list in this manner is that it's unnecessarily of a higher computational complexity (in memory). It can scale up unnecessarily poorly in consequence.

A technique of instead adding new items into a table, and then producing a list from that as a final step, is often better. (Note that was done in the OP's earlier version.)

Using another data-structure might be faster, and doing less additional work (than the original) might be faster. But that revision by janhardo likely doesn't scale up as well as it might do, because it's falling afoul of O(n^2) instead of O(n) memory use for assembling some results (and ensuing memory management cost) unneceassarily.

@sand15 Using Maple 2015.2, I was able to "see" that the conditional block of the original code, in,

   showstat(`plots/sparsematrixplot`,4..8);

corresponded to this portion of the inert representation of that original procedure,

   op([5,3],ToInert(eval(`plots/sparsematrixplot`)));

I also noticed that the fragment of that inert code,

   _Inert_STATSEQ(_Inert_ASSIGN(_Inert_LOCAL(8), _Inert_FLOAT(_Inert_INTPOS(0), _Inert_INTPOS(0))))

corresponds to the assignment,

    delta := 0.

So I used subsop to replace the conditional block with itself followed by that simple assignment.

I suppose that it'd be simpler merely to replace that conditional block altogether, by that simple assignment. Ie, in Maple 2015.2,

restart;
`plots/sparsematrixplot` := 
   FromInert(subsop([5,3]=_Inert_STATSEQ(_Inert_ASSIGN(_Inert_LOCAL(8),
                                                      _Inert_FLOAT(_Inert_INTPOS(0),
                                                                   _Inert_INTPOS(0)))),
                    ToInert(eval(`plots/sparsematrixplot`)))):

showstat(`plots/sparsematrixplot`);

One can ascertain that delta is the 8th local of the original procedure. The following is one way to figure out that the inert form of the desired new simple assignment has to be done to _Inert_LOCAL(8) .

restart;

# local sequence of the procedure
op([2],ToInert(eval(`plots/sparsematrixplot`)));

      _Inert_LOCALSEQ(_Inert_NAME("i"), _Inert_NAME("j"), 
      _Inert_NAME("m"), _Inert_NAME("n"), 
      _Inert_NAME("otherOptions"), _Inert_NAME("poly"), 
      _Inert_NAME("polygons"), _Inert_NAME("delta"))

# delta is the 8th local
op(8,%);
                      _Inert_NAME("delta")

The fact that the assignment delta:=0. happens to already occur in the original code means that there's an even terser variant (which can only be done once per session, because it clobbers the original procedure). It's luck that one could also do it this way:

restart;
`plots/sparsematrixplot` := 
   FromInert(subsop([5,3]=op([5,3,4],ToInert(eval(`plots/sparsematrixplot`))),
                    ToInert(eval(`plots/sparsematrixplot`)))):
showstat(`plots/sparsematrixplot`);

These kinds of manipulation looks complicated, but if you manipulate inert representations of procedures often enough it gets easier.

Naturally, it gets much more involved if one wants to check that a target-to-replace exists as one expects. One can guard against mistakes by wrapping in checks of the Maple version number, or testing that the block to replace is of the expected form, etc.

I could mention that if this kind of manipulation gets too involved then (as long as the proc doesn't reference some module's non-exports) one might as well print the proc, copy, and edit the new version as desired.

@awass 

It might be faster if you also changed your h like so,

h:=proc(k,x) option remember;
  local z,w;
  nans(parameters=[k,x]);
  z:=nans(0.98);
  w:=rhs(z[2]),rhs(z[3])
end proc;

 

I'm glad that you got at least one way to work.

note that, in my other way,
    plot([ x->h(7.7,x)[1], x->h(7.7,x)[2], 0..30 ]);
the last entry of the list is 0..30 ie. a plain range, and no longer x=0..30 ie. of type name=range. I'm not sure if you copied literally or whether that might have tripped up you attempt. If not, it's be difficult to say more without the full code of your example.

@dharr Note that the OP is attempting to use a parametric calling sequence of the plot command. And that expects two items in the list, before the name=range.

If you merely delay the evaluation of h(7.7,x)  (eg. by your way, or by simply wrapping in single-quotes) then plot won't see a pair, and it won't receive valid input.

It'd be similar to this:

f:=proc(x,y)
 if not [x,y]::[numeric,numeric] then return 'procname'(_passed) end if;
 cos(x),sin(x*y) 
end proc:

plot([f(x,7.7),x=0..3]);
Error, (in plot) incorrect first argument [f(x,7.7), x = 0 .. 3]

I will submit a bug report against this regression in behavior in Maple 2026.

It seems to me that the problem might occur only if the filled (or filledregion) option is supplied.

In the meantime I'll note that the orientation can be otherwise specified in a wrapping call to plots:-display. Apologies if you were already aware of that, and naturally there's still a bug. Eg.

with(plots):

display(contourplot3d(-x/(x^2+y^2+1),x=-3..3,y=-3..3,filled),
        orientation=[90,45,45]);
1 2 3 4 5 6 7 Last Page 1 of 609