The integral is indeed 2*Pi for any complex a. This follows from the Residue theorem.

Maple finds the antiderivative Ei(1, exp(I*x))  and using Newton-Leibniz concludes that the integral is 0.
Unfortunately the antiderivative is discontinuous but Maple cannot check this (at least for a symbolic a).
Such situations are common and will not disappear soon (most of them).

The system can be reduced to 2 equations eliminating e.g. x3.
It is easy to obtain  60 < x1 < 85,  50 < x2 < 70.

S12:=[floor((1/2)*x1)+floor((1/5)*x2)+floor(-(2/9)*x1-(7/18)*x2+1/3)-8,
      floor((1/5)*x1)+floor((1/2)*x2)+floor(-(1/6)*x1-(7/24)*x2+3/4)-14]:
f:=(i,j) -> eval( S12,[x1=60+i/500*25,x2=50+j/500*20]): 
M:=Matrix(500, (i,j) -> `if`(f(i,j)=[0,0],1,0)):
plots:-sparsematrixplot(M, size=[500,500],
  labels=["",""],tickmarks=[[1=60,500=85], [1=50,500=70]]);

It is possible to write a code and obtain all these regions, but probably a motivation would be necessary
(without it the problem does not seem to be interesting enough).

Increasing Digits:=15 (which is always a good idea for numerics) is enough in this case.

(Note that the summation order is not necessarily the same in I1 and I2 and `+` is not associative for floats.)

See also here.

Such 5-fold integrals containing discontinuous/nondifferentiable functions  are hard.
You must know to help Maple for them.
For example, use the method from my previous comment ==>

restart;
dd:= u=w*y..1, v=w*x..1, x=0..1, y=0..1, w=0..1:
int(abs((w*y-u)*(w*x-v)-w^2*x*y), dd);
evalf(%);

       -91/1080+4*ln(2)*(1/15)
       0.1005799889

Now, a numerical check is mandatory.

int(abs((w*y-u)*(w*x-v)-w^2*x*y)*piecewise(0 < -w*x*y+u*x and -w*x*y+u*x < 1, 1)*
piecewise(0 < -w*x*y+v*y and -w*x*y+v*y < 1, 1),
 [u = 0 .. 1, v = 0 .. 1, w = 0 .. 1, x = 0 .. 1, y = 0 .. 1],
 numeric, method=_CubaDivonne, epsilon=1e-5);

        0.100580642789240

We may conclude that the symbolic computation is correct [ 99.9% say :-) ].

restart:
with(Statistics):
#N :=3;# ... Some integer value >= 2;  
X := RandomVariable(Binomial(N, 1/2)): 
F:=CDF(X, s):
N:=3:
plot(F, s=-1..N+1, gridlines=true, axis[1]=[gridlines=N+1], discont, thickness=3);

A symmetric real matrix has orthogonal eigenvectors but the eigenvectors found by Eigenvectors are not necessarily orthogonal unless the eigenvalues are distinct or the matrix contains floats.
To obtain the orthogonal ones, use GramSchmidt.
In your case:

#  A = U . DiagonalMatrix(v) . U^+;

f := sin(x) / ( ( sin(a*x) )^a * ( sin((1-a)*x) )^(1-a) ):
MultiSeries:-limit(f, x=0,right) assuming a>0;

Same for a<0.

Maple is correct.

pde.mw

You have copied the expressions from the web page. Don't do this, because (unfortunately)  they are not designed to be copied.
Just open the help page (type ?Jacobian in Maple) and copy from there.
 

@Klausklabauter 

In this case all you need are the definitions and some Maple syntax.
Example

Maple sorts polynomials using monomial orders.
There is no such order with x^2 > y^2 > x*y.

Some tricks can be used e.g.

sort(x^2+y^2+x*y*``, order=plex(``,y),ascending);
  

but I don't like them because I prefer to use Maple for computations (and LaTeX for typesetting).

You should not expect 1.
Take the simpler A:=<1,2;3,-4>  to see why.

Maple computes the classical Cholesky decomposition A = L . L^*  (for a positive definite matrix A).  U is L^*.

It seems that you want the generalized one (the LDL decomposition)  A = LL . DD . LL^*

LL and DD can be obtained from L using:

S:=DiagonalMatrix(Diagonal(L));
DD:=S^2;
LL:=L.S^(-1);

Using limit for expressions containg floats is inherently problematic.
Compare:

limit(sin(1/3-x/3)/(3.*x-3.), x=1) ;
                         -0.1111111111
limit(sin(1/3.-x/3)/(3.*x-3.), x=1) ;
                         Float(undefined)