Applications, Examples and Libraries

Share your work here

CoolProp is an open source C++ library of thermophysical properties for pure fluids, pseudo-pure fluids, and humid air. Ian Bell has recently developed a wrapper for Maple (get the wrapper and library at Github). Compiling CoolProp gives a library (a DLL on Windows) you can call in Maple via define_external().

I started exploring CoolProp a few days ago, and here's a few simple examples of what you can do

The saturation pressure (in kPa) of the refrigerant R134a at 253 K

The pressure (in kPa) of the refrigerant R22 that produces a two-phase mixture of quality 0.3 with an enthalpy of 300 kJ/kg

And since I'm a fan of engineering visualization, here's a refrigeration cycle on a P-h-T chart, generated in Maple with CoolProp.

Here's a Maple application that uses CoolProp to analyze a refrigeration cycle (together with a CoolProp DLL for 64-bit Windows).

Analysis_of_a_Refrig.zip

I'd like to encourage anyone with an interest in thermophysical modeling to download CoolProp and explore its functionality. It's certainly opened up a new field of applications for me.

Samir

Event?

January 2014: Pages of oldest Russian Maple Application Center have been opened 10000000 times.

http://webmath.exponenta.ru/

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

 

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.

 

Collatz:= proc(b,a)  #Axes switched
local z:= -a+b*I, k;  #real part negated
     for k to 31 while abs(Im(z)) < 1 do
          z:= (1+4*z-(1+2*z)*cos(Pi*z))/4
     end do;
     k #escape time
end proc:

#Test evalhf'ability:

evalhf(Collatz(0,1));

32.

plotsetup(
     jpeg, plotoutput= "C:/Users/Carl/desktop/Collatz.jpg",
     plotoptions="height= 1024, width= 1024, quality= 95"
);

 

CodeTools:-Usage(
     plots:-densityplot(
          Collatz,
          -1..1, # imaginary range
          -0.5..4.5, #negative of real range
          colorstyle= HUE, grid= [1024, 1024], style= patchnogrid,
          labels= [Im,-Re], labelfont= [TIMES, BOLD, 14],
          axes= boxed,
          caption= cat("      Happy New Year ",                  

                StringTools:-FormatTime("%Y")),
          captionfont= [HELVETICA, BOLDOBLIQUE, 18]
     )
);

memory used=24.08MiB, alloc change=24.00MiB, cpu time=7.78s, real time=7.79s

 

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.

Greetings to all.

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

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

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.

 

The Gross-Pitaevskii equation and Bogoliubov spectrum
  

Pascal Szriftgiser1 and Edgardo S. Cheb-Terrab2 

(1) Laboratoire PhLAM, UMR CNRS 8523, Université Lille 1, F-59655, France

(2) Maplesoft, Canada

 

Departing from the equation for a quantum system of identical boson particles, i.e.the Gross-Pitaevskii equation (GPE), the dispersion relation for plane-wave solutions are derived, as well as the Bogoliubov equations and dispersion relations for small perturbations `&delta;&varphi;` around the GPE stationary solutions.

Stationary and plane-wave solutions to the Gross-Pitaevskii equation

 


Problem: Given the Gross-Pitaevskii equation, NULL

I*`&hbar;`*psi[t] = (G*abs(psi)^2+V)*psi-`&hbar;`^2*%Laplacian(psi)/(2*m)

  

a) Derive a relationship between the chemical potential mu entering the phase of stationary, uniform solutions, the atom-atom interaction constant G and the particle density n = abs(psi)^2 in the absence of an external field (V = 0).

  

b) Derive the dispersion relation for plane-wave solutions as a function of G and n. 

  

 

Background: The Gross-Pitaevskii equation is particularly useful to describe Bose Einstein condensates (BEC) of cold atomic gases [3, 4, 5]. The temperature of these cold atomic gases is typically in the w100 nano-Kelvin range. The atom-atom interaction are repulsive for G > 0 and attractive for G < 0 , where G is the interaction constant. The GPE is also widely used in non-linear optics to model the propagation of light in optical fibers.

Solution

   

The Bogoliubov equations and dispersion relations

 

 

Problem: Given the Gross-Pitaevskii equation,

  

a) Derive the Bogoliubov equations, that is, equations for elementary excitations `&delta;&varphi;` and conjugate(`&delta;&varphi;`)around a GPE stationary solution `&varphi;`(x, y, z), NULL

 

