Carl Love

Carl Love

28065 Reputation

25 Badges

13 years, 25 days
Himself
Wayland, Massachusetts, United States
My name was formerly Carl Devore.

MaplePrimes Activity


These are answers submitted by Carl Love

S:= {[-2, 1, 1, -2, 1, 1], [-2, 1, 1, 2, -1, -1], [0, -1, 1, 0, 1, -1], [0, 1, -1, 0, -1, 1]}:

Vectors u and v are multiples iff (u.v)*u = (u.u)*v. Proof of this is left as an exercise.

IsMultiple:= (u::~Vector, v::~Vector)-> ArrayTools:-IsZero(u.v*u - u.u*v):

If we exclude the zero vector, which is a multiple of everything, IsMultiple is an equivalence relation on the remaining vectors. ListTools:-Categorize will partition a list into equivalence classes.

S1:= S minus {[0$nops(S[1])]};   
L:= ListTools:-Categorize(IsMultiple, [S1[]]);
       [[-2, 1, 1, -2, 1, 1]], [[-2, 1, 1, 2, -1, -1]],
         [[0, -1, 1, 0, 1, -1], [0, 1, -1, 0, -1, 1]]

Select a representative (the first element) of each equivalence class.

map(`?[]`, {L}, [1]);
{[-2, 1, 1, -2, 1, 1], [-2, 1, 1, 2, -1, -1], [0, -1, 1, 0, 1, -1]}

The answer is, briefly, no, Maple cannot work with matrices and vectors of arbitrary, symbolic dimension. (Considering the number of questions that I see about it here, this could be one of the next frontiers for symbolic computation.) If you can deduce (with mindpower) a recurrence along a single dimension, you might be able to solve it with rsolve. Multi-dimensional recurrences are another frontier.

This is of no consequence in your program currently, but the statement

r:= Vector[(N1+1)*N2+1];

in your procedure main doesn't do anything useful. If you did want to pre-allocate a Vector, the statement should be with ( ) instead of [ ].

subsindets(
     test1r2,
     `^`,
     x-> `if`(op(2,x) = 1/2,`if`(op(1,x)::`+`, `Math.Sqrt`(op(1,x)), x), `Math.Pow`(op(x)))
);

Answering your second question, there is the command interface(prompt= ...). For example, I have the following line in my initialization file:

interface(prompt= "(**)");

Erik wrote:

I am working on a problem, which basicly is dealing with putting n identical marbles in k different boxes. I need to run through all possible combinations and do something with it.

Your "combinations" are a combinatorial structure called "compositions".

I have heard that a recursive procedure sometimes can be transformed into an iterative procedure, which will maybe make it execute faster.

That is for a small but important subset of recursions called "tail recursion". Your procedure is not that type. An example of tail recursion is

Factorial:= (n::nonnegint)-> `if`(n=0, 1, n*Factorial(n-1));

The issue with recursion is usually memory more than time.

I can't find a way to do it iteratively....

I found two pre-coded iterative ways in Maple. The first uses a very old package called combstruct. The second uses Joe Riel's (excellent) Iterator package available for free download at the Maple Applications Center here. The problem is discussed in Donald Knuth's (excellent and long-awaited) The Art of Computer Programming: Volume 4. But these methods are much slower than your recursive method. The reason to use an iterative method is that with a recursive method you essentially have to hold all combinations in memory at the same time. This is no problem for the binomial(45,5) ~ 1.3 M combinations in your example. But if you just increase your number of boxes to 10, then you have binomial(49,9) ~ 2 G, each one being an Array with 10 elements, and it's not possible in the memory of most computers.

With Joe's Iterator package, the code for your procedure is simply:

Combi:= proc(n::posint, k::posint)
local boxContent;
     for boxContent in Iterator:-BoundedComposition([n$k], n) do
          # Do something here
     end do
end proc;

With combstruct, it is

Combi:= proc(n::posint, k::posint)
local
     Iter:= combstruct[iterstructs](Composition(n+k), size= k),
     NV:= Iter[nextvalue],
     boxContent
