acer

32864 Reputation

29 Badges

20 years, 141 days
Ontario, Canada

Social Networks and Content at Maplesoft.com

MaplePrimes Activity


These are answers submitted by acer

restart;

kernelopts(version);

`Maple 2024.2, X86 64 LINUX, Oct 29 2024, Build ID 1872373`

t := (5/9)*Pi:

e := tan(t) + 4*sin(t);

-tan((4/9)*Pi)+4*sin((4/9)*Pi)

radnormal(expand(convert(e,exp)));

-3^(1/2)

Download simp_one_way.mw

ps. Note that simplification-to-zero is quite often easier than is simplification to a particular form. Here, simplify(e+sqrt(3)) produces exact zero (0), directly, and that may well be related to why is can also handle it.

I agree with the OP that a strong (but potentially costly) simplifier would be useful some times. My oracle implementation of such gave me the incantation above which attained the result.

However, that particular sequence of commands makes reasonable sense here, and is something I'd try if left to my own devices: expanding, after conversion from trig to exp form, can often lead to a form containing things like -1 raised to rational powers (ie. powers of roots of -1), and radnormal or evala are natural ways to deal with simplifying such to nicer radicals.

For fun, (and noting that some special case handlers are also possible).

nb. The eval of calls by uneval-quoted `[]` is there to side-step inadventent arithmetic-on-lists(!) at intermediate replacement stages.

You can of course remove function or `<=` or any other type from the following code. I'm guessing at what you want covered.

restart;

H := ee -> eval(subsindets(ee, {`+`,`*`,`^`,`.`,
                                `=`,`<`,`<=`,function},
                           u -> '`[]`'(op(0,u),op(u)))):

 

H( (a+b*c)^2 );

[`^`, [`+`, [`*`, b, c], a], 2]

H( a*b(x)*c*d );

[`*`, a, [b, x], c, d]

H( sin(Pi*x*f(t))=Q );

[`=`, [sin, [`*`, Pi, x, [f, t]]], Q]

# automatic simplification by kernel
H( 2*(x+y) );
H( q >= p );

[`+`, [`*`, 2, x], [`*`, 2, y]]

[`<=`, p, q]

H( Int(f(x),x=1..s) );

[Int, [f, x], [`=`, x, 1 .. s]]

H( b/c );

[`*`, b, [`^`, c, -1]]

# What other operators are you interested in?
lprint(whattype(1..s));

`..`

Download op_fun02.mw

Keep in mind that this is doing what's asked of it, working with the expression to which its input evaluates.  Eg. b/c vs the expression b*c^(-1), or q<=p vs the expression p<=q, etc.

It seems to me that having the transformer (3rd argument) of the subsindets call be something simple that appears like your simple original definition/requirement, is a natural kind of approach. I mean that in a hand-wavy kind of manner. If you start complicating the definition/requirement then writing a more complicated process might be order.

The description of what you want, and the details of the requirements, are key here. You've only provided a single example and a short description in your original Question. So extending that to other kinds of examples is dubious because of your thin specification.

Note that your Question's title mentions "of an expression", as if all expressions were to be covered. That doesn't explain how to deal with expressions that evaluate to something other than themselves. (Look at the operands of the result of a numtheory[cfrac] call...)

Also, your Question title just says listlist, as if your single example's output's list structure involving op(0,..) and op(...) were the only possible consequent form. But it's not, right? I mean you don't even have to flatten everything into a single sublist.

[edit] I'm a bit cautious about the recursion. The examples I've given can also be covered by,
    ee->eval(subsindets(ee,Not(atomic),u->'`[]`'(op(0,u),op(u))))
where the only difference is in replacing a specially curated choice of types by the blanket Not(atomic). It's not that I'm so worried about a runaway, as I am that you might not want so many things descended into.

Not everything that is not of type atomic is going to behave so nicely under op (see table example in my Reply to janhardo). So I would rather see a clear and careful collection of mentioned types, even if it's longer code.

There are some cases where the recursion by this particular subsindets approach will vary from that done by the recursive construction in the code given by janhardo. Eg, Vector([3,f(x)]) , but I'm not sure whether that matters to you. You didn't mention that you needed, say, to be able to reconstruct expressions from the resulting listlist.

Your full and precise goal is unclear.

I find that sometimes Maple does a bit better with simplifying in ln form, for arctrig.

restart;

kernelopts(version);

   Maple 2025.1, X86 64 LINUX, Jun 12 2025, Build ID 1932578

f:=subs(r=(sqrt(5)+1)/2,7*arctan(r)^2
        +2*arctan(r^3)^2-arctan(r^5)^2):

simplify(combine(convert(f,expln)));

   7/8*Pi^2

But even there, combine as intermediate step allowed it to work. (Sometimes, a factor call at just the right intermediate point helps...) It's not always obvious what intermediate form is simpler or will lead to further improvement.

