acer

32333 Reputation

29 Badges

19 years, 325 days
Ontario, Canada

Social Networks and Content at Maplesoft.com

MaplePrimes Activity


These are answers submitted by acer

I'd be interested to know whether this procedure `UpdatingRead` accomplishes what you want. It seems to work, in my 64bit Maple 18.02 and 2015.1 on Windows 7.

proc()
  assign(':-UpdatingRead',
         proc(fn::string)
           local oldanames;
           oldanames := {anames(':-user')};
           read fn;
           if IsWorksheetInterface(':-Standard') then
             try
               map(:-_UIUtils:-VariableManager:-AddToVariableManager,
                   {anames(':-user')} minus oldanames);
             catch:
               WARNING("variable manager update interrupted");
             end try;
           end if;
           NULL;
         end proc);
  NULL;
end proc();

The idea is to call it with the (string) name of the .m file.

UpdatingRead( cat(kernelopts(homedir),"/Downloads/foo.m") );

acer

Here are two ways, for lists of numeric entries (like your example L).

The first of these seems faster than the second, for long L.

ListTools:-SelectFirst(nops(L),`>`,L,0,output=indices);

[ListTools:-SearchAll(1,map(signum,L))];

But the first of those may have stack limit problems if tried repeatedly for long L (10^6 entries, say).

If you also want L to be able to contain nonnumeric entries then adjust the predicate in the SelectFirst approach above.

acer

Do you need something fancier than this?

G := (L::list(name),n::posint)->op((`*`@op)~(combinat[choose](L,n))):

G( [a,b,c], 2 );
                                a b, a c, b c

G( [q,r,s,t], 3 );
                         q r s, q r t, q s t, r s t

G( [q,r,s,t], 2 );
                        q r, q s, q t, r s, r t, s t

It wasn't clear to me that your input must be names. If not, then change the operator's first parameter from L::list(name) to just L::list, say.

acer

This particular problem can be solved using eliminate and without having to use solve or fsolve. In my opinion this is clean and simple. Others are free to hold their own opinions.

eq1:=convert( X=cos(theta) + 0.8e-1*cos(3.*theta), rational):
eq2:=convert( Y=-sin(theta)+ 0.8e-1*sin(3.*theta), rational):

sols:=[eliminate({eq1,eq2},{Y,theta})]:

S:=seq(eval(Y,sols[i][1]),i=1..nops(sols)):

Now S[1] is a formula for Y in terms of X alone.

There is another solution. I give just one way of obtaining it -- there are others. Observe that the formula for X is even (in theta), and formula for Y is odd. Hence we note that for a given X we can look to Y and -Y. No need to take my word for it: substitute -theta for theta, and re-eliminate.

osols:=[eliminate(subs(theta=-theta,{eq1,eq2}),{Y,theta})]:

oS:=seq(eval(Y,osols[i][1]),i=1..nops(osols)):

The other solution is oS[1] which is also a formula for Y in terms of X alone.

plot([S[1],oS[1]],X=-1-0.8e-1..1+0.8e-1,color=red);

Another thing that the eliminate command returns is further restrictions (second operand) on its main result (first operand). In this case there are no qualifying restrictions involving the eliminated variable(s).

seq(sols[i][2], i=1..nops(sols));  # restrictions on X (NULL)

                           {}, {}, {}

seq(osols[i][2], i=1..nops(osols));  # restrictions on X (NULL)

                           {}, {}, {}

No exact, symbolic method will work for all such problems in general. Approximate floating-point back-solving (eg. fsolve) is sometimes slower (like in this case), and it is prone to solution-tracking issues, especially for multivalued solutions. So if you can obtain a reasonably short explicit formula in a short amount of time then it's useful. Of course it is also easy to find example problems for which exact symbolic solving takes an enormous amount of resources to obtain an impractically un-useful result.

I'll probably ignore irrelevant comments.

acer

