The procedure P solves your problem:

P:=proc(k, i)

local L, s, M, K;

if not (k>=2*i and 2*i>=4) then

error `should be k>=2*i>=4` fi;

if not type(k, posint) or not type(i, posint) then

error `should be k is integer and i is integer` fi;

L:=[];

for s from k-2*i to k-i do

M:=combinat[composition](s+i, i);

K:=seq([seq(M[j,l]-1, l=1..i)], j=1..nops(M));

L:=[op(L), K];

od;

L;

end proc;

Example:

P(8, 3);

Matrix А isn't singular:

LinearAlgebra[Determinant](A);

                 0.03893229110      

with(geom3d):

Cd:=coordinates(M):

subs(x=Cd[1], y=Cd[2], z=Cd[3], Equation(P));

Объём_части_эллипсои.mws

Read more on the Russian site  http://forum.exponenta.ru/viewtopic.php?t=12179  and in the attached file 

For the functional coloring all surfaces should be given by explicit equations, rather than using graphical primitives.

A := plot3d([u, v, 1], u = 0 .. 1, v = 0 .. 1, axes = normal, numpoints = 10000, style = surface, color = sqrt(u^2+v^2+1)):

B := plot3d([u, 0, v], u = 0 .. 1, v = 0 .. 1, axes = normal, numpoints = 10000, style = surface, color = sqrt(u^2+v^2)):

C := plot3d([0, u, v], u = 0 .. 1, v = 0 .. 1, axes = normal, numpoints = 10000, style = surface, color = sqrt(u^2+v^2)):

E := plot3d([1, u, v], u = 0 .. 1, v = 0 .. 1, axes = normal, numpoints = 10000, style = surface, color = sqrt(u^2+v^2)):

F := plot3d([u, 1, v], u = 0 .. 1, v = 0 .. 1, axes = normal, numpoints = 10000, style = surface, color = sqrt(u^2+v^2+1)):

G := plot3d([u, v, 0], u = 0 .. 1, v = 0 .. 1, axes = normal, numpoints = 10000, style = surface, color = sqrt(u^2+v^2)):

plots[display](A, B, C, E, F, G, view = [-.27 .. 1.27, -.27 .. 1.27, -.27 .. 1.27], lightmodel = light4, orientation = [20, 70]);

It is clear that points equidistant from the origin have the same color. Detailization of coloring depends on numpoints option.

plot3d([ sqrt(5+4*sin(t))*sin(s), sqrt(5+4*sin(t))*cos(s),cos(t)], t=0..2*Pi, s=0..2*Pi, scaling=constrained, color=(5+4*sin(t))^2+cos(t)^2, axes=normal, style=surface,  view=[-3.7..3.7, -3.7..3.7, -1.7..1.7], numpoints=10000, lightmodel=light4, orientation=[20, 60]);

Firstly, are created a list of all pairs of diagonals for such squares, and then in a loop filled with the remaining entries. The total number of such  Latin squares is equal to 960. For illustration  3 squares are displayed .

restart:

L:=combinat[permute](5):

ND:=0: Diag:=[]:

for i in L do

A:=convert(i, set) minus {i[3]}:

M:=combinat[permute](A):

for j in M do

if i[1]<>j[1] and i[2]<>j[2] and i[4]<>j[3] and i[5]<>j[4] and i[1]<>j[4] and i[2]<>j[3] and i[4]<>j[2] and i[5]<>j[1] then ND:=ND+1: Diag:=[op(Diag),[[seq(i[k], k=1..5)], [j[1], j[2], i[3], j[3], j[4]]]]: fi:

od: od:

ND: # The number of pairs of diagonals

Diag: # The list of pairs of diagonals

B:={1, 2, 3, 4, 5}:

N:=0: List:=[]:

for i from 1 to nops(Diag) do

for a in B minus {Diag[i,1,1], Diag[i,2,5], Diag[i,1,2], Diag[i,2,2]} do

for b in B minus {Diag[i,1,1], Diag[i,2,5], Diag[i,1,3]} do

for c in B minus {Diag[i,1,1], Diag[i,2,5], Diag[i,2,4], Diag[i,1,4]} do

for d in B minus {Diag[i,1,2], Diag[i,2,4], Diag[i,1,1], Diag[i,2,1]} do

for e in B minus {Diag[i,1,2], Diag[i,2,4], Diag[i,1,3]} do

for f in B minus {Diag[i,1,2], Diag[i,2,4], Diag[i,2,5], Diag[i,1,5]} do

