acer

32385 Reputation

29 Badges

19 years, 335 days
Ontario, Canada

Social Networks and Content at Maplesoft.com

MaplePrimes Activity


These are answers submitted by acer

[Following your later edit to the Question:] It's easy enough to get from poly1 to poly3.

restart;

poly1:=(A*B/(A+B)+X)/(X+Y*X/(Y+X));

                            A B     
                           ----- + X
                           A + B    
                           ---------
                                Y X 
                           X + -----
                               Y + X

applyrule( (a::name*b::name)/(a::name+b::name)=f(a,b), poly1 );

                          f(A, B) + X
                          -----------
                          X + f(X, Y)

So that could leave getting from poly2 to poly1.

Hopefully I've interpreted your question properly.

If you can only obtain an explicit representation in terms of one of the two variables (and it's not the one you want...) then you could use the implicitplot command, or you could reflect the result from the simpler plot command about the line x=y.

The first of those ways is more straightforward to remember, but in (numerically) difficult cases the second approach can often produce a cleaner curve while taking less computational resources.

Here's an example. Note that special options like gridrefine are often needed to get implicitplot to produce as smooth a curve, in my experience.

impeq := x = sin(y) + y:

with(plots): with(plottools):

implicitplot(impeq,x=-10..10,y=-10..10,gridrefine=2);

display(transform((a,b)->[b,a])(plot(eval(x,impeq),y=-10..10)),labels=[x,y]);

Download refl.mw

acer

Please spend some more time thinking about all the responses you've received to your earlier, related Questions about these kinds of expressions involving very large exponents and round-off error.

forget(evalf);
Digits:=20:
evalf(y(99.6));
                                 270
                            -1 10   

forget(evalf);
Digits:=300:
evalf(y(99.6));
                          0.0184074142

acer

The vertical bar and the name sn are the result of using the extended typesetting mode. The value is the same in either case. It's the particular mode of prettyprinting (typesetting) of the output that differs.

restart;

interface(typesetting=standard):

JacobiSN((1/2)*sqrt(2)*x+EllipticK(I), I)

JacobiSN((1/2)*2^(1/2)*x+EllipticK(I), I)

interface(typesetting=extended):

JacobiSN((1/2)*sqrt(2)*x+EllipticK(I), I)

JacobiSN((1/2)*2^(1/2)*x+EllipticK(I), I)

kernelopts(version);

`Maple 2015.2, X86 64 WINDOWS, Nov 13 2015, Build ID 1087698`

 

Download sn.mw

It is possible to disable the typesetting of output on a function by function basis, using the RuleAssistant. Using the dropmenu in that popup assistant one can retain extended typesetting mode in general while changing the typeset look of function calls of particular names such as JacobiSN. I cannot recall right now whether I ever figured out how to do this entirely programmatically without using the popup. But without the means to do it programmatically then this functionality of the assistant seems not very useful to me, since I don't see how to enforce the per-function disabling of output typesetting upon reopening a worksheet.

acer

The result you see for n=5 is correct.

Your simplification of c1 is not strong enough to show that c1=120 when n=5. The result is larger and more complicated by default than that returned by linalg[cond], however. (See my followup comment below.)

I show that the exact LinearAlgebra:-ConditonNumber result is correct for n=5 below, by using an alternative ("stronger", for this example) simplification.

restart;

with( LinearAlgebra ):

Digits := 50:

n := 5;

M := [ seq( cos( Pi*j/n ), j=0..n ) ]:

V := VandermondeMatrix( M ):

5

c1 := ConditionNumber(V):

c1 := simplify(convert(simplify(c1),radical));

evalf(c1);

120

120.

c2 := linalg[cond](convert(V,matrix)):

c2 := simplify(%);

evalf(c2);

120

120.

kernelopts(version);

`Maple 2015.2, X86 64 LINUX, Nov 13 2015, Build ID 1087698`

 

 

Download trigsimp.mw

acer

There is a collspased (hidden) Execution Group in the Document Block in the Table Cell in the 4th row and 3rd column.

That's the same row whose first column Cell contains "inclination angle of the swash plate".

The problematic Cell curretly shows a TextArea component whose current entry is "1.2".

To reveal it, use the mouse pointer and select the content of that Cell (or an larger selection that contains that Cell). Then, while it is selected, use the main menubar item View->Group and Block Management->Expand Document Block. After that the offending Execution Group will be visible, below the Execution Group that contains the TextArea.

Delete the problematic Execution Group, which contains a call to the Tutor. You can use Ctl-Delete for that. The Select the upper Execution Group's content (the TextArea), and use the menubar item View->Group and Block Management->Collapse Document Block.

That corrects the problem with the Tutor launching every time the sheet is executed in its entirety.

There is another mistake. The Action code on that TextArea component in the Cell in the 4th row and 3rd column contains some code. Right-click on it invoke the "Edit Content Changed Action..." popup menu item, so as to edit that code. The mistake is that it is missing the final end use; line.

acer

Notice how the mL and g in your axis label are in italic. But if you used :-min in the label then it would be typeset in upright (roman).

It would look better if all these unit symbols and abbreviations (prefixed or not) were typeset consistently. And one simple way to get that effect is to create the axis labels using the Unit() command.

For example, using your previous code,

plot(fprod(t), t = -2..2,
     labels = [typeset("Time (", Unit(min), ")"),
               typeset("Rate (", Unit(mL/(min*g)), ")")],
     color = "blue");

acer

Is this what you mean?

restart:

ode:=diff(x(t), t, t)+(1000*(x(t)^2-1))*(diff(x(t), t))+x(t) = 0:
sollistA:=dsolve({ode,x(0)=-1,D(x)(0)=-0.8},numeric,output=listprocedure):

# method 1, numerical differentiation, after dsolve. Not generally robust.
X:=eval(x(t),sollistA):
plot(t->fdiff(X(tt),[tt,tt],{tt=t}), 0..1);

# method 2, isolate higher derivative from the original ode, if you can,
# and evaluate each function on the right using the listprocedure.

rhs(isolate( ode, diff(x(t), t, t) ));
DDXA:=eval(%, sollistA):
plot(DDXA, 0..1);

# method 3, augment the original DE with an extra equation,
# and then use that function from the listprocedure.

sollistB:=dsolve({ode,diff(x(t),t,t)=ddx(t),x(0)=-1,D(x)(0)=-0.8},
                 numeric,output=listprocedure):

DDXB:=eval(ddx(t),sollistB):
plot(DDXB, 0..1);

acer

The first of these just happens to work (but you can also force it with a few `collect` calls and a mapped simplify...).

It might be possible to harangue the terms so that the "Nu" name appears at the left of each product. I didn't go there.

restartNULL

ee := 1-(1/4)*k1*`Νu`*(1+k2)+(1/64*(`Νu`*k1*k2^2+2*`Νu`*k1*k2+`Νu`*k1+4))*k1*`Νu`-(1/2304)*`Νu`*k1*(`Νu`^2*k1^2*k2^3+3*`Νu`^2*k1^2*k2^2+3*`Νu`^2*k1^2*k2+`Νu`^2*k1^2+20*`Νu`*k1*k2+20*`Νu`*k1-64*k2)+(1/147456*(`Νu`^2*k1^2*k2^4+4*`Νu`^2*k1^2*k2^3+6*`Νu`^2*k1^2*k2^2+4*`Νu`^2*k1^2*k2+`Νu`^2*k1^2+56*`Νu`*k1*k2^2+112*`Νu`*k1*k2+56*`Νu`*k1-640*k2^2-640*k2+144))*`Νu`^2*k1^2-(1/14745600)*k1^2*`Νu`^2*(`Νu`^3*k1^3*k2^5+5*`Νu`^3*k1^3*k2^4+10*`Νu`^3*k1^3*k2^3+10*`Νu`^3*k1^3*k2^2+5*`Νu`^3*k1^3*k2+120*`Νu`^2*k1^2*k2^3+`Νu`^3*k1^3+360*`Νu`^2*k1^2*k2^2+360*`Νu`^2*k1^2*k2-2944*`Νu`*k1*k2^3+120*`Νu`^2*k1^2-5888*`Νu`*k1*k2^2-1520*`Νu`*k1*k2+1424*`Νu`*k1-13312*k2)+(1/2123366400*(`Νu`^4*k1^4*k2^6+6*`Νu`^4*k1^4*k2^5+15*`Νu`^4*k1^4*k2^4+20*`Νu`^4*k1^4*k2^3+15*`Νu`^4*k1^4*k2^2+220*`Νu`^3*k1^3*k2^4+6*`Νu`^4*k1^4*k2+880*`Νu`^3*k1^3*k2^3+`Νu`^4*k1^4+1320*`Νu`^3*k1^3*k2^2-9344*`Νu`^2*k1^2*k2^4+880*`Νu`^3*k1^3*k2-28032*`Νu`^2*k1^2*k2^3+220*`Νu`^3*k1^3-21008*`Νu`^2*k1^2*k2^2+4704*`Νu`^2*k1^2*k2+7024*`Νu`^2*k1^2-205312*`Νu`*k1*k2^2-205312*`Νu`*k1*k2+14400*`Νu`*k1+409600*k2^2))*`Νu`^2*k1^2-(1/416179814400)*k1^3*`Νu`^3*(`Νu`^4*k1^4*k2^7+7*`Νu`^4*k1^4*k2^6+21*`Νu`^4*k1^4*k2^5+35*`Νu`^4*k1^4*k2^4+35*`Νu`^4*k1^4*k2^3+364*`Νu`^3*k1^3*k2^5+21*`Νu`^4*k1^4*k2^2+1820*`Νu`^3*k1^3*k2^4+7*`Νu`^4*k1^4*k2+3640*`Νu`^3*k1^3*k2^3-23744*`Νu`^2*k1^2*k2^5+`Νu`^4*k1^4+3640*`Νu`^3*k1^3*k2^2-94976*`Νu`^2*k1^2*k2^4+1820*`Νu`^3*k1^3*k2-118160*`Νu`^2*k1^2*k2^3+364*`Νu`^3*k1^3-22064*`Νu`^2*k1^2*k2^2+49168*`Νu`^2*k1^2*k2-1435648*`Νu`*k1*k2^3+24304*`Νu`^2*k1^2-2871296*`Νu`*k1*k2^2-1216192*`Νu`*k1*k2+9625600*k2^3+219456*`Νu`*k1+9625600*k2^2-3990528*k2):
NULL

ans1 := sort(simplify(ee, size))

-(1/416179814400)*(k2+1)^7*`Νu`^7*k1^7+(1/2123366400)*(k2-6/7)*(k2+1)^5*`Νu`^6*k1^6-(1/92897280)*(k2^2-(93/40)*k2+21/10)*(k2+1)^3*`Νu`^5*k1^5+(79/33177600)*(k2^3+(11145/3871)*k2^2+(41787/15484)*k2+4631/7742)*(k2+1)*`Νu`^4*k1^4-(6541/25401600)*(k2^3+(65315/26164)*k2^2+(107003/52328)*k2+234171/418624)*`Νu`^3*k1^3+(119/10368)*(k2^2+(4959/2975)*k2+657/952)*`Νu`^2*k1^2-(2/9)*(k2+27/32)*`Νu`*k1+1

var := (`minus`(indets(ee, name), {k1, k2}))[1]:

-(1/416179814400)*(k2+1)^7*`Νu`^7*k1^7+(1/14863564800)*(7*k2-6)*(k2+1)^5*`Νu`^6*k1^6-(1/3715891200)*(40*k2^2-93*k2+84)*(k2+1)^3*`Νu`^5*k1^5+(1/6502809600)*(15484*k2^4+60064*k2^3+86367*k2^2+51049*k2+9262)*`Νu`^4*k1^4-(1/1625702400)*(418624*k2^3+1045040*k2^2+856024*k2+234171)*`Νu`^3*k1^3+(1/2073600)*(23800*k2^2+39672*k2+16425)*`Νu`^2*k1^2+(-(2/9)*k1*k2-(3/16)*k1)*`Νu`+1

simplify(ans1-ans2);

0

length(ans1), length(ans2);

499, 451

``


Download Polysimp.mw

acer

Without attempting to manipulate your procedure's body in order to deal with the numeric difficulties more gracefully...

restart:

y:=x->-9.8455282400*10^9142*exp(-(2/3)*x^(3/2))/(x^(1/4)*sqrt(Pi))
      +3.3889331940*10^(-9169)*exp((2/3)*x^(3/2))/(x^(1/4)*sqrt(Pi))
      +(16/153)*x^(7/6)*sqrt(Pi)*exp((2/3)*x^(3/2))
      +Pi*((1/2)*exp(-(2/3)*x^(3/2))*(-1+exp((2/3)*x^(2/3)))/(x^(1/4)*Pi)
      -(16/153)*x^(7/6)*exp((2/3)*x^(3/2))/sqrt(Pi)):

forget(evalf): Digits:=20:
P1:=plottools:-transform((x,y)->[x+1000,y])(plot(expand(y(x+1000)),x=0..1)):

forget(evalf): Digits:=9200:
P2:=plot(y, 1000..1001, style=point, adaptive=false, numpoints=25):
Digits:=20:

plots:-display(P1,P2);

Download trighplot.mw

acer

Just assign the result of plots:-display to P.

  P := plots:-display( P1, P2 ):

Then P is a structure that contains both. So then either of the following will once more show it.

  plots:-display( P );
  print( P );

Also, in an earlier question you mentioned animating, if I recall. If so then try,

  plots:-display( P1, P2,... Pn, insequence=true );

which again you can assign to a new name.

acer

See section 6.5 of the Programming Manual, under the subheading "The remember, cache, and system Options".

acer

You didn't provide values for a, D1, and D2. So I just made some up.

If you wanted you could also treat those as parameters (in the dsolve/numeric sense).

Below, I'll treat q as a parameter for dsolve. See here for details on that.

restart:

with(plots):

odes := D2/(3*a*D1+4*D2)*1/P(x)*diff(P(x),x)
        -q/(q*S(x)-1)*diff(S(x),x)
        -3*a*D1/(3*a*D1+4*D2)*cot(x)=0,
        diff(S(x),x$2)+diff(S(x),x)*cot(x)
        +4*Pi*(3*a*D1+4*D2)/((q*S(x)-1)*D2)=0:

ics := P(Pi/2)=1, S(Pi/2)=-1, D(S)(Pi/2)=0:

vals := [a=-0.2222, D1=0.027, D2=.00416]:

sol := dsolve({eval(odes,vals),ics}, numeric, parameters=[q]):

F:=proc(qi)
     sol( parameters = [ q = qvals[qi] ] );
     odeplot( sol, eval( [x, (3*D1*a+4*D2)*P(x)/((1-q*S(x))*D2)], [vals[],q=qvals[qi]] ),
              0.87 .. Pi/2,
              color=qcolors[qi], linestyle=qlinestyles[qi],
              tickmarks = [ [ seq((1/10)*i*Pi = (180*i*(1/10))*`°`, i = 1 .. 8) ],
                            default ]
              );
   end proc:

qvals := [0.0, 0.6, 0.7, 0.8, 0.9]:
qcolors := [gold,red,blue,green,magenta,cyan]:
qlinestyles := [dot,solid,dash,dashdot,longdash,spacedash]:

plots:-display( seq(F(i),i=1..nops(qvals)),
                legendstyle=[location=left], size=[600,400] );



Download odecurves.mw

acer

The curve is almost flat in the given range of x=-4..4. So to get an idea of the qualitative aspects of the curve one can increase working precision (to deal with accuracy issues), subtract off the ostensibly constant part, and manually create tickmarks.

 

restart:

ee := 1+18*(sinh(9*x-9)-sinh(3*x-477))^2/(9*cosh(9*x-9)+cosh(3*x-477))^2;

1+18*(sinh(9*x-9)-sinh(3*x-477))^2/(9*cosh(9*x-9)+cosh(3*x-477))^2

evalf[1000](convert(series(ee,x,6),polynom)): evalf[5](%);

19.000+0.12152e-199*x-0.36455e-199*x^2+0.72910e-199*x^3-0.10937e-198*x^4+0.13124e-198*x^5

[evalf[1000](eval(ee,x=-4)-19), evalf[1000](eval(ee,x=4)-19)]:
evalf[5](%);

[-0.53648e-190, -0.17314e-187]

yoffscal := 1e-192:
yticks := [ seq( yadj*yoffscal=sprintf("19%+.1e",yadj*yoffscal), yadj=-10..10 ) ]:
Digits := 300:

plot( ee-19, x=-4..4, ytickmarks=yticks, view=-yoffscal*10..0) ;

plot( ee-19, x=-4..4, axis[2]=[mode=log]) ;

 

Download flat.mw

The method=_d01amc seems to work fast here, in Maple 18.02. That is a compiled NAG Library function designed for semi-infinite integrands.

This is a hint that the problem was in something like an attempt at either expanding the integrand or poking at it for discontinuities in evalf/Int's internal control mechanism. For similar problems seen in the past a workaround has often been to force a suitable method or to unapply the integrand to get an operator rather than an expression.

If you raise Digits to 15 in your original, and use 5 steps, then the results appear to concur with those from forcing this method.

I was also able to simplify the integrand more compactly.

I also changed your initial constants to exact values (ie. 1 rather than 1.0, etc). It's quite often better to keep floats out of things... as much as possible until after you've finished simplifying expressions.

restart

Ts := 1:

0

``

Digits := 15:

 

 

Download int_question.mw

 

The printed results A[i,2] of the above, using epsilon=1e-8, were,

                       -1.06754466261525
                       -1.05325302242161
                       -1.03558981773781
                       -1.01386179523186
                       -0.987289540024826
                       -0.955029158917161
                       -0.916216297715231
                       -0.870041095327763
                       -0.815862074159109
                       -0.753362760097712

That ran very fast, in my Maple 18.02 on 64bit Linux.

And (only) raising Digits to 15 in your original with epsilon=1e-5, and using 5 steps, I got,

                       -1.05325301997235
                       -1.01386179283364
                       -0.955029156513430
                       -0.870041093087944
                       -0.753362758062382

acer

First 216 217 218 219 220 221 222 Last Page 218 of 336