acer

32587 Reputation

29 Badges

20 years, 37 days
Ontario, Canada

Social Networks and Content at Maplesoft.com

MaplePrimes Activity


These are answers submitted by acer

For this example you could set Normalizer:=simplify which appears to avoid so-called hidden zeroes while pivoting as well as keep the reduced row entries in check.

The default of Normalizer = normal is well-suited to entries all as rational polynomials. Using simplify is a big hammer, but figuring out the smallest hammer that will suffice is difficult in general.

acer

The "alloc" means allocated by the kernel. In Maple 18 the kernel can actually free memory to the OS when doing garbage collection, so that number can now go both up and down. In recent versions the kernel will also free memory to the OS on `restart`, at which point this number reverts to the size of memory allocated at kernel launch and initialization.

The "used" figure is an indicator of how much memory has being processed by the garbage collection mechanism. It's similar to a running total of size of memory of objects (no longer referenced) and reclaimed for re-use by the kernel. Since the same bits of memory may be reclaimed and re-used over and over, that running total can be a number much higher than the total memory of the machine. Also, this value is not reset to zero upon `restart`.

I suppose that you know that you can turn these messages off in the CLI, using kernelopts(printbytes=false) .

acer

S:= Array( 1..4, -6..3, 1..5 );

And use whatever ranges for the indices that you'd like. 

acer

In 2D Math input mode try typing it in with the keystrokes for & then \ and then ^ in order to escape the caret.

[edited] Oh, it seems there is a mention of this as an item lower down in the table for 2DMathShortcutKeys.

acer

Are you looking for purely real solutions?

Unless I transcribed something wrong when changing 45.5 to exact 91/2 then I don't quite see how there can be any purely real solution to the system.

Problem_7_Instructor_4_rat.mw

acer

Could it be converted from set to logic notation? I wonder... how. (How best, I suppose. `subs` might suffice.)

restart:                                                                   

with(Logic):

Equivalent( (A &and &not B) &and (A &and &not C), A &and &not (B &or C ) );

                                     true

So, if that is acceptable, how could this be improved,

restart:                                                        
ee := (A minus B) intersect (A minus C) = A minus (B union C):  
ff := subsindets( subs( `union`=`&or`, `intersect`=`&and`, ee ),
                  specfunc(`minus`),                            
                  z->op(1,z) &and &not op(2,z) ):               
evalindets( ff, `=`, z->Logic:-Equivalent(op(z)) );             

                                     true

I have a hazy recollection that this sort of question might have been asked before on this site.

acer

For Matrix M, you can do it with,

map( convert, M, exp);

acer

I don't really understand how you have (or have not) defined M[n+1], as expressions or operators, etc. But just using the expression you gave, here are some ideas for a multidimensional Array and plotting slices of it.

I used Maple 18, FWIW.

restart:

for n from 1 to 10 do
  R[n]:=n^2;
end do:

A:=Array( 1..10,1..10,1..5,1..2,
          (n,b,theta,m) -> (b*(theta^m)*R[n])/((theta^m)+R[n]^m)
          , datatype=float[8]  # only if it evaluates to floats
         ):

with(plots):

opts := transparency=0.3, style=surface:
all_fixed_m_plots := [ seq( surfdata(A[n,..,..,1], 1..10, 1..5,
                               labels=["b", "theta_m", "n"],
                               opts, color=COLOR(HUE,n/10)),
                      n=1..10 ) ]:

display( all_fixed_m_plots );
all_b_m_plots := [ seq( [ seq( surfdata( A[..,..,b,m],
                                         caption=sprintf("b=%a m=%a",b,m) ),
                               m=1..2 ) ], b=1..5 ) ]:

display(Array( all_b_m_plots ));

display(Array([
    surfdata(A[..,..,3,1]),
    surfdata(A[..,7,..,2]),
    surfdata(A[..,3,5,..])
              ]));
