Close. I see two changes that you need to make.

First, you want each sum to start at 0 - so move the sm:=0 to between the two do statments. Likewise, you want to print the sum only after all n terms have been added - so move the print statement to between the two end do statements. Like this:

restart:
n := N*20:
for N from 1 to 30 do
  sm := 0;
  for k from 1 to n do
    sm := sm+(1/(k^2+3)):
  end do:
  print(N,n,evalf(sm)):
end do:

You do have to be a little careful about the assignment n:=N*20. If this is done when N has a value (as it will have at the end of the do loop) then n will have a specific value. Your code does work correctly assuming the restart is executed immediately before the assignment to n. Personally, I would probably move this assignment inside the first loop (together with the resetting of sm to 0).

Doug

---------------------------------------------------------------------
Douglas B. Meade  <><
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu
Phone:  (803) 777-6183         URL:    http://www.math.sc.edu

Thatnks, Preben, for making this observation.

I think the real expression that the OP wants to work with is

f := Sum( exp(-beta*(2*x0-1)-alpha*(2*x0-1)), x0=0..1 )
   / Sum( Sum( exp(-beta*(4*x0*x1-2*x0-2*x1+1)-alpha*(2*x0-1)),
               x1=0..1 ),
          x0=0..1 );
                  1                                                       
               -----                                                      
                \                                                         
                 )                                                        
                /    exp(-beta (2 x0 - 1) - alpha (2 x0 - 1))             
               -----                                                      
               x0 = 0                                                     
   -----------------------------------------------------------------------
     1   /  1                                                            \
   ----- |-----                                                          |
    \    | \                                                             |
     )   |  )                                                            |
    /    | /    exp(-beta (4 x0 x1 - 2 x0 - 2 x1 + 1) - alpha (2 x0 - 1))|
   ----- |-----                                                          |
   x0 = 0\x1 = 0                                                         /

Note that the sigma is not clearly in the numerator, not in front of the fraction. This also has no problems when it is printed.

While the previous usage was ambiguous, shouldn't Maple have responded in a different way? Maybe noting that the summation index is being used in two sums? An error? A warning?

And, as I think more carefully about this, when you evaluate the expression with alpha=0 and beta=0 before evaluating the sums the exponents lose there dependence on x0 and x1 (and x2). So, when the sum is evaluated there is no problem with the duplicated index. But, when you evaluate the sums before giving a value to alpha and beta, the indices remain in the terms of the sum. This is where Maple gets confused.

Doug

---------------------------------------------------------------------
Douglas B. Meade  <><
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu
Phone:  (803) 777-6183         URL:    http://www.math.sc.edu

How, exactly, have you entered your function?

Here is how would reproduce your expression.

f := Sum( exp(-beta*(2*x0-1)-alpha*(2*x0-1))
   / Sum( Sum( exp(-beta*(4*x0*x1-2*x0-2*x1+1)-alpha*(2*x0-1)), x1=0..1 ),
               x0=0..1 ),
          x0=0..1 );
  1                                                                          
-----                                                                        
 \                                                                           
  )                  exp(-beta (2 x0 - 1) - alpha (2 x0 - 1))                
 /    -----------------------------------------------------------------------
-----   1   /  1                                                            \
x0 = 0----- |-----                                                          |
       \    | \                                                             |
        )   |  )                                                            |
       /    | /    exp(-beta (4 x0 x1 - 2 x0 - 2 x1 + 1) - alpha (2 x0 - 1))|
      ----- |-----                                                          |
      x0 = 0\x1 = 0                                                         /

Note that I am showing the 1D input to be sure there are no issues with implied multiplication or other subtleties that do not come across in your typeset version of the expression. Also, I used the inert Sum to keep Maple from automatically simplifying the expression.

To see the value of this expression:

value(f);
(exp(beta + alpha))/(exp(-beta + alpha) + exp(beta + alpha)

   + exp(beta - alpha) + exp(-beta - alpha)) + (exp(-beta - alpha))/(exp(-beta

   + alpha) + exp(beta + alpha) + exp(beta - alpha) + exp(-beta - alpha))