;
     while not Iter[finished] do
          boxContent:= NV() -~ 1;
          # Do something here
     end do
end proc;

If we need to stay with this recursive procedure, then is there a way to optimize the code?

Like I said, your procedure is already quite efficient. I can make a factor-of-five improvement by using evalhf, but this will only be possible if your "# Do something here" is doable within evalhf, which is a very restricted environment. Here is your procedure done with evalhf:

Combi:= module()
export boxContent;  #A vector whose i'th entry holds the number of marbles in the i'th box

local
    Combi:= proc(k::float, boxContent::Vector, activeBox::float, marblesRemaining::float)       
    local i::integer;
        
        if activeBox < k then
            for i from 0 to marblesRemaining do
                boxContent[activeBox]:= i;
                thisproc(k, boxContent, activeBox+1, marblesRemaining-i)
            end do
        else
            boxContent[activeBox]:= marblesRemaining;
            #Do something here        
        end if;
    end proc,

    ModuleApply:= proc(n::posint, k::posint)
        boxContent:= Vector(k, datatype= float[8]);
        evalhf(Combi(k,boxContent,1,n))
    end proc
;
end module;

The call for all three procedures is simply

Combi(n, k)

I was able to do it in under a second by using simplify with side relations (see ?simplify,siderels ). I refer to it here by its syntax skeleton simplify(..., [...]). Unlike regular simplify---which may use heuristic ideas of what is simpler---simplify(..., [...]) uses a pure polynomial algebraic algorithmic approach. 

Replacing a common subexpression with a new variable is a type of side relation. Here are the substitutions that you wanted to apply:

There are a number of factors that can make a vast difference in the amount of time and memory that the algorithm uses. These factors are referred in the help page as "monomial orderings", but the help page fails to mention just how important these are. The naive approach to your problem would be

simplify(kappa, [newpar]);

I let this run several minutes until I was nearly out of memory (it was using about 4Gig) before I killed it. Then I tried applying the side relations individually from most complicated to least (this last point is very important), like this:

newkappa:= kappa:
for k from nops([newpar]) by -1 to 1 do
     newkappa:= simplify(newkappa, [newpar[k]])
end do:
newkappa;

This works nearly instantaneously (0.140 secs) and achieves a great simplification: The only variables in the result are the s's. I also did the C=4, K=7 case in 0.656 secs. So, there you go: I took it from several days of computation to under a second.

Here's the modified worksheet:

sC.mw

In the line before the else, there appears #Iota.... The # comments out the rest of the line, so it unbalances the parentheses.

Your erroneous result seems to be related to a bad simplification of factorials. After evaluating at x=1 and r=1, your summand has the subexpression

ex:= (1+k)!/k!^2/(2-k)!:

which simplifies as

simplify(ex);

Obviously this is wrong for any integer k. A workaround is assume(k::nonnegint).

simplify(ex) assuming k::nonnegint;

To get rid of all of the RootOfs at once, you could do

evalindets(newpar, RootOf, x-> allvalues(x)[1]);

First method (my preference):  second:= eval(t, test[2]);

Second method: second:= rhs(test[2][]) or rhs(test[2][1]) or op([2,1,2], test)

Unfortunately, there is no analog to the first method for selecting the first value; it only works on expressions of the form (variable = value).

The trick is to first use command combine to bring the A into the sum, and then use algsubs(A*W[n](t)=5, ...). Finally, use expand to bring constants back out of the sums.

A*B*Sum(W[n](t)*q[n](z), n= 1..infinity);
                     /infinity               \
                     | -----                 |
                     |  \                    |
                     |   )                   |
                 A B |  /     W[n](t) q[n](z)|
                     | -----                 |
                     \ n = 1                 /
combine(%);
                  infinity                   
                   -----                     
                    \                        
                     )                       
                    /     A B W[n](t) q[n](z)
                   -----                     
                   n = 1                     
algsubs(A*W[n](t)=5, %);
                      infinity           
                       -----             
                        \                
                         )               
                        /     5 B q[n](z)
                       -----             
                       n = 1             
