acer

9 years, 18 days


These are answers submitted by acer

ideas

October 17 2014 acer 10356
0 3

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

plots:-display

October 16 2014 acer 10356
0 3

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

acer

premature evaluation

October 16 2014 acer 10356
2 1

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

unapply

October 16 2014 acer 10356
1 0

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

row replacement

October 15 2014 acer 10356
0 3
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

Maple 15

October 15 2014 acer 10356
1 1

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

convert

October 15 2014 acer 10356
1 3
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

iterated map

October 14 2014 acer 10356
3 8

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

Ok, so part of the instruction for this question is that you use the LinearAlgebra:-LeastSquares command. (Perhaps your instructor wants you to know how to set up such problems by hand...)

restart:

r:=Vector[row]([1,cos(x),sin(x),cos(2*x),sin(2*x)]);

                r := [1, cos(x), sin(x), cos(2 x), sin(2 x)]

varvec:=Vector([a,b,c,d,e]);

                                          [a]
                                          [ ]
                                          [b]
                                          [ ]
                                varvec := [c]
                                          [ ]
                                          [d]
                                          [ ]
                                          [e]

# Here is the formula we wish to use to approximate y(x).
# Notice that it has 5 unknown parameters. I'll refer to those a the 5 unknowns, below.

r . varvec;

              a + cos(x) b + sin(x) c + cos(2 x) d + sin(2 x) e

# Here are the given x and y data values.
# Notice that there are 9 data pairs.

xdata:=Vector([1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0,9.0]):
ydata:=Vector([2.33,0.0626,-2.16,-2.45,-0.357,2.21,2.75,0.636,-2.45]):

# If we plug a particular xdata value into the formula, and set that equal
# to the corresponding ydata value, then we have an equation that is linear
# in the 5 unknowns.

eval( r, x=xdata[1] ) . varvec = ydata[1];

       a + 0.5403023059 b + 0.8414709848 c - 0.4161468365 d + 0.9092974268 e = 2.33

# We could do the same thing for all 9 data pairs. That would give us 9
# equations that are linear in the 5 unknowns. That would be an overdetermined
# linear system, which could also be expressed in Matrix-Vector form and be solved
# using the LeastSquares command.

#
# First, here is the left-hand side of system    M . varvec = ydata
#

M := Matrix(9,5, (i,j)-> eval(r[j], x=xdata[i]) );

         [1   0.5403023059   0.8414709848  -0.4161468365   0.9092974268]
         [                                                             ]
         [1  -0.4161468365   0.9092974268  -0.6536436209  -0.7568024953]
         [                                                             ]
         [1  -0.9899924966   0.1411200081   0.9601702867  -0.2794154982]
         [                                                             ]
         [1  -0.6536436209  -0.7568024953  -0.1455000338   0.9893582466]
         [                                                             ]
    M := [1   0.2836621855  -0.9589242747  -0.8390715291  -0.5440211109]
         [                                                             ]
         [1   0.9601702867  -0.2794154982   0.8438539587  -0.5365729180]
         [                                                             ]
         [1   0.7539022543   0.6569865987   0.1367372182   0.9906073557]
         [                                                             ]
         [1  -0.1455000338   0.9893582466  -0.9576594803  -0.2879033167]
         [                                                             ]
         [1  -0.9111302619   0.4121184852   0.6603167082  -0.7509872468]

# Now find the least squares solution, the result of which are approximate
# values for the 5 unknowns in `varvec`.
#

C := LinearAlgebra:-LeastSquares( M, ydata );

                             [0.0000897022769103875]
                             [                     ]
                             [     2.61796597219174]
                             [                     ]
                        C := [     1.07445085124322]
                             [                     ]
                             [  -0.0131699592951893]
                             [                     ]
                             [   0.0796678388801930]

# Using those 5 solution values in place of the 5 unknown parameters in `varvec` gives us
# the solution: a formula for y as a function of x.
#

sol := r . C;
         sol := 0.0000897022769103875 + 2.61796597219174 cos(x)

            + 1.07445085124322 sin(x) - 0.0131699592951893 cos(2 x)

            + 0.0796678388801930 sin(2 x)

plots:-display(
  plots:-pointplot( <xdata|ydata>, symbolsize=15, color=blue ),
  plot( sol, x=min(xdata)..max(xdata) )
              );

acer

results

October 11 2014 acer 10356
1 1

How about just passing `values` from `MCN` to `Display`?

restart:

Display:= proc(fin::list(indexed=nonnegint),values)
  local S,e;
  S := `%+`(seq(`if`(rhs(e)=0, [][],
                     `if`(rhs(e)=1, values[op(lhs(e))],
                          `%*`(values[op(lhs(e))], rhs(e)))
                     ),
                e = fin));
  InertForm:-Display(InertForm:-Value(S) = S)
end proc:

MCN := proc(values::list)
  local a, N, const, tval, res;
  N := nops(values); const := add(x[i], i = 1 .. N) <= 30;
  tval := add(x[i]*values[i], i = 1 .. N); 
  for a from 100 by -1 to 1 do 
    try
      res := Optimization:-LPSolve(tval, {tval = a, const},
                                   assume = {nonnegint}, depthlimit = 200) 
    catch "no feasible integer point found": 
      return a,Display(res[2],values);
    catch "no feasible point found for LP subproblem":
      return a,Display(res[2],values);
    catch:
      error;
    end try;
  end do;