When alpha=0 and beta=0, this reduces to:

eval( f, [alpha=0,beta=0] );
                              1                  
                            -----                
                             \                   
                              )          1       
                             /    ---------------
                            -----   1   /  1    \
                            x0 = 0----- |-----  |
                                   \    | \     |
                                    )   |  )    |
                                   /    | /    1|
                                  ----- |-----  |
                                  x0 = 0\x1 = 0 /
value( % );
                                      1
                                      -
                                      2

As you report, the 3d plot of f shows the value at the origin is 1/4:

plot3d( f, alpha=0..5, beta=0..5, axes=boxed );

But, if you change this slightly to force Maple to evaluate the sums, you get a plot with the correct value when alpha=0 and beta=0:

plot3d( value(f), alpha=0..5, beta=0..5, axes=boxed );

This difference does not appear at the top level:

eval( value(f), [alpha=0,beta=0] );
                                      1
                                      -
                                      2
value( eval( f, [alpha=0,beta=0] ) );
                                      1
                                      -
                                      2

This suggests to me that there might be a problem somewhere inside plot3d.

Very interesting!

Doug

---------------------------------------------------------------------
Douglas B. Meade  <><
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu
Phone:  (803) 777-6183         URL:    http://www.math.sc.edu

This would be easier if I knew exactly what your system of equations looked like, but I'll try regardless.

I think what you want to do is:

[ seq( subs( fsolve(sys), P ), a2=90..95 ) ];

For example:

sys := { P^2+Q^2=a2^2, P+Q=1 };
                         /            2    2     2\ 
                        { P + Q = 1, P  + Q  = a2  }
                         \                        / 
v2 := [ seq( subs( fsolve(sys), P ), a2=90..95 ) ];
[64.13764609, 64.84477446, 65.55190235, 66.25902980, 66.96615680, 67.67328338]

In your code, your second example refers to "solutions". Each time through the loop you refer to the same object (solutions). That is why you got the same value for v2 for each value of a2. In your first loop you solved a different system (because I assume your definition of "sys" includes the parameter a2).

I hope this is helpful.

Doug

---------------------------------------------------------------------
Douglas B. Meade  <><
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu
Phone:  (803) 777-6183         URL:    http://www.math.sc.edu

My guess is that this is a Lagrange multiplier problem with l as the Lagrange multiplier.

If I'm correct, I'll guess the problem is to maximize f(x,z)=x^2*z subject to the constraint 4*x*z-x^2=12. If this is correct, then here's one way to work with this problem:

F := x^2*z;
                                     2  
                                    x  z
G := 12-4*x*z+x^2;
                                     2        
                               12 + x  - 4 x z
eq1 := diff(F,x)=lambda*diff(G,x);
                         2 x z = lambda (2 x - 4 z)
eq2 := diff(F,z)=lambda*diff(G,z);
                               2              
                              x  = -4 lambda x
solve( {eq1,eq2,G=0}, {x,z,lambda} );
   /            /  2                 \             /  2                 \  
  { x = 2 RootOf\_Z  + 1, label = _L1/, z = -RootOf\_Z  + 1, label = _L1/, 
   \                                                                       

               1       /  2                 \\ 
    lambda = - - RootOf\_Z  + 1, label = _L1/ }
               2                             / 
allvalues( % );
    /                            1  \    /                          1  \ 
   { x = 2 I, z = -I, lambda = - - I }, { x = -2 I, z = I, lambda = - I }
    \                            2  /    \                          2  / 

Of course, as others have already shown, adding Explicit to the solve command can avoid the intermediate step:

solve( {eq1,eq2,G=0}, {x,z,lambda}, Explicit );
    /                            1  \    /                          1  \ 
   { x = 2 I, z = -I, lambda = - - I }, { x = -2 I, z = I, lambda = - I }
    \                            2  /    \                          2  / 

You'll have to decide what to do with the complex-valued solutions.

I hope this is on target and helpful,

Doug

---------------------------------------------------------------------
Douglas B. Meade  <><
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu
Phone:  (803) 777-6183         URL:    http://www.math.sc.edu