expand(%);
                         /infinity       \
                         | -----         |
                         |  \            |
                         |   )           |
                     5 B |  /     q[n](z)|
                         | -----         |
                         \ n = 1         /

Let me know if that's enough explanation, or if you need a more-detailed example.

The default solution method rkf45 fails miserably at this task. (I used x(4)=10, u(4)=10, z(4)=10.) Using method= dverk78 gives about 5 digits of accuracy. I am trying some other methods.

restart:
sys:= {diff(u(t),t) = 4*u(t), diff(x(t),t) = 6*x(t)+2, diff(z(t),t) = x(t)+u(t)}:
Sol1:= dsolve(sys union {u(4)=10, x(4)=10, z(4)=10}, numeric):
Sol1(0);
         [t = 0., u(t) = HFloat(1.11364764942213e-6),
           x(t) = HFloat(-0.3333333330052356),
           z(t) = HFloat(7.1111113895777045)]

Sol2:= dsolve(sys union eval({u(0)=u(t), x(0)=x(t), z(0)=z(t)}, Sol1(0)), numeric):
Sol2(4);
          [t = 4., u(t) = HFloat(9.742963826802734),
            x(t) = HFloat(7.58924709847216),
            z(t) = HFloat(9.533948806446041)]

Sol3:= dsolve(sys union {u(4)=10, x(4)=10, z(4)=10}, numeric, method= dverk78):
Sol3(0);
     [t = HFloat(0.0), u(t) = HFloat(1.12536549286537e-6),
       x(t) = HFloat(-0.33333333294284945),
       z(t) = HFloat(7.111111392517568)]

Sol4:= dsolve(
     sys union eval({u(0)=u(t), x(0)=x(t), z(0)=z(t)}, Sol3(0)),
     numeric, method= dverk78
):
Sol4(4);
      [t = HFloat(4.0), u(t) = HFloat(9.999698295626974),
        x(t) = HFloat(10.000134201687256),
        z(t) = HFloat(9.999946940854626)]

In terms of elegance and simplicity of the algorithm used, it'd be hard to beat rem:

f:= x^5+a*x^4+b*x^3+c*x^2+d*x+e:
rem(f, x^3-1, x);
                          2                    
                 (c + 1) x  + (a + d) x + b + e


@casperyc 

It really surprises me that this would happen because I simplified it.

Here's a simple example that shows why not to simplify. These are the same type of expressions as occur in your Matrices (the ones in your Matrices are longer).

simplify(1+(a+b+c+d)*(e+f+g+h)*(i+j+k+l));

 a e i + a e j + a e k + a e l + a f i + a f j + a f k + a f l
 + a g i + a g j + a g k + a g l + a h i + a h j + a h k
 + a h l + b e i + b e j + b e k + b e l + b f i + b f j
 + b f k + b f l + b g i + b g j + b g k + b g l + b h i
 + b h j + b h k + b h l + c e i + c e j + c e k + c e l
 + c f i + c f j + c f k + c f l + c g i + c g j + c g k
 + c g l + c h i + c h j + c h k + c h l + d e i + d e j
 + d e k + d e l + d f i + d f j + d f k + d f l + d g i
 + d g j + d g k + d g l + d h i + d h j + d h k + d h l + 1

Before simplification, there are 2 multiplications and 10 additions; after, there are 128 multiplications and 64 additions. So the complexity increased by a factor of 16. One way to guard against this is to include the size modifier to simplify.

simplify(1+(a+b+c+d)*(e+f+g+h)*(i+j+k+l), size);
      1 + (a + b + c + d) (e + f + g + h) (i + j + k + l)

But, based on an extremely brief viewing of some of the entries in your Matrices, I guess that simplify will provide no benefit in your case, even with the size option.

The functional way:

f:= x-> exp(-exp(x^2)):
plot(z-> Int(f, 0..z), -3..3);

Note the capital I in Int and the lack of a variable appearing before the ranges.

First 356 357 358 359 360 361 362 Last Page 358 of 395