"{[[i `&hbar;` (&PartialD;)/(&PartialD;t) `delta&varphi;`=-(`&hbar;`^2 (&nabla;)^2`delta&varphi;`)/(2 m)+(2 G |`&varphi;`|^2+V-mu) `delta&varphi;`+G `&varphi;`^2 (`delta&varphi;`),,],[i `&hbar;` (&PartialD;)/(&PartialD;t)( `delta&varphi;`)=+(`&hbar;`^2 (&nabla;)^2(`delta&varphi;`))/(2 m)-(2 G |`&varphi;`|^2+V-mu) (`delta&varphi;`)-G `delta&varphi;` ((`&varphi;`))^(2),,]]"

  


b) Show that the dispersion relations of these equations, known as the Bogoliubov spectrum, are given by

  

epsilon[k] = `&hbar;`*omega[k] and `&hbar;`*omega[k] = `&+-`(sqrt(`&hbar;`^4*k^4/(4*m^2)+`&hbar;`^2*k^2*G*n/m)),

  


where k is the wave number of the considered elementary excitation, epsilon[k] its energy or, equivalently, omega[k] its frequency.

Solution

   

NULL

References

NULL

[1] Gross-Pitaevskii equation (wiki)

[2] Continuity equation (wiki)
[3] Bose–Einstein condensate (wiki)

[4] Dispersion relations (wiki)

[5] Advances In Atomic Physics: An Overview, Claude Cohen-Tannoudji and David Guery-Odelin, World Scientific (2011), ISBN-10: 9812774963.

[6] Nonlinear Fiber Optics, Fifth Edition (Optics and Photonics), Govind Agrawal, Academic Press (2012), ISBN-13: 978-0123970237.

 

 

QuantumMechanics.pdf     Download QuantumMechanics2.mw

Edgardo S. Cheb-Terrab
Physics, Maplesoft

Hi,
Relevant developments in Physics happened during the last month and a 1/2, some of them of interest beyond the use of this package. Among the most exciting things I can mention:

  1. The redefinition of the derivative rule for the complex components (abs, argument, conjugate, Im, Re, signum) together with the introduction of Wirtinger calculus, as an user-option to work with complex variables. In other words: it is now possible to compute taking z and its conjugate as independent variables and in equal footing.
  2. Introduction of textbook mathematical display for all the inert functions of the mathematical language, also for unknown functions f(x).
  3. New options in Physics:-Setup to indicate that some mathematical objects are real (different from assume(x, real), while integrated with `is` and `coulditbe`).
  4. A rather large number of micro developments advancing the integration of Physics with simplify, expand and combine.
  5. Another large number of micro developments for quantum mechanics.
  6. New options in Physics:-Setup to redefine sum as Physics:-Library:-Add, and with that have access to multiindex summation directly from sum, for instance as in sum(f(i, j), i + j <= n), including related typesetting copy & paste.

As usual the latest version of the package is available for download in the Maplesoft Physics: Research & Development webpage  and in the zip there is a worksheet illustrating all these developments. Below I'm copying the section related to the new redefinesum option of Physics:-Setup and multiindex summation.

Thanks to everyone who provided feedback, it has been of great value and at the root of this new round of developments.

December 4

 
• 

New option in Setup: redefinesum, so that the sum command is redefined in such a way that
    a) the sum arguments are processed in a way avoiding premature evaluation and related unexpected results or error interruptions
    b) the sum command includes new functionality present in Physics:-Library:-Add to perform sum over integer values of many indices, as in

"(&sum;)S(i,j)"     or  "(&sum;)S(i,j)" 

restart; with(Physics); Setup(notation = true)

`* Partial match of  'notation' against keyword 'mathematicalnotation'`

 

[mathematicalnotation = true]

(1.1)

New option: redefine sum so that its arguments are processed by the more modern Physics:-Library:-Add and so that it can perform multiindice summation.

 

Example:

By default, the sum command is not redefined, so the value of redefinesum is

Setup(redefinesum)

[redefinesum = false]

(1.2)

Consider this multiindex summation functionality of the Physics:-Library:-Add command

Library:-Add(f[i, j], 1 <= i+j and i+j <= n)

Physics:-Library:-Add(f[i, j], i+j <= n, lowerbound = 1)

(1.3)

For instance, for n = 2,