The simplify command has to deal with alternate forms and paths, while avoiding getting caught in loops or inadvertantly trying something that can "blow up" (like some expand). It's heuristic, and often tweaked or improved. That's one reason why it's useful to report such weaknesses.

ps. Here we could convert to either ln or expln, to change the arctrig.

No this type check is not supposed to expand/simplify/collect-wrt-y. My claim is supported in a few places by documentation.

It's supported by parts of the 3rd bullet point of the help-page for topic type,linear . While discussing a different but closely situation it mentions, "[would otherwise require]...that the polynomial is expanded (collected) in x, and the type test is only syntactic". The upshot there is that this check does not expand or otherwise simplify, and is more syntactic in nature.

Also, degree(A,y) returns 2, while degree applied to collect(A,y) or expand(A) or simplify(A) returns 1. That kind of behaviour is explicitly documented on the degree help-page (later example, where "...normalization is necessary").

And the type,linear help-page explicitly defines its check as including a check that the degree is 1, in the 2nd bullet-point of its Description. The code also reflects that "definition":

showstat(`type/linear`);


`type/linear` := proc(f, x)
   1   if _npassed = 1 then
   2       type(f,polynom) and degree(f) = 1
       elif type(x,list) then
   3       type(f,linear({op(x)}))
       else
   4       type(f,polynom(anything,x)) and degree(f,x) = 1
       end if
end proc
 

Download type_linear.mw

Having said that, it'd be reasonable (I think) for the is command to be able to figure this out, mathematically rather than syntactically. And the type linear is also considered a property, in this context. However, in Maple 2024.2 the call is(A,linear(y)) unfortunately returns false, for your unexpanded A example. You might submit a bug report against that.

ps. quotes in italics above denote quotations from the documentation.

[edit] One might not always want to expand or simplify everything (eg. trigs/radicals/spec.-func./etc may be present). A practical alternative might be to do the check on collect(A,y) (like the last ?degree example). Or possibly frontend an expand.

How would you want to handle an example like, say, (y^2-1)/(y+1)+7 ?

Some examples,

A := (x+y+1)^2-(x+y-1)^2;

(x+y+1)^2-(x+y-1)^2

collect(A,y);

4*y+(x+1)^2-(x-1)^2

type(%,linear(y));

true

frontend(expand,[A]);

4*x+4*y

type(%,linear(y));

true

B := (y^2-1)/(y+1)+7;

(y^2-1)/(y+1)+7

collect(B,y);

(y^2-1)/(y+1)+7

type(%,linear(y));

false

frontend(expand,[B]);

y^2/(y+1)-1/(y+1)+7

type(%,linear(y));

false

simplify(B); # or normal(B); # loss of knowledge

y+6

type(%,linear(y));

true

Download type_lin_ex1.mw

I was able to get the following. If you are willing to accept a condition like x>0 for the odetest then it seems reasonable to try also passing that to dsolve.

The dsolve call happens to use the extra condition x<1, but the odetest check doesn't. (Without it, the dsolve call goes away for quite a while...)

Memory use is high, and sometimes it misbehaves, especially if doing further computations in the session. (Once I got a "bad id" error from the kernel.)

restart;

kernelopts(version);

`Maple 2024.2, X86 64 LINUX, Oct 29 2024, Build ID 1872373`

ode:=-x*sqrt((1 - x)/(x + 1))*(x + 1)*arcsech(x)*diff(y(x), x)*exp(y(x)/arcsech(x) + exp(y(x)/arcsech(x))) - y(x)*exp(y(x)/arcsech(x) + exp(y(x)/arcsech(x))) + 2*x*sqrt((1 - x)/(x + 1))*(x + 1)*arcsech(x)^2 = 0;

-x*((1-x)/(x+1))^(1/2)*(x+1)*arcsech(x)*(diff(y(x), x))*exp(y(x)/arcsech(x)+exp(y(x)/arcsech(x)))-y(x)*exp(y(x)/arcsech(x)+exp(y(x)/arcsech(x)))+2*x*((1-x)/(x+1))^(1/2)*(x+1)*arcsech(x)^2 = 0

new := simplify(convert(ode,expln)) assuming x>0;

(ln(x)*(diff(y(x), x))*(-x^2+1)^(1/2)*x-ln((-x^2+1)^(1/2)+1)*(diff(y(x), x))*(-x^2+1)^(1/2)*x-y(x))*exp(((ln(x)-ln((-x^2+1)^(1/2)+1))*exp(y(x)/(-ln(x)+ln((-x^2+1)^(1/2)+1)))-y(x))/(ln(x)-ln((-x^2+1)^(1/2)+1)))+2*x*(-x^2+1)^(1/2)*(ln(x)-ln((-x^2+1)^(1/2)+1))^2 = 0

sol := CodeTools:-Usage( dsolve( new ) ) assuming x>0, x<1;

memory used=3.78GiB, alloc change=81.56MiB, cpu time=52.98s, real time=49.35s, gc time=11.25s