# The following are animations. Click on them and then use the menubar to play.
plots:-animate(plots:-surfdata,[ 'A[i,..,..,1]' ], i=[$1..10]);
plots:-animate(plots:-surfdata,[ 'A[..,..,i,1]' ], i=[$1..5]);
plots:-animate(plots:-surfdata,[ 'A[..,i,..,2]' ], i=[$1..10]);
plots:-animate(plots:-surfdata,[ 'A[..,3,..,i]' ], i=[$1..2]);

acer

Use the plots:-display command to combine multiple plots.

acer

This is a FAQ, and it is and example of what is sometimes called premature evaluation. Under Maple's normal evaluation rules the first argument to sum is evaluated up front, before it does it's computation.

So in your usage what sum sees as the first argument is just the result of Comm(j), which is 0 because the test j=0 does not return true. So your call just does sum( 0, j=0..0 ).

You can instead use delay quotes (unevaluation quotes) on the call Comm(j), or you can use add which has special evaluation rules, or you can use the Physics package to redefine sum as Edgardo has mentioned.

restart:

Comm := n -> `if`( n=0, 1, 0 ):

Comm(j);

                                       0

sum( Comm(j), j=0..0 );

                                       0

sum( 'Comm'(j), j=0..0 );

                                       1

add( Comm(j), j=0..0 );

                                       1

acer

We'll have to take you on your word, that you really do want operators rather than just expressions, for the phi[j,k].

In that case you can do this using `unapply`.

N_x:=5;
    N_y:=4;
    N_elx:=N_x-1;
    N_ely:=N_y-1;
    h_x:=(x_e-x_s)/N_elx;
    h_y:=(y_e-y_s)/N_ely;
    x_n:=[seq(x_s+j*h_x,j=0..N_elx)];
    y_n:=[seq(y_s+j*h_y,j=0..N_ely)];

    # Partition of unity.
    for j from 2 by 1 to N_x-1 do
        for k from 2 by 1 to N_y-1 do
            phi[j,k]:=unapply( (x-x_n[j-1])*(x-x_n[j+1])
                               *(y-y_n[k-1])*(y-y_n[k+1])
                               /((x_n[j]-x_n[j-1])*(x_n[j]-x_n[j+1])
                                *(y_n[k]-y_n[k-1])*(y_n[k]-y_n[k+1])),
                               [x,y]);
        od;
    od;

acer

restart:

with(LinearAlgebra):

lign1:=x-y+2*z=1:
lign2:=(-2)*x+y+z=0:
lign3:=(-4)*x+y+7*z=2:
lign4:=3*x-2*y+z=1:

T := GenerateMatrix([lign1, lign2, lign3, lign4],[x, y, z], augmented);

lign5:=x+y+z=6:
R := GenerateMatrix([lign5],[x, y, z], augmented);

# One way
Q:=<R,T[2..,..]>;

# Another way
Q:=copy(T):
Q[1,..]:=R:
Q;

# A third way, overwriting T's first row
T[1,..]:=R:
T;

acer

