acer

29899 Reputation

29 Badges

18 years, 203 days
Ontario, Canada

Social Networks and Content at Maplesoft.com

MaplePrimes Activity


These are answers submitted by acer

Here's one way to handle an example like that.

The steps are to select the indices in which you're interested, and then to map table-lookup across those indices, and then take the min of the resulting values.

f := rand(1..100):

T := table([seq([R||i=f(),a||i=f(),b||i=f(),b||i=f()][], i=1..3)]);

table( [( R2 ) = 98, ( a1 ) = 45, ( R3 ) = 38, ( R1 ) = 93, ( a2 ) = 59, ( b3 ) = 96, ( a3 ) = 69, ( b2 ) = 100, ( b1 ) = 6 ] )

 

min(map[2](`?[]`,T,select(hastype,[indices(T)],suffixed(a))));

45

 

 

The steps of that are,

 

select(hastype,[indices(T)],suffixed(a));

[[a1], [a2], [a3]]

map[2](`?[]`,T,%)

[45, 59, 69]

min(%);

45


Download anthr_tab_sel2.mw

It's not clear whether you're asking about how to programmatically determine the version of Maple in which code gets run, or whether you're asking about how to determine the version in which a worksheet was last saved.

I'm guessing it's the former.

The command kernelopts(version) returns a name that contains the version number. That could be converted to a string.

You can chop up that string to take only a particular portion, eg.

str:=convert(kernelopts(version),string);

   str := "Maple 2018.2, X86 64 LINUX, Nov 16 2018, Build ID 1362973"

str[7 .. StringTools:-Search(".",str)-1];

            "2018"

You could use either InertForm:-Display or (undocumented) Typesetting:-Typeset for this.

restart;

 

R := 10:

 

plots:-textplot([600,-Pi*floor(R)*ln(R)-0.01,
                 InertForm:-Display(-Pi*floor(r)*ln(r)),
                 font=[Helvitica,bold,14]],labels=[T," "]);

plots:-textplot([600,-Pi*floor(R)*ln(R)-0.01,
                 Typesetting:-Typeset(-Pi*floor(r)*ln(r)),
                 font=[Helvitica,bold,14]],labels=[T," "]);


Download textplot_ts_ex1.mw

Your given example can be handled by the define command.

restart;

 

define(T,T(a::nonunit(algebraic)+b::nonunit(algebraic))=T(a)+T(b),
         T(a::nonunit(algebraic)*x^n::integer)=a*T(x^n),
         'conditional'(T(a::anything)=a*T(1),
                       _type(a,And(Not(identical(1)),freeof(x)))));

 

T(alpha__1*x^2*y+alpha__2*x^4*t+alpha__3*t*y);

t*alpha__2*T(x^4)+y*alpha__1*T(x^2)+alpha__3*t*y*T(1)


Download define_ex.mw

Is this the kind of effect you're after?

Note that there's no need for any isolation/resubstitution to attain the first target form: collect is enough.

Note a sign difference in the second term of the target form.

The same algsubs result can also be obtained directly from the original expression (expanded), without the intermediate target form.

restart;

expr := 3*G*(`Δγ`*H - `σy`(`Δγ`) + q)/(-q + `σy`(`Δγ`))^2;

3*G*(`Δγ`*H-`σy`(`Δγ`)+q)/(-q+`σy`(`Δγ`))^2

targ := collect(expr,H,simplify);

3*G*`Δγ`*H/(-q+`σy`(`Δγ`))^2-3*G/(-q+`σy`(`Δγ`))

RR := G*`Δγ`=y/3*(q - `σy`(`Δγ`));

G*`Δγ` = (1/3)*y*(-`σy`(`Δγ`)+q)

algsubs(RR, targ);

-3*G/(-q+`σy`(`Δγ`))-H*y/(-q+`σy`(`Δγ`))

ans := simplify(%);

(H*y+3*G)/(-`σy`(`Δγ`)+q)

algsubs(RR, expand(expr));

-3*G/(-q+`σy`(`Δγ`))-H*y/(-q+`σy`(`Δγ`))

