@Markiyan Hirnyk No, poles alone don't imply the divergence of an integral, of course. But in my experience poles in multi-dimension integrals almost always confound the numerical integration routines. Consider this relatively simple case:

evalf(Int(1/(x^2+y^2)^(1/4), [x= -1..1, y= -1..1]));

When will you no longer be a "Maple newbie"? You've been a member of MaplePrimes for nearly five years (indeed, 86 days longer than I've been a member), and you've been a regular contributor for over three years. You have a reputation of 1340, which places you at #30 on the all-time list of highest-reputation members.

There are numerous points where the denominator is zero.

@John Fredsted Thanks. It's a reconciliation of the dichotomous nature of Boolean logic with the trichotomous nature of total ordering. The third case, a[k] = b[k], is the unstated else clause of the if statement, which forces either comparison of the next pair (if there is a next pair) or comparison of the list (or Vector) lengths (if there's no next pair).

@magiblot 

To answer your first question "How does the computation change....": If you don't specify an initial point, the algorithm uses some method to guess at an initial point, or perhaps just picks one at random. But there must be an initial point for the iteration whether you supply it or not.

To answer you second question "Would any point....": No, not necessarily. The relationship between the initial point and point which is converged upon can be extremely complex, even for simple polynomial objective functions. For any given extremum, the set of initial points for which the algorithm converges to that extremum is called the extremum's basin of attraction (for that algorithm). These basins can be fractals of unimaginable intricacy and complexity, in which case they're called strange attractors. There are Wikipedia articles on both of these italicized terms.

But given an extremum, there is some open set containing that extremum that is a subset of its basin of attraction. Were it not so, there'd be no practical value in having an initial point option. In other words, if you start close enough to the extremum, then you'll find it. There's no a priori way of knowing how close "close enough" is.

@magiblot Yes, you may solve for the corner points directly rather than considering those as separate Lagrange multiplier problems. This is certainly the way to go if you're solving these problems by hand. If you're using Maple, I have a technique, which I'll post shortly, which does the interior and all the Lagrange multiplier problems for all the boundaries, in all dimensions (including the 0D corner points), with a single solve command.

@Earl I'm not sure if this is relevant, but is there any significant difference between the edges of the Catalan surface and the Mobius strip, which you called single-edged? The edges are topologically homeomorphic, although the surfaces are not.

Indeed, it is true that solve ignores real assumptions, even if you use option useassumptions. This is mentioned in the help. Usually it is better to formulate assumptions as inequalities and include them in the set of equations passed to solve. I don't see how that can help in this case. 

@mj1751 Using the paper, you need to change F[k+1]:= solve(eq) to F[k+1]:= solve(F[k+1] = eq).

@mj1751 Looking at the paper's equation now, I see that you've entered it incorrectly. Try entering it the way that it appears in the paper; do not attempt to manually solve it for F[k+1]. I don't have time right now, but I'll get back to it in about 3 - 4 hours.

@mj1751 But don't you see that the results we are getting differ from the desired results by a factor of exactly  3? Somewhere in the process of the paper being written, typeset, or read, a simple factor of 3 is being missed.

@nikkarkian If the equation is degree-4 in v, and you use option explicit with the solve command, then it'll work.

solve(eq, {v}, explicit);

@kuwait1 In the example that you present, the eval does nothing. The same thing would happen had you used simply Re(f(3,1)).

Quick advice (I don't have much time right now): Use Sum not sum. Use evalf. Use a lower value of Digits. What's the point of using 50 digits when a is only accurate to 4 digits?

@awass You misunderstand what Re is intended to do. It doesn't extract the real part of an expression, except in the trivial case that the expression is in the form a + b*I with a and b real (* footnote). Rather, Re is a placeholder that represents the real part of an expression, and that real part may be extracted by some later command, very likely evalc. So since Re((2-I*x)^4) returns unevaluated (even under assuming x::real), that expression is what is passed to expand. Then the procedure `expand/Re` tells expand that the argument to Re can be expanded. Once expanded, the expression is in the form a + b*I, the Re is evaluated again, and this time it works.

*footnote: If a procedure `Re/f` exists, then the call Re(f(x)) calls that procedure if f(x) returns unevaluated.

@magiblot Each of your problems is specifed over a region in the plane, i.e., a 2-dimensional space. Each constraint, when converted to an equality, is a line, a 1-dimensional space. Each pair of these equalities when solved together give a finite set of points (or no solution at all), which is a 0-dimensional space.

If the original problem were specifed in a 3D region, you'd need to check the interior of the region for critical points, then check each 2D surface that forms the boundaries, then each 1D curve that forms the boundaries of those, then finally the 0D corner points.

Do you understand?