Slides of the presentation at the VII Workshop Fast Computational and Applied Mathematics developed in graduate school at the National University of Trujillo. January 8, 2014.

Visualización_Geomét.pdf

L. Araujo C.

An isometry of the Euclidean plane is a distance-preserving transformation of the plane.  There are four types: translations,  rotations,  reflections,  and  glide reflections. See http://en.wikipedia.org/wiki/Euclidean_plane_isometry#Classification_of_Euclidean_plane_isometries

Procedure  AreIsometric  checks two plane sets  Set1  and  Set2  and  if there is an isometry of the plane mapping the first set to the second one, then the procedure returns true and false otherwise. Global variable  T  saves the type of isometry and all its parameters. For example, it returns  the rotation center and rotation angle, etc.

Each of the sets  Set1  and  Set2   is the set (or list) consisting of the following objects in any combinations:

1) The points that are defined as a list of coordinates  [х, y] .

2) Segments, which are defined as a set (or list)  of two points  {[x1, y1], [x2, y2]}  or  [[x1, y1], [x2, y2]]  .

3) Curves, which should  be defined as a list of points   [[x1, y1], [x2, y2], ..., [xn, yn]].

Of course, if  n = 2, then the curve is identical to the segment.

Code of the procedure:

AreIsometric:=proc(Set1::{set,list}, Set2::{set,list}) 

local n1, n2, n3, n4,s1, S, s, l1, l2, S11, f, x0, y0, phi, Sol, x, y, M1, M2, A1, A2, A3, A4, B1, B2, B3, B4, line1, line2, line3, line4, u, v, Sign, g, M, Line1, Line2, Line3, A, B, C, h, AB, CD, Eq, Eq1, T1, T2, i, S1, S2, T11; 

global T; 

uses combinat;     

S1:={};  S2:={};  T1:={}; T2:={}; 

for i in Set1 do 

if i[1]::realcons  then S1:={op(S1),i} else 

S1:={op(i), op(S1)};  T1:={op(T1), seq({i[k],i[k+1]}, k=1..nops(i)-1)} fi;  

od; 

for i in Set2 do 

if i[1]::realcons  then S2:={op(S2),i} else 

S2:={op(i), op(S2)};  T2:={op(T2), seq({i[k],i[k+1]}, k=1..nops(i)-1)} fi; 

od;

n1:=nops(S1);  n2:=nops(S2);  n3:=nops(T1); n4:=nops(T2); 

if is(S1=S2) and is(T1=T2) then T:=identity;  return true fi; 

if n1<>n2 or n3<>n4 then return false fi; 

if n1=1 then T:=[translation, <S2[1,1]-S1[1,1], S2[1,2]-S1[1,2]>];  return true fi;

f:=(x,y,phi)->[(x-x0)*cos(phi)-(y-y0)*sin(phi)+x0, (x-x0)*sin(phi)+(y-y0)*cos(phi)+y0];  g:=(x,y)->[(B^2*x-A^2*x-2*A*B*y-2*A*C)/(A^2+B^2), (A^2*y-B^2*y-2*A*B*x-2*B*C)/(A^2+B^2)]; 

_Envsignum0 := 1;

s1:=[S1[1], S1[2]];  S:=select(s->is((s1[2,1]-s1[1,1])^2+(s1[2,2]-s1[1,2])^2=(s[2,1]-s[1,1])^2+(s[2,2]-s[1,2])^2),permute(S2, 2));    

for s in S do   

# Checking for translation    

l1:=s[1]-s1[1]; l2:=s[2]-s1[2]; 

if is(l1=l2) then S11:=map(x->x+l1, S1); 

if n3<>0 then T11:={seq(map(x->x+l1, T1[i]), i=1..nops(T1))}; fi; 

if n3=0 then  if is(S11=S2) then T:=[translation, convert(l1, Vector)]; return true fi;  else 

if is(S11=S2) and is(T11=T2) then T:=[translation, convert(l1, Vector)]; return true fi; fi; 