y(x) = -ln(ln(2*x-c__1))*ln(x)+ln(ln(2*x-c__1))*ln((-x^2+1)^(1/2)+1)

simplify(convert( odetest(sol,ode), expln)) assuming x>0;

0

Download de_hard02.mw

You could simply use Not(atomic), which imo is still pretty legible.

Or, if you'd like you can set that up as a new named type (and, optionally, put it in your init file if wanted).

restart;

 

TypeTools:-AddType( nonatomic, Not(atomic) );

 

type( sin, atomic );

true

type( sin(x)+cos(x), atomic );

false

type( sin(x)+cos(x), nonatomic );

true

Download AddType.mw

Is this the kind of thing you're after, making a distribution and random-variables based upon that PDF, and then, say, combining such RVs?

with(Statistics)

dz := proc (t) options operator, arrow; piecewise(t <= 0, 0, 0 < t and t < 2, -(3/80)*(-2+t)^3*(4+6*t+t^2), 2 <= t, 0) end proc

proc (t) options operator, arrow; piecewise(t <= 0, 0, 0 < t and t < 2, -(3/80)*(-2+t)^3*(4+6*t+t^2), 2 <= t, 0) end proc

int(dz(t), t = 0 .. 2)

1

zdist := Distribution(PDF = dz)

Z1 := RandomVariable(zdist)

Z2 := RandomVariable(zdist)

PDF(Z1, s)

piecewise(s <= 0, 0, 0 < s and s < 2, -(3/80)*(s-2)^3*(s^2+6*s+4), 2 <= s, 0)

z12pdf := PDF(Z1+Z2, s)

(1/1971200)*piecewise(s <= 0, 0, s <= 2, s^11-220*s^9+1320*s^8+7920*s^7-103488*s^6+147840*s^5+887040*s^4-2365440*s^3+2838528*s, s <= 4, -(s^4+28*s^3+228*s^2+536*s+80)*(-4+s)^7, 4 < s, 0)

int(z12pdf, s = -infinity .. infinity)

1

CDF(Z1+Z2, s)

piecewise(s <= 0, 0, s <= 2, (1/23654400)*s^12-(1/89600)*s^10+(1/13440)*s^9+(9/17920)*s^8-(3/400)*s^7+(1/80)*s^6+(9/100)*s^5-(3/10)*s^4+(18/25)*s^2, s <= 4, -(1/23654400)*s^12+(1/89600)*s^10-(1/13440)*s^9-(9/17920)*s^8+(3/400)*s^7-(1/40)*s^6-(9/175)*s^5+(3/5)*s^4-(176/105)*s^3+(288/175)*s^2+(256/385)*s-53/75, 4 < s, 1)

Mean(Z1+Z2)

36/35

Mean(Z1)

18/35

PDF(Z1+Z2, 3/2); evalf(%)

1781183979/4037017600

.4412128347

PDF(Z1+Z2, 3/2, numeric)

.4412128347

Download SphereFinal_ac.mw

The Help-page for topic procedure is describing technical aspects of a fundamental programming component of the Maple language, and the use of the word variable there is intentional jargon related to the Computer Science context of programming languages.

In contrast, in many mathematical contexts (for example, differential equations, integration, optimization, etc) the word variable has quite another set of connotations. And several other Maple's Help-pages reflect that, so as to be more broadly and readily understandable.

So, a word like "variable" can sometimes be used in a technical C.S. way when describing part of the Maple programming language, and can sometimes be used in a mathematical or general way.

One way is to use assumptions on t.
 

with(Statistics)

x := RandomVariable(Uniform(1, 2))

y := RandomVariable(Uniform(2, 3))

`assuming`([PDF(x*y, t)], [t >= 2, t <= 6])

piecewise(t <= 3, -ln(2)+ln(t), t <= 4, ln(3)-ln(2), 4 < t, ln(3)-ln(t)+ln(2))

`assuming`([PDF(x*y, t)], [t >= 2])

piecewise(t <= 3, -ln(2)+ln(t), t <= 4, ln(3)-ln(2), t <= 6, ln(3)-ln(t)+ln(2), 6 < t, 0)

`assuming`([PDF(x*y, t)], [t >= 0])

piecewise(t <= 2, 0, t <= 3, -ln(2)+ln(t), t <= 4, ln(3)-ln(2), t <= 6, ln(3)-ln(t)+ln(2), 6 < t, 0)

Download MultiplyPDFFunctions_ac.mw

So, you must be using 2D Input mode (the default), for which (if done first),

   f(x) := foo

gets parsed as an operator assignment, ie, just as if it had been entered as,

   f := x-> foo

In 1D Maple notation the syntax f(x):=foo is interpreted always as an
assignment to the remember table of the operator/procedure.

