acer

32622 Reputation

29 Badges

20 years, 43 days
Ontario, Canada

Social Networks and Content at Maplesoft.com

MaplePrimes Activity


These are answers submitted by 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


Z := (N^2*sin(N*q+A)*h*f[n]+y[n]*N^2*sin(N*q+A)
      +N*cos(N*q+A)*h*g[n]+g[n]*sin(N*q+A)-g[n]
      *sin(N*h+N*q+A))/(N^2*sin(N*q+A)):

 

applyrule(cos(a::anything)/sin(a::anything)=cot(a),
          frontend(expand,[Z]));

h*f[n]+y[n]+cot(N*q+A)*h*g[n]/N+g[n]/N^2-g[n]*sin(N*h+N*q+A)/(N^2*sin(N*q+A))

Download cotrule.mw

Now, did you also prefer to have cos(N*h)+cot(N*q+A)*sin(N*h) instead of sin(N*h+N*q+A)/sin(N*q+A) ?

acer

This is a very common mistake. The formal parameters A, omega, and k of your original procedure SW are not the same as the (in this example, global) A, omega, and k in the expressions W1, W2, and W.

Here are two ways to make it work.

restart;
with(plots):

W1:=A*cos(omega*t-k*x):
W2:=A*cos(omega*t+k*x):
W:=W1+W2:

SW:=(A,omega,k)->animate(plot,[eval([W1,W2,W],[:-A=A,:-omega=omega,:-k=k]),
                               x=-4..4,color=[red,green,blue],
                               scaling=constrained],t=0..5,frames=10):

display(SW(2,2*Pi,5),insequence);

AltSW:=unapply('animate'(plot,[[W1,W2,W],
                               x=-4..4,color=[red,green,blue],
                               scaling=constrained],t=0..5,frames=10),
               [A,omega,k]):

display(AltSW(2,2*Pi,5),insequence);

acer

If you have your heart set on using a RealRange, then you can test the condition using the is command. Note that if...then will utilize evalb in your original, which is insufficient.

evalb( 0.5 in RealRange(0,1) );

                                                0.5  in  RealRange(0, 1)

is( 0.5 in RealRange(0,1) );   

                                                          true

If you don't need to mess about with Open() end-points then you might find it simpler to do it as follows,

0.5>=0 and 0.5<=1;

                                                          true

If the RealRange is being passed to you as a result of another computation you could also try it like so. Again, this doesn't account yet for the Open endpoints cases -- if you have such then using is might be simpler,

r := RealRange(0,1):  # given elsewhere...

0.5>=op(1,r) and 0.5<=op(2,r);

                                                           true

acer


 

We could of course convert the results (in power) to MW by by an explicit call to the convert command.

 

restart:

p := 43.5e5 * Unit(kg*m^2/s^3):

simplify(p);

0.435e7*Units:-Unit('W')

convert(p, units, MW);

4.350000000*Units:-Unit('MW')

with(Units:-Standard):

p + 11.2e8 * Unit(g*m^2/s^3);

5470000.000*Units:-Unit('W')

convert(%, units, MW);

5.470000000*Units:-Unit('MW')

 

In Maple 18 or 2015 we could use the UseUnit command to set the preference for the dimension of power.

 

restart:

Units:-UseUnit(MW):  # needs Maple 2015

p := 43.5e5 * Unit(kg*m^2/s^3):

simplify(p);

4.350000000*Units:-Unit('MW')

with(Units:-Standard):

p + 11.2e8 * Unit(g*m^2/s^3);

5.470000000*Units:-Unit('MW')

 

An older way: use a new systemn that is like SI, but with MegaWatts.

 

restart:

Units:-AddSystem(MySys,Units:-GetSystem(SI),MW);

Units:-UseSystem(MySys);

p := 43.5e5 * Unit(kg*m^2/s^3):

simplify(p);

4.350000000*Units:-Unit('MW')

