acer

32405 Reputation

29 Badges

19 years, 346 days
Ontario, Canada

Social Networks and Content at Maplesoft.com

MaplePrimes Activity


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

The ImportMatrix command should be able to handle this.

acer

There are several ways to do this. Here are two.

restart:                         

f:=RandomTools:-Generate(list(integer(range=-30..50), 5), makeproc):

seq( f(), i=1..4 );                                                 

    [14, -25, 28, 13, 7], [38, -4, -14, 3, -13], [21, 25, 12, -6, 29],

    [-17, 19, 16, -23, 15]

restart:                                                            

f:=rand(-30..50):                                                   

seq( [seq(f(),j=1..5)], i=1..4 );                                   

    [14, -25, 28, 13, 7], [38, -4, -14, 3, -13], [21, 25, 12, -6, 29],

    [-17, 19, 16, -23, 15]

I expect that there are several faster ways.

I expect that you could also use Statistics:-Sample and a discrete distribution.

Here's another one,

restart:

f:=(n,k)->convert(LinearAlgebra:-RandomMatrix(k,n,generator=-30..50),
                  listlist)[]:

f(5,4);

    [38, 43, 9, 23, 36], [-4, -18, -5, -9, -12], [-22, -24, -5, 14, 21],

    [28, 21, 46, -11, 15]

Do you need it very, very fast? Do you need to create a large number of these lists? Do you need to create them all up front, before you work on them?

acer

Change that last line to this,

plots:-display(p1,p2);

The display comamnd you're trying to use is part of the plots package. Either reference it as plots:-display or first load the package using the command with(plots): and then call it as just display.

acer

Perhaps I don't understand what your trying to accomplish. But why not use exact values for x?

for x from -1 by 0.5 to 1 do print(x, evalf(0^(x^2), 20)) end do; 

                                    -1, 0.

                                   -0.5, 0.

                             0., Float(undefined)

                                    0.5, 0.

                                    1.0, 0.

for x from -1 by 1/2 to 1 do print(x, evalf(0^(x^2), 20)) end do;

                                    -1, 0.

                                   -1/2, 0.

                                     0, 1.

                                    1/2, 0.

                                     1, 0.

for x from -1 by 1/2 to 1 do print(x, evalf(0.0^(x^2), 20)) end do;

                                    -1, 0.

                                   -1/2, 0.

                                     0, 1.

                                    1/2, 0.

                                     1, 0.

And of course there are lots of ways that you can make the behavior conditional upon the value of x. Examples include,

for x from -1 by 0.5 to 1 do
  print(x,`if`(x=0,1.0,evalf(0^(x^2), 20)));
end do;

                                    -1, 0.

                                   -0.5, 0.

                                    0., 1.0

                                    0.5, 0.

                                    1.0, 0.

for x from -1 by 0.5 to 1 do                
  print(x,evalf(0^(`if`(x=0,0,x^2)), 20));
end do;

                                    -1, 0.

                                   -0.5, 0.

                                    0., 1.

                                    0.5, 0.

                                    1.0, 0.

acer

Change the square brackets in the assignment to Eqns (inside plotsol1s) to round brackets.

Round brackets are expression delimiters. Square brackets create a list.

acer

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