Internally, the system does subsequent 2D Input interpretations according
to whichever syntax you've used first (after restart). If you first use the
syntax f:=x->foo in 2D Input mode then subsequently the syntax
f(x):=foo gets interpreted as a remember table assignment.

Many releases ago the syntax f(x):=foo would elicit a disambiguation
popup dialog box, asking the user to specify which interpretation was
intended. Nowadays it tries to to distinguish according to which syntax
was already/first used.

restart

"f(x):=2*x"

proc (x) options operator, arrow, function_assign; 2*x end proc

eval(f)

proc (x) options operator, arrow, function_assign; 2*x end proc


The 2D parser now "knows" that you are using the syntax to mean
operator-assignment, and will also interpret additonal instances to
be operator-assignment (unless you restart).

"f(x):=sin(x) "

proc (x) options operator, arrow, function_assign; sin(x) end proc

eval(f)

proc (x) options operator, arrow, function_assign; sin(x) end proc

f(16)

sin(16)

f(x)

sin(x)

restart

f := proc (x) options operator, arrow; 2*x end proc

proc (x) options operator, arrow; 2*x end proc


The 2D parser will now interpret the following as an assignment
to the remember table of the operator assigned to f. (It'll stop doing
that if you restart.)

"f(x):=sin(x) "

sin(x)


But the body of the operator remains the same.

eval(f)

proc (x) options operator, arrow; table( [( x ) = sin(x) ] ) 2*x end proc


If you pass something else (that is not exactly the name x) then
the formula in the operator will get used to compute the result.

f(y)

2*y

f(8)

16

The following is done by looking up the entry for the name x in the
remember table of f. The result is obtained via lookup, not by using
the formula in the operator assigned to f.

f(x)

sin(x)


We can even examine the remember table of the operator/procedure.

op(4, eval(f))

table( [( x ) = sin(x) ] )

Here, plot receives sin(x) after f(x) gets evaluated
plot(f(x), size = [300, 200])

Here, plot receives 2*y after f(y) gets evaluated
plot(f(y), size = [300, 200])

 

 

Download 2d_oper_assignment.mw

restart;

T := table([([-.687272727272727, -.687272727272727, 1.], [2.25818181818182, 2.258181\
81818182, 1.])=COLOR(RGB,.63529412,.60392157,.60000000),([-2.47, -1.47, 2.], [.\
53, 1.53, 2.])=COLOR(RGB,.60784314,.47450980,.46274510),([.245, -2.755, -2.], [
2.745, -.255, -2.])=COLOR(RGB,.44705882,.52549020,.67058824),([-2.6, -.6, 3.],
[.379, 2.379, 3.])=COLOR(RGB,.57647059,.34509804,.32549020),([2.659, -1.341, -3\
.], [4.652, .652000000000001, -3.])=COLOR(RGB,.34901961,.46666667,.70980392),([
-3.178, .822000000000001, 5.], [-1.185, 2.815, 5.])=COLOR(RGB,.47058824,0.,.\
54901961e-1),([.672727272727271e-1, -.932727272727273, 0.], [3.06181818181818,
2.06181818181818, 0.])=COLOR(RGB,.62352941,.63529412,.65490196),([2.8895, -2.11\
05, -4.], [4.3895, -.6105, -4.])=COLOR(RGB,.22745098,.40392157,.78039216),([1.0\
1, -.99, -1.], [4.02571428571429, 2.02571428571429, -1.])=COLOR(RGB,.53725490,.\
57647059,.65490196),([-3.63, -.63, 4.], [-1.13875, 1.86125, 4.])=COLOR(RGB,.537\
25490,.20392157,.18823529)]);

table( [( [-2.47, -1.47, 2.], [.53, 1.53, 2.] ) = COLOR(RGB, .60784314, .47450980, .46274510), ( [.245, -2.755, -2.], [2.745, -.255, -2.] ) = COLOR(RGB, .44705882, .52549020, .67058824), ( [-.687272727272727, -.687272727272727, 1.], [2.25818181818182, 2.25818181818182, 1.] ) = COLOR(RGB, .63529412, .60392157, .60000000), ( [-3.178, .822000000000001, 5.], [-1.185, 2.815, 5.] ) = COLOR(RGB, .47058824, 0., 0.54901961e-1), ( [-2.6, -.6, 3.], [.379, 2.379, 3.] ) = COLOR(RGB, .57647059, .34509804, .32549020), ( [2.659, -1.341, -3.], [4.652, .652000000000001, -3.] ) = COLOR(RGB, .34901961, .46666667, .70980392), ( [-3.63, -.63, 4.], [-1.13875, 1.86125, 4.] ) = COLOR(RGB, .53725490, .20392157, .18823529), ( [1.01, -.99, -1.], [4.02571428571429, 2.02571428571429, -1.] ) = COLOR(RGB, .53725490, .57647059, .65490196), ( [2.8895, -2.1105, -4.], [4.3895, -.6105, -4.] ) = COLOR(RGB, .22745098, .40392157, .78039216), ( [0.672727272727271e-1, -.932727272727273, 0.], [3.06181818181818, 2.06181818181818, 0.] ) = COLOR(RGB, .62352941, .63529412, .65490196) ] )

