Magnet lattices for particle accelerators (the sequence of focusing and bending magnets and drift sections making up a beam line) are often designed numerically using computer codes like MAD that model each beam-line element using either a matrix description or numeric integration, or some other algorithm. The Lattice package for Maple—recently published in Maplesoft's Application Center—allows modelling such beam lines using the full algebraic power of Maple. In this way, analytic solutions to beam-optics problems can be found in order to establish feasibility of certain solutions or study parameter dependencies. Beam lines are constructed in an intuitive way from standard beam-optical elements (drifts, bends, quadrupoles etc.) using Maple's object-oriented features, in particular Records which represent individual elements or whole beam lines. These Records hold properties like the first-order transport matrices, element length and some other properties plus, for beam lines, the elements the line is composed of. Operations like finding matched Twiss (beam-envelope) functions and dispersion are implemented and the results can be plotted. The Lattice package knows about beam matrices and can track such matrices as well as particles in a beam. The tracking function (map) is implemented separately from the first-order matrices thus allowing nonlinear or even scattering simulations; at present sextupole elements, compensating-wire elements and scattering-foil elements take advantage of this feature. Standard textbook problems are programmed and solved easily, and more complicated ones are readily solved as well €”within the limits set by Maple's capabilities, memory and computing time. Many investigations are possible by using standard Maple operations on the relevant properties of the beam lines defined. An interesting possibility is, e.g., using Maple's mtaylor command to truncate the transfer map of a beam line to a desired order, making it a more manageable function.

The package can output a beam-line description in MAD8 format for further refinement of the solution, cross-checking the results and possibly more detailed tracking.

Developed initially for my own use I have found the package useful for a number of accelerator design problems, teaching at the US Particle Accelerator School as well as in modelling beam lines for experiments I have been involved in. Presently at Version 1.0 the package still has limits; esp. the higher-order descriptions are not yet as complete as desirable. Yet already many practical design problems can be tackled in great detail. The package is backwards compatible to Maple 15. It can be found in the Application center together with help database files (old and new style) and a Users Guide.

An example showing the flavor of working with the package is attached (FODO_example.mw). It analyzes a FODO cell, the basic cell used in many ring accelerators.

Uli Wienands, aka Mac Dude

The method of solving underdetermined systems of equations, and universal method for calculating link mechanisms. It is based on the Draghilev’s method for solving systems of nonlinear equations. 
When calculating link mechanisms we can use geometrical relationships to produce their mathematical models without specifying the “input link”. The new method allows us to specify the “input link”, any link of mechanism.

Example.
Three-bar mechanism.  The system of equations linkages in this mechanism is as follows:

f1 := x1^2+(x2+1)^2+(x3-.5)^2-R^2;
f2 := x1-.5*x2+.5*x3;
f3 := (x1-x4)^2+(x2-x5)^2+(x3-x6)^2-19;
f4 := sin(x4)-x5;
f5 := sin(2*x4)-x6;

Coordinates green point x'i', i = 1..3, the coordinates of red point x'i', i = 4..6.
Set of x0'i', i = 1..6 searched arbitrarily, is the solution of the system of equations and is the initial point for the solution of the ODE system. The solution of ODE system is the solution of system of equations linkages for concrete assembly linkage.
Two texts of the program for one mechanism. In one case, the “input link” is the red-green, other case the “input link” is the green-blue.
After the calculation trajectories of points, we can always find the values of other variables, for example, the angles.
Animation displays the kinematics of the mechanism.
MECAN_3_GR_P_bar.mw 
MECAN_3_Red_P_bar.mw

(if to use another color instead of color = "Niagara Dark Orchid", the version of Maple <17)

Method_Mechan_PDF.pdf

The method of solving underdetermined systems of equations, and universal method for calculating link mechanisms. It is based on the Draghilev’s method for solving systems of nonlinear equations. 
When calculating link mechanisms we can use geometrical relationships to produce their mathematical models without specifying the “input link”. The new method allows us to specify the “input link”, any link of mechanism.

Example.
Three-bar mechanism.  The system of equations linkages in this mechanism is as follows:

f1 := x1^2+(x2+1)^2+(x3-.5)^2-R^2;
f2 := x1-.5*x2+.5*x3;
f3 := (x1-x4)^2+(x2-x5)^2+(x3-x6)^2-19;
f4 := sin(x4)-x5;
f5 := sin(2*x4)-x6;