end proc:

MCN([6, 9, 20]);

acer

LinearAlgebra

October 11 2014 acer 10356
0 3

You have overlooked that you are trying to access a command from a package.

Either call it as LinearAlgebra:-RowSpace(...) or load the package before using it. For example,

with(LinearAlgebra):
RowSpace(...);

acer

backslash to escape

October 09 2014 acer 10356
1 0

In 2D math input mode, type \ (ie. the backslash key) before typing _ (ie. underscore, usually Shift-minus keystrokes) in order to have the underscore symbols be displayed literally.

That should work before and after Maple 17. In Maple 16 and earlier typing a single underscore would lower the prompt to the subscript position. In Maple 17 and later that happens when two underscores are typed in succession.

You can often use the resulting name as intended, even when it's inadvertantly displayed as a subscript in the input. But it looks weird afterwards, and may be confusing if you come back to the Document at a later date. Using the backslash helps it look as intended, with underscores displayed literally.

See the help topic tasks/building/annotate/mathModeShortcut , or the table entry on escaping characters in the help topic worksheet/documenting/2DMathShortcutKeys .

acer

AllSolutions option

October 09 2014 acer 10356
1 0

The `int` command provides an option `AllSolutions` which allows it to attempt and return a piecewise solution with conditions on parameters in the range of integration.

For this example this produces results equivalent to what one gets under assumptions. The result for the case x<0 is the same as what I get under assumptions if I also call `simplify`(...under the same assumption on x).

If you set up the call to `int` or `Int` yourself then this can be done by adding the option `AllSolutions`. The alternate spelling as `allsolutions` is also ok.

restart:

ee := Int((x-tau)^(0.5-1)*(tau^2+(tau^(1.5))*(8/(3))-tau-(2)*tau^0.5),
                         tau=0..x, allsolutions):

value(ee);

value( convert(ee,rational) );

If the integral call was produced by some other command then you can also enable this behavior with an environment variable (so that you don't have to mess with `op` in order to shoehorn that option into the `int` or `Int` call).

restart:

ee := Int((x-tau)^(0.5-1)*(tau^2+(tau^(1.5))*(8/(3))-tau-(2)*tau^0.5),
                         tau=0..x):

_EnvAllSolutions:=true:

value(ee);

value( convert(ee,rational) );

That all works in my Maple 18.01.

The result in the piecewise solution for condition x<=0 evaluates to zero for x=0 (the empty integration range). When I tried this example using assumptions on x instead of the `AllSolutions` option I could get a direct result for the assumption x<0 but not for x≤0.

acer

Expand Execution Group

October 07 2014 acer 10356
1 3

In a Document, apply a context-menu operation to either 2D Math input of output. This produces a black right-arrow and a result.

Now select the right-arrow and the output beside it with the mouse cursor. Then from the GUI's main menu select View , and Expand Execution Group

The hidden command that was applied, as the context-menu action, should be revealed.

On my Linux, those main menu choice are available as the keyboard acceleration, Alt-v then x

acer

errors

October 06 2014 acer 10356
2 2

Welcome to the site. Welcome to Maple! It's not as hard to use as it might seem at first.

There were several errors in your code, as you guessed.

There were mismatched start- and endbraces in both lambda1 and lambda2.

The assignments to R[0] and R[1] were just equations with `=`, not actual assignments with `:=`.

You were assigning to `R` after assigning to R[0] and R[1]. Assigning to `R` will wipe out those assignments to R[0] and R[1].

You were mixing up expressions and operators. You had `R` defined as an operator, but then were using it with `subs` as if it were an expression.

Here is a way, with RR as an expression.

restart:
f := .25; g := 1; R[0] := 5; R[1] := 4.8; nDays := 30;
lambda1 := (1 - f + sqrt((1 - f) + 4*g*f))*1/2;
lambda2 := (1 - f - sqrt((1 - f) + 4*g*f))*1/2;
k1 := (-lambda2*R[0] + R[1] )/(lambda1 - lambda2);
k2 := ( lambda2*R[0] - R[1])/(lambda1 - lambda2);
RR := k1*lambda1^N + k2*lambda2^N;
for n from 2 to nDays do
  evalf(subs(N=n,RR));
end do;

Here is a way, with RR as an operator.

restart:
f := .25; g := 1; R[0] := 5; R[1] := 4.8; nDays := 30;
lambda1 := (1 - f + sqrt((1 - f) + 4*g*f))*1/2;
lambda2 := (1 - f - sqrt((1 - f) + 4*g*f))*1/2;
k1 := (-lambda2*R[0] + R[1] )/(lambda1 - lambda2);
k2 := ( lambda2*R[0] - R[1])/(lambda1 - lambda2);
RR := N -> k1*lambda1^N + k2*lambda2^N;
for n from 2 to nDays do
  evalf(RR(n));
end do;

I suggest you study Robert Israel's great post of Top Ten Maple Errors.

acer

1 2 3 4 5 6 7 Last Page 1 of 100