eval(Physics[Library]:-Add(f[i, j], i+j <= n, lowerbound = 1), n = 2)

f[0, 1]+f[1, 0]+f[0, 2]+f[1, 1]+f[2, 0]

(1.4)

This functionality can be plugged directly into the sum command. For that purpose, set redefinesum to true

Setup(redefinesum = true)

[redefinesum = true]

(1.5)

You can now compute directly with sum. The left-hand side is inert while the right-hand side is computed

(%sum = sum)(f[i, j], i+j <= 2)

%sum(f[i, j], i+j <= 2) = f[0, 0]+f[0, 1]+f[1, 0]+f[0, 2]+f[1, 1]+f[2, 0]

(1.6)

(%sum = sum)(f[i, j], 1 <= i+j and i+j <= 2)

%sum(f[i, j], 1 <= i+j and i+j <= 2) = f[0, 1]+f[1, 0]+f[0, 2]+f[1, 1]+f[2, 0]

(1.7)

value(%sum(f[i, j], 1 <= i+j and i+j <= 2) = f[0, 1]+f[1, 0]+f[0, 2]+f[1, 1]+f[2, 0])

f[0, 1]+f[1, 0]+f[0, 2]+f[1, 1]+f[2, 0] = f[0, 1]+f[1, 0]+f[0, 2]+f[1, 1]+f[2, 0]

(1.8)

The formula for the integer power of a sum

(a+b+c)^n = sum(factorial(n)*a^p*b^q*c^r/(factorial(p)*factorial(q)*factorial(r)), p+q+r = n)

(a+b+c)^n = sum(Physics:-`*`(Physics:-`*`(Physics:-`*`(Physics:-`*`(factorial(n), Physics:-`*`(1, Physics:-`^`(Physics:-`*`(Physics:-`*`(factorial(p), factorial(q)), factorial(r)), -1))), Physics:-`^`(a, p)), Physics:-`^`(b, q)), Physics:-`^`(c, r)), p+q+r = n)

(1.9)

eval((a+b+c)^n = sum(Physics[`*`](Physics[`*`](Physics[`*`](Physics[`*`](factorial(n), Physics[`*`](1, Physics[`^`](Physics[`*`](Physics[`*`](factorial(p), factorial(q)), factorial(r)), -1))), Physics[`^`](a, p)), Physics[`^`](b, q)), Physics[`^`](c, r)), p+q+r = n), n = 2)

(a+b+c)^2 = a^2+2*a*b+2*a*c+b^2+2*b*c+c^2

(1.10)

eval((a+b+c)^n = sum(Physics[`*`](Physics[`*`](Physics[`*`](Physics[`*`](factorial(n), Physics[`*`](1, Physics[`^`](Physics[`*`](Physics[`*`](factorial(p), factorial(q)), factorial(r)), -1))), Physics[`^`](a, p)), Physics[`^`](b, q)), Physics[`^`](c, r)), p+q+r = n), n = 3)

(a+b+c)^3 = a^3+3*a^2*b+3*a^2*c+3*a*b^2+6*a*b*c+3*a*c^2+b^3+3*b^2*c+3*b*c^2+c^3

(1.11)

Verify whether this equation is true

(`@`(evalb, expand))((a+b+c)^3 = a^3+3*a^2*b+3*a^2*c+3*a*b^2+6*a*b*c+3*a*c^2+b^3+3*b^2*c+3*b*c^2+c^3)

true

(1.12)

Besides this new functionality, the redefined sum does a more modern handling of its arguments, consider a typical problem posted in Maple primes

a := 1; b := 2; j := 3

1

 

2

 

3

(1.13)

In the following summation, j is a dummy summation index, so the value just assigned, j := 3, is not expected to interfer with the summation. This is the case with the redefined sum

sum(f(j), j = a .. b)

f(1)+f(2)

(1.14)

while without redefining sum the input above is interrupted with an error message. Likely, in this other case also reported in Mapleprimes

g := proc (j) options operator, arrow; if j::odd then G[j] else 0 end if end proc

proc (j) options operator, arrow; if j::odd then G[j] else 0 end if end proc

(1.15)

the following two summations can be performed after having redefining sum:

sum(g(i), i = 1 .. f)

sum(g(i), i = 1 .. f)

(1.16)

For the summation above, without redefining sum, it returns 0 instead of unevaluated, because of a premature evaluation of the function g(i) with an unassigned index i before performing the summation. Returning unevaluated as (1.16) permits evaluate the sum at a latter moment, for instance attributing a value to f