Coordinates green point x'i', i = 1..3, the coordinates of red point x'i', i = 4..6.
Set of x0'i', i = 1..6 searched arbitrarily, is the solution of the system of equations and is the initial point for the solution of the ODE system. The solution of ODE system is the solution of system of equations linkages for concrete assembly linkage.
Two texts of the program for one mechanism. In one case, the “input link” is the red-green, other case the “input link” is the green-blue.
After the calculation trajectories of points, we can always find the values of other variables, for example, the angles.
Animation displays the kinematics of the mechanism.
MECAN_3_GR_P_bar.mw 
MECAN_3_Red_P_bar.mw

(if to use another color instead of color = "Niagara Dark Orchid", the version of Maple <17)

Method_Mechan_PDF.pdf

The method of solving underdetermined systems of equations, and universal method for calculating link mechanisms. It is based on the Draghilev’s method for solving systems of nonlinear equations. 
When calculating link mechanisms we can use geometrical relationships to produce their mathematical models without specifying the “input link”. The new method allows us to specify the “input link”, any link of mechanism.

Example.
Three-bar mechanism.  The system of equations linkages in this mechanism is as follows:

f1 := x1^2+(x2+1)^2+(x3-.5)^2-R^2;
f2 := x1-.5*x2+.5*x3;
f3 := (x1-x4)^2+(x2-x5)^2+(x3-x6)^2-19;
f4 := sin(x4)-x5;
f5 := sin(2*x4)-x6;

Coordinates green point x'i', i = 1..3, the coordinates of red point x'i', i = 4..6.
Set of x0'i', i = 1..6 searched arbitrarily, is the solution of the system of equations and is the initial point for the solution of the ODE system. The solution of ODE system is the solution of system of equations linkages for concrete assembly linkage.
Two texts of the program for one mechanism. In one case, the “input link” is the red-green, other case the “input link” is the green-blue.
After the calculation trajectories of points, we can always find the values of other variables, for example, the angles.
Animation displays the kinematics of the mechanism.
MECAN_3_GR_P_bar.mw 
MECAN_3_Red_P_bar.mw

(if to use another color instead of color = "Niagara Dark Orchid", the version of Maple <17)

Method_Mechan_PDF.pdf

In this course you will learn automatically using Maple course Statics applied to civil engineering especially noting the use of components properly. Let us see the use of Maple to Engineering.

Static_for_Engineering.mw

(in spanish)

Atte.

Lenin Araujo Castillo

Ambassador of Maple

#Most  dediction of  depth of field of optical lens  involves various simplification,  hence cannot be used  for  close up photography.  With Maple, it is easy to obtain  precise   Depth of field  formuar for optical lens  without  any simplification

> restart; h := H = F^2/(N*coc)+F; E0 := 1/d+1/D0 = 1/F; E1 := 1/(d+e)+1/D2 = 1/F; E2 := a/(d+e) = coc/delta; E3 := a = F/N; eq := {E0, E1, E2, E3, h}; var := {D2, N, coc, d, delta}; e := -delta; sol1 := solve(eq, var); t1 := op(sol1)[1]; Dfar := op(t1)[2]; e := delta; sol2 := solve(eq, var); t2 := op(sol2)[1]; Dnear := op(t2)[2];

>

>
>

>
>
>

IntegerPoints2  procedure generalizes  IntegerPoints1  procedure and finds all the integer points inside a bounded curved region of arbitrary dimension.  We also use a brute force method, but to find the ranges for each variable  Optimization[Minimize]  and   Optimization[Maximize]  is used instead of  simplex[minimize]  or  simplex[minimize] .

Required parameters of the procedure: SN is a set or a list of  inequalities and/or equations with any number of variables, the Var is the list of variables. Bound   is an optional parameter - list of ranges for each variable in the event, if  Optimization[Minimize/Maximize]  fails. By default  Bound  is NULL.

If all constraints are linear, then in this case it is recommended to use  IntegerPoints1  procedure, as it is better to monitor specific cases (no solutions or an infinite number of solutions for an unbounded region).

Code of the procedure:

IntegerPoints2 := proc (SN::{list, set}, Var::(list(symbol)), Bound::(list(range)) := NULL)

local SN1, sn, n, i, p, q, xl, xr, Xl, Xr, X, T, k, t, S;

uses Optimization, combinat;

n := nops(Var);

if Bound = NULL then

SN1 := SN;