with(Units:-Standard):

p + 11.2e8 * Unit(g*m^2/s^3);

5.470000000*Units:-Unit('MW')

 


Download megawatts.mw

acer

I chose the height for the 2D part to be minz+(minz-maxz)*0.1 using the min and max of the 3D part. It has not meaning, in terms of the expression `ee` that is plotted as the surface. Change that if you wish.

restart:

with(plots):
with(plottools):

ee := sin(x)*y^2:

P3:=plot3d(ee, x=-Pi..Pi, y=-1..1, shading=zhue):
minz,maxz := op(plottools:-getdata(P3)[2][3]);

-1., 1.

P2:=densityplot(-ee, x=-Pi..Pi, y=-1..1, style=surface, colorstyle=HUE):

# The shadings are not an exact match. This is not satisfactory. Read on...

display( P3,
         transform((x,y)->[x,y,minz+(minz-maxz)*0.1])(P2) );

# Here, the shadings match.
# The style of the surface is inherited from that of the densityplot.

display( transform((X,Y)->[X,Y,eval(ee,[x=X,y=Y])])(P2),
         transform((X,Y)->[X,Y,minz+(minz-maxz)*0.1])(P2) );

# Or we could correct the raised surface to have a patch style.

display( subsindets(transform((X,Y)->[X,Y,eval(ee,[x=X,y=Y])])(P2),
                    specfunc(anything,STYLE), u->STYLE(PATCH)),
         transform((X,Y)->[X,Y,minz+(minz-maxz)*0.1])(P2) );

# Or we could create the densityplot with the `patch` style.

P2b:=densityplot(-ee, x=-Pi..Pi, y=-1..1, style=patch, colorstyle=HUE):
display( transform((X,Y)->[X,Y,eval(ee,[x=X,y=Y])])(P2b),
         transform((X,Y)->[X,Y,minz+(minz-maxz)*0.1])(P2b) );

PC:=contourplot(-ee, x=-Pi..Pi, y=-1..1, color=black):

display( subsindets(transform((X,Y)->[X,Y,eval(ee,[x=X,y=Y])])(P2),
                    specfunc(anything,STYLE), u->STYLE(PATCH)),
         transform((X,Y)->[X,Y,minz+(minz-maxz)*0.1])(display(P2,PC)) );

display( transform((X,Y)->[X,Y,eval(ee,[x=X,y=Y])])(P2),
         transform((X,Y)->[X,Y,minz+(minz-maxz)*0.1])(display(P2,PC)) );

display( subsindets(transform((X,Y)->[X,Y,eval(ee,[x=X,y=Y])])(P2),
                    specfunc(anything,STYLE), u->STYLE(PATCHCONTOUR)),
         transform((X,Y)->[X,Y,minz+(minz-maxz)*0.1])(display(P2,PC)) );

 


Download surfcont.mw

acer

Inside your procedure which calls Color(), try either replacing that with ColorTools:-Color() , or add the line uses ColorTools; at the start of the procedure, or put that use...end use inside your procedure.

Does your procedure color_my_rgb call Color?

restart;

f:=proc(nm)
  local c;
  c:=Color(nm);
  c[];
end proc:

# This does not work.
use ColorTools in
  f("Orange");
end use;
                              Color("Orange")[]

restart;

f:=proc(nm)
  local c;
  use ColorTools in
    c:=Color(nm);
  end use;
  c[];
end proc:

f("Orange");
                         1.00000000, 0.64705882, 0.

restart;

f:=proc(nm)
  local c;
    c:=ColorTools:-Color(nm);
  c[];
end proc:

f("Orange");
                         1.00000000, 0.64705882, 0.

restart;

f:=proc(nm)
  uses ColorTools;
  local c;
    c:=Color(nm);
  c[];
end proc:

f("Orange");
                         1.00000000, 0.64705882, 0.

acer

First 219 220 221 222 223 224 225 Last Page 221 of 339