15082 Reputation

24 Badges

12 years, 30 days

MaplePrimes Activity

These are Posts that have been published by Kitonum

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:


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.


Edit. The attached file has been corrected.

knight's tour is a sequence of moves of a knight on a chessboard such that the knight visits every square only once. This problem mentioned in  page of tasks still without  Maple implementation. 

The post presents the implementation of this task in Maple. Required parameter of the procedure (named  KnightTour)  is  address  - the address of the initial square in the algebraic notation. The second parameter  opt  is optional parameter:

1) if  opt is sequence  (by default) then  the procedure returns the sequence of moves of the knight in the usual algebraic notation,

2)  if  opt is diagram   then  the procedure returns the plot of moves of the knight and  sequentially numbers all the visited squares,

3) if  opt is animation  then  the procedure returns an animation of moves of the knight.

In the procedure is used a solution with maximum symmetry by George Jelliss,


Code of the procedure:

KnightTour := proc(address::symbol, opt::symbol := sequence)

local L, n, L1, k, i, j, squares, border, chessboard, letters, digits, L2, L3, Tour, T, F;

uses ListTools, plottools, plots;

L := [a1, b3, d2, c4, a5, b7, d8, e6, d4, b5, c7, a8, b6, c8, a7, c6, b8, a6, b4, d5, e3, d1, b2, a4, c5, d7, f8, h7, f6, g8, h6, f7, h8, g6, e7, f5, h4, g2, e1, d3, e5, g4, f2, h1, g3, f1, h2, f3, g1, h3, g5, e4, d6, e8, g7, h5, f4, e2, c1, a2, c3, b1, a3, c2];

n := Search(address, L);

L1 := [L[n .. 64][], L[1 .. n-1][]];

if opt = sequence then return L1[] fi;

k := 0;

for i to 8 do

for j from `if`(type(i, odd), 1, 2) by 2 to 8 do

k := k+1;

squares[k] := polygon([[i-1/2, j-1/2], [i-1/2, j+1/2], [i+1/2, j+1/2], [i+1/2, j-1/2]], color = grey);

od;  od;

squares := convert(squares, list);

border := curve([[1/2, 1/2], [1/2, 17/2], [17/2, 17/2], [17/2, 1/2], [1/2, 1/2]], color = black, thickness = 4);

chessboard := display(squares, border);

letters := [a, b, c, d, e, f, g, h];

digits := [$ 1 .. 8];

L2 := convert~(L1, string);

L3 := subs(letters=~digits, map(t->[parse(t[1]), parse(t[2])], L2));

Tour := curve(L3, color = red, thickness = 3);

T := textplot([seq([op(L3[i]), i], i = 1 .. 64)], align = above, font = [times, bold, 16]);

if opt = diagram then return display(chessboard, Tour, T, axes = none) fi;

F := seq(display(chessboard, curve(L3[1 .. s], color = red, thickness = 3), textplot([seq([op(L3[i]), i], i = 1 .. s)], align = above, font = [times, bold, 16])), s = 1 .. 64);

display(seq(F[i]$5, i = 1 .. 64), insequence = true, axes = none);

end proc:


 Examples of use:


KnightTour(f3, diagram);



KnightTour(f3, animation);


@Markiyan Hirnyk   I try not to use this package, as I think the results are not reliable enough. Here is the example, where instead of the three real roots it finds only one, despite the hint to look for the three roots:


DirectSearch:-SolveEquations(100^x=x^100, AllSolutions, solutions=3);


There are many other examples, particularly in discrete optimization in which it returns false results. Here is one example (well-known to you).

This post is related to the this thread

The recursive procedure  PosIntSolve  finds the number of non-negative or positive solutions of any linear Diophantine equation

a1*x1+a2*x2+ ... +aN*xN = n  with positive coefficients a1, a2, ... aN .  

Formal parameters: L is the list of coefficients of the left part, n  is the right part,  s (optional) is nonneg (by default) for nonnegint solutions and  pos  for positive solutions.

The basic ideas:

1) If you make shifts of the all unknowns by the formulas  x1'=x1-1,  x2'=x2-1, ... , xN'=xN-1  then  the number of positive solutions of the first equation equals the number of non-negative solutions of the second equation.

2) The recurrence formula (penultimate line of the procedure) can easily be proved by induction.


The code of the procedure:


PosIntSolve:=proc(L::list(posint), n::nonnegint, s::symbol:=nonneg)

local N, n0;

option remember;

if s=pos then n0:=n-`+`(op(L)) else n0:=n fi;


if N=1 then if irem(n0,L[1])=0 then return 1 else return 0 fi; fi;

add(PosIntSolve(subsop(1=NULL,L),n0-k*L[1]), k=0..floor(n0/L[1]));

end proc:


Examples of use.


Finding of the all positive solutions of equation 30*a+75*b+110*c+85*d+255*e+160*f+15*g+12*h+120*i=8000:


PosIntSolve([30,75,110,85,255,160,15,12,120], 8000, pos);





To test the procedure, solve (separately for non-negative and positive solutions) the simple equation  2*x1+7*x2+3*x3=2000  in two ways (by the  procedure and brute force method):


PosIntSolve([2,7,3], 2000);

PosIntSolve([2,7,3], 2000, pos);








for x from 0 to 2000/2 do

for y from 0 to floor((2000-2*x)/7) do

for z from 0 to floor((2000-2*x-7*y)/3) do

if 2*x+7*y+3*z=2000 then k:=k+1 fi;

od: od: od:



for x from 1 to 2000/2 do

for y from 1 to floor((2000-2*x)/7) do

for z from 1 to floor((2000-2*x-7*y)/3) do

if 2*x+7*y+3*z=2000 then k:=k+1 fi;

od: od: od:







Another example - the solution of the famous problem: how many ways can be exchanged $ 1 using the coins of smaller denomination.




 Edit.  The code has been slightly edited 

This post is related to the question. There were  proposed two ways of finding the volume of the cutted part of a sphere in the form of a wedge.  Here the procedure is presented that shows the rest of the sphere. Parameters procedure: R - radius of the sphere, H1 - the distance the first cutting plane to the plane  xOy,  H2 -  the distance the second cutting plane to the plane  zOy. Necessary conditions:  R>0,  H1>=0,  H2>=0,  H1^2+H2^2<R^2 . For clarity, different surfaces are painted in different colors.


Pic := proc (R::positive, H1::nonnegative, H2::nonnegative)

local A, B, C, E, F;

if R^2 <= H1^2+H2^2 then error "Should be H1^(2)+H2^(2)<R^(2)" end if;

A := plot3d([R*sin(theta)*cos(phi), R*sin(theta)*sin(phi), R*cos(theta)], phi = arctan(sqrt(-H1^2-H2^2+R^2), H2) .. 2*Pi-arctan(sqrt(-H1^2-H2^2+R^2), H2), theta = 0 .. Pi, color = green);

B := plot3d([R*sin(theta)*cos(phi), R*sin(theta)*sin(phi), R*cos(theta)], phi = -arctan(sqrt(-H1^2-H2^2+R^2), H2) .. arctan(sqrt(-H1^2-H2^2+R^2), H2), theta = 0 .. arccos(sqrt(R^2-H2^2-H2^2*tan(phi)^2)/R), color = green);

C := plot3d([R*sin(theta)*cos(phi), R*sin(theta)*sin(phi), R*cos(theta)], phi = -arctan(sqrt(-H1^2-H2^2+R^2), H2) .. arctan(sqrt(-H1^2-H2^2+R^2), H2), theta = arccos(H1/R) .. Pi, color = green);

E := plot3d([r*cos(phi), r*sin(phi), H1], phi = -arccos(H2/sqrt(R^2-H1^2)) .. arccos(H2/sqrt(R^2-H1^2)), r = H2/cos(phi) .. sqrt(R^2-H1^2), color = blue);

F := plot3d([H2, r*cos(phi), r*sin(phi)], phi = arccos(sqrt(-H1^2-H2^2+R^2)/sqrt(R^2-H2^2)) .. Pi-arccos(sqrt(-H1^2-H2^2+R^2)/sqrt(R^2-H2^2)), r = H1/sin(phi) .. sqrt(R^2-H2^2), color = gold);

plots[display](A, B, C, E, F, axes = none, view = [-1.5 .. 1.5, -1.5 .. 1.5, -1.5 .. 1.5], scaling = constrained, lightmodel = light4, orientation = [60, 80]);

end proc:


Example of use:

Pic(1,  0.5,  0.3);




3 4 5 6 7 8 9 Page 5 of 11