for sn in SN1 do

if type(sn, `<`) then

SN1 := subs(sn = (`<=`(op(sn))), SN1) fi od;

for i to n do

p := Minimize(Var[i], SN1); q := Maximize(Var[i], SN1);

xl[i] := eval(Var[i], p[2]); xr[i] := eval(Var[i], q[2]) od else

assign(seq(xl[i] = lhs(Bound[i]), i = 1 .. n));

assign(seq(xr[i] = rhs(Bound[i]), i = 1 .. n)) fi;

Xl := map(floor, convert(xl, list)); Xr := map(ceil, convert(xr, list));

X := [seq([$ Xl[i] .. Xr[i]], i = 1 .. n)];

T := cartprod(X); S := table();

for k while not T[finished] do

t := T[nextvalue]();

if convert(eval(SN, zip(`=`, Var, t)), `and`) then

S[k] := t fi od;

convert(S, set);

end proc:

In the first example, we find all the integer points in the four-dimensional ball of radius 10:

Ball := IntegerPoints2({x1^2+x2^2+x3^2+x4^2 < 10^2}, [x1, x2, x3, x4]):  # All the integer points

nops(Ball);  # The total number of the integer points

seq(Ball[1000*n], n = 1 .. 10);  # Some points

                                                                    48945

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

                                [-5, 4, -6, -2], [-4, -5, 1, 5], [-4, -1, 6, 1], [-4, 3, 5, 6]

In the second example, with the visualization we find all the integer points in the inside intersection of  a cone and a cylinder:

A := <1, 0, 0; 0, (1/2)*sqrt(3), -1/2; 0, 1/2, (1/2)*sqrt(3)>:  # Matrix of rotation around x-axis at Pi/6 radians

f := unapply(A^(-1) . <x, y, z-4>, x, y, z):  

S0 := {4*x^2+4*y^2 < z^2}:  # The inner of the cone

S1 := {x^2+z^2 < 4}:  # The inner of the cylinder

S2 := evalf(eval(S1, {x = f(x, y, z)[1], y = f(x, y, z)[2], z = f(x, y, z)[3]})):

S := IntegerPoints2(`union`(S0, S2), [x, y, z]);  # The integer points inside of the intersection of the cone and the rotated cylinder

Points := plots[pointplot3d](S, color = red, symbol = solidsphere, symbolsize = 8):

Sp := plot3d([r*cos(phi), r*sin(phi), 2*r], phi = 0 .. 2*Pi, r = 0 .. 5, style = surface, color = "LightBlue", transparency = 0.7):

F := plottools[transform]((x, y, z)->convert(A . <x, y, z>+<0, 0, 4>, list)):

S11 := plot3d([2*cos(t), y, 2*sin(t)], t = 0 .. 2*Pi, y = -4 .. 7, style = surface, color = "LightBlue", transparency = 0.7):

plots[display]([F(S11), Sp, Points], scaling = constrained, orientation = [25, 75], axes = normal);

In the third example, we are looking for the integer points in a non-convex area between two parabolas. Here we have to specify ourselves the ranges to enumeration (Optimization[Minimize] command fails for this example):

P := IntegerPoints2([y > (-x^2)*(1/2)+2, y < -x^2+8], [x, y], [-4 .. 4, -4 .. 8]);

A := plots[pointplot](P, color = red, symbol = solidcircle, symbolsize = 10):

B := plot([(-x^2)*(1/2)+2, -x^2+8], x = -4 .. 4, -5 .. 9, color = blue):

plots[display](A, B, scaling = constrained);

 IntegerPoints2.mw

This post is my attempt to answer the question from   here : how to find all integer points (all points with integer coordinates) in the intersection of two cubes. The following procedure  IntegerPoints  solves a more general problem: it finds all the integer points of a bounded polyhedral region of arbitrary dimension, defined by a system of linear inequalities and / or equations.

Required parameters of the procedure: SN is a set or a list of linear inequalities and/or equations with any number of variables, the Var is the list of variables. The procedure returns the set of all integer points, satisfying the conditions  SN .

Code of the procedure:

restart;

IntegerPoints := proc (SN::{list, set}, Var::list)

local SN1, sn, n, Sol, k, i, s, S, R;

uses PolyhedralSets, SolveTools[Inequality];

SN1 := convert(evalf(SN), fraction);

for sn in SN1 do