fi;   

# Checking for rotation   

x0:='x0'; y0:='y0'; phi:='phi'; u:='u'; v:='v'; Sign:='Sign';    

if  is(s1[1]-s[1]<>s1[2]-s[2]) then  

M1:=[(s1[1,1]+s[1,1])/2, (s1[1,2]+s[1,2])/2]; M2:=[(s1[2,1]+s[2,1])/2, (s1[2,2]+s[2,2])/2]; A1:=s1[1,1]-s[1,1]; B1:=s1[1,2]-s[1,2]; A2:=s1[2,1]-s[2,1]; B2:=s1[2,2]-s[2,2];    line1:=A1*(x-M1[1])+B1*(y-M1[2])=0; line2:=A2*(x-M2[1])+B2*(y-M2[2])=0;  

if is(A1*B2-A2*B1<>0) then Sol:=solve({line1, line2}); x0:=simplify(rhs(Sol[1]));   y0:=simplify(rhs(Sol[2])); u:=[s1[1,1]-x0,s1[1,2]-y0]; v:=[s[1,1]-x0,s[1,2]-y0];    else   

if is(s[2]-s1[1]=s[1]-s1[2])  then   x0:=(s1[1,1]+s[1,1])/2;  y0:=(s1[1,2]+s[1,2])/2; 

if is([x0,y0]<>s1[1]) then  u:=[s1[1,1]-x0,s1[1,2]-y0]; v:=[s[1,1]-x0,s[1,2]-y0]; else 

u:=[s1[2,1]-x0,s1[2,2]-y0]; v:=[s[2,1]-x0,s[2,2]-y0]; fi;

else  A3:=s1[2,1]-s1[1,1];  B3:=s1[2,2]-s1[1,2]; A4:=s[2,1]-s[1,1];  B4:=s[2,2]-s[1,2];  line3:=B3*(x-s1[1,1])-A3*(y-s1[1,2])=0;  line4:=B4*(x-s[1,1])-A4*(y-s[1,2])=0;Sol:=solve({line3, line4}); x0:=simplify(rhs(Sol[1])); y0:=simplify(rhs(Sol[2]));   

if is(s1[1]=s[1]) then    u:=s1[2]-[x0,y0]; v:=s[2]-[x0,y0]; else   

u:=s1[1]-[x0,y0]; v:=s[1]-[x0,y0];  fi;  fi;  fi;   

Sign:=signum(u[1]*v[2]-u[2]*v[1]);   phi:=Sign*arccos(expand(rationalize(simplify((u[1]*v[1]+u[2]*v[2])/sqrt(u[1]^2+u[2]^2)/sqrt(v[1]^2+v[2]^2)))));      S11:=expand(rationalize(simplify(map(x->f(op(x), phi), S1))));   

if n3<>0 then T11:={seq(expand(rationalize(simplify(map(x->f(op(x), phi), T1[i])))), i=1..nops(T1))}; fi; 

if n3=0 then  if is(S11=expand(rationalize(simplify(S2))))  then T:=[rotation, [x0,y0], phi]; return true fi;  else 

if is(S11=expand(rationalize(simplify(S2)))) and  is(T11=expand(rationalize(simplify(T2)))) then  

T:=[rotation, [x0,y0], phi]; return true fi;  fi; 

fi; 

od;   

# Checking for reflection or glide reflection   

for s in S do    

AB:=s1[2]-s1[1]; CD:=s[2]-s[1];  

if is(AB[1]*CD[2]-AB[2]*CD[1]=0) then  M:=(s1[2]+s[1])/2;

if  is(AB[1]*CD[1]+ AB[2]*CD[2]>0) then  A:=AB[2]; B:=-AB[1];    Line1:=A*(x-M[1])+B*(y-M[2])=0;  else 

A:=AB[1]; B:=AB[2];  Line2:=A*(x-M[1])+B*(y-M[2])=0; fi;  