new := table([map(tt->tt[1]=T[tt[]],[indices(T)])[]]);

table( [( [-3.178, .822000000000001, 5.] ) = COLOR(RGB, .47058824, 0., 0.54901961e-1), ( [0.672727272727271e-1, -.932727272727273, 0.] ) = COLOR(RGB, .62352941, .63529412, .65490196), ( [2.8895, -2.1105, -4.] ) = COLOR(RGB, .22745098, .40392157, .78039216), ( [1.01, -.99, -1.] ) = COLOR(RGB, .53725490, .57647059, .65490196), ( [-.687272727272727, -.687272727272727, 1.] ) = COLOR(RGB, .63529412, .60392157, .60000000), ( [-2.47, -1.47, 2.] ) = COLOR(RGB, .60784314, .47450980, .46274510), ( [.245, -2.755, -2.] ) = COLOR(RGB, .44705882, .52549020, .67058824), ( [-3.63, -.63, 4.] ) = COLOR(RGB, .53725490, .20392157, .18823529), ( [-2.6, -.6, 3.] ) = COLOR(RGB, .57647059, .34509804, .32549020), ( [2.659, -1.341, -3.] ) = COLOR(RGB, .34901961, .46666667, .70980392) ] )

map(print,sort([indices(T)])):

[[-3.63, -.63, 4.], [-1.13875, 1.86125, 4.]]

[[-3.178, .822000000000001, 5.], [-1.185, 2.815, 5.]]

[[-2.6, -.6, 3.], [.379, 2.379, 3.]]

[[-2.47, -1.47, 2.], [.53, 1.53, 2.]]

[[-.687272727272727, -.687272727272727, 1.], [2.25818181818182, 2.25818181818182, 1.]]

[[0.672727272727271e-1, -.932727272727273, 0.], [3.06181818181818, 2.06181818181818, 0.]]

[[.245, -2.755, -2.], [2.745, -.255, -2.]]

[[1.01, -.99, -1.], [4.02571428571429, 2.02571428571429, -1.]]

[[2.659, -1.341, -3.], [4.652, .652000000000001, -3.]]

[[2.8895, -2.1105, -4.], [4.3895, -.6105, -4.]]

map(print,sort([indices(new)])):

[[-3.63, -.63, 4.]]

[[-3.178, .822000000000001, 5.]]

[[-2.6, -.6, 3.]]

[[-2.47, -1.47, 2.]]

[[-.687272727272727, -.687272727272727, 1.]]

[[0.672727272727271e-1, -.932727272727273, 0.]]

[[.245, -2.755, -2.]]

[[1.01, -.99, -1.]]

[[2.659, -1.341, -3.]]

[[2.8895, -2.1105, -4.]]


Download nm_table_q1.mw

Maple uses a extension mechanism for some kernel builtin commands like expand. This recognizes procedures with a name in a pattern like `expand/f` to be the thing to use when asked to expand a call to f.

There is a default Library procedure assigned to `expand/Sum`. That is a somewhat complicated procedure that tries a bit to respect the notion that expanding (splitting) a Sum is not always mathematically valid.

