MaplePrimes Posts

MaplePrimes Posts are for sharing your experiences, techniques and opinions about Maple, MapleSim and related products, as well as general interests in math and computing.

Latest Post
  • Latest Posts Feed
  • 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/

    Dear friends. 

    This is to alert you to a minor problem with your website. I use a variety of operating systems and browsers to access your site. With Firefox 22 and Firefox 26 and OpenSuse 12.2 when I click on a question with a lot of commentary like "Perl vs. Maple (MPF)" from a couple of days ago or "collect x/2-y/2 1/2" the question appears to load for a while but then an empty page appears, containing only the string "2014" and nothing else. Maybe you want to look into this.

    Best regards,

    Marko Riedel

     

    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.

     

    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.

    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.

    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

    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.

    First 79 80 81 82 83 84 85 Last Page 81 of 297