You write that "I obtained k from non dimensionalising my wave problem" and "we solve for eigenvalues,k". These seem at odds with one another.

Based on the equations you solve, the values of k are the roots of an characteristic equation for the 3x3 matrix you call eigen.

The eigenvalues of eigen will be the numbers alpha with the property that det( eigen - alpha*IdentityMatrix(3) ) = 0. Because the entries of eigen depend on k and lambda, the 3 eigenvalues of eigen will also depend on k and lambda.

Your analysis based upon the sign of the real part of the eigenvalues makes sense, but first you have to get the eigenvalues of the matrix. I'm almost definite they are not the values for k.

I think you need to work out some specific (simpler) cases first. Once your problem is well-defined, I'm pretty confident Maple will be a useful tool to complete the analysis.

Doug

---------------------------------------------------------------------
Douglas B. Meade  <><
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu
Phone:  (803) 777-6183         URL:    http://www.math.sc.edu

There is still a mismatch in dimensions.

It's not really true that the k values are eigenvalues. These are the values of the parameter that ensure the matrix eigen is singular (i.e., that 0 is an eigenvalue of eigen).

For each value of lambda, you find 4 values of k and 4 corresponding 3x3 matrices. I still do not see how you have 4 eigenvectors.

Doug

---------------------------------------------------------------------
Douglas B. Meade  <><
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu
Phone:  (803) 777-6183         URL:    http://www.math.sc.edu

You wrote:

Desired result:  For each list, produce a list of permutations of the list (as you would receive from combinat[permute]) such that in each permutation, the nonnegative integers in the list appear in ascending order from left to right and no two nonnegative integers are adjacent to one another.

But this is not what your code does. First, your list containsonly nonnegative integers:  the integers 1..n and then k copies of 0. So, you MUST vuolate the condition that "no two nonnegative integers are adjacent to one another".