When you load a command from the old student package (using with, like you've done) Maple does some initialization of the package and redefines some things. In particular, it redefines `expand/Sum`.

It actually redefines it as a call to readlib of the new, replacement procedure. That's not key. What's key is that the initialization assigns another, simpler procedure to `expand/Sum`. That new procedure turns a Sum call of an addition of terms into an addition of Sum calls, regardless of the fact that is not always right. I suspect that the rationale was the students would find the functionality useful, and were less likely to be working with divergent sums (infinite series), etc, or that the benefit of the functionality outweighed the danger.

You can compare these:

restart;
interface(verboseproc=3):
with(student,changevar);
showstat(eval(`expand/Sum`));

restart;
interface(verboseproc=3):
showstat(`expand/Sum`);

Here are various corrections and adjustments, including:

- Utilize the nolines option to plots:-inequal, to get rid of dotted line demarking part of one strict inequality.
- Utilize InertForm:-Display for (right tick) unevaluation-quoted expressions used in the plots. Otherwise, subsequent re-use or evaluation of the created plots will cause those nicely formatted names in inequalities to evaluate to lengthy expressions. NB. the typeset option does not(!) prevent that. You can test this worked by just assigning the plot to a name, and then printing that after the fact; it's working properly if the inlined math in the plot survives.
- Get rid of some unnecessary quotes.
- Pass gridrefine=4 on to implicitplot, called from plots:-inequal, to get smoother edges of polygon regions.

restart

with(Optimization); with(plots); with(Student[VectorCalculus]); with(LinearAlgebra)

DATA := [delta = .9, a = 0.1e-1, g = .25, c = 0.5e-1, sigma = 0.8e-1, Cn = .4, Crm = .1, Cr = 0.1e-1, Pr = .35, t = 0.5e-1, alpha = .95, s = 0.1e-1, upsilon = .95, rho0 = .4, b = .5]

`&pi;_WD` := (((alpha*((-g*i2 + a)*Cr + 2*Crm + 2*c - 2*Pr)*Cr*b + ((g*i2 - a)*Cr - 2*Crm - 2*c + 2*Pr)*alpha - (-g*i2 + a)*Cr)*d + alpha*((-g*i2 + a)*Cr + 2*Crm + 2*c - 2*Pr)*b + 2*g*i2 - 2*a)*rho0 + 2*((Cr*b - 1)*d + b)*(delta + Cn - Pr - 1))^2*b*d/(8*(Cr*b*d + b - d)^2*(((Cr*alpha*b - alpha + 1)*d + alpha*b)*rho0^2 - 2*b*d*(delta - 1))):

`&pi;_D` := b*(Cr*rho0*t*(Cr*alpha*b-alpha-1)*d^2+((alpha*((-g*i2+a)*Cr+2*Crm+2*c+3*t-2*Pr)*Cr*b+((g*i2-a)*Cr-2*Crm-2*c-2*t+2*Pr)*alpha+(g*i2-a)*Cr-2*t)*rho0+(2*(Cr*b-1))*(sigma*t+Cn-Pr+delta-1))*d+(alpha*((-g*i2+a)*Cr+2*Crm+2*c+2*t-2*Pr)*b+2*g*i2-2*a)*rho0+2*b*(sigma*t+Cn-Pr+delta-1))^2*d/(8*((Cr*b-1)*d+b)^2*(((Cr*alpha*b-alpha+1)*rho0^2-2*b*(delta-1))*d+alpha*b*rho0^2))

`&pi;R` := subs(DATA, `&pi;_WD`)

`&pi;D` := subs(DATA, `&pi;_D`)

Di := `&pi;R`-`&pi;D`

# Step 1: try to figure out the scales of the domain where to use plots:-inequal
(*
simplify(`&pi;R`-`&pi;D`);
S := solve(%=0, Ce):

eval([S], i2=0.01);
eval(`&pi;R`-`&pi;D`, {i2=0.01, d=0.8});
eval(`&pi;R`-`&pi;D`, {i2=0.01, d=0.3});

eval([S], i2=0.09);
eval(`&pi;R`-`&pi;D`, {i2=0.09, d=0.8});
eval(`&pi;R`-`&pi;D`, {i2=0.09, d=0.3});
*)

LX := 0.09: LY := 0.8:

# Plot region where πR > πD
p := plots:-display(
  plots:-inequal(`&pi;R` > `&pi;D`, i2 = 0 .. LX,
                 d = 0 .. LY, ':-color' = "Chartreuse",
                 ':-nolines', ':-optionsimplicit'=[':-gridrefine'=4]),
  axes=boxed):

display(p,plot(x, x=0..0.06, color="Chartreuse", thickness=10,
               legend=InertForm:-Display('`&pi;R`' > '`&pi;D`'),
        view = [0 .. LX, 0 .. LY],labels = ['i2', 'd']));

eps  := 0.001:

OneHatch := display(plottools:-polygon(
                [[-2*LX, -2*LX-eps], [2*LX, 2*LX-eps], [2*LX, 2*LX+eps],
                 [-2*LX, -2*LX-eps], [-2*LX, -2*LX-eps]],
                color=gray, style=line)):

Tr := plottools:-transform((x, y) -> [x, y*LY/LX]):

back := display(seq(plottools:-translate(Tr(OneHatch), 0, h),
                    h=-2*LY..2*LY, LY/10),
                view=[-LX..LX, -LY..LY], axes=framed):

display(p
  , back
  , plot(x, x = -LX .. -(1/2)*LX, color = "Chartreuse", thickness = 10,
         legend = InertForm:-Display('`&pi;R` > `&pi;D`'))
  , title=typeset("Hatched region =   ", InertForm:-Display('`&pi;R`' < '`&pi;D`'))
  , view = [0 .. LX, 0.3 .. LY]);

display(plottools:-rectangle([0.01, 0.65], [0.03, 0.7], color=white)
  , textplot([0.02, 0.675, Pi[m]^WD > Pi[m]^:-D], font=[Times, bold, 10])
  , p
  , back
  , plottools:-rectangle([0.01, 0.4], [0.03, 0.45], color=white),
    textplot([0.02, 0.435, InertForm:-Display('`&pi;WD` < `&pi;D`')]
  , font=[Times, bold, 10])
  , view = [0 .. LX, 0.3 .. LY]);

 

 

Download Question_1_regional_ac.mw

Here is a further revision of Carl's revision of Kitonum's procedure.

The main point of this revision is to allow the same kinds of color functionality (for the curves and the text contour labels) as contourplot has for its curves.

(Aside: If you have an old version which does not support the colorbar option -- on densityplot or contourplot -- then you could just remove those from the code. The code below references it just to disable it.)

[edit] The code immediately below has been re-edited and corrected/improved, related to Replies below this Answer.

restart;

kernelopts(version);

`Maple 2024.2, X86 64 LINUX, Oct 29 2024, Build ID 1872373`

ContoursWithLabels := proc(
     Expr::algebraic,
     Range1::(name=range(realcons)), Range2::(name=range(realcons)),
     {contours::{posint, {list}(realcons)}:=8},
     {GraphicOptions::{list,set}({name, name= anything}):= NULL},
     {ImplicitplotOptions::{list,set}({name, name=anything}):=NULL},
     {TextOptions::{list,set}({name, name=anything}):=NULL},
     {ContourOptions::{list,set}({name, name=anything}):=NULL},
     {DensityOptions::{list,set}({name, name=anything}):=NULL} )
   local r1,r2,f,L1,h,S1,P,r,M,C,T,p,p1,m,n,i,G,GL,mm,PO,data,
         x:=lhs(Range1),y:=lhs(Range2);
     f := unapply(Expr,(x,y));
     if contours::posint then
          r1 := rand(convert(rhs(Range1), float));
          r2 := rand(convert(rhs(Range2), float));
          L1 := select(type,(f@op)~({seq([r1,r2](),i=1..205)}),realcons);
          h := (L1[-6]-L1[1])/contours;
          S1 := [seq(L1[1]+h/2+h*(n-1), n=1..contours)]
     else #contours::{set,list}(realcons)
          S1 := [contours[]]
     end if;
     userinfo(1, ContoursWithLabels, print('Contours' = evalf[2](S1)), `\n`);
     r := k-> rand(20..k-20);
     G := plots:-contourplot(Expr,Range1,Range2,':-contours'=S1, ':-grid'=[5,5],
                             ':-contourlabels'=false,':-colorbar'=false,ContourOptions[]);
     GL := ListTools:-Reverse(select(u->type(u,specfunc(CURVES)) and nops(u)>0,[op(G)]));
     for C from 1 to nops(S1) do
          P := plots:-implicitplot(Expr = S1[C], Range1, Range2,
                                   gridrefine=3, ImplicitplotOptions[]);
          PO := select(type,[op(GL[C])],specfunc({COLOR,COLOUR}))[-1];
          data := [plottools:-getdata(P)];
          if data = [] then next; end if;
          for p in data do
               p1 := convert(p[3], listlist); n := nops(p1);
               if n < 500 then
                    m := `if`(40 < n, r(n)(), round(n/2));
                    M[`if`(40 < n, [p1[1..m-11], p1[m+11..n]], [p1])[]] := PO;
                    T[[p1[m][], evalf[2](S1[C])]] := PO;
               else
                    h := trunc(n/2); m := r(h)();
                    M[p1[1..m-11], p1[m+11..m+h-11], p1[m+h+11..n]] := PO;
                    T[[p1[m][], evalf[2](S1[C])], [p1[m+h][], evalf[2](S1[C])]] := PO;
               end if;
          end do;
     end do;
     `if`( ( (not M::table) or nops([indices(M)])=0 )
           and ( (not T::table) or nops([indices(T)])=0 )
           and DensityOptions=NULL,
          PLOT(),
          plots:-display([
          `if`(DensityOptions = NULL,NULL,
               plots:-densityplot(Expr, Range1, Range2, ':-colorbar'=false,
                                  ':-style'=':-surface', DensityOptions[])),
           `if`(M::table,seq(plot(mm,':-color'=M[mm[]], GraphicOptions[]),mm=[indices(M)]),NULL),
           `if`(T::table,seq(plots:-textplot(mm,':-color'=T[mm[]], TextOptions[]),mm=[indices(T)]),NULL)
                    ], 'axes'=':-box', 'gridlines'=false, _rest)
        );
end proc:

 

ff := (11-x-y)^2+(1+x+10*y-x*y)^2;

(11-x-y)^2+(-x*y+x+10*y+1)^2

ContoursWithLabels(
     ff, x=0..20, y=0..15, contours=11,
     TextOptions = [font=[HELVETICA,BOLD,8], color="DarkBlue"],
     ContourOptions = [colorscheme="DivergeRainbow"],
     ImplicitplotOptions = [gridrefine=4],
     size = [900,500], labels=[x,y]);

ContoursWithLabels(
     ff, x=0..20, y=0..15, contours=[seq(50+i*50,i=0..10),seq(1000+i*1100,i=0..8)],
     TextOptions = [font=[HELVETICA,BOLD,8]],
     ContourOptions = [colorscheme="DivergeRainbow"],
     ImplicitplotOptions = [gridrefine=3],
     size = [900,500], labels=[x,y]);

PP := 0.3800179925e-3*exp(-0.6065722618e-3*(x-29.51704536)^2
      +(0.6650093594e-3*(x-29.51704536))*(a-12.94061928)-0.1106850312e-2*(a-12.94061928)^2);

0.3800179925e-3*exp(-0.6065722618e-3*(x-29.51704536)^2+0.6650093594e-3*(x-29.51704536)*(a-12.94061928)-0.1106850312e-2*(a-12.94061928)^2)

ContoursWithLabels(
     PP, x=-20..20, a=-20..20, contours=[seq(1e-4..4e-4,5e-5)],
     TextOptions = [font=[HELVETICA,BOLD,8]],
     ContourOptions = [colorscheme="DivergeRainbow"],
     labelfont = [TIMES,BOLDITALIC,12], axesfont = [HELVETICA,10],
     size = [500,500]);

ContoursWithLabels(
     PP, x=-20..20, a=-20..20, contours=[seq(1e-4..4e-4,5e-5)],
     TextOptions = [font=[HELVETICA,8]],
     ContourOptions = [color=black],
     GraphicOptions = [thickness=0.5],
     labelfont = [TIMES,BOLDITALIC,12], axesfont = [HELVETICA,10],
     size = [500,500]);

ContoursWithLabels(
     PP, x=-20..20, a=-20..20, contours=[seq(1e-4..4e-4,5e-5)]);

ContoursWithLabels(
     PP, x=-20..20, a=-20..20, contours=[seq(1e-4..4e-4,5e-5)],
     DensityOptions = [colorscheme="DivergeRainbow", style=surface],
     TextOptions = [font=[HELVETICA,BOLD,8]],
     ImplicitplotOptions = [gridrefine=2],
     ContourOptions = [color=black],
     labelfont = [TIMES,BOLDITALIC,12], axesfont = [HELVETICA,10],
     size = [500,500]);

ContoursWithLabels(
     PP, x=-20..20, a=-20..20, contours=[seq(1e-4..4e-4,5e-5)],
     DensityOptions = [colorscheme=["zgradient",["Purple","LightBlue","Orange"]],colorbar=true],
     TextOptions = [font=[HELVETICA,BOLD,9],color="#303030"],
     ContourOptions = [color="#303030"],
     labelfont = [TIMES,BOLDITALIC,12], axesfont = [HELVETICA,10],
     size = [700,400]);

ContoursWithLabels(
     PP, x=-20..20, a=-20..20, contours=[seq(1e-4..4e-4,5e-5)],
     DensityOptions = [colorscheme=[black,white]],
     TextOptions = [font=[HELVETICA,BOLD,10]],
     ContourOptions = [colorscheme=[white,"LightGray",black]],
     GraphicOptions = [thickness=0.5],
     labelfont = [TIMES,BOLDITALIC,12], axesfont = [HELVETICA,10],
     size = [700,400]);

ContoursWithLabels(
     y^2+x-1, x=-10..10, y=-10..10, contours=[$-1 .. 5],
     TextOptions = [font=[HELVETICA,BOLD,7], color="DarkBlue"],
     ContourOptions = [colorscheme="DivergeRainbow"]);

ContoursWithLabels(
     y^2+x-1, x=-10..10, y=-10..10, contours=[$-1 .. 5],
     TextOptions = [font=[HELVETICA,BOLD,10]],
     ContourOptions = [colorscheme="DivergeRainbow"]);

Download CPL05a.mw

ps. I notice that my M2024.2 has trouble rendering the minus sign in a textplot numeral/float, when the font size is 8 or 9. Not central to this code, but I'll report it.

I'm not suggesting that this is a good choice of template function. It's merely an illustration, a bit like you had.

Not also that NonlinearFit might do local optimization, so quality of the fit can depend on the choice of initialvalues. See also the DirectSearch v.2 package on the Application Center, which can do similar kinds of fitting but doing global optimization over supplied parameter ranges.

notes:
- The randomize() call makes it get a different L after each restart.
- The float range -5.0..5.0 makes it produce floats in the range. A range -5..5 would just produce integer values in the range.

restart;

randomize():

with(plots):

setoptions(size=[500,200]);

N := 10:

L:=LinearAlgebra:-RandomVector(N,generator=-5.0..5.0):

listplot(L);

F1 := Statistics:-NonlinearFit( a*cos(e*x+d)+b*x+c, <$1..N>, L, x
                                , parameterranges=[a=-5..5, d=-Pi..Pi]
                                , initialvalues=[e=ceil(1/2*N/Pi)]
                              );

HFloat(2.946500373957529)*cos(HFloat(1.660098536783793)*x+HFloat(1.972035777846103))+HFloat(0.03533616937042544)*x-HFloat(0.4678613378149076)

display(listplot(L), plot(F1,x=1..N));

Download Mma_trans_01.mw

1 2 3 4 5 6 7 Last Page 1 of 342