simplify(%);

(H*y+3*G)/(-`σy`(`Δγ`)+q)

sort(ans, q);

(H*y+3*G)/(q-`σy`(`Δγ`))


Download ceeeb_ex.mw

[edit] Even though your intermediate target form is not necessary to achieve the algsubs result, it could also be obtained with,

   convert(expr,parfrac,q);
   sort(%,order=plex(H,q));   # optional

For some values of parameter Cl the dsolve command may here return a result containing an inert integral, rather than a fully explicit answer.

In that case you could apply the value command. (This aspect of dsolve is mentioned in a bullet point a section of the ?dsolve,details Help page.)

Eg,

     plot(value(sol(t)), t = 0 .. 48)

dsH_ex.mw

note: There are other ways to deal with it (changing the Int integral to use some specific numeric methods, or doing a numeric dsolve, etc.) Let us know if you'd like to see such.

@jrive First I'll show a form similar to that in nm's Answer, and discuss more on how that is not equivalent to the original even if all parameters are real. (See also nm's pertinent note, "sqrt(a)/b is same as sqrt(a/b^2) assuming b>0", by which he qualifies his answer.)

restart;

sol := (-v + sqrt(-4*a^2*R^2 + v^2))/(2*a*omega*L);

(1/2)*(-v+(-4*R^2*a^2+v^2)^(1/2))/(a*omega*L)

alt := evalindets(expand(sol),`*`,
                  u->`if`(membertype(anything^(1/2),u),sqrt(u^2),u));

-(1/2)*v/(a*omega*L)+(1/2)*((-4*R^2*a^2+v^2)/(a^2*omega^2*L^2))^(1/2)


Example for which they're not equal.

eval([sol,alt], [R=1,v=1,a=1,L=1,omega=-1]);

[1/2-(1/2)*(-3)^(1/2), 1/2+(1/2)*(-3)^(1/2)]

combine(simplify( sol-alt )) assuming real, a*omega*L>0;

0

combine(simplify( sol-alt )) assuming real, a*omega*L<0;

(-4*R^2*a^2+v^2)^(1/2)/(a*omega*L)

Download rad_comb_ex.mw

Regarding your comment about bringing a term into a radical, consider this simpler example,

ee := sqrt(x^2)/x;

(x^2)^(1/2)/x

simplify(ee) assuming x>0;

1

simplify(ee) assuming x<0;

-1


If we absorb that x^1 into the radical (by squaring it) and make no other changes then we've mistakenly lost the multiplicative sign of x. If x is positive then it's ok. But if x is negative then we've mistakenly lost the negative sign.

With that in mind, here is a variant on your original query, with a form which retains the signum of a*omega*L outside the radical. Notice the broader equivalence, compared to your original target.

restart;

sol := (-v + sqrt(-4*a^2*R^2 + v^2))/(2*a*omega*L);

(1/2)*(-v+(-4*R^2*a^2+v^2)^(1/2))/(a*omega*L)

alt2 := -v/(2*a*omega*L) + sqrt((-4*R^2*a^2 + v^2)/(a^2*omega^2*L^2))/(2*signum(a*omega*L));

-(1/2)*v/(a*omega*L)+(1/2)*((-4*R^2*a^2+v^2)/(a^2*omega^2*L^2))^(1/2)/signum(a*omega*L)

combine(simplify(sol - alt2)) assuming real;

0

alt2 assuming a*omega*L>0;

-(1/2)*v/(a*omega*L)+(1/2)*((-4*R^2*a^2+v^2)/(a^2*omega^2*L^2))^(1/2)

alt2 assuming a*omega*L<0;

-(1/2)*v/(a*omega*L)-(1/2)*((-4*R^2*a^2+v^2)/(a^2*omega^2*L^2))^(1/2)

Download rad_comb_ex3.mw


Note: If, say, x=-3 then by convention sqrt(x^2) taken to be 3, not -3. That positive answer is known as the principal square root. See also the third paragraph here. Or see notes on the square-root case on principal value or n-th roots.

Is this the kind of thing you're trying to accomplish?

(I made no effort to optimize, or figure out what you're doing with tau at the end.)

I'll find a value for parameter beta that attains the goal g(10)=1. I'll do this in two ways, first using my own solving procedure and second using Doug Meade's Shoot package (see URL below).

restart

Digits:=15:

with(plots):

alias( U=u(x,y), V=v(x,y) ):

PDE:={ U*diff(U,x)+V*diff(U,y)-nu*diff(U,y$2)=0,

       diff(U,x)+diff(V,y)=0 }:

simsubs := eta(x, y) = y*sqrt((1/2)*u[0]/(nu*x)):

stream := psi(x, y) = sqrt(2*nu*x*u[0])*f(eta(x, y)):

Usubs := U = diff(rhs(stream), y):

Vsubs := V = -(diff(rhs(stream), x)):

ODE := simplify(subs(Usubs, Vsubs, simsubs, PDE)):

simsubs2 := solve(subs(eta(x, y) = eta, simsubs), {y}):

ODE := simplify(subs(simsubs2, ODE), symbolic):

FNS := {f(eta), g(eta), h(eta)};

{f(eta), g(eta), h(eta)}

ODE:={ diff(f(eta),eta)=g(eta),

       diff(g(eta),eta)=h(eta),

       diff(h(eta),eta)=-f(eta)*h(eta) }:

IC:={ f(0)=0, g(0)=0, h(0)=beta }:

BC:=g(10)=1:

 

findbeta := proc(b) local tt;
  if not b::numeric then return 'procname'(b); end if;
  tt := dsolve(eval(ODE union IC, beta=b),{f(eta),g(eta),h(eta)},numeric):
  eval(g(eta),tt(10));
end proc:

betaval := beta = fsolve(findbeta(beta)-1, beta=0..1);

beta = .469600090249286

Sf := dsolve(eval(ODE union IC, betaval),{f(eta),g(eta),h(eta)},numeric,output=listprocedure):


Check that g(10)=1 for this solution, computed using betaval.

eval(g(eta), Sf)(10);

HFloat(1.0000000000000007)

odeplot( Sf, [ [eta,f(eta)],

            [eta,g(eta)], [eta,h(eta)] ],0..4,
            labels=[`eta`,``], legend=["f","f '","f ''"],
            title="Blasius solution for flat-plate boundary layer" );

fp0:=subs( Sf, g(eta) ):

fp := proc(x)

  if not type(evalf(x),`numeric`) then

    'procname'(x);

  else

    fp0(x);

  end if;

end proc:

eta[`99%`] = fsolve(fp(eta) = .99, eta = 3 .. 4);

eta[`99%`] = 3.47188559058681

delta[`99%`] = subs(eta = rhs(%), combine(rhs(op(simsubs2))));

delta[`99%`] = 3.47188559058681*2^(1/2)/(u[0]/(nu*x))^(1/2)

tau := mu*simplify(diff(subs(simsubs, rhs(Usubs)), y));

(1/2)*mu*(nu*x*u[0])^(1/2)*((D@@2)(f))((1/2)*y*2^(1/2)*(u[0]/(nu*x))^(1/2))*2^(1/2)*u[0]/(nu*x)

tau[w] := subs(((D@@2)(f))(0) = (subs(Sf, h(eta)))(0), subs(y = 0, tau));

HFloat(0.234800045124643)*mu*(nu*x*u[0])^(1/2)*2^(1/2)*u[0]/(nu*x)


Now do the same thing using,
    http://www.math.sc.edu/~meade/maple/Shoot9/Shoot9.zip

libname:=cat(kernelopts(homedir),"/mapleprimes/Shoot/Shoot9"),libname:
with(Shoot);

[shoot]

infolevel[shoot] := 1:

S:=shoot( ODE, IC, BC, FNS, beta=1,

          abserr=Float(5,-7), output=listprocedure, method=taylorseries ):

shoot: Step #  1

shoot: Parameter values :  beta = 1
shoot: Step #  2
shoot: Parameter values :  beta = HFloat(0.4062401739572171)
shoot: Step #  3
shoot: Parameter values :  beta = HFloat(0.46805773973620646)

shoot: Step #  4
shoot: Parameter values :  beta = HFloat(0.4695991426488529)
shoot: Step #  5
shoot: Parameter values :  beta = HFloat(0.46959998836128514)


Check that g(10)=1 for this solution (approximately), computed using betaval.
(You could vary the options to shoot.)

eval(g(eta), S)(10);

.999999999972365

 

 

odeplot( S, [ [eta,f(eta)],

            [eta,g(eta)], [eta,h(eta)] ],0..4,
            labels=[`eta`,``], legend=["f","f '","f ''"],
            title="Blasius solution for flat-plate boundary layer" );

fp0:=subs( S, g(eta) ):

fp := proc(x)

  if not type(evalf(x),`numeric`) then

    'procname'(x);

  else

    fp0(x);

  end if;

end proc:

eta[`99%`] = fsolve(fp(eta) = .99, eta = 3 .. 4);

eta[`99%`] = 3.47188688127711

delta[`99%`] = subs(eta = rhs(%), combine(rhs(op(simsubs2))));

delta[`99%`] = 3.47188688127711*2^(1/2)/(u[0]/(nu*x))^(1/2)

tau := mu*simplify(diff(subs(simsubs, rhs(Usubs)), y));

(1/2)*mu*(nu*x*u[0])^(1/2)*((D@@2)(f))((1/2)*y*2^(1/2)*(u[0]/(nu*x))^(1/2))*2^(1/2)*u[0]/(nu*x)

tau[w] := subs(((D@@2)(f))(0) = (subs(S, h(eta)))(0), subs(y = 0, tau));

HFloat(0.23479999418064257)*mu*(nu*x*u[0])^(1/2)*2^(1/2)*u[0]/(nu*x)

NULL

Shoot_Blasius_solution_ac.mw

Is this the kind of effect you're after?

[edit:The solution(s) in terms of explicit radicals can be obtained by using the explicit option of the solve command. The simplification to your target form can be achieved via the evala command -- though there are a few longer incantations that also get it.]

restart

NULLeq1 := (R2+Rhi)*R1/(R1+R2+Rhi) = RloNULL

(R2+Rhi)*R1/(R1+R2+Rhi) = Rlo

 

NULLeq2 := R1*Rlo/(R1+Rlo)+R2 = RhiNULL

R1*Rlo/(R1+Rlo)+R2 = Rhi

evala({solve({eq1, eq2}, {R1, R2}, explicit)})

{{R1 = ((Rhi-Rlo)*Rhi)^(1/2)*Rlo/(Rhi-Rlo), R2 = ((Rhi-Rlo)*Rhi)^(1/2)}, {R1 = -((Rhi-Rlo)*Rhi)^(1/2)*Rlo/(Rhi-Rlo), R2 = -((Rhi-Rlo)*Rhi)^(1/2)}}


Download for_help_jrive_ac.mw

Maybe there is something nicer which "understands" the inert `%.` (as well as recognizes the distinct Matrix instances with equal entries as being equivalent), but it's less ad hoc.

ras5_v1_acc.mw

[edit] Here is very a minor revision to that, splitting to show the simplify step separately,

restart; Digits := 14: with(PDEtools): with(LinearAlgebra): with(plots):

pde:= diff(u(t,x),t$2) = diff(u(t,x),x$2):

 Subs:= {r1(t,x) = diff(u(t,x),t), r2(t,x)= diff(u(t,x),x)}:

eq1 := map(z->diff(z,t),<r1(t,x),r2(t,x)>) = Matrix(2,[0,1,1,0]) %.map(z->diff(z,x),<r1(t,x),r2(t,x)>):

#value(subs(Subs,eq1)):

GenerateStencil := proc(F,N,{orientation:=center,stepsize:=s,showorder:=true, showerror:=false}) loc...

eq2:=GenerateStencil(diff(r1(t,x),t),2,orientation=right,stepsize=s,showorder=false):
eq3:=GenerateStencil(diff(r2(t,x),t),2,orientation=right,stepsize=s,showorder=false):
eq4:=GenerateStencil(diff(r1(t,x),x),2,orientation=center,stepsize=h,showorder=false):
eq5:=GenerateStencil(diff(r2(t,x),x),2,orientation=center,stepsize=h,showorder=false):
 

eq6 := <rhs(eq2),rhs(eq3)> = Matrix(2,[0,1,1,0]) %. <rhs(eq4),rhs(eq5)>:

eq7 := [seq(seq(r1(t+i*s,x+j*h)= r1[n+i,m+j],i=0..1),j=-1..1)]:
eq8 :=[seq(seq(r2(t+i*s,x+j*h)= r2 [n+i,m+j],i=0..1),j=-1..1)]:

eq9 := map(z->subs(eq7,eq8,z),eq6):

eq10 := [seq(seq(r1(t+i*s,x+j*h)= r1[n+i,m+j]-E1[n+i,m+j],i=0..1),j=-1..1)]:

eq11 :=[seq(seq(r2(t+i*s,x+j*h)= r2 [n+i,m+j]-E2[n+i,m+j],i=0..1),j=-1..1)]:

eq12 := map(z->subs(eq10,eq11,z),eq6):

eq13 := simplify(eq12-eq9,'mindeg');

(Vector(2, {(1) = (E1[n, m]-E1[n+1, m])/s, (2) = (E2[n, m]-E2[n+1, m])/s})) = `%.`(Matrix(2, 2, {(1, 1) = 0, (1, 2) = 1, (2, 1) = 1, (2, 2) = 0}), Vector(2, {(1) = (1/2)*(-r1[n, m-1]+E1[n, m-1]+r1[n, m+1]-E1[n, m+1])/h, (2) = (1/2)*(-r2[n, m-1]+E2[n, m-1]+r2[n, m+1]-E2[n, m+1])/h}))-`%.`(Matrix(2, 2, {(1, 1) = 0, (1, 2) = 1, (2, 1) = 1, (2, 2) = 0}), Vector(2, {(1) = (1/2)*(-r1[n, m-1]+r1[n, m+1])/h, (2) = (1/2)*(-r2[n, m-1]+r2[n, m+1])/h}))

factor(subs[eval](`%.`=`*`,subsindets(subsindets(eq13,Matrix,
                                                 u->freeze(convert(u,listlist))),
                                      Vector,u->freeze(convert(u,list))))):
eq14 := subsindets(subsindets(subsindets(thaw(%),And(`*`,satisfies(u->ormap(type,u,list))),
                                         `%.`@op),listlist,Matrix),list,Vector);

(Vector(2, {(1) = (E1[n, m]-E1[n+1, m])/s, (2) = (E2[n, m]-E2[n+1, m])/s})) = `%.`(Matrix(2, 2, {(1, 1) = 0, (1, 2) = 1, (2, 1) = 1, (2, 2) = 0}), Vector(2, {(1) = -(1/2)*(-r1[n, m-1]+r1[n, m+1])/h+(1/2)*(-r1[n, m-1]+E1[n, m-1]+r1[n, m+1]-E1[n, m+1])/h, (2) = -(1/2)*(-r2[n, m-1]+r2[n, m+1])/h+(1/2)*(-r2[n, m-1]+E2[n, m-1]+r2[n, m+1]-E2[n, m+1])/h}))

simplify(%);

(Vector(2, {(1) = (E1[n, m]-E1[n+1, m])/s, (2) = (E2[n, m]-E2[n+1, m])/s})) = `%.`(Matrix(2, 2, {(1, 1) = 0, (1, 2) = 1, (2, 1) = 1, (2, 2) = 0}), Vector(2, {(1) = (1/2)*(E1[n, m-1]-E1[n, m+1])/h, (2) = (1/2)*(E2[n, m-1]-E2[n, m+1])/h}))

Download ras5_v1_accc.mw

The technique might be adapted to handle several kinds of more involved Matrix/vector expressions with inert `%.`. That's related to a major reason why I made the code so involved -- for generality.

If results is a list, set, Vector/Array, or expression sequence (of, say, 7 elements) then you could do any of,

    return  results[2];

    return results[3..6];

    return results[2],results[4],results[5];

    return results[1],results[7];

Here are a few edits. I've added comments as to why I've made the changes.

Note that I've also allowed your utility formatters (Code-Edit region) to continue to function like the original.

(In the actual Maple GUI the units render more nicely. The double-braces appearing around units below are an artefact of this website's old backend.)

restart

kernelopts(version)

`Maple 2024.0, X86 64 LINUX, Mar 01 2024, Build ID 1794891`

conv2ftin:=proc(arg)

NULL

 

NULLNULL

totvol := proc (ht) options operator, arrow; Pi*.4^2*ht^2*ht+(1/3)*Pi*.4^2*ht^2*((2/3)*ht) end proc

proc (ht) options operator, arrow; .6143558967*ht^3 end proc

The procedures defined before this package is loaded
(ie. your utilities, as well as `totvol`) get the global
bindings for arithmetic.
That allows the utilities to work as intended, and
it allows the body of `totvol` to automatically simplify
at creation time, as well as to print as you wanted.

with(Units)

Automatically loading the Units[Simple] subpackage
 

Utilize the foot-pound-second system instead of SI.

UseSystem(FPS)

totvol(2.369133117)

8.169367285

fsolve(totvol(ht) = 8.169367284, ht)

2.369133117

volReq := 130*Unit('gallon'['US_liquid'])

130*Units:-Unit(gal)

volReq := 1.2*volReq

156.0*Units:-Unit(gal)

Alternate to next pair of commands below.
This produces result in ft, without having to first
convert volReq from gal to ft^3. For the FPS system
length dimensions computations simplify to ft.

volReq := convert(volReq, 'units', ft^3)

20.85416667*Units:-Unit(ft^3)

By default fsolve computes real roots.
It's suitable for your floating-point example.
Note that the result is in ft, not meters.

Also, by assigning directly to ht__main_tank we
avoid assigning to name `ht` which is used as
dummy name in a few places. This allows these
parts of your worksheet to be recomputed out-of-order,
since name ht is not clobbered.

ht__main_tank := fsolve(totvol(ht) = volReq, ht)

3.237856544*Units:-Unit(ft)

 

With the computed height in ft, we don't need
to convert all of these results separately.

ht__pan := (2/3)*ht__main_tank; rad := .4*ht__main_tank; total_height := ht__main_tank+ht__pan

2.158571029*Units:-Unit(ft)

1.295142618*Units:-Unit(ft)

5.396427573*Units:-Unit(ft)

Note that your utilities still function as intended.
(They don't continue to work properly if you just follow the
Answers suggesting that all you need to do it utilize `unapply`...)

conv2ftin(ht__main_tank); conv2ftin(ht__pan); conv2ftin(rad); conv2ftin(total_height)

3.*Units:-Unit(ft)+2.854278528*Units:-Unit(`in`)

2.*Units:-Unit(ft)+1.902852348*Units:-Unit(`in`)

1.*Units:-Unit(ft)+3.541711416*Units:-Unit(`in`)

5.*Units:-Unit(ft)+4.757130900*Units:-Unit(`in`)

convert(totvol(ht__main_tank), 'units', gallons)

156.0000001*Units:-Unit(gal)

Download Tank_Design_Calculation_-_Units_Question_(v00)_ac.mw

You can utilize the assuming facility to place separate assumptions during individual computations.

That allows you to compute under the assumption m::integer, but with no assumption on alpha.

For example, using Maple 2022 because the OP did,

restart;

kernelopts(version);

`Maple 2022.2, X86 64 LINUX, Oct 23 2022, Build ID 1657361`

expr2 := 4*cos(Pi*m+2*alpha*Pi)/(Pi*(2*m + 1));

4*cos(2*Pi*alpha+Pi*m)/(Pi*(2*m+1))

combine(expand(expr2)) assuming m::integer;

4*(-1)^m*cos(2*alpha*Pi)/(Pi*(2*m+1))

Download assuming_trig_ex.mw

It sounds as if you want nice pretty-printed ways to both enter and print those (2D Input and Output), and not just the means to compute.

For a more involved attempt you could look at this old Answer thread (and also some later revisions in worksheets attached to replies in the thread.)

For a preview, here's some simple functionality along those lines, again with 2D Input/Output.

restart;
`convert/degpolar`:=proc(x)
   degpolar(abs(x),180/Pi*argument(x));
end proc:
`print/degpolar`:=proc(ma,ar)
   uses Typesetting;
   mrow(Typeset(`&angle;`(ma,ar)),mi("&deg;"));
end proc:
`evalc/degpolar`:=proc(ma,ar)
   ma*(cos(ar*Pi/180)+I*sin(ar*Pi/180));
end proc:

`&angle;` := degpolar

degpolar

 

 

 

a := 1/2+sqrt(3)/2*I;

1/2+((1/2)*I)*3^(1/2)

ans := convert(a, degpolar);

degpolar(1, 60)

evalc(ans);

1/2+((1/2)*I)*3^(1/2)

 

For the next I've used the symbol from the Operators Palette.

 

q := `&angle;`(230, 0)/`&angle;`(92, -53.13)

degpolar(230, 0)/degpolar(92, -53.13)

convert(evalc(q), degpolar)

degpolar(2.500000001, 53.12999999)

Download degpolar_ch.mw

You specifically inquired about implicitplot3d, so in coding below I will address that in particular. It's documented [123] and not too difficult to color/shade a surface constructed by plot3d, etc. Hence if you could break up your implicit x-y-z region into explicitly defined portions then you could use another 3d plotting command. So I'm supposing that you wouldn't be asking about implicitplot3d if you had only an example that was subject to such recourse.

I'm just going to address coloring (arbitrarily) with implicitplot3d, since dharr's already treated your particular surface examples with a split into explicitly formulated portions.

You can in fact utilize a float[8] Array to color the triangles by which the GUI renders the surface, when it displays an ISOSURFACE plotting structure generated by the implicitplot3d command. You can use such an Array to specify coloring by either RGB or HSV values.

But the big problem is that those triangles don't match the set of x-y-z points stored as data in the ISOSURFACE. The ISOSURFACE values are just for the regularly spaced points in the ranges. The location of the vertices of those triangles are a result of the interpolating process by which the GUI's plot driver figures out where the surface might be, and we are not privy to that. So we lack programmatic access to the x-y-z values for the vertices of the GUI-rendered triangles.

Even the number of the mesh points of stored function values doesn't match the number of triangles in the rendered surface. And that number of triangles is also a function of the proportion of the plotting window which the surface occupies.

Here's an example. This code does not form those triangles. The GUI does that. This code merely adds a COLOR substructure to a usual implicitplot3d result, to specifiy hue values for its rendered triangles.

subsindets(plots:-implicitplot3d(x^2+y^2+z^2=3.9,x=-2..2,y=-2..2,z=-2..2,grid=[6,6,6],
                                 lightmodel=none,orientation=[60,75,0]),
           specfunc(ISOSURFACE),
           Iso->ISOSURFACE(op(Iso),
                           COLOUR(HSV,Array(1..1528,1..3,
                                            (ijk,hsv)->
                                              `if`(hsv=1,
                                                   piecewise[1.0](irem(ijk,3)=0,0.1,
                                                                  irem(ijk,3)=1,0.4,
                                                                  irem(ijk,3)=2,0.7),
                                                   1.0),
                                            datatype=float[8]))));

 

 

 

Download ISO_col_0.mw

note: If you changed that to utilize grid=[6,6,7] instead then in Maple 2016.2 and later you could see an example of the GUI message spit out when the coloring Array doesn't have the right number of entries.

As it happens that COLOR substructure's Array can match either the number of segments of the number of vertices (whichever is closest, I suspect). And the GUI spits out a message when there's a bad mismatch. But in neither case do we know the x-y-z coordinates of those vertices, and so we don't know how to populate that COLOR structure's Array. (If I had such a means I'd ditch ISOSURFACE myself, and do the math and draw those portions by my own construction, with full control...)

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