else     u:=[AB[1]+CD[1], AB[2]+CD[2]];  A:=u[2]; B:=-u[1];     M:=[(s1[1,1]+s[1,1])/2, (s1[1,2]+s[1,2])/2]; Line3:=A*(x-M[1])+B*(y-M[2])=0;   fi;    C:=-A*M[1]-B*M[2];  h:= simplify(expand(rationalize(s[1]-g(op(s1[1])))));    S11:=expand(rationalize(simplify(map(x->g(op(x))+h, S1))));  

if n3<>0 then T11:={seq(expand(rationalize(simplify(map(x->g(op(x))+h, T1[i])))), i=1..nops(T1))}; fi;    

if n3=0 then   if is(S11=expand(rationalize(S2))) then 

Eq:=A*x+B*y+C=0; Eq1:=`if`(is(coeff(lhs(Eq), y)<>0), y=solve(Eq, y),  x=solve(Eq, x)); 

if h=[0,0] then  T:=[reflection, Eq1] else T:=[glide_reflection,Eq1,convert(h, Vector)] fi; return true fi; else  

if is(S11=expand(rationalize(S2))) and is(T11=expand(rationalize(T2))) then 

Eq:=A*x+B*y+C=0; Eq1:=`if`(is(coeff(lhs(Eq), y)<>0), y=solve(Eq, y),  x=solve(Eq, x)); 

if h=[0,0] then T:=[reflection, Eq1] else

T:=[glide_reflection,Eq1,convert(h, Vector)] fi; return true fi;  fi;   

od;      

T:='T';   false;  

end proc:

Three simple examples:

AreIsometric({[4, 0], [7, 4], [14, 0]}, {[4, 14], [9, 14], [10, 6]});  T;

AreIsometric({[2, 0], [2, 2], [5, 0]}, {[3, 3], [3, 6], [5, 3]});  T;

S1 := {[[5, 5], [5, 20], [10, 15], [15, 20], [15, 5]]}: 
S2 := {[[21, 11], [30, 23], [31, 16], [38, 17], [29, 5]]}: 
S3 := {[[50, 23], [41, 11], [51, 16], [49, 5], [58, 17]]}: 
AreIsometric(S1, S2); T; AreIsometric(S1, S3);

plots[display](plottools[curve](op(S1), thickness = 2, color = green), plottools[curve](op(S2), thickness = 2, color = green), plottools[curve](op(S3), thickness = 2, color = red), scaling = constrained, view = [-1 .. 60, -1 .. 25]);  # Green sets are isometric,  green and red sets aren't

Example with animation:

S1:={[0,0],[-1,2],[2,4],[4,2]}:  S2:={[8,3],[6,6],[8,8],[10,4]}:

AreIsometric(S1, S2);  T;  

with(plots): with(plottools): 

#  For clarity, instead of points polygons depicted with vertices at these points 

N:=50:   

A:=seq(rotate(polygon([[0,0],[-1,2],[2,4],[4,2]], color=blue), (k*Pi)/(2*N), [3,7]), k=0..N):  B:=polygon([[8,3],[6,6],[8,8],[10,4]], color=green):  E:=line([3,7], [6,6], color=black, linestyle=2): 

C:=seq(rotate(line([3,7], [2,4], color=black, linestyle=2), (k*Pi)/(2*N), [3,7]), k=0..N):  L:=curve([[3,7],[2,4], [-1,2], [0,0],[4,2], [2,4]], color=black, linestyle=2):  T:=textplot([3, 7.2, "Center of rotation"]): 

Frames:=seq(display(A[k], B, E,T,L, C[k]), k=1..N+1):   

display(seq(Frames[k],k=1..N+1), insequence=true, scaling=constrained);

Finding unique solutions to the problem of Queens (m chess queens on an n×n chessboard not attacking one another). Used the procedures  Queens  and  QueenPic  . See  http://www.mapleprimes.com/posts/200049-Queens-Problem-And-Its-Visualization-With-Maple

Queens(8, 8);  M := [ListTools[Categorize](AreIsometric, S)]:

nops(M);