if type(sn, `<`) then SN1 := subs(sn = (`<=`(op(sn))), SN1)

end if; end do;

if IsBounded(PolyhedralSet(SN1)) = false then error "The region should be bounded" end if;

n := nops(Var);

Sol := LinearMultivariateSystem(SN, Var);

if Sol = {} then return {} else

k := 0;

for s in Sol do if nops(indets(s[1])) = 1 then

S[0] := [[]];

for i to n do

S[i] := [seq(seq([op(j1), op(j2)], j2 = [isolve(eval(s[i], j1))]), j1 = S[i-1])] end do;

k := k+1; R[k] := op(S[n]);

end if; end do;

convert(R, set);

map(t->rhs~(t), %);

end if;

end proc:

Examples of use:

IntegerPoints({x > 0, y > 0, z > 0, 2*x+3*y+z < 12}, [x, y, z]);

  {[1, 1, 1], [1, 1, 2], [1, 1, 3], [1, 1, 4], [1, 1, 5], [1, 1, 6], [1, 2, 1], [1, 2, 2], [1, 2, 3], [2, 1, 1], [2, 1, 2],

                                   [2, 1, 3], [2, 1, 4], [2, 2, 1], [3, 1, 1], [3, 1, 2]}

IntegerPoints({x > 0, y > 0, z > 0, 2*x+3*y+z = 12}, [x, y, z]);

                                    {[1, 1, 7], [1, 2, 4], [1, 3, 1], [2, 1, 5], [2, 2, 2], [3, 1, 3], [4, 1, 1]}

IntegerPoints([x > 0, y > 0, z > 0, 2*x+3*y+z = 12, x+y+z <= 6], [x, y, z]);

                                                           {[1, 3, 1], [2, 2, 2], [4, 1, 1]}

isolve({x > 0, y > 0, z > 0, 2*x+3*y+z < 12});  #  isolve fails with these examples

              Warning, solutions may have been lost

isolve({x > 0, y > 0, z > 0, 2*x+3*y+z = 12});

              Warning, solutions may have been lost

In the following example (with a visualization) we find all integer point in the intersection of a square and a triangle:

S1 := {x > 0, y > 0, x < 13/2, y < 13/2}:

S2 := {y > (1/4)*x+1, y < 2*x, y+x < 12}:

S := IntegerPoints(`union`(S1, S2), [x, y]):

Region := plots[inequal](`union`(S1, S2), x = 0 .. 7, y = 0 .. 7, color = "LightGreen", nolines):

Points := plot([op(S)], style = point, color = red, symbol = solidcircle):

Square := plottools[curve]([[0, 0], [13/2, 0], [13/2, 13/2], [0, 13/2], [0, 0]], color = blue, thickness = 3):

Triangle := plottools[curve]([[4/7, 8/7], [4, 8], [44/5, 16/5], [4/7, 8/7]], color = blue, thickness = 3):

plots[display](Square, Triangle, Points, Region, scaling = constrained);

In the following example (with a visualization) we find all integer point in the intersection of two cubes. The second cube is obtained from the first cube by rotation with orthogonal matrix  A  and by a translation:

A := <1/3, 2/3, 2/3; -2/3, 2/3, -1/3; -2/3, -1/3, 2/3>:

f := unapply(A^(-1).<x+5, y-4, z-7>, x, y, z):

S1 := {x > 0, y > 0, z > 0, x < 6, y < 6, z < 6}:

S2 := eval(S1, {x = f(x, y, z)[1], y = f(x, y, z)[2], z = f(x, y, z)[3]}):

S := IntegerPoints(`union`(S1, S2), [x, y, z]);

Points := plots[pointplot3d](S, color = red, symbol = box):

Cube := plottools[cuboid]([0, 0, 0], [6, 6, 6], color = blue, linestyle = solid):

F := plottools[transform]((x, y, z)->convert(A.<x, y, z>+<-5, 4, 7>, list)):

plots[display](Cube,  F(Cube), Points, scaling = constrained, linestyle = solid, transparency = 0.7, orientation = [25, 75], axes = normal);

In the example below, all the ways to exchange $ 1 coins of 1, 5, 10, 25 and 50 cents, if the number of coins no more than 8, there is no pennies and there is at least one 50-cent coin:

IntegerPoints({x1 = 0, x2 >= 0, x3 >= 0, x4 >= 0, x5 >= 1,  x1+5*x2+10*x3+25*x4+50*x5 = 100, x1+x2+x3+x4+x5 <= 8}, [x1, x2, x3, x4, x5]);

