Joe Riel

9530 Reputation

23 Badges

20 years, 24 days

MaplePrimes Activity


These are answers submitted by Joe Riel

A more efficient technique is to use the Iterator package available in the Maplesoft Application Center.  For example

with(Iterator):
m := 2:
n := 3:
P := MixedRadixTuples([3 $ m*n]):
(h,g) := ModuleIterator(P):
# Create a Matrix alias to the Vector used by the iterator
M := ArrayTools:-Alias(g(),[m,n]):
in P do printf("%{}d\n\n",M); end do:
0 0 0
0 0 0
0 0 0
0 0 1
0 0 0
0 0 2
0 0 1
0 0 0
0 0 1
0 0 1
...

In the following code I used the BinaryGrayCode export of the Iterator package, available in the Maple Application Center, to implement the Ryser algorithm. 

 

permanentRyser4 := proc(M::Matrix)
local G, V, aj, g, h, i, j, k, m, n, s, sj, sig;
m, n := op(1, M);
if m <> n then
error "expecting a square Matrix, got dimensions %1, %2", m, n
end if;
   G := Iterator:-BinaryGrayCode(n,'ret_delta');
V := Vector(n);
s := 0;
(h,g) := ModuleIterator(G);
j := g();
h();
sig := -1;
while h() do
sj := sign(j[1]);
aj := abs(j[1]);
for i to n do
V[i] := V[i] + sj * M[i,aj];
end do;
s := s + sig*mul(V[k], k=1..n);
sig := -sig;
end do;
expand((-1)^n*s);
end proc:
M := Matrix(20,rand(10)):
CodeTools:-Usage(permanentRyser4(M));
memory used=1.73GiB, alloc change=32.00MiB, cpu time=30.68s, real time=30.77s, gc time=1.64s
15220180698446432599207679453939
CodeTools:-Usage(LinearAlgebra:-Permanent(M));
memory used=3.18GiB, alloc change=0.73GiB, cpu time=49.49s, real time=44.12s, gc time=11.44s
15220180698446432599207679453939

I prefer an analytical (Lagrangian) approach. 

# Consider a ramp of mass M with angle from horizontal alpha
# that can slide along the x-axis and a block of mass m
# that slides along the ramp. Coefficient of friction
# between block and ramp is mu1, between ramp and ground
#
# M = mass of ramp
# m = mass of block
# alpha = ramp angle
# mu1 = coefficient of friction between ramp and block
# mu2 = coefficient of friction between ramp and ground
# g = gravitational acceleration
# Fx = x-component of external force applied to ramp
# Fy = y-component of external force applied to ramp

# Choose generalized coordinates u and v
#
# u = the x-distance from origin to tip of ramp,
# v = the distance of block along ramp, from tip.
# Assign an alias and some shortcuts
alias(u=u(t), v=v(t));
du := diff(u,t):
dv := diff(v,t):
# The kinetic energy is
T := 1/2*M*du^2 + 1/2*m*((du+dv*cos(alpha))^2 + (dv*sin(alpha))^2);
Tu := frontend(diff, [T,du]);
# normal force between ramp and ground
Fn2 := M*g + m*(g + diff(v,t,t)*sin(alpha)) - Fy;
eq2 := simplify(diff(Tu,t)) = Fx - signum(du)*mu2*Fn2;
Tv := frontend(diff, [T,dv]);
# normal force between block and ramp
Fn1 := m*g*cos(alpha) - m*diff(u,t,t)*sin(alpha);
eq1 := simplify(diff(Tv,t)) = -m*g*sin(alpha) - signum(dv)*mu1*Fn1;
eqs := [eq1,eq2];
eqs := map(lhs-rhs, eqs);
eqs := map(collect, eqs, [diff(u,t,t),diff(v,t,t)]);

 Here I substitute some arbitrary values for the parameters and initial conditions and plot u and v

deqs := subs([m=1,M=10,g=1,alpha=Pi/4,Fy=0,Fx=0,mu1=0,mu2=0.1]
, {seq(eq=0,eq=eqs), eval(u,t=0)=0, eval(v,t=0)=0, D(op(0,u))(0)=1, D(op(0,v))(0)=1});
dsol := dsolve(deqs,numeric);
plots:-odeplot(dsol, [[t,u],[t,v]], 0..10);

 