seq(op(M[k])[1], k = 1 .. %);   # 12 unique solutions from total 92 solutions

QueensPic([%], 4);  #  Visualization of obtained solutions

Finding unique solutions to the problem "Polygons of matches"  (all polygons with specified perimeter and area). See http://www.mapleprimes.com/posts/142884-Polygons-Of-The-Matches

N := 12: S := 6: Polygons(N, S);

M := [ListTools[Categorize](AreIsometric, T)]:

n := nops(M);

seq([op(M[k][1])], k = 1 .. n);  #  7 unique solutions from total 35 solutions with perimeter 12 and area 6

Visual([%], 4);  #  Visualization of obtained solutions

 Isometry_of_plane_se.mw

I work entitled Point Exeter made ​​for Fast VII workshop on applied and computational mathematics 2014 Trujillo Peru.

  Punto_de_Exeter.mw   (version in spanish)

Atte.

Lenin Araujo Castillo

Physics Pure

Computer Science

It is a relatively recent innovation that complex-number computations can be done in the evalhf environment. When combined with plots:-densityplot, this makes escape-time fractals in the complex domain very easy to plot. This fractal is based on the Collatz problem. This Wikipedia article has a high-resolution picture of this fractal. I've switched the real and imaginary axes and reversed the direction of the real axis purely for asthetic reasons.

Download Collatz_fractal.mw

The Stone-Weierstass theorem  in its simplest form asserts that every continuous function defined on a closed interval [a,b] can be uniformly approximated as closely as desired by a polynomial function. Let us consider a concrete function (say, arcsin(sqrt(x))) on a concrete interval (for example,[0,1]) and a concrete rate (for instance, 0.01). The question arises: what can be  the degree of an approximating polynomial?
Looking in the constructive proof of the Weierstrass theorem (for example, see
W. Rudin, Principles of mathematical analysis. Third Ed. McGraw-Hill Inc. New York-...-Toronto. 1976, pp. 159-160 SWT.docx), we find the inequality for degree n in terms of the modulus of the  continuity delta and the maximum of the modulus M of a function f on [0,1]: 4*M*sqrt(n)*(1-delta^2)^n < epsilon/2.
Next, we find the modulus of the continuity of arcsin(sqrt(x)) with help of Maple (namely, the DirectSearch package):
>restart;
>CM := proc (delta) DirectSearch:-Search(abs(arcsin((x+delta)^(1/2))-arcsin(x^(1/2))),
 {0 <= x, 0 <= x+delta, x <= 1, x+delta <= 1}, maximize)
end proc
. Now delta is fitting to satisfy CM(delta) < 0.01:
>Digits := 15: CM(0.9999640e-4);

[0.999995686126010e-2, [x = .999900003599999], 18].
At last, we find the required degree, taking into account M=Pi/2 for arcsin(sqrt(x)) on [0,1]:
>DirectSearch:-SolveEquations((4*Pi*(1/2))*sqrt(n)*(1-0.9999640e-4^2)^n = (1/2)*10^(-2), {n >= 10^9}, tolerances = 10^(-8));

[3.68635028417869*10^(-35), Vector(1, {(1) = -0.607153216591882e-17}),[n = 1.77870508105403*10^9], 74]
The obtained result is unexpected and impressive. However, this is only an estimate of the degree for the chosen construction. There are different ways to construct an approximating polynomial. For example, let us take the interpolating polynomial.
>with(CurveFitting): Digits := 200: P := PolynomialInterpolation([seq([(1/200)*j,
evalf(arcsin(sqrt((1/200)*j)), 180)], j = 0 .. 200)], x);

8.57260524574724504043891781488113281218267308627010084942700641\
2116721658995225354525109649870447266086431479184935898860221001\
6810627259201248204607733508370522655937863029427984169024474693\
605019813*10^(-24)*x^200+
3.4102471291087052576144785068387656673244314487588\
37173451046570851636655790486463697061695256004409457030\
661587523327337363549630285194598139656219506035056874382\
5412929520214254752642899246978334986*10^83*x^199+...
The whole long output of sort(P) can be seen in the attached file.
>DirectSearch:-Search(abs(arcsin(sqrt(x))-P), {x >= 0, x <= 1}, maximize, tolerances = 10^(-10));