eval(sum(g(i), i = 1 .. f), f = 3)

G[1]+G[3]

(1.17)

And this other sum where f is given from the begining also returns 0 without redefining sum

sum(g(i), i = 1 .. 3)

G[1]+G[3]

(1.18)

Problems like this other one reported in Mapleprimes here also get resolved with this redefinition of sum.

 

 

Download sum_in_physics.mw

Edgardo S. Cheb-Terrab
Physics, Maplesoft

 

The attached presentation is the first one of a sequence of three that we wanted to do on Quantum Mechanics using Computer Algebra. The level is that of an advanced undergraduate QM course. Tackling this topic within a computer algebra worksheet in the way it's done below, however, is an entire novelty, and illustrates well the kind of computations that can be done today with Maple & Physics.

Ground state of a quantum system of identical boson particles
  

Pascal Szriftgiser1 and Edgardo S. Cheb-Terrab2 

(1) Laboratoire PhLAM, UMR CNRS 8523, Université Lille 1, F-59655, France

(2) Maplesoft

 

Departing from the Energy of a quantum system of identical boson particles, the field equation is derived. This is the Gross-Pitaevskii equation (GPE). A continuity equation for this system is also derived, showing that the velocity flow satisfies `&x`(VectorCalculus[Nabla], `#mover(mi("v"),mo("&rarr;"))`) = 0, i.e.: is irrotational.  

The Gross-Pitaevskii equation

 

NULL


Problem: derive the field equation describing the ground state of a quantum system of identical particles (bosons), that is, the Gross-Pitaevskii equation (GPE).

 

Background: The Gross-Pitaevskii equation is particularly useful to describe Bose Einstein condensates (BEC) of cold atomic gases [3, 4, 5], that is, an ensemble of identical quantum boson particles that interact with each other with an interaction constant G. The temperature of these cold atomic gases is typically in the w100 nano-Kelvin range. The atom-atom interactions are repulsive for G > 0 and attractive for G < 0  (which could lead to some instabilities). The GPE is also widely used in non-linear optics to model the propagation of light in optical fibers. In this area, GPE is known as "non-linear Schrödinger equation", and the non-linearity comes from the Kerr effect [6].

Solution

   

Continuity equation for a quantum system of identical particles

   

References

``

[1] Gross-Pitaevskii equation (wiki)

[2] Continuity equation (wiki)
[3] Bose–Einstein condensate (wiki)

[4] Bose-Einstein Condensation in Dilute Gases, C. J. Pethick and H. Smith, Second Edition, Cambridge (2008), ISBN-13: 978-0521846516.

[5] Advances In Atomic Physics: An Overview, Claude Cohen-Tannoudji and David Guery-Odelin, World Scientific (2011), ISBN-10: 9812774963.

[6] Nonlinear Fiber Optics, Fifth Edition (Optics and Photonics), Govind Agrawal, Academic Press (2012), ISBN-13: 978-0123970237.

 


Downlioad: QuantumMechanics1.mw,    QuantumMechanics1.pdf

Edgardo S. Cheb-Terrab
Physics, Maplesoft

The DirectSearch package is a powerful Maple  tool. However, every soft has its advantages and disadvantages. In particular, the DS has problems in the case of a thin feasible set in higher dimensions. Recently a serious bug in the DS was detected by me. Solving an optimization problem, the DirectSearch produces the error communication

Warning, initial point [x1 = 1., x2 = 1., x4 = 2., y1 = 2., y2 = 3., y4 = 2.] does not satisfy the inequality constraints; trying to find a feasible initial point
Error, (in DirectSearch:-Search) cannot find feasible initial point; specify a new one
 while that initial point satisfies the constraints.

 

restart

DirectSearch:-Search(((x2-x1)^2+(y2-y1)^2)*((x4-x1)^2+(y4-y1)^2), {seq(parse(y || j) >= -(2/3)*parse(x || j)+2, j = 1 .. 4), seq(parse(y || j) >= (1/2)*parse(x || j)-3/2, j = 1 .. 4), seq(parse(y || j) <= 4, j = 1 .. 4), seq(parse(y || j) <= -3*parse(x || j)+16, j = 1 .. 4), seq(parse(y || j) <= 2*parse(x || j)+2, j = 1 .. 4), (x2-x1)*(x4-x1)+(y2-y1)*(y4-y1) = 0, (x3-x2)*(x2-x1)+(y3-y2)*(y2-y1) = 0, (x4-x1)*(x4-x3)+(y4-y1)*(y4-y3) = 0, (x4-x3)*(x3-x2)+(y4-y3)*(y3-y2) = 0}, maximize, initialpoint = [x1 = 1, x2 = 1, x3 = 2, x4 = 2, y1 = 2, y2 = 3, y3 = 3, y4 = 2])