for g in B minus {Diag[i,1,3], Diag[i,1,1], Diag[i,2,1]} do

for h in B minus {Diag[i,1,3], Diag[i,1,2], Diag[i,2,2]} do

for l in B minus {Diag[i,1,3], Diag[i,2,4], Diag[i,1,4]} do

for m in B minus {Diag[i,1,3], Diag[i,2,5], Diag[i,1,5]} do

for p in B minus {Diag[i,2,2], Diag[i,1,4], Diag[i,1,1], Diag[i,2,1]} do

for q in B minus {Diag[i,2,2], Diag[i,1,4], Diag[i,1,3]} do

for r in B minus {Diag[i,2,2], Diag[i,1,4], Diag[i,2,5], Diag[i,1,5]} do

for u in B minus {Diag[i,2,1], Diag[i,1,5], Diag[i,1,2], Diag[i,2,2]} do

for v in B minus {Diag[i,2,1], Diag[i,1,5], Diag[i,1,3]} do

for w in B minus {Diag[i,2,1], Diag[i,1,5], Diag[i,2,4], Diag[i,1,4]} do

if nops({Diag[i,1,1],a,b,c,Diag[i,2,5]})=5 and nops({d,Diag[i,1,2],e,Diag[i,2,4],f})=5 and nops({g,h,Diag[i,1,3],l,m})=5 and nops({p,Diag[i,2,2],q,Diag[i,1,4],r})=5 and nops({Diag[i,2,1],u,v,w,Diag[i,1,5]})=5 and nops({Diag[i,1,1],d,g,p,Diag[i,2,1]})=5 and nops({a,Diag[i,1,2],h,Diag[i,2,2],u})=5 and nops({b,e,Diag[i,1,3],q,v})=5 and nops({c,Diag[i,2,4],l,Diag[i,1,4],w})=5 and nops({Diag[i,2,5],f,m,r,Diag[i,1,5]})=5 then N:=N+1: List:=[op(List), Matrix([[Diag[i,1,1],a,b,c,Diag[i,2,5]], [d,Diag[i,1,2],e,Diag[i,2,4],f], [g,h,Diag[i,1,3],l,m], [p,Diag[i,2,2],q,Diag[i,1,4],r], [Diag[i,2,1],u,v,w,Diag[i,1,5]]])]: fi:

od: od: od: od: od: od: od: od: od: od: od: od: od: od: od: od: od:

N;  # N -  the total number of squares,  List - the list of all squares

List[1], List[100], List[960];

A possible way :

with(LinearAlgebra):

M:=RandomMatrix(6);

G:=DiagonalMatrix(RandomVector(6));

M - DiagonalMatrix([seq(M[i,i], i=1..6)])+DiagonalMatrix([a$6]);

M - DiagonalMatrix([seq(M[i,i], i=1..6)])+G;

a:=2: b:=4: c:=7: d:=2:

if b>3 then z:=a:  y:=c elif b<3 then z:=d:  y:=a end if:

z;  y;

2

7

Add  numpoints  option:

plot3d(sin(sqrt(x^2+y^2))/sqrt(x^2+y^2), x=-3..3, y=-3..3, axes=normal, numpoints=3000);

Use

Optimization[Minimize]  command.

The same as in Mathematica (based on Markiyan's idea):

A := seq(plots[tubeplot]([x, y, x^2-y^2, y = -1 .. 1, radius = 0.02], scaling = constrained, color = grey, style = surface), x = -1 .. 1, .2):

B := seq(plots[tubeplot]([x, y, x^2-y^2, x = -1 .. 1, radius = 0.02], scaling = constrained, color = grey, style = surface), y = -1 .. 1, .2):

plots[display](A, B, orientation = [120, 60], lightmodel = light4);

It is easy to write a procedure that makes the same for any surface.

See this link  http://www.maplesoft.com/applications/view.aspx?SID=33406! The installation package includes detailed instructions for installation!

Example:

A:=Matrix([[1, 0.00009], [2, 1]]):

map(fnormal,A, 10, 1e-4);    

       [ 1   0. ]

       [ 2   1  ]

Of course, the constant a should not be an integer, otherwise the denominator would be 0. Unfortunately, Maple does not convert the sum of trigonometric functions in the product, so the denominator should be simplified by hand:

simplify(algsubs(1-cos(2*Pi*a)=2*sin(Pi*a)^2,-2*Pi*sin(Pi*a)/(-1+cos(2*Pi*a))));

                                                             Pi/sin(Pi*a)