[0.7028873160870935332477114389520278374486329450431055674880288416078\

033259753063018233397798614e-2, [x = .999760629733897552108099038488344\

76319065496787157065017717228830101\

791752323133523143936216508553686883680060439608736578363\

678796478147136266075441732651036025656505033942652374763794644368578081487], 22]
See SWtheorem.mw

The following limit does not return a value. Then the evalf gives a wrong answer.

The answer should be "undefined" or -infinity .. infinity.

limit(exp(n)/(-1)^n, n = infinity) assuming n::posint; evalf(%);

                       /exp(n)              \
                  limit|------, n = infinity|
                       |    n               |
                       \(-1)                /

                               0.

The same happens if you delete the assumption.

A similar problem occurs with

limit(sin(Pi/2+2*Pi*n), n = infinity) assuming n::posint;
                            -1 .. 1
without the assumption this would be appropriate.

A little late in the game, with no Maple 18 wishlist yet I'd like to get it started

granted all wishes for previous versions not yet applied should still, hopefully, be under consideration.  

1 - The abitlity to angle x-axis labels.

2 - Textplot option for text to rotate with graph as it is rotated.  Currently all textplot texts are held to the horizontal

Thanks to the community through Maplesoft Mapleprimes that could develop in Computational Mathematics Achievement Day at our institution.

libro de geometrias:

http://www.slideshare.net/kaizzerz/geometria-analitica-charles-h-lehmann

Greetings to all.

This past year I have on occasion shared mathematical adventures with cycle index computations and Maple, e.g. at these links:

• Fascinating cycle indices with Maple
• The enumeration of non-isomorphic graphs
• Automorphisms of the Petersen Graph.

Befitting the season I am sending another post to continue this series of cycle index computations. I present two Maple implementations of Power Group Enumeration as described by Harary and Palmer in their book "Graphical Enumeration" and by Fripertinger in his paper "Enumeration in Musical Theory." It was a real joy working with Maple to implement the computational aspects of their work, i.e. the Power Group Enumeration Theorem. Moreover the resulting software is easy to read, simple and powerful and has a straightforward interface, taking advantage of many different capabilities present in Maple.

The problem I am treating is readily described. Consider a cube in 3 space and its symmetries under rotation, i.e. rigid motions. We ask in how many different ways we may color the edges of the cube with at most N colors where all colors are completely interchangable, i.e. have the symmetric group acting on them in addition to the edge permutation group of the cube. At the following Math Stackexchange Link  I have posted the Maple code to implement the algorithms / formulas of Harary / Palmer / Fripertinger to solve this problem. The reader is invited to study and test these algorithms. It seems to me an excellent instance of computational combinatorics fun.

To conclude I would like to point out that these algorithms might be candidates for a Polya Enumeration Theorem (PET) package that I have been suggesting for a future Maple release at the above posts, the algorithms being of remarkable simplicity while at the same time providing surprisingly sophisticated combinatorics and enumeration methods.

Season's greetings!

Marko Riedel

I find this very strange. I see questions posted at Maple prime using a language other than English.

I do not know if this is a Mapleprime bug, or not. If this is not a bug, why is it allowed to post questions in a language other than English?
 

It would be nice if Maple had a procedure which could turn a procedurelist into a listprocedure.
I use these words in the sense they are used in dsolve/numeric.
Thus by a procedurelist I mean a procedure returning lists (of numbers).
By a listprocedure I mean a list of procedures all having the same formal parameters and each returning one number.
Thus a `convert/listprocedure`should accept a procedurelist p as input and give as output the corresponding listprocedure [p1,p2, ... pN] where N is the number of elements in the output from p.

I tried making a `convert/listprocedure` myself and in doing so found that it was not totally trivial.
I had lots of problems but did end up with something that seems to work.