As I understand your code, what you want are all of the permutations in which the integers 1..n appear in order and the k zeros are inserted in such a way that (I'm guessing now) no two zeros are adjacent to one another.

If I am understanding your problem correctly, you are correct that generating the full set of permutations of the n+k entries in the list is very wasteful. The way I would try to approach this is as follows: you know the positive integers 1..n will appear in order. All you really have to choose is where to put the k zeros. Even if you create all of these possible locations that's only C(n+k,k) instead of (n+k)!/k! - a savings of n! already.

But, we're not yet done. You still have to decide which of these lists of locations meet your second condition. For that, here's my proposal (assuming I correctly guessed that you don't want to have adjacent zeros):

Q := (n,k) -> remove( L -> has( L[2..-1]-L[1..-2], 1 ),
                      select( ListTools:-Sorted,
                              combinat:-permute(n+k,k)
                             )
                     ):
Q(5,3);

[[1, 3, 5], [1, 3, 6], [1, 3, 7], [1, 3, 8], [1, 4, 6], [1, 4, 7], [1, 4, 8], 
  [1, 5, 7], [1, 5, 8], [1, 6, 8], [2, 4, 6], [2, 4, 7], [2, 4, 8], [2, 5, 7], 
  [2, 5, 8], [2, 6, 8], [3, 5, 7], [3, 5, 8], [3, 6, 8], [4, 6, 8]]

What this does is first construct the k-element permutations of 1..n+k, keep only those that are sorted, and then removes any list that contains consecutive integers (as detected by a 1 in the list of first differences). Even this takes quite a while when n=8 and k=6, but it's doable (and 8!=40320 times faster. I can imagine other ways to construct this list of locations, but I'll let you decide exactly how you want to proceed.

I'll also leave it to you to work out the code to take each list of zero locations and insert the integers 1..n to obtain the final permutation of  your original list.

I hope this is on target, and useful,

Doug

---------------------------------------------------------------------
Douglas B. Meade  <><
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu
Phone:  (803) 777-6183         URL:    http://www.math.sc.edu

Now we are at the same point of (mis)understanding. Once you figure out exactly what you are trying to do, I'm very confident we'll be able to point you towards a workable solution.

Doug

---------------------------------------------------------------------
Douglas B. Meade  <><
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu
Phone:  (803) 777-6183         URL:    http://www.math.sc.edu

In your problem, for a given value of lambda, there are 4 values of k.

For each (lambda,k) pair there is a 3x3 matrix, with 3 eigenvalues and 3 eigenvectors. So, for each lambda, there are 12 eigenvalues and 12 eigenvectors.

But, you refer to u, v, and w as the 3 components of the "four eigenvectors".

I must be misinterpreting something. I could understand, and implement, all of this if you told me you were only interested in one of the three eigenvalues for each value of k. Should I always look at the largest eigenvalue (in absolute value, real part, ...)? the smallest eigenvalue? or something else?

Why are you surprised that the command to construct K "only displays all 'k' values for one value of lambda"? That's what I thought you wanted. Aren't these what you are referring to as k[1], k[2], k[3], and k[4]?

I hope you are understanding where I'm getting confused. Once we get past this roadblock, I'm sure all will be fine.

Doug

---------------------------------------------------------------------
Douglas B. Meade  <><
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu
Phone:  (803) 777-6183         URL:    http://www.math.sc.edu

Identifying zero can be a challenge, even for Maple.

When Georgios' ideas doesn't work, I've been known to resort to a plot. You can try plotting the two expressions, but it's often more effective to plot their difference. You might need to move the axes out of the way, with axes=boxed or axes=none as an optional argument to plot.

Once you are convinced graphically that the two expressions are equivalent, then you still have to show it algebraically.

If you have a specific situation in mind, a few more details might give us some more ideas about how you can proceed.

Doug

---------------------------------------------------------------------
Douglas B. Meade  <><
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu
Phone:  (803) 777-6183         URL:    http://www.math.sc.edu

You are mistaken here.

The plot command generates a collection of points and then connects consecutive points (ordered by the horizontally). The vertical asymptotes that appear when the (line) plot of the tangent function are actually the segments that connect the last point on the left of the discontinuity with the first poin to the right of the discontinuity.

When you ask for a point plot, the same points are produced, but they are not connected - so you don't see the "asymptotes".

When you specify discontinuities=true for a line plot, Maple tries to identify discontinuities, and separates the graph into disjoint pieces by not connecting the points on either side of a discontinuity.

This is why I suggested in my previous post that asking for an implicitplot with style=points might be a way to get a pseudo-plot of a function. But, in this case, Maple does come up with quite a few points along the asymptotes. This was a surprise to me, and one that I hope to better understand soon.

Doug

---------------------------------------------------------------------
Douglas B. Meade  <><
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu
Phone:  (803) 777-6183         URL:    http://www.math.sc.edu

Thank you for posting your worksheet. That is very helpful - but it does leave me with a few questions.

The key part of your code is:

for l from 1 to 2 do # by 0.01 do
  fsol := [fsolve(eval(detE, lambda = l), {k})];
  Ml := map(subs, fsol, lambda = l, eigen2);
  EM := map([Eigenvectors], Ml);
end do;

detE is a polynomial in k (degree 4) and lambda (degree 3) with highest total degree 4 as well. For a given value of lambda, Maple's fsolve command does not have any trouble finding the 4 roots; these are stored in fsol. These are the 4 values of k that you refer to; to get them in a more convenient form, I suggest adding:

K := map( rhs@op, fsol );

The entries of eigen2 depend on lambda and k ( sqrt(k), k, k^(3/2), and k^2 ). In fact, I notice that the lambda dependence is only on the main diagonal and always with the same term (91/100*lambda, with your parameters). You could subtract this part of the main diagonal from the matrix without changing the eigenvectors of the matrix. This, and related observations, might improve your ultimate efficiency.

Next, you substitute each of the 4 values of k and the corresponding value of lambda into the 3x3 matrix, eigen2. What you call Ml is a list of 4 3x3 matrices. Then, the associated EM is a list of 4 eigenvalue decompositions. Each of these decompositions is list containing a 3x1 matrix with the eigenvalues and a 3x3 matrix with the eigenvectors.

So, here's how I see all of your results:

K = [ k1, k2, k3, k4 ]

Ml = [ [3x3 matrix], [3x3 matrix], [3x3 matrix], [3x3 matrix] ]

EM = [ [ [3x1 vector], [3x3 matrix] ], [ [3x1 vector], [3x3 matrix] ],
       [ [3x1 vector], [3x3 matrix] ], [ [3x1 vector], [3x3 matrix] ] ]

You describe wanting to construct a 4x4 matrix from the 4 k values and the 12 elements of 4 eigenvectors.

This is where I get confused. For each value of k there are three eigenvalues and 3 eigenvectors. What, exactly, am I supposed to use for the u[1], v[1], and w[1]? for the u[2], v[2], and w[2]? for the u[3], v[3], and w[3]? for the u[4], v[4], and w[4]? Am I correct in observing that you don't use the eigenvalues for anything in your matrix?

I assume you are trying to get one such matrix for each value of lambda.

Let me conclude with some final remarks about your final matrix. All 4 columns have the same structure:

V := (k,u,v,w) -> < 3*I/10*u + k*v-w, k*u+I*v, -k*v+(3/10-k^2)*w, (k^2-7/5)*k*v+(k^2-17/10)*w >:
V(k,u,v,w);
                        [     3                    ]
                        [     -- I u + k v - w     ]
                        [     10                   ]
                        [                          ]
                        [        k u + I v         ]
                        [                          ]
                        [           /3     2\      ]
                        [    -k v + |-- - k | w    ]
                        [           \10     /      ]
                        [                          ]
                        [/ 2   7\       / 2   17\  ]
                        [|k  - -| k v + |k  - --| w]
                        [\     5/       \     10/  ]

So, once we figure out what information goes into a particular set of k[i], u[i], v[i], w[i], we can construct the final matrix as:

M := < V(k[1],u[1],v[1],w[1]) | V(k[2],u[2],v[2],w[2]) |
       V(k[3],u[3],v[3],w[3]) | V(k[4],u[4],v[4],w[4]) >;

Note that it appears to me that there is a misplaced [1] in the fourth row of the first column of your matrix.

I hope my attempts at explaining my understanding of your problem either help you find a solution to your problem, or show you where you need to correct my understanding so that I, or others, can help you.

Doug

---------------------------------------------------------------------
Douglas B. Meade  <><
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu
Phone:  (803) 777-6183         URL:    http://www.math.sc.edu

How interesting.

This exact question was asked, and answered, earlier this week on MaplePrimes.

This question is unusual enough that I have to wonder if this is part of an assignment. And, if so, if asking us for help is allowed.

Doug

---------------------------------------------------------------------
Douglas B. Meade  <><
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu
Phone:  (803) 777-6183         URL:    http://www.math.sc.edu

The different plot commands are used in very different ways. Some of them have completely different sets of options and features. It does not make sense to have one interface with every possible customizable feature (and many disabled because they are not supported for a particular command). It's probably more work to create the separate GUIs, each with only the supported options, but this does make life simpler, overall.

The real question is why some features are supported in some commands and not others.

The particular question behind this discussion deals with discontinuities. This is, in general, a hard problem. Given an equation f(x,y) = g(x,y), how - exactly - do you expect Maple to identify all singularities? Implicit plots are challenging. Maple samples a collection of points on a regular grid and tries to find the ones that come closest to satisfying the equation. That's not so bad. What becomes a problem is how to connect these dots. Oftentimes you can circumvent some discontinuities by changing the plot style to point (style=point) and increasing the number of points:

plot( tan(x), x=-5..5, y=-5..5, style=point, numpoints=200 );

This is not unreasonable. But, I have to admit I am surprised by the corresponding plot from implicitplot:

plots:-implicitplot( y=tan(x), x=-5..5, y=-5..5, style=point, grid=[50,50] );

Maybe someone else can provide an explanation for implicitplot's selection of points on the discontinuity? Is this a consequence of Maple's use of the complex-valued tangent function?

Doug

---------------------------------------------------------------------
Douglas B. Meade  <><
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu
Phone:  (803) 777-6183         URL:    http://www.math.sc.edu