Error, (in DirectSearch:-Search) cannot find feasible initial point; specify a new one

 

eval({seq(parse(y || j) >= -(2/3)*parse(x || j)+2, j = 1 .. 4), seq(parse(y || j) >= (1/2)*parse(x || j)-3/2, j = 1 .. 4), seq(parse(y || j) <= 4, j = 1 .. 4), seq(parse(y || j) <= -3*parse(x || j)+16, j = 1 .. 4), seq(parse(y || j) <= 2*parse(x || j)+2, j = 1 .. 4), (x2-x1)*(x4-x1)+(y2-y1)*(y4-y1) = 0, (x3-x2)*(x2-x1)+(y3-y2)*(y2-y1) = 0, (x4-x1)*(x4-x3)+(y4-y1)*(y4-y3) = 0, (x4-x3)*(x3-x2)+(y4-y3)*(y3-y2) = 0}, [x1 = 1, x2 = 1, x3 = 2, x4 = 2, y1 = 2, y2 = 3, y3 = 3, y4 = 2])

{0 = 0, -1 <= 2, -1 <= 3, 2 <= 4, 2 <= 6, 2 <= 10, 2 <= 13, 3 <= 4, 3 <= 6, 3 <= 10, 3 <= 13, -1/2 <= 2, -1/2 <= 3, 2/3 <= 2, 2/3 <= 3, 4/3 <= 2, 4/3 <= 3}

(1)

``

 

Download opti.mw

Greetings to all.

As some of you may remember I made several personal story type posts concerning my progress in solving enumeration problems with the Polya Enumeration Theorem (PET). This theorem would seem to be among the most exciting in mathematics and it is of an amazing simplicity so that I never cease to suggest to mathematics teachers to present it to gifted students even before university. My previous efforts are documented at your site, namely at this MaplePrimes link I and this MaplePrimes link II.

I have been able to do another wonderful cycle index computation using Maple recently and I would like to share the Maple code for this problem, which is posted at Math StackExchange.com (this post includes the code) This time we are trying to compute the cycle index of the automorphism group of the 3-by-3-by-3 cube under rotations and reflections. I suggest you try this problem yourself before you look at my solution. Enjoy!

I mentioned in some of my other posts concerning PET that Maple really ought to come with a library of cycle indices and the functions to manipulate them. I hope progress has been made on this issue. I had positive feedback on this at the time here at your website. Do observe that you have an opportinuity here to do very attractive mathematics if you prepare a worksheet documenting cycle index facilities that you may eventually provide. This is good publicity owing to the fact that you can include images of the many geometric objects that appear which all look quite enticing and moreover potential readers get rewarded quickly as they discover that it takes little effort to master this theorem and proceed to work with symmetries themselves and investigate them. This sort of thing also makes nice slides.

With best wishes for happy combinatorics computing,

Marko Riedel

Hi
In connection with recent developments in the Physics package, we now have mathematical typesetting for all the inert functions of the mathematical language. Hey! This is within the Physics update available on the Maplesoft Physics: Research & Development webpage

I think this is an interesting development that will concretely change the computational experience with these functions: it is not the same to compute with something you see displayed as %exp(x) instead of the same computation but flowing with it nicely displayed as an exponential function with the e in grey, reflecting that Maple understands this object as the exponential inert function, with known properties (all those of the active exp function), and so Maple can compute with the inert one taking these properties into account while not executing the function itself - and this is the essence of the inert function behaviour.

Introducing mathematical display, copy and paste for all these inert functions of the mathematical language concretely increases the mathematical expressiveness of the system, for teaching, working and also for presenting ideas.

Attached is a brief illustration.

Edgardo S. Cheb-Terrab
Physics, Maplesoft

InertMathematicalFun.mw  InertMathematicalFun.pdf

First 45 46 47 48 49 50 51 Last Page 47 of 71