nops(%);

                              {[0, 0, 0, 0, 2], [0, 0, 0, 2, 1], [0, 0, 5, 0, 1], [0, 1, 2, 1, 1], [0, 2, 4, 0, 1],

                                                 [0, 3, 1, 1, 1], [0, 4, 3, 0, 1], [0, 5, 0, 1, 1]}

                                                                                    8

Integer_points.mw

Addition: Below in my comments another procedure  IntegerPoints1  is presented that solves the same problem.

Some Maple 18 short (and I believe elegant) code for doing gravitational simulations with N bodies in space:

N_body_problem.mw

Initial velocities have been tweaked to keep the system stable for the duration of the animation.

Please feel free to fiddle with its parameters, velocities and positions and/or N itself, to produce more interesting animations or re-use the code therein (You can safely ignore the (c), it's there just for archiving purposes).

The following are animations from three runs with N=4, N=3 and N=2, no other parameters changed.

There has been a spate of Questions posted in the past week about computing eigenvalues. Invariably, the Questioners have computed some eigenvalues by applying fsolve to a characteristic polynomial obtained from a floating-point matrix via LinearAlgebra:-Determinant. They are then surprised when various tests show that these eigenvalues are not correct. In the following worksheet, I show that the eigenvalues computed by the fsolve@Determinant method (when applied to a floating-point matrix) are 100% garbage for dense matrices larger than about Digits x Digits. The reason for this is that computing the determinant introduces too much round-off error into the coefficients of the characteristic polynomial. The best way to compute the eigenvalues is to use LinearAlgebra:-Eigenvalues or LinearAlgebra:-Eigenvectors. Furthermore, very accurate results can be obtained without increasing Digits.

Download Eigenvalues.mw

A new Maple e-book, Multivariate Calculus Study Guide, is now available. Part of the Clickable Calculus collection of interactive Maple e-books, this guide takes full advantage of Maple’s Clickable Math approach. It has over 600 worked examples, the vast majority of which are solved using interactive, Clickable Math techniques. 

Deisgned to help students taking this course, instructors may also find this e-book useful as a guide to using Clickable Math to teach Multivariate Calculus.

See Multivariate Calculus Study Guide for more information.

eithne

A new Maple e-book, Multivariate Calculus Study Guide, is now available. Part of the Clickable Calculus collection of interactive Maple e-books, this guide takes full advantage of Maple’s Clickable Math approach. It has over 600 worked examples, the vast majority of which are solved using interactive, Clickable Math techniques. 

Deisgned to help students taking this course, instructors may also find this e-book useful as a guide to using Clickable Math to teach Multivariate Calculus.

See Multivariate Calculus Study Guide for more information.

eithne

Consider the well-known Euler's formula  

 e^ix = cos x + i sin x   

When we calculate that for  x = π  we get:

e^iπ = cos π + i sin π   or

e^iπ = −1 + i × 0   (because cos π = −1 and sin π = 0)  or

e^iπ = −1  or  e^iπ + 1 = 0

It seems absolutely magical that such a neat equation combines  5  fundamental constants: e ,  i ,  π , 1 , 0

The purpose of this post - to give a simple visualization of equality  e^iπ = −1  (statical and animated) by expanding  e^iπ  in a series of complex numbers. These numbers we represent as vectors in the plane. We will see that the partial sums of this series are broken lines like a spiral, twisting around the point -1 steadily approaching to it.

Euler procedure has one required parameter  n is positive integer - the number of displayed terms of the series  for  e^iπ  

Optional parameter  a  is any symbol (by default  a=NULL). We use this option if  instead of a static spiral want to see an animated spiral. 

Procedure code can be found in the attached file  Euler.mw

Examples of use.

The first example shows  8 terms of the series (broken line of 8 units):

Euler(8);

The terms of the series where  n> = 10  on the same plot can not be seen as very small. In this case, we use  the second plot with magnification of  100 : 1 .  

The second example:

Euler(14);

In the third example, we see an animated broken line. It's  first 9 units represented  on the left plot, and then for n> = 10 on the right plot:

Euler(13, a);

Euler.mw

In this paper we will demonstrate the importance of using simple to complex algorithms applied to complex systems in civil and mechanical engineering. In order to develop solutions that developers need to be involved in issues of advanced dynamic computer science. We show how is that with the Maple scientific program and through component-based algorithms can generate power then then be inserted into specific algorithms. Will form patterns with movements of rotation and revolution of their axes, in each case to model and analyze the curves thereof comprising. With these modelalos and curve analysis we can predict manufacturing costs, freight, inter alia estrcturas which they can be used with the correct use of Maplesoft.

IX_Fast_2016.pdf

Solid_Algorithms_applied_in_complex_3D_structures_for_Civil_Engineering_with_Maplesoft.mw

(in spanish)

Lenin Araujo Castillo

This post is my attempt to answer the question from here .  

The procedure  ContoursWithLabels  has 2 required parameters: Expr  is an expression in  x  and  y  variables,  Range1  and  Range2  are ranges for  x  and  y . In this case, the output is the list of floats for the contours and 8 black contours (with labels) (the axis of coordinates as a box). 

The optional parameters: Number is positive integer - the number of contours (by default Number=8),  S is a set of real numbers  C  for contours (for which Expr=C) (by default  S={}),  GraphicOptions  is a list of graphic options for plotting (by default  GraphicOptions=[color = black, axes = box]),  Coloring  is an equality  Coloring=list of color options for  plots[dencityplot]  command (by default Coloring=NULL). 

The code of the procedure:

restart;

ContoursWithLabels := proc (Expr, Range1::(range(realcons)), Range2::(range(realcons)), Number::posint := 8, S::(set(realcons)) := {}, GraphicOptions::list := [color = black, axes = box], Coloring::`=` := NULL)

local r1, r2, L, f, L1, h, S1, P, P1, r, M, C, T, p, p1, m, n, A, B, E;

uses plots, plottools;

f := unapply(Expr, x, y);

if S = {} then r1 := rand(convert(Range1, float)); r2 := rand(convert(Range2, float));

L := [seq([r1(), r2()], i = 1 .. 205)];

L1 := convert(sort(select(a->type(a, realcons), [seq(f(op(t)), t = L)]), (a, b) ->is(abs(a) < abs(b))), set);

h := (L1[-6]-L1[1])/Number;

S1 := [seq(L1[1]+(1/2)*h+h*(n-1), n = 1 .. Number)] else

S1 := convert(S, list)  fi;

print(Contours = evalf[2](S1));

r := k->rand(20 .. k-20); M := []; T := [];

for C in S1 do

P := implicitplot(Expr = C, x = Range1, y = Range2, op(GraphicOptions), gridrefine = 3);

P1 := [getdata(P)];

for p in P1 do

p1 := convert(p[3], listlist); n := nops(p1);

if n < 500 then m := `if`(40 < n, (r(n))(), round((1/2)*n)); M := `if`(40 < n, [op(M), p1[1 .. m-11], p1[m+11 .. n]], [op(M), p1]); T := [op(T), [op(p1[m]), evalf[2](C)]] else

if 500 <= n then h := floor((1/2)*n); m := (r(h))(); M := [op(M), p1[1 .. m-11], p1[m+11 .. m+h-11], p1[m+h+11 .. n]]; T := [op(T), [op(p1[m]), evalf[2](C)], [op(p1[m+h]), evalf[2](C)]]

fi; fi; od; od;

A := plot(M, op(GraphicOptions));

B := plots:-textplot(T);

if Coloring = NULL then E := NULL else E := ([plots:-densityplot])(Expr, x = Range1, y = Range2, op(rhs(Coloring)))  fi;

display(E, A, B);

end proc:

Examples of use:

ContoursWithLabels(x^2+y^2, -3 .. 3, -3 .. 3);

ContoursWithLabels(x^2-y^2, -5 .. 5, -5 .. 5, {-20, -15, -10, -5, 0, 5, 10, 15, 20}, [color = black, thickness = 2, axes = box], Coloring = [colorstyle = HUE, colorscheme = ["White", "Red"], style = surface]);

The next example, I took from here:

ContoursWithLabels(sin(1.3*x)*cos(.9*y)+cos(.8*x)*sin(1.9*y)+cos(.2*x*y), -5 .. 0, 2 .. 5, {seq(-2 .. 2, 0.5)}, [color = black, axes = box], Coloring = [colorstyle = HUE, colorscheme = ["Cyan", "Red"], style = surface]);

There are many more examples can be found in the attached file. 

ContoursWithLabels1.mw

Edit. The attached file has been corrected.