MapleSim currently does not handle complex values.  At best you might get sections of a plot where a signal is undefined when it is complex. If the complex signal is not used by the model, you might be able to extract the data from the model and then post process it.

@Bendesarts How do you want that function to effect the amplitude?  As a factor?  You could use a multiply block to multiply the output of the Chirp with your signal.

One way to generate the amplitude signal is via a Modelica code. Attach the Modelica Custom Component template, replace the content of the code block with the following (I've modified your equations slightly, but you can change this to whatever is desired):

model Amplitude
    import SI = Modelica.SIunits;
    extends Modelica.Blocks.Interfaces.SO;
    parameter SI.Duration Tstart = 1;
    parameter SI.Duration Tstop = 10;
protected
    parameter SI.Duration Trise = Tstop - Tstart;
equation
    y =  if time < Tstart then 0
        elseif time < Tstart + Trise/2 then
            2*((time - Tstart)/Trise)^2
        elseif time < Tstop then
            -2*((time - Tstart - Trise/2)/Trise)^2 + (time - Tstart - Trise/2)/Trise + 1/2
        else
            1;
end Amplitude;

Click the Generate Component from Source button and drag the component into the worksheet. A better way to do this may be to incorporate the Chirp block in this Modelica code, however, I'm not sure what you eventually want to achieve.

In Maple, a list (similarly a set and lots of other things) is not a mutable structure. When a list is modified an entirely new list is created. The original list continues to exist in memory, until no longer referenced and garbage collected.  The subsop procedure isn't necessarily anymore efficient at creating a new list. While direct assignment is syntactically allowed for lists with less than 100 elements, it probably should not be. I strongly discourage using the direct assignment method. If you need a mutable structure, use one: a one-dimensional Array or Vector is much more efficient for this.

As for your question, yes, Maple will raise an error if you attempt to modify an element of a list with more than 100 elements:

L := [0$101]:L[1] := 1;

Error, assigning to a long list, please use Arrays

You may be wondering why this is the case. It has to do with the fundamental design of Maple, particularly with its Simplification table. Expressions in Maple are stored in the Simplication table. For immutable objects, only one expression is stored for a given entity.  For example, there will only ever be one list identical to [1,2,3] stored in the table; similarly the expression a+3. If you later create another list identical to [1,2,3], Maple will use the original list (more accurately, its address). This is why immutable objects are immutable; if you could modify one, you'd modify every variable that was assigned that object.

This also explains why building an expression sequence (or list or set) an element at a time is O(n^2). Consider the loop

L := NULL:
for i to 100 do
  L := L,i;
end do:

When complete, L is assigned the expression sequence 1,2,...,100. However, to get there, the intermediate expression sequences were created:
L := 1
L := 1,2
L := 1,2,3
....

Each such expression sequence is stored in the Simplification table.  The storage for the first assignment is, say, 1 unit. The storage for the second is 2 units.  The storage for the entire lot is 1+2+...+100 units = 5050, which is substantially larger than the 100 units needed by the final sequence. The memory allocated for the intermediate expressions may eventually be reclaimed by the garbage collector, but there is an immediate cost of allocating, as there is a cost to lookup/insert each intermediate expression into the Simplification table.

For this trivial example one can avoid the additional expense by doing

L := seq(1..100):

Depending on what is being computed, that isn't always feasible, at least not easily. An alternative is to allocate an Array, either a fixed size, if the total number of slots is known, or dynamically increase its size as needed.  Then one could do

A := Array(1..100):
for i to A do
   A[i] := i;
end do:
L := seq(A[i],i=1..100):

The memory/time cost of that approach is O(n) [here n = 100].

I agree that rtable_scanblock is poorly documented.  The call you want is

rtable_scanblock(A,[],(val,ind,res)->(`if`(val>=3,[ind,op(res)],res)),[]);

The procedure appends the index to the list.  The empty list at the end is to initialize 'res'.

This is an inefficient way to operate; it builds up a list an element at a time, which is O(n^2), where n is the number of elements in the final list. A more efficient approach is

V := Array(1..0): k := 0:
rtable_scanblock(A,[],proc(val,ind) global k,V; if val>=3 then k:=k+1; V(k):=ind; fi; end);
                               ''
V;
             [[1,3],[2,1],[2,2]]

 

Simpler is to use unix tools such as uniq.  However, here's a Maple approach (I haven't actually tested it):