[edit: I initially misread your question as being about how to plot the intersection of the surfaces, and I initially focused on your Question's title. Sorry. Later down I get to the region, and to restricting the surface faces.]

There is a command for this in the plots package, called intersectplot. It produces the red spacecurve below.

If you don't want the gray plane as well, then just omit it from the call to display.

I deliberately used x^2=4-y^2 , for a circular cylinder below (though the axes scaled are displayed unconstrained). Change that as you wish.

restart:

with(plots):

display(intersectplot( x^2=4-y^2, z=2-y, x=0..4, y=0..2, z=0..2, thickness=4),
        plot3d(2-y, x=0..4, y=0..2, style=surface, color=gray, transparency=0.4),
        implicitplot3d(x^2=4-y^2, x=0..4, y=0..2, z=0..2, style=surface,
                       color=gold, transparency=0.4, grid=[30,30,10]),
        lightmodel=none, axes=box);

 

There are several other ways to accomplish this example, of course.

And, using your x=4-y^2 one could obtain the following. (I constrain the scaling of the axes here, for fun).

display(intersectplot( x=4-y^2, z=2-y, x=0..4, y=0..2, z=0..2, thickness=4),
        plot3d(2-y, x=0..4, y=0..2, style=surface, color=gray, transparency=0.4),
        implicitplot3d(x=4-y^2, x=0..4, y=0..2, z=0..2, style=surface,
                       color=gold, transparency=0.4, grid=[30,30,10]),
        lightmodel=none, axes=box, scaling=constrained);

As I mentioned, apoligies for focusing on the intersecting curve, while you really mentioned wanting to plot the enclosed region.

Note that the range arguments of the plot3d command can be variable. That's how I restrict the plane z=2-y in one of the plot3d calls below.

Also, a coordinate-transformed surface from plot3d can sometimes have less jagged rendered edges than a surface from implicitplot3d. Changing coordinate systems can be another way to deal with that. Below I take the portion of the parabolic sheet (what you called a cylinder) as being x=f(z,y) a mathematical function. And then I transformed to flip x and z coordinates.

restart:
with(plots):

display(intersectplot( x=4-y^2, z=2-y, x=0..4, y=0..2, z=0..2, thickness=4),
        plot3d(2-y, x=0..4-y^2, y=0..2, style=surface,
               color=gold, transparency=0.4),
        plottools:-transform((z,y,x)->[x,y,z])(plot3d(4-y^2, z=0..2-y, y=0..2,
                                               style=surface, color=gold,
                                               transparency=0.4)),
        lightmodel=Light4, glossiness=1.0, axes=box, scaling=constrained,
        view=[0..4, 0..2, 0..2]);

I'll leave the bottom and back plane portions to you. The variable range is the key.

acer


restart:

kernelopts(version);

`Maple 18.02, X86 64 LINUX, Oct 20 2014, Build ID 991181`

# -0.15356213  says wolframalpha.com for ee
# -0.15356212695352467  says wolframalpha for F5

MK := 1/96*(3^(3/5) - 9)*sqrt(10-2*sqrt(5))*GAMMA(7/5):

evalf[500](MK): evalf[15](%);

-.153562126953523

ee := Int(sin((x - 1/x)^5)^3/x^3, x=0..infinity);

Int(sin((x-1/x)^5)^3/x^3, x = 0 .. infinity)

with(IntegrationTools):

Split(ee, 1);

Int(sin((x-1/x)^5)^3/x^3, x = 0 .. 1)+Int(sin((x-1/x)^5)^3/x^3, x = 1 .. infinity)

# ee = f1 + f2

f1,f2 := op(%);

Int(sin((x-1/x)^5)^3/x^3, x = 0 .. 1), Int(sin((x-1/x)^5)^3/x^3, x = 1 .. infinity)

# ee = F1 + f2

F1 := student[changevar](t=x, student[changevar](t=1/x, f1, t), x);

Int(sin((1/x-x)^5)^3*x, x = 1 .. infinity)

# ee = F1a + f2

F1a := Int(subs(x=-x,op(1,F1)), op(2,F1));

Int(-sin((x-1/x)^5)^3*x, x = 1 .. infinity)

# Last was ok (even)
simplify(F1 - F1a);

0

# ee = EE = F1a + f2

EE := frontend(combine, [F1a+f2], [{`+`,specfunc(Int)},{}]);

Int(-sin((x-1/x)^5)^3*x+sin((x-1/x)^5)^3/x^3, x = 1 .. infinity)

# ee equals EEa

EEa := simplify(student[changevar](u=x-1/x, EE, u));

Int(-sin(u^5)^3*u, u = 0 .. infinity)

# ee equals EEb

EEb := student[changevar](w=u^5, EEa, w);

Int(-(1/5)*sin(w)^3/w^(3/5), w = 0 .. infinity)

func := Int(op(1,EEb), w=2*k*Pi..(2*k+2)*Pi);

Int(-(1/5)*sin(w)^3/w^(3/5), w = 2*k*Pi .. (2*k+2)*Pi)

# ee equals add(func, k=0..infinity))

add(func, k=0..2); # you get the idea...

Int(-(1/5)*sin(w)^3/w^(3/5), w = 0 .. 2*Pi)+Int(-(1/5)*sin(w)^3/w^(3/5), w = 2*Pi .. 4*Pi)+Int(-(1/5)*sin(w)^3/w^(3/5), w = 4*Pi .. 6*Pi)

# evalf(ee) equals evalf(Sum(func, k=0..infinity)))

evalf[15](Sum(func, k=0..infinity));

-.153562126953523

 


Download mk1.mw

acer

SFloat@op seems to work even if UseHardwareFloats=true.

> restart:

> t := HFloat(1.00000000000002):

> lprint(t);
HFloat(1.00000000000001998)

> (SFloat@op)(t);

                             1.00000000000001998

> lprint(%);
1.00000000000001998

> evalf((SFloat@op)(t));

                                 1.000000000
> restart:

> UseHardwareFloats:=true:

> t := HFloat(1.00000000000002):

> lprint(t);
HFloat(1.00000000000001998)

> (SFloat@op)(t);

                             1.00000000000001998

> lprint(%);
1.00000000000001998

> evalf((SFloat@op)(t));

                                 1.000000000

acer

AFAIK the text in a TextArea component doesn't have any Character or Paragraph style associated with it.

Combined with it not being a property of the component, this all indicates to me that it probably can't be done.

I recall resorting to marked up text at some time in the past, to get colored, bold text inside a MathContainer component.

acer

Is the list L below the kind of thing that you want?

restart:

f := (x) -> (310*(x+0.5)^0.2)+70:

z := [ -0.5, 0.5 ]:

P1 := plot( f, z[1]..z[2], view=-100..400 ):

P1;

L := [ seq( z[2]-(z[2]-z[1])/(50/(51-i))^(3), i=1..50 ) ];

# The entries of L are getting successively closer together.
#
seq( L[i+1]-L[i], i=1..49 );

P2 := plots:-pointplot( [ seq( [t,f(t)], t=L ) ], view=0..400, symbolsize=15 ):

plots:-display( P1, P2 );

# Or, do you need to ensure that
#    abs(f(L[i+1])-f(L[i])) < abs(f(L[i])-f(L[i-1]))
# ?
# That happens to be the case for your f. The successive
# abs(f(L[i+1])-f(L[i])) are decreasing.
#
plots:-listplot( [ seq( abs(f(L[i+1])-f(L[i])), i=1..49 ) ] );

Or do you want the reverse, with the L[i] or f(L[i]) becoming closer toward the left?

Or do you just want the L[i] taken uniformly over -0.5 .. 0.5, including end-points, where for your particular f it happens that abs(f(L[i+1])-f(L[i])) is decreasing?

restart:

f := (x) -> (310*(x+0.5)^0.2)+70:

z := [ -0.5, 0.5 ]:

L := [ seq( z[1] + (i-1)*(z[2]-z[1])/(50-1), i=1..50 ) ];

[ seq( abs(f(L[i+1])-f(L[i])), i=1..49 ) ];

plots:-listplot( [ seq( abs(f(L[i+1])-f(L[i])), i=1..49 ) ] );

plots:-display(
  plots:-pointplot( [ seq( [t,f(t)], t=L ) ], view=0..400, symbolsize=15 ),
  plot( f, z[1]..z[2], view=-100..400 ) );

Sorry that I'm finding your explanation difficult to read.

acer

I'm not sure whether any of this might be useful to you.

restart:

`diff/E`:=proc(t) K(t); end proc:

diff(E(y),y);

                                    K(y)

`D/E`:=t->K:

D(E);

                                      K

D[1](E)(y);

                                    K(y)

(D@@2)(E);

                                    D(K)

(D@@2)(E)(y);

                                   D(K)(y)

acer

This looks like a regression bug to me. I've submitted a bug report (SCR).

I'm guessing that you already know that you can use the strings ["1,"10","100"] and parse, to get it to work in both versions. But that's just a workaround; it should still work with numeric values.

acer

As Carl mentions, you are being bitten by premature evaluation. Note that sum and Sum adhere to Maple's usual evaluation rules for procedure calls: first the supplied arguments are evaluated, and only after that are they passed to the called procedure.

You can search this site, or the web, and find lots of similar examples.

The problem is one of usage -- that the new user has not yet learned how the system normally works.

And it comes up for lots of other commands too (again, especially for new users, until they learn how Maple works). You can find examples for plot, fsolve, Optimization:-Minimize, and so on.

Carl mentions two ways to make things work, but in fact there are sometimes more than just those two. Some commands like plot, fsolve, and Minimize also accept a so-called operator form for the first argument. For example plot(f,0..1) instead of plot(f(x),x=0..1), where there is no dummy variable and hence no function call to be accidentally prematurely evaluated.

And for some other commands, such as sum and product there are related commands add and mul with special evaluation rules. The basic idea is that add follows a special rule: it receives its first argument without that having being evaluated. For this and a few other reasons using add is often suggested for finite summation, over using sum.

If you want to see what the call sum(pdf4(i),i=3..3) is actually sending to the sum procedure then call pdf4(i) alone, outside of a call to sum or anything else, and see what result is returned.

I won't repeat the two ways that Carl mentioned in his good Answer. Both those approaches often suffice with both the commands like plot as mentioned above, and also with the commands like sum as mentioned above. However I often prefer using the approaches I've mentioned above, if possible, rather than having to rewrite my custom procs to return unevaluated on nonnumeric input.

What you've shown is a (mostly new user related) usage problem. But there are indeed flaws. 2D Math Input allows for the typeset summation (Sigma) symbol only for sum and Sum. There is no good hook between the typeset summation symbol and add, unless you make your code convoluted by utilizing special exports of the Physics package. Yet add would often serve better than sum. And 2D Math Input is the default for new users. And the documentation pushes palette use big time. And combined that is all quite unfortunate because it make this issue arise far more than ought to be necessary.

Note that new users to Mathematica often have equivalent or similar problems, until they learn to use NumericQ.

acer

Matrix(5,5,(i,j)->i+j)

Also if you intend on indexing into this Matrix (let's suppose you assigned it to the name A) later on then use the syntax such as A[i,j] and not A[i][j] which can be quite inefficient.

In this particular case that you've asked about you can also do simply Matrix(5,5,`+`) , but the way above shows how you can construct more generally.

acer

I use Maple 18.02 below, as this Question is marked as related to Maple 18.

Here are a few ways that work with your approach of using double pairs of unevaluation quotes around each of the plot calls.

restart;
with(plots):
with(plottools):
Parab1 := ''plot(a*x^2, x = -1 .. 1, y = -3 .. 3)'':
Parab2 := ''plot(a*x^2+1, x = -1 .. 1, y = -3 .. 3)'':
L:=Array(1..2):
L[1]:=Parab1:
L[2]:=Parab2:

Explore((unapply('display'(L),a))(A), parameters = [[A = -1.0 .. 1.0, label=a]]);

And here is another way that also works with that approach,

restart;
Parab1 := ''plot(a*x^2, x = -1 .. 1, y = -3 .. 3)'':
Parab2 := ''plot(a*x^2+1, x = -1 .. 1, y = -3 .. 3)'':
L:=Array(1..2):
L[1]:=Parab1:
L[2]:=Parab2:

F:=proc(A)
     plots:-display(eval(L,a=A));
end proc:
Explore(F(a), parameters=[a=-1.0 .. 1.0]);

And here is another way that also works with that approach,

restart;
with(plots):
with(plottools):
Parab1 := ''plot(a*x^2, x = -1 .. 1, y = -3 .. 3)'':
Parab2 := ''plot(a*x^2+1, x = -1 .. 1, y = -3 .. 3)'':
L:=Array(1..2):
L[1]:=Parab1:
L[2]:=Parab2:

eval('Explore'('display'(L), parameters = [[a = -1.0 .. 1.0]]));

Here is a way which doesn't use unevaluation quotes.

restart:
H:=proc(a)
     plots:-display(Array([plot(a*x^2, x = -1 .. 1, y = -3 .. 3),
                           plot(a*x^2+1, x = -1 .. 1, y = -3 .. 3)]));
end proc:  
Explore(H(a), parameters=[a=-1.0 .. 1.0]);

acer

[This answer has been edited: see comments below for modifications, obtaining 16 solutions with all of t,t1,t2,t3 between 0 and 2*Pi]

Here I'm getting 128 solutions with all of t,t1,t2,t3 between -2*Pi and 2*Pi.

(I used Maple 18.02 on 64bit Linux.)

restart:
f1 := -8100+(-30+70*cos(t1)-40*cos(t2))^2
      +(-70*sin(t1)+40*sin(t2))^2:
f2 := (-20-80*cos(t3))^2+(-15+70*cos(t1)+10*cos(t1+t))^2
      +(-70*sin(t1)-10*sin(t1+t)+80*sin(t3))^2-5625:
f3 := (-20-80*cos(t3))^2+(15+40*cos(t2)+10*cos(t1+t))^2
      +(-40*sin(t2)-10*sin(t1+t)+80*sin(t3))^2-5625:
f4 := 10*cos(t1+t)*(30-70*cos(t1)+40*cos(t2))
      -10*sin(t1+t)*(70*sin(t1)-40*sin(t2)):
origsys := {f1,f2,f3,f4}:
sys := map(expand,origsys):
T := [indets(sys,function)[]]:
Tsub := [seq(x=freeze(x),x=T)]:
extra := {seq(sin(r)^2+cos(r)^2-1,r=[t,t1,t2,t3])}:
Digits,oldDigits := 20*Digits,Digits:
raw := thaw(RootFinding:-Isolate(eval(sys union extra,Tsub),eval(T,Tsub))):
Digits := oldDigits:
gen := proc(cand)
  local c, v;
  c[t]:=select(has,cand,t);
  c[t1]:=select(has,cand,t1);
  c[t2]:=select(has,cand,t2);
  c[t3]:=select(has,cand,t3);
  [seq(seq(seq(seq([e,f,g,h],e=c[t3]),f=c[t2]),g=c[t1]),h=c[t])];
end proc:
sols := {seq(seq(
      `if`(evalf[2*Digits](max(map(abs,eval(origsys,eq))))<10^(-Digits+1),
           eq, NULL),
      eq=evalf[2*Digits](gen(
map(eq->select(u->is(abs(rhs(u))<=2*Pi),
    [seq(lhs(eq)=rhs(eq)+k*Pi,k=-4..4)])[],
    (map(eq->fsolve(eq,indets(eq,name))[],raw[i])))
))),
    i=1..nops(raw))}:

sols := evalf[15](sols):
nops(sols);
                                     128

max(seq(max(map(abs,evalf[2*Digits](eval(origsys,s)))), s in sols));
                                          -11
                                6.02316 10   

#seq(print(evalf[10](s)), s=sort(sols));
First 223 224 225 226 227 228 229 Last Page 225 of 336