On Maple 15, you might do something like this (using a reworking of Carl's Display).

restart;

PF := module() local ModuleApply, Fact, handler, rgb2hex;
  ModuleApply := proc(N::posint, {color::list(nonnegint):=[0,100,0]}, $)
    local e, h, res;
    uses ColorTools;
    if not ( nops(color)=3 and andmap(t->t>=0 and t<=255,color) ) then
      error "expecting color option as a list of 3 integers between 0 and 255";
    end if;
    h := rgb2hex(color);
    res := Fact(N)[2];
    N=`+`(seq(`if`(rhs(e)=0, [][],
                   handler(ithprime(op(lhs(e))),h)
                   * `if`(rhs(e)=1, 1,
                          nprintf(`#mn("%ld")`,rhs(e)))
                   ), e = res ));
  end proc;
  Fact := proc(a::posint)
    local cst, obj, res;
    uses numtheory, Optimization;
    cst := add(x[i], i = 1 .. pi(prevprime(a))) <= 6;
    obj := add(x[i]*ithprime(i), i=1 .. pi(prevprime(a)))-a;
    res := LPSolve(obj, {cst ,obj>=0},
                   assume={nonnegative,integer});
  end proc;
  handler := proc(n::posint, s::string)
    nprintf(`#mn("%ld", mathcolor=%a)`,n,s);
  end proc;
  rgb2hex:=proc(rgb::list(integer))
    local h,i,r;
    h:="#";
    for i to 3 do
      h:=h,iquo(rgb[i],16,'r'),r;
    end do;
    cat(op(subs([0="0",1="1",2="2",3="3",4="4",5="5",6="6",7="7",8="8",
                 9="9",10="A",11="B",12="C",13="D",14="E",15="F"],[h])));
  end proc;
end module:
  
PF(30);

seq(print(PF(2*i,color=[255,0,0])), i=2..20);

 

In Maple 18 the color argument can be made nicer, as follows.

restart:
PF := module() local ModuleApply, Fact, handler;
  Fact := proc(a::posint)
    local cst, obj, res;
    uses numtheory, Optimization;
    cst := add(x[i], i = 1 .. pi(prevprime(a))) <= 6; 
    obj := add(x[i]*ithprime(i), i=1 .. pi(prevprime(a)))-a; 
    res := LPSolve(obj, {cst ,obj>=0},
                   assume={nonnegative,integer});
  end proc:
  ModuleApply := proc(N::posint, {color::name:=':-Green'}, $)
    local e, h, res;
    uses ColorTools;
    h := RGB24ToHex(NameToRGB24(sprintf("%s",color)));
    res := Fact(N)[2];
    N=`+`(seq(`if`(rhs(e)=0, [][],
                   handler(ithprime(op(lhs(e))),h)
                   * `if`(rhs(e)=1, 1,
                          nprintf(`#mn("%ld")`,rhs(e)))
                  ), e = res ));
  end proc:
  handler := proc(n::posint, s::string)
    nprintf(`#mn("%ld", mathcolor=%a)`,n,s);
  end proc:
end module:
PF(30);
seq(print(PF(2*i,color=red)), i=2..20);

acer

f := int(cosh(a*x)*cos(b*x), x):

simplify( convert( f, trigh) );


                 cosh(a x) sin(b x) b + sinh(a x) a cos(b x)
                 -------------------------------------------
                                    2    2                  
                                   a  + b                   

acer

This is about 33 times faster on my 64bit Linux system. How fast do you need it?

I expect that it could be made at least 10-20 times faster again, if the iterated action were burned into a evalhf'able or Compilable procedure that acted inplace on an Array of `n` values. And that could be split with Threads:-Task.

It seems to work for the example parameter values you gave. I would rather know that the approach works for other parameter values in which you are interested, before going to town on further optimization.

I put your original procedure at the end.

restart:

expr:=tan(theta_n) - C_a*Z_0*(-v^2*(Pi*n-theta_n)^2/x_l^2
             +omega_a^2)*x_l/(v*(Pi*n-theta_n)):

Equate([omega_a,v,x_l,C_a,Z_0],[100e9, 1e8, 0.3, 20e-14, 50.0]):
ee:=simplify(eval(expr,%));

                                                         2              
                          0.003333333333 (Pi n - theta_n)  - 300.0000000
     ee := tan(theta_n) + ----------------------------------------------
                                          Pi n - theta_n                

# We hope that this provides a contraction mapping which
# converges quickly:
ff:=simplify(map(arctan,op(1,ee)=-op(2..-1,ee)))
      assuming theta_n>-Pi/2, theta_n<Pi/2;

                          /                               2              \
                          |0.003333333333 (Pi n - theta_n)  - 300.0000000|
   ff := theta_n = -arctan|----------------------------------------------|
                          \                Pi n - theta_n                /

F:=N->unapply(eval(evalf(rhs(ff)),n=N),theta_n,proc_options=[hfloat]):

# Applying F to n produces a procedure, which we can
# iterate at initial values of theta_n.
F(117);

          proc(theta_n)                                           
          option hfloat;                                          
                                                                  
            -1.*arctan((0.003333333333*(%1)^2 - 300.0000000)/(%1))
          end proc;                                               
                                                                  
                                                                  
          %1 := 367.5663405 - 1.*theta_n                          

# For most all initial theta_n, F(117) converges in a
# small number of iterations. It might be that 10 iterations
# is enough, in general.
seq( (F(117)@@i)(11.0), i=0..10 );

     11.0, -0.334174918606374, -0.389865087453356, -0.390129467304736,

       -0.390130722191024, -0.390130728147372, -0.390130728175643,

       -0.390130728175778, -0.390130728175778, -0.390130728175778,

       -0.390130728175778

# Let's see whether many initial guesses converge to the same value
# after 10 iterations.
plot( (F(117)@@10), 0.1..11.0 );

# Let's always use 1.0 as the initial theta_n guess, and produce a
# plot for various values of n.
forget(evalf):
CodeTools:-Usage( plot( '(F(n)@@10)(1.0)', n=1..1000) );

memory used=6.09MiB, alloc change=0 bytes, cpu time=151.00ms, real time=146.00ms, gc time=0ns

restart:
# Now lets put it into a procedure.

new_get_theta_n_array:=proc(max_n::integer, omega_a::float, v::float,
                            x_l::float, C_a::float, Z_0::float)

    local ff, F, theta_n_array;
    ff:=simplify(expand(arctan(C_a*Z_0*(-v^2*(Pi*n-theta_n)^2
                                /x_l^2+omega_a^2)*x_l/(v*(Pi*n-theta_n)))));
    F:=N->unapply(eval(evalf(ff),n=N),'theta_n',proc_options=[hfloat]);
    # Here we always start iterating from guess theta_n=1.0, but
    # it might be even faster to start iterating from the previous
    # n's final iterate. If we keep it as is, though, we could Thread:-Task
    # this action and split across all the values of n to be computed.
    # Also, if convergence were ever slow then we could consider making
    # the number of iterations ba according to some additional proc argument
    # rather than fixed at 10 iterations.
    theta_n_array:=Array([seq( '(F(n)@@10)(1.0)', n=1..1000 )],datatype=float[8]);
    if ArrayNumElems(theta_n_array) <> max_n then

        printf("Bad Array Dimensions! Got too many or not enough solutions.");

        theta_n_array:="CHECK: get_theta_n_array()":

    end if;

    theta_n_array;

end:

CodeTools:-Usage( new_get_theta_n_array(1000,100e9,1e8,0.3,20e-14,50.0) ):

memory used=9.21MiB, alloc change=32.00MiB, cpu time=151.00ms, real time=158.00ms, gc time=20.00ms

plots:-listplot(%);

restart:

get_theta_n_array:=proc(max_n::integer, omega_a::float, v::float,
                        x_l::float, C_a::float, Z_0::float)

    local theta_n_array:=Array(

        select(x->Re(x)<evalf(Pi/2) and Re(x)>=-evalf(Pi/2),
               [seq(

        #NOTE: careful with fsolve - in some cases returns
        #      unevaluated equation

                    fsolve(tan(theta_n) - C_a*Z_0*(-v^2*(Pi*n-theta_n)^2
                             /x_l^2+omega_a^2)*x_l/(v*(Pi*n-theta_n)),
                           theta_n) , n=1..max_n)]

               , real)

                               , datatype=float[8]);

    if ArrayNumElems(theta_n_array) <> max_n then

        printf("Bad Array Dimensions! Got too many or not enough solutions.");

        theta_n_array:="CHECK: get_theta_n_array()":

    end if;

    theta_n_array;

end:

ans:=CodeTools:-Usage( get_theta_n_array(1000,100e9,1e8,0.3,20e-14,50.0) ):

memory used=308.13MiB, alloc change=117.85MiB, cpu time=5.28s, real time=5.28s, gc time=130.00ms

plots:-listplot(ans);

 


Download itsme0.mw

acer

First 237 238 239 240 241 242 243 Last Page 239 of 338