count := proc(file :: string)
local cnt_consec, cnt_uniq, line, max_consec, prev;
    prev := NULL;
    cnt_consec := 0;
    cnt_uniq := 0;
    max_consec := 0;
    do
        line := readline(file);
        if line = prev then
            cnt_consec := cnt_consec+1;
        else
            prev := line;
            if cnt_consec > max_consec then
                max_consec := cnt_consec;
            end if;
            cnt_consec := 0;
            cnt_uniq := cnt_uniq+1;
            if line = 0 then break; end if;
        end if;
    end do;
    (cnt_uniq, max_consec);
end proc:

For the particular example you chose, all three use LinearAlgebra:-Diagonal, so you might as well call it directly. MTM:-diag is equal to ArrayTools:-Diagonal, which is a somewhat more general version of LinearAlgebra:-Diagonal, i.e. it applies to Arrays as well as Matrices and can be used to convert a Vector to a diagonal Matrix. MTM is a package for calling Matlab from Maple, some of the commands work without Matlab.

The OP asks when it is appropriate to use seq vs map.  For typical usage this is more a stylistic choice than one of efficiency.  Sometimes the choice depends on what you want to do with the result: if it has the same form as the "input" then map is natural.

Note that if you are using seq inside a procedure (my normal usage), then you should declare the index variable as a local. That isn't strictly required, however, if you don't there can occasionally be issues.  For example

(**) protect(k);
(**) proc() seq(k, k=1..3) end();
Error, (in unknown) attempting to assign to `k` which is protected.  Try declaring `local k`; see ?protect for details.

 

 

The most likely cause of this is that a parameter value was entered incorrectly. Try viewing the Modelica source for the model, it may be easy to spot it in the listing.

The way to do that is to convert the parameters to variables in the source.  A symbol declared as a parameter cannot change during a simulation, though its initial value can be computable.

Use collect with the factor transformation applied to the coefficients:

collect(a, mu(F,H), factor);

 

The real desire here, I think, is to keep intact 1/tau - tau.  One approach is to freeze, factor, and thaw:

ex := 1/(1/tau-tau):
thaw(map(factor, subs(ex = freeze(ex), Expr)));

I modified this; reread the post and realized the user wanted to print the output to a file.  This does that.

Compute := proc( data :: string, cols :: posint, ofile :: string )
local A, B, cnt, fd1, fd2, fmt, i, j, line, p, row, x;

    fd1 := fopen(data,'READ');
    fd2 := fopen(ofile, 'WRITE');
    fmt := StringTools:-Repeat("%f",cols);
    A := Matrix(0,1..cols); # working Matrix
    cnt := 0;               # number of Matrices read
    i := 0;                 # number of rows in current Matrix
    do
        line := readline(fd1);
        if line = 0 or line = "" or line[1] = "#" then
            if i > 0 then
                # Matrix was started; assume it is complete.
                cnt := cnt+1;
                # Specify the (possible) submatrix of A;
                B := A[1..i,..];
                # Compute charpoly and extract coefficients as Vector
                p := LinearAlgebra:-CharacteristicPolynomial(B,x);
                p := PolynomialTools:-CoefficientVector(p,x,'termorder=reverse');
                fprintf(fd2,"%{c,}a\n",p); # see ?rtable_printf for details
                i := 0;
            end if;
            if line = 0 then
                # end of file, exit loop
                break;
            end if;
        else
            # Parse row and append to Matrix
            i := i+1;
            row := sscanf(line, fmt);
            for j to cols do
                A(i,j) := row[j];
            end do;
        end if;
    end do;
    fclose(fd1,fd2);
    cnt; # return number of Matrices
end proc:
Compute("read-data-compute-charpoly.dat",3,"read-data-compute-charpoly.poly");

The resulting file contains

1,-15.,-18.,0
1,-6.,13.,-28.
1,-3.,1.,3.

 

First 26 27 28 29 30 31 32 Last Page 28 of 114