convert-listprocedur.mw

But my main point is that Maple ought to have some such facility either as described available to the user or by changing fsolve, complexplot or what have you, so that procedurelists are accepted.

Hi

It's been 3+ months since we launched this new, experimental, Maple Physics: Research & Development webpage, containing fixes and new developments around the clock made available to everybody. Today we are extending this experience to Differential Equations and Mathematical functions, launching the Maple Differential Equations and Mathematical Functions: Research & Development Maplesoft webpage. Hey!

With these pages we intend to move the focus of developments directly into the topics people are actually working on. The experience so far has been really good, putting our development at high RPM, an exciting roller-coast of exchange and activity.

As with the Research version of Physics, when suggestions about DEs or Mathematical Functions are implemented or issues are fixed, typically within a couple of days when that is possible, the changes will be made available to everybody directly in this new Maplesoft webpage. One word of clarification: for now, these updates will not include numerical ODE or numerical PDE solutions nor their numerical plotting. Sorry guys. One step at a time :)

This first update today concerns Differential Equations: dsolve and pdsolve can now handle linear systems of equations also when entered in Vector notation (Matrices and Vectors), related to a post in Mapleprimes from October/29. Attached is a demo illustrating the idea.

Everybody is welcome to bring suggestions and post issues. You can do that directly in Mapleprimes or writing to physics@maplesoft.com. While Differential Equations and Mathematical Functions are two areas where the Maple system is currently more mature than in Physics, these two areas cover so many subjects, including that there are the Research and the Education perspectives, that the number of possible topics is immense. 

DEsAndMathematicalFu.pdf   DEsAndMathematicalFu.mw

Edgardo S. Cheb-Terrab
Physics, DEs and Mathematical Functions, Maplesoft

Thу post is the answer to  http://www.mapleprimes.com/questions/200355-Clever-Cutter-3?reply=reply

This message contains two procedures . The first procedure called  LargestIncircle  symbolically looks for the largest circle inscribed in a convex polygon . Polygon must be specified as a list of lists of its vertices . If the polygon does not have parallel sides , the problem has a unique solution , otherwise may be several solutions . It is easy to prove that there is always a maximum circle that touches three sides of the polygon . Procedure code based on this property of the optimal solution.

The second procedure called  MinCoveringCircle  symbolically seeking smallest circle that encloses a given set of points. The set of points can be anyone not necessarily being vertices of a convex polygon. Procedure code based on the following geometric properties of the optimal solution , which can easily be proved :

1) The minimum covering circle is unique.

2) The minimum covering circle of a set Set can be determined by at most three points in  Set  which lie on the boundary of the circle. If it is determined by only two points, then the line segment joining those two points must be a diameter of the minimum circle. If it is determined by three points, then the triangle consisting of those three points is not obtuse.

Code of the 1st procedure:

restart;

LargestIncircle:=proc(L::listlist)

local n, x, y, P, S, T, Eqs, IsInside, OC, Pt, E, R, t, Sol, P0, R1, Circle;

uses combinat, geometry;

n:=nops(L); _EnvHorizontalName := x; _EnvVerticalName := y;

P:=[seq(point(A||i, L[i]), i=1..n)];

if nops(convexhull(P))<n then error "L should be a convex polygon" fi;

S:={seq(line(l||i,[A||i, A||(i+1)]), i=1..n-1), line(l||n,[A||n, A||1])};

Eqs:=map(lhs,{seq(Equation(l||i), i=1..n)});

T:=choose({seq(l||i, i=1..n)}, 3);

IsInside:=proc(Point, OC, Eqs)

convert([seq(is(eval(Eqs[i], {x=Point[1], y=Point[2]})*eval(Eqs[i], {x=OC[1], y=OC[2]})>0), i=1..n)], `and`);

end proc;

OC:=[add(L[i,1], i=1..n)/n, add(L[i,2], i=1..n)/n]; 

R:=0; point(Pt,[x,y]);

   for t in T do

     solve([distance(Pt,t[1])=distance(Pt,t[2]),      distance(Pt,t[2])=distance(Pt,t[3])],[x,y]);

     Sol:=map(a->[rhs(a[1]), rhs(a[2])], %);

     P0:=op(select(b->IsInside(b, OC, Eqs), Sol)); if P0<>NULL then R1:=distance(point(E,P0),      t[1]);

     if convert([seq(is(R1<=distance(E, i)), i=S minus t)], `and`) and is(R1>R) then R:=R1;      Circle:=[P0, R] fi; fi;

   od; 

Circle; 

end proc:

Code the 2nd procedure:

MinCoveringCircle:=proc(Set::{set, list}(list))

local S, T, t, d, A, B, C, P, R, R1, O1, Coord, L;

uses combinat, geometry; 

if type(Set, list) then S:=convert(Set, set) else S:=Set fi;

T:=choose(S, 2);

   for t in T do

      d:=simplify(distance(point(A, t[1]), point(B, t[2]))); midpoint(C, A, B);

      if convert([seq(is(distance(point(P, p), C)<=d/2), p in S minus {t[1], t[2]})], `and`)   then return            [coordinates(C), d/2] fi;

   od;

T:=choose(S, 3);

R:=infinity;

   for t in T do

      point(A, t[1]); point(B, t[2]); point(C, t[3]);

      if is(distance(A, B)^2<=distance(B, C)^2+distance(A, C)^2 and        distance(A,C)^2<=distance(B,        C)^2+distance(A, B)^2 and distance(B, C)^2 <= distance(A, C)^2 + distance(A, B)^2) then

      circle(c, [A,B,C],'centername'=O1); R1:=radius(c); Coord:=coordinates(O1);

      if convert([seq(is(distance(point(P, p), O1)<=R1), p in S minus {t[1], t[2], t[3]})], `and`) and is(R1<R) then    R:=R1; L:=[Coord, R] fi; fi;

   od; 

L; 

end proc: 

Examples of using the first procedure:

L:=[[0,2],[1,4],[4,4],[5,1],[3,0]]:

C:=LargestIncircle(L);

A:=plottools[disk](op(C), color=green):

B:=plottools[polygon](L, style=line, thickness=2):

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

The procedure can be used to find the largest circle inscribed in a convex shape, bounded by the curves, if these curves are approximated by broken lines:

L:=[seq([k,k^2], k=-1..2, 0.1)]:

C:=LargestIncircle(L);

A:=plottools[disk](op(C), color=green):

B:=plottools[polygon](L, style=line, thickness=2):

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

Examples of using the second procedure:

L:=[[0,2],[1,4],[4,4],[5,1],[3,0]]:

C:=MinCoveringCircle(L);

A:=plottools[circle](op(C), color=red, thickness=2):

B:=plottools[polygon](L, color=cyan):

plots[display](A,B);

The smallest circle covering 30 random points:

> roll:=rand(1..10):

L:=[seq([roll(), roll()], i=1..30)]:

C:=MinCoveringCircle(L);

A:=plottools[circle](op(C), color=red, thickness=2):

B:=seq(plottools[disk](L[i],0.07,color=blue), i=1..nops(L)):

C1:=plottools[disk](C[1], 0.07, color=red):

T:=plots[textplot]([C[1,1]-0.2, C[1,2]-0.2, "C"], font=[TIMES,ROMAN,16]):

plots[display](A, B, C1, T);

The attached presentation is the second one of a sequence of three that we want to do on Quantum Mechanics using Computer Algebra. The first presentation was about the equation for a quantum system of identical particles, the Gross-Pitaevskii equation (GPE). This second presentation is about the spectrum of its solutions. The level is that of an advanced undergraduate QM course. The novelty is again in the way these problems can be tackled nowadays in a computer algebra worksheet with Physics.

QuantumMechanics.pdf     Download QuantumMechanics2.mw

Edgardo S. Cheb-Terrab
Physics, Maplesoft