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
  • Last week, one of our Maple Learn developers, Valerie McKay-Crites, published a Maple Learn document, based on the very popular Maple application by Highschool Teacher, Jason Schattman called "Just Move It Over There, Dear!".

    In the Maple application, Schattman explains the math behind moving a rectangular sofa down a hallway with a 90-degree turn. In the 3D Moving Sofa Problem Estimate, Valerie uses Schattman’s math to determine the largest rectangular sofa that can be taken down a flight of stairs and down a hallway with a 90-degree turn. Both applications reminded me of how interesting the Moving Sofa Problem is, which inspired me to write a blog post about it!

    If you’ve ever been tasked with moving a rectangular sofa around a 90-degree turn, you might wonder:

    What is the largest sofa that can make the move?

     

     Icon

Description automatically generated with medium confidence

     

    Following these steps as outlined in Schattman’s "Just Move It Over There, Dear!", will guarantee that the sofa will make the turn:

    1. Measure the width of the hallway (h)
    2. Measure the length (L) and width (w) of the sofa.
    3. If L + 2w is comfortably less than triple the width of the hall, you'll make it!

    When we work out the math exactly, we see that if the sofa's length plus double its width is less than 2*h*sqrt(2), the sofa will make the turn!

     

    Chart, line chart

Description automatically generated

     

    This problem is easy if we only consider rectangular sofas, however, the problem becomes significantly more complex if we consider sofas of different shapes and areas. In mathematics, this problem is known as the Moving Sofa Problem, and it is unsolved. If we look at a hallway with a 90-degree turn and legs of width 1 m (i.e. h = 1 above), the largest known sofa that can make the turn is Gerver’s Sofa which has an area of 2.2195 m2, this area is known as the Sofa Constant. Gerver’s Sofa, created in 1992, was constructed with 18 curve sections:

    Icon

Description automatically generated

     

    Check out this GIF of the sofa moving through the turn. It provides some insight into why Gerver’s sofa is such an interesting shape:

    What is fascinating is that no mathematician has yet to prove that Gerver’s sofa is the sofa with the largest area capable of making the 90-degree turn.

    The Moving Sofa Problem, is a great example of how math is embedded in our everyday lives. So, don’t stop being curious about the math around you as it can be fascinating and sometimes unproven!

    If you are curious to learn more about the moving sofa problem check out this video by Numberphile, featuring Dan Romik from UC Davis: https://www.youtube.com/watch?v=rXfKWIZQIo4&t=1s

    Hello,

      2022 was a wonderful year of progress in using Maple/MapleSim for almost everything my mathematical world.     

    I just wanted to wish all the Maplesoft user community a very productive and Happy New Year for 2023.   I look forward to continue to find great nuggets of capability and insider techniques for using Maple in my endeavors.

    Kindest Regards to ALL.
    Happy New Year - 2023.
    Bill


    This post is inspired by a serie of questions from @JAMET.
    I wondered if it was possible to prove plane geometry theorems with the geometry package.

    Here is an illustration for the Poncelet's theorem for the triangle (French designation), see for instance
    https://en.wikipedia.org/wiki/Poncelet%27s_closure_theorem

    Are any of you interested in challenging the geometry package to proof other plane geometry theorems?
     

    restart:


    Poncelet's theorem for the triangle

    Let ABC a triangle, c its incircle (center omega, radius r) and C its circumcircle (center Omega, radius R).
    Let D the distance between omega and Omega.

    then R^2 - D^2 - 2*r*R = 0


    Proof

    Without loss of generality one sets :

        x(A) = y(A) = 0
       and  y(B) = 0

    ABC is a non degerated triangle provided x(B) <> 0 and y(C) <> 0
     

    with(geometry):

    kernelopts(version);

    `Maple 2020.2, X86 64 WINDOWS, Nov 11 2020, Build ID 1502365`

    (1)

    assume(x__B <> 0):
    assume(y__C <> 0):

    point(A, 0, 0);
    point(B, x__B, 0);
    point(C, x__C, y__C);

    A

     

    B

     

    C

    (2)

    triangle(T, [A, B, C])

    T

    (3)

    bisector(bA, A, T);
    bisector(bB, B, T);

    eA := isolate(Equation(bA, [x, y]), y):
    eB := isolate(Equation(bB, [x, y]), y):

    xc := solve(rhs(eA)=rhs(eB), x):
    yc := eval(rhs(eA), x=xc):

    point(omega, xc, yc):
    r := distance(omega, line(lAB, [A, B]))

    bA

     

    bB

     

    (abs(x__B)^2)^(1/2)*abs(x__B*y__C/((x__B^2-2*x__B*x__C+x__C^2+y__C^2)^(1/2)+(x__C^2+y__C^2)^(1/2)+(x__B^2)^(1/2)))/(x__B^2)^(1/2)

    (4)

    circumcircle(C, T, 'centername' = Omega);
    R := radius(C);

    C

     

    ((1/4)*x__B^2+(1/4)*(-x__B*(x__C^2+y__C^2)+x__C*x__B^2)^2/(x__B^2*y__C^2))^(1/2)

    (5)

    Oo := distance(Omega, omega)

    (((1/2)*x__B-(x__B*(x__C^2+y__C^2)^(1/2)+x__C*(x__B^2)^(1/2))/((x__B^2-2*x__B*x__C+x__C^2+y__C^2)^(1/2)+(x__C^2+y__C^2)^(1/2)+(x__B^2)^(1/2)))^2+(-(1/2)*(-x__B*(x__C^2+y__C^2)+x__C*x__B^2)/(x__B*y__C)-y__C*(x__B^2)^(1/2)/((x__B^2-2*x__B*x__C+x__C^2+y__C^2)^(1/2)+(x__C^2+y__C^2)^(1/2)+(x__B^2)^(1/2)))^2)^(1/2)

    (6)

    S := simplify(R^2 - Oo^2 - 2*r*R)  assuming x__B::real, x__C::real, y__C::real

    ((x__B^2*(-abs(y__C)*signum(y__C)+y__C)*(x__C^2+y__C^2)^(1/2)+x__C^2*(y__C*abs(x__B)-signum(y__C)*abs(x__B*y__C)))*(x__B^2-2*x__B*x__C+x__C^2+y__C^2)^(1/2)+(-signum(y__C)*(x__C-x__B)^2*abs(x__B*y__C)+(abs(x__B)^2+x__C*(x__C-2*x__B))*abs(x__B)*y__C)*(x__C^2+y__C^2)^(1/2))/(y__C*((x__B^2-2*x__B*x__C+x__C^2+y__C^2)^(1/2)+(x__C^2+y__C^2)^(1/2)+abs(x__B))^2)

    (7)

    simplify(S) assuming x__B > 0, y__C > 0;
    simplify(S) assuming x__B > 0, y__C < 0;
    simplify(S) assuming x__B < 0, y__C > 0;
    simplify(S) assuming x__B < 0, y__C < 0;

    0

     

    0

     

    0

     

    0

    (8)

     

     

    Download PoTh_proof.mw


    Improvements of the geometry package                                                                                           

    It already appears that (some) assumptions are not (always) correctly taken into account. This is a weak point which requires, as in the attached mw, to use an indirect approache to construct the incircle.

    As a matter of fact, the procedure incircle, whose first lines are

    showstat(incircle)
    
    geometry:-incircle := proc(inci, T)
    local cname, d, A, B, l1, l2, dis, x, y, tmp, msg;
       1   if nargs < 2 or 3 < nargs then
       2     error "wrong number of arguments"
           end if;
       3   if geometry:-form(T) <> ('triangle2d') then
       4     error "wrong type of arguments"
           end if;
       5   if nargs = 3 and op(1,args[3]) = ('centername') and type(op(2,args[3]),'name') then
       6     cname := op(2,args[3])
           else
       7     cname := cat('center_',inci)
           end if;
       8   if geometry:-method(T) = (':-points') then
       9     d := geometry:-DefinedAs(T);
      10     A := op(1,d);
      11     B := op(2,d);
      12     msg := sprintf("find the bisector of %a at vertex %a",T,A);
      13     userinfo(2,geometry,msg);
      14     geometry:-bisector('l1',A,T);
      15     msg := sprintf("find the bisector of %a at vertex %a",T,B);
      16     userinfo(2,geometry,msg);
      17     geometry:-bisector('l2',B,T);
      18     msg := sprintf("find the intersection of the two bisectors");
      19     userinfo(2,geometry,msg);
      20     geometry:-intersection(cname,l1,l2);
    

    requires that the two bissectors are not parallel(call to geometry:-intersection).

    Since the non-parallelism of bisectors is trivial for all non-degenerate triangles, why doesn't incircle inherit this property rather than not being able to decide if the bisectors are parallel or not?)

    Here is a detail of what happens and the endless loop in which incircle seems to be caught

    kernelopts(version);
                      Maple 2015.2, APPLE UNIVERSAL OSX, Dec 20 2015, Build ID 1097895
    
    AreParallel(bA, bB, 'condition'):
            AreParallel: hint: cannot determine if -y__C*(x__B^2)^(1/2)*(x__C*(x__B^2)^(1/2)-x__B*(x__B^2)^(1/2)-x__B*((x__B-x__C)^2+y__C^2)^(1/2))+y__C*(x__B^2)^(1/2)*(x__B*(x__C^2+y__C^2)^(1/2)+x__C*(x__B^2)^(1/2)) is zero
    
    assume(lhs(condition) <> 0);
    AreParallel(bA, bB);
                                 false
    intersection(J, bA, bB);
                                   J
    
    
    infolevel[geometry] := 4:
    incircle(inc, T);
    incircle: find the bisector of T at vertex A
    incircle: find the bisector of T at vertex B
    incircle: find the intersection of the two bisectors
    intersection: find the intersection between two lines l1 and l2
    intersection: one point of intersection
    incircle: find the radius of the incircle
    line: define the line from two points
    circle: define the circle from its center and radius
    circle: hint: abs(x__B^2*y__C/(csgn(x__B)*x__B+(x__B^2-2*x__B*x__C+x__C^2+y__C^2)^(1/2)+(x__C^2+y__C^2)^(1/2)))/(x__B^2)^(1/2) > 0
    Error, (in geometry:-circle) not enough information: the radius might not be positive
    assume(abs(x__B^2*y__C/(csgn(x__B)*x__B+(x__B^2-2*x__B*x__C+x__C^2+y__C^2)^(1/2)+(x__C^2+y__C^2)^(1/2)))/(x__B^2)^(1/2) > 0):
    
    assume(abs(x__B^2*y__C/(csgn(x__B)*x__B+(x__B^2-2*x__B*x__C+x__C^2+y__C^2)^(1/2)+(x__C^2+y__C^2)^(1/2)))/(x__B^2)^(1/2) > 0): 
    incircle(inc, T);
    incircle: find the bisector of T at vertex A
    incircle: find the bisector of T at vertex B
    incircle: find the intersection of the two bisectors
    intersection: find the intersection between two lines l1 and l2
    AreParallel: hint: cannot determine if -y__C*(x__B^2)^(1/2)*(x__C*(x__B^2)^(1/2)-x__B*(x__B^2)^(1/2)-x__B*((x__B-x__C)^2+y__C^2)^(1/2))+y__C*(x__B^2)^(1/2)*(x__B*(x__C^2+y__C^2)^(1/2)+x__C*(x__B^2)^(1/2)) is zero
    intersection: two given lines intersect each other if -y__C*(x__B^2)^(1/2)*(x__C*(x__B^2)^(1/2)-x__B*(x__B^2)^(1/2)-x__B*((x__B-x__C)^2+y__C^2)^(1/2))+y__C*(x__B^2)^(1/2)*(x__B*(x__C^2+y__C^2)^(1/2)+x__C*(x__B^2)^(1/2)) <> 0
    Error, (in geometry:-intersection) not enough information
    
    assume(-y__C*(x__B^2)^(1/2)*(x__C*(x__B^2)^(1/2)-x__B*(x__B^2)^(1/2)-x__B*((x__B-x__C)^2+y__C^2)^(1/2))+y__C*(x__B^2)^(1/2)*(x__B*(x__C^2+y__C^2)^(1/2)+x__C*(x__B^2)^(1/2)) <> 0):
    incircle(inc, T);
    incircle: find the bisector of T at vertex A
    incircle: find the bisector of T at vertex B
    incircle: find the intersection of the two bisectors
    intersection: find the intersection between two lines l1 and l2
    intersection: one point of intersection
    incircle: find the radius of the incircle
    line: define the line from two points
    circle: define the circle from its center and radius
    circle: hint: abs(x__B^2*y__C/(csgn(x__B)*x__B+(x__C^2+y__C^2)^(1/2)+(x__B^2-2*x__B*x__C+x__C^2+y__C^2)^(1/2)))/(x__B^2)^(1/2) > 0
    Error, (in geometry:-circle) not enough information: the radius might not be positive
    
    
    

     

     

    Recently here Mapleprimes member @Icz asked about generating all minimal edge cuts of a graph. I gave a brute force algorithm based on testng all bipartitions of the vertices. I then looked into improving the method by using fundamental cutsets to generate the cutsets. A description of the method and the ideas behind it is given here, together with a comparison of the two methods.

    [Edit - see below for a newer version.]

    Find all minimal edge cuts of a connected undirected graph using fundamental cutsets. (Assume no self-loops or multiple edges.)
    The MinimalEdgeCuts procedure is in the startup code edit region (See Edit -> Startup Code.) (and is reproduced at the end of the worksheet).

    restart

    with(GraphTheory)

    Choose a graph.

    G := Graph({{1, 2}, {1, 3}, {1, 4}, {2, 3}, {3, 4}}); DrawGraph(G, size = [200, 200])

    GRAPHLN(undirected, unweighted, [1, 2, 3, 4], Array(1..4, {(1) = {2, 3, 4}, (2) = {1, 3}, (3) = {1, 2, 4}, (4) = {1, 3}}), `GRAPHLN/table/1`, 0)

    The objective is to find all minimal edge cuts. If we partition the vertices of a graph in two non-empty sets (parts), then a cut is defined as a set of edges of the graph that have one end in the first part and the other end in the second part. Removal of these edges breaks the graph into two (or more) disconnected components.

    For example the cut with edges {1,2} and {2,3} leaves two components with vertices [1,3,4] in one component and vertex [2] in the other component.

    cut := {{1, 2}, {2, 3}}; IsCutSet(G, cut); DrawGraph(DeleteEdge(G, cut, inplace = false), size = [200, 200])

    {{1, 2}, {2, 3}}

    true

    A cut is minimal if adding back any of the edges makes it connected again, i.e, the cut leave only two components. The above cut is minimal. Sometimes the term cutset is used to mean a minimal cut, but sometimes, as in Maple's IsCutSet, cutset means any removal of edges that increases the number of components.

    This graph has 6 minimal cuts:

    mincuts := MinimalEdgeCuts(G)

    Vector[column](%id = 36893490386482949164)

    This partitioning process suggests a brute force algorithm to find all minimal cuts. Find all partitions of the vertices into two (non-empty) sets. The two induced subgraphs G__1 and G__2 have some of the edges of the graph. All edges of the graph that are not in the two subgraphs form a cut. For example for the above partition we have:

    partition := [[1, 3, 4], [2]]; G__1 := InducedSubgraph(G, partition[1]); 'Edges(G__1)' = Edges(G__1); G__2 := InducedSubgraph(G, partition[2]); 'Edges(G__2)' = Edges(G__2); cut := `minus`(`minus`(Edges(G), Edges(G__1)), Edges(G__2))

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

    GraphTheory:-Edges(G__1) = {{1, 3}, {1, 4}, {3, 4}}

    GraphTheory:-Edges(G__2) = {}

    {{1, 2}, {2, 3}}

    Now we have to check that each component is connected, or that there are a total of two components. Since there are in this case, the cut is minimal.

    IsConnected(G__1) and IsConnected(G__2); nops(ConnectedComponents(DeleteEdge(G, cut, inplace = false)))

    true

    2

    Note that we do have to test that there are only two components. For example for the path graph 1--2--3, the partition [[2],[1,3]] has three components (the three individual vertices). The associated cut {{1,2},{2,3}} is not minimal because we could put edge {1,2} back and the graph would still be disconnected.

     

    The number of partitions to be tested is 2^(n-1)-1, where n is the number of vertices. Any subset of the vertices could be in the first set (except the empty set and the set of all vertices), and then the other vertices must be in the second set. However that will give double the possibilities, since swapping set 1 and set 2 makes no difference. In the present case there are 2^3-1 = 7 possibilities, as seen below, where 0 indicates that the vertex in that position is in set 1 and 1 indicates that it is in set 2. For example, possibility 4 is the above partition [[1, 3, 4], [2]].

    lprint("   1 2 3 4"); Print(Iterator:-SetPartitions(4, parts = 2), showrank)

       1 2 3 4
    1: 0 0 0 1
    2: 0 0 1 0
    3: 0 0 1 1
    4: 0 1 0 0
    5: 0 1 0 1
    6: 0 1 1 0
    7: 0 1 1 1

    The algorithm described is implemented as procedure MinimalEdgeCutsPartition in the startup code. (The minimal cuts are not produced in the same order as above.)

     

    However, a better approach is by building up the cuts fron a set of fundamental cuts (or cutsets, using the word in the restrictive sense). These are associated with a spanning tree. We choose any spanning tree, which we highlight with red edges. This is a tree graph (connected with no cycles) that has the same vertices as the graph.

    treeG := SpanningTree(G); HighlightEdges(G, treeG); DrawGraph(G, size = [200, 200])

    GRAPHLN(undirected, unweighted, [1, 2, 3, 4], Array(1..4, {(1) = {2, 3, 4}, (2) = {1}, (3) = {1}, (4) = {1}}), `GRAPHLN/table/23`, 0)

    Any spanning tree graph has n-1 edges. This is also the number of vertices minus the number of connected components (=1), which is the rank of the graph, which we will call r.
    We will find one fundamental cutset for each edge of the tree, so there will be r fundamental cutsets.

    NumberOfVertices(G)-1, NumberOfEdges(treeG), GraphRank(G)

    3, 3, 3

    Select one of the tree edges, say {1,2}. Removal of this edge partitions the tree into two components. We find the vertex partition and the corresponding cut by the procedure above. (We do not have to test for only two components because that follows from the tree structure.). This is the example we did above. The first three cuts in mincuts are the fundamental ones. Each one has one of the tree edges.

    f[1], f[2], f[3] := entries(mincuts[1 .. 3], nolist)

    {{1, 2}, {2, 3}}, {{1, 3}, {2, 3}, {3, 4}}, {{1, 4}, {3, 4}}

    It is possible to generate all the cuts of the graph involving the tree edges by taking all 2^r = 8 subsets of the set of fundamental cutsets (this includes the empty set), and adding the elements of each subset with a ring sum operation.

    A ring sum of two edge sets includes all edges in both sets except those that are common to both (the symmetric difference or disjunctive union):

    `&oplus;` := proc (edges__1, edges__2) options operator, arrow; `minus`(`union`(edges__1, edges__2), `intersect`(edges__1, edges__2)) end proc

    proc (edges__1, edges__2) options operator, arrow; `minus`(`union`(edges__1, edges__2), `intersect`(edges__1, edges__2)) end proc

    (Maple has this as the symmdiff command, but not in infix form.)
    We can also represent a set of edges by a vector, with entry 1 if the corresponding edge is in the set and 0 otherwise. The edge order we choose here is the order in Edges(G). Then the ring sum operation is the same as addition modulo 2 of the corresponding vectors, though we won't make much use of this in implementing the algorithm. Set up to convert to vectors.

    edges := Edges(G); e := NumberOfEdges(G); edgetable := table(`~`[`=`](edges, [`$`(1 .. e)])); tovec := proc (edges) options operator, arrow; Vector(e, `~`[`=`](map(proc (x) options operator, arrow; edgetable[x] end proc, edges), 1)) end proc

    {{1, 2}, {1, 3}, {1, 4}, {2, 3}, {3, 4}}

    Let's try the ring sum on f[1] and f[2]. It is a cut, and since it now has two tree edges it is not a fundamental cut. (It happens to be minimal because it splits the graph into two components.)

    f[1]*`&oplus;`*f[2] = `&oplus;`(f[1], f[2]); IsCutSet(G, rhs(%))

    {{1, 2}, {2, 3}}*`&oplus;`*{{1, 3}, {2, 3}, {3, 4}} = {{1, 2}, {1, 3}, {3, 4}}

    true

    Do the same thing in the vector representation. The same logic about the tree edges (first three entries) shows the fundamental vectors are basis vectors and the three vectors V__1, V__2, `mod`(V__1+V__2, 2)are linearly independent (with respect to the addition mod 2 operation).

    V[1], V[2], V[3] := tovec(f[1]), tovec(f[2]), tovec(f[3]); '`mod`(V[1]+V[2], 2)' = `mod`(V[1]+V[2], 2)

    V[1], V[2], V[3] := Vector(5, {(1) = 1, (2) = 0, (3) = 0, (4) = 1, (5) = 0}), Vector(5, {(1) = 0, (2) = 1, (3) = 0, (4) = 1, (5) = 1}), Vector(5, {(1) = 0, (2) = 0, (3) = 1, (4) = 0, (5) = 1})

    `mod`(V[1]+V[2], 2) = Vector[column](%id = 36893490386597633668)

    The ring sum of any two cuts is another cut, which may be minimal (a cutset) or a disjoint union of cutsets (= two cutsets together). Combining the elements of the 8 subsets gives the eight cuts for this graph associated with the chosen tree graph, i.e., it has all combinations of tree edges. In vector form they are:

    cutvecs := seq(seq(seq(`mod`(i*V[1]+j*V[2]+k*V[3], 2), `in`(i, [0, 1])), `in`(j, [0, 1])), `in`(k, [0, 1]))

    Vector[column](%id = 36893490386597646332), Vector[column](%id = 36893490386597647052), Vector[column](%id = 36893490386597647772), Vector[column](%id = 36893490386597648252), Vector[column](%id = 36893490386597648972), Vector[column](%id = 36893490386597649452), Vector[column](%id = 36893490386597649932), Vector[column](%id = 36893490386597650172)

    Although we constructed these cuts relative to a specific spanning tree, they can be shown to be all the cuts of the graph. Consider the set of all edges of the graph. It might seem that this is a cut that has been missed, because it separates the graph into two or more (actually 4) components. However, looking at the edges

    edges

    {{1, 2}, {1, 3}, {1, 4}, {2, 3}, {3, 4}}

    and recalling the definition of a cut, we see that this is not actually a cut. We need the edges in a cut all to have one end in one set of vertices and the other end in the set of other vertices, and this is not possible here. (Maple's IsCutSet is not this strict and returns true.) We can show it is not possible with IsBipartite.

    IsBipartite(Graph(edges))

    false

    Although we have found all cuts, we now need to test them to find which ones are minimal. The first case with no edges does not cut the graph into two components and we reject it. (It is usually included as a cut to make the cuts into a vector space.)  The 6th one is not minimal. Its first and third entries are 1, meaning it is the ring sum of the first and third fundamental cuts. It splits the graph into three components:

    cut := `&oplus;`(f[1], f[3]); DrawGraph(DeleteEdge(G, cut, inplace = false), size = [200, 200])

    {{1, 2}, {1, 4}, {2, 3}, {3, 4}}

    This is a case where the ring sum is a disjoint union (not minimal), and arises because f[1] and f[3] have no edges in common:

    `intersect`(f[1], f[3])

    {}

    For the present graph, all but the empty set and `&oplus;`(f[1], f[3]) are minimal, and so there are 6 minimal cuts. In vector form:

    cutvecs := convert(map(tovec, mincuts), set)[]

    Vector[column](%id = 36893490386597690532), Vector[column](%id = 36893490386597690652), Vector[column](%id = 36893490386597690772), Vector[column](%id = 36893490386597690892), Vector[column](%id = 36893490386597691012), Vector[column](%id = 36893490386597691132)

    So the algorithm is to generate all cuts and test them for minimality. If we neglect the empty cut, then we have to generate 2^r-1 = 2^(n-1)-1 cuts for testing. This is the same as the number of partitions that we generated in the partition algorithm. However, for the fundamental cuts we do not need to test that the components are connected, and more importantly, the ring sum is a more efficient way of finding the edges than partitioning and then using InducedSubgraphs.

     

    The tests for minimality can be made more efficient by finding cases that can be classified without having to check the number of components. Two simple cases with modest savings (implemented here) are:

    1. The fundamental cuts are minimal by construction.

    2. If two cuts have no edges in common, then the ring sum is a disjoint union and is not minimal. This is the case for example for f[1] and f[3].

     

    There may be other ways to classify without doing the more difficult minimality test, for example variations of (2.) for ring sums with more than two cuts. Hovever, the cutset algorithm is already significantly more efficient than the partition algorithm. Perhaps working with the vectors would be more efficient. Other more efficient algorithms may be known that I am not aware of.

     

    The MinEdgeCuts algorithm in the startup code is reproduced here:

     

    MinimalEdgeCuts:=proc(G::GRAPHLN)
      local mincutsets,n,edges,treeedges,treeG,
        i,j,r,vertexpartition,partitionedges,nf,fc,
        fsets,ref,fcutsets,cutsetedges;
      uses GraphTheory;
      if not IsConnected(G) then return Vector([]) end if;
      n:=NumberOfVertices(G);
      edges:=Edges(G);
      # choose a spanning tree and find corresponding fundamental cutsets
      treeG:=SpanningTree(G);
      treeedges:=Edges(treeG);
      r:=n-1;    # =nops(treeedges) = GraphRank(G)
      mincutsets:=table();  # first r ones are the fundamental ones
      for j to r do
        vertexpartition:=ConnectedComponents(DeleteEdge(treeG,treeedges[j],'inplace'=false));
        partitionedges:=`union`(map(Edges,map2(InducedSubgraph,G,vertexpartition))[]);
        mincutsets[j]:=edges minus partitionedges;
      end do;
      j:=r;
      # now find ringsums of all subsets of the set of fundamental cutsets
      fsets := Iterator:-BinaryGrayCode(r,'rank'=2); #skip empty set
      for ref in fsets do;
        fcutsets:=seq(ifelse(ref[i]=1,mincutsets[i],NULL),i=1..r);
        nf:=nops([fcutsets]);
        if nf=1 then              
          next      # fundamental cutsets already in mincutsets
        elif nf=2 and (fcutsets[1] intersect fcutsets[2] = {}) then
          next      # pair of disjoint cutsets not minimal
        # todo other detections of non-minimal cutsets
        else        # rest the hard way
          cutsetedges:=fcutsets[1];
          for i from 2 to nf do
                  fc:=fcutsets[i];
                  cutsetedges:=(cutsetedges union fc) minus (cutsetedges intersect fc);
          end do;
          vertexpartition:=ConnectedComponents(DeleteEdge(G,cutsetedges,'inplace'=false));
          if nops(vertexpartition)=2 then # cutset is minimal
            j:=j+1;
            mincutsets[j]:=cutsetedges
          end if
        end if;
      end do;
      Vector(j,mincutsets);
    end proc:

    NULL

    Download CutSets.mw

    Does everyone have a good idea of ​​the work of the Draghilev method? For example, in this answer https://www.mapleprimes.com/questions/235407-The-Second-Example-Of-Finding-All-Solutions#answer291268
    ( https://www.mapleprimes.com/questions/235407-The-Second-Example-Of-Finding-All-Solutions )
    there was a very successful attempt by Rouben Rostamian to calculate the line of intersection of surfaces without applying the Draghilev method.
    Let now not 3d, but 8d. And how will the solve command work in this case? Imagine that aij ((i=1..7,j=1..8)) are partial derivatives, and xj (,j=1..8) are derivatives, as in the above example. f8 is responsible for the parametrization condition.

     

    restart;
     f1 := a11*x1+a12*x2+a13*x3+a14*x4+a15*x5+a16*x6+a17*x7+a18*x8; 
     f2 := a21*x1+a22*x2+a23*x3+a24*x4+a25*x5+a26*x6+a27*x7+a28*x8; 
     f3 := a31*x1+a32*x2+a33*x3+a34*x4+a35*x5+a36*x6+a37*x7+a38*x8; 
     f4 := a41*x1+a42*x2+a43*x3+a44*x4+a45*x5+a46*x6+a47*x7+a48*x8;
     f5 := a51*x1+a52*x2+a53*x3+a54*x4+a55*x5+a56*x6+a57*x7+a58*x8; 
     f6 := a61*x1+a62*x2+a63*x3+a64*x4+a65*x5+a66*x6+a67*x7+a68*x8; 
     f7 := a71*x1+a72*x2+a73*x3+a74*x4+a75*x5+a76*x6+a77*x7+a78*x8;
     f8 := x1^2+x2^2+x3^2+x4^2+x5^2+x6^2+x7^2+x8^2-1; 
    allvalues(solve({f1, f2, f3, f4, f5, f6, f7, f8}, {x1, x2, x3, x4, x5, x6, x7, x8}))


    And this is how the Draghilev method works in this case.
     

    restart; with(LinearAlgebra):
    f1 := a11*x1+a12*x2+a13*x3+a14*x4+a15*x5+a16*x6+a17*x7+a18*x8; 
    f2 := a21*x1+a22*x2+a23*x3+a24*x4+a25*x5+a26*x6+a27*x7+a28*x8; 
    f3 := a31*x1+a32*x2+a33*x3+a34*x4+a35*x5+a36*x6+a37*x7+a38*x8; 
    f4 := a41*x1+a42*x2+a43*x3+a44*x4+a45*x5+a46*x6+a47*x7+a48*x8; 
    f5 := a51*x1+a52*x2+a53*x3+a54*x4+a55*x5+a56*x6+a57*x7+a58*x8;
    f6 := a61*x1+a62*x2+a63*x3+a64*x4+a65*x5+a66*x6+a67*x7+a68*x8; 
    f7 := a71*x1+a72*x2+a73*x3+a74*x4+a75*x5+a76*x6+a77*x7+a78*x8;
    
    n := 7; 
    x := seq(eval(cat('x', i)), i = 1 .. n+1):
     F := [seq(eval(cat('f', i)), i = 1 .. n)]: 
    A := Matrix(nops(F), nops(F)+1):
     for j to nops(F) do for i to nops(F)+1 do A[j, i] := op(1, op(i, op(j, F))) 
    end do:
               end do: 
    
    # b[i] and b[n+1] are solutions of a linear homogeneous system and at the 
    # same time they serve as the right-hand sides of an autonomous ODE.
    for i to n do
    
     b[i] := Determinant(DeleteColumn(ColumnOperation(A, [i, n+1]), n+1)) 
                                                                         end do:
     b[n+1] := -Determinant(DeleteColumn(A, n+1)):
    


    Only the original seven linear homogeneous equations with eight variables are needed. We solve them according to Cramer's rule, and in order to have uniqueness when solving the ODE, we use a point on the curve (according to the theory). (This point is sought in any convenient way.)
    If we want to get a parameterization, then additionally, directly in dsolve, we can add the following:
     

    for i to n do 
    b[i] := simplify(Determinant(DeleteColumn(ColumnOperation(A, [i, n+1]), n+1))) 
    end do:
    b[n+1] := simplify(-Determinant(DeleteColumn(A, n+1))); 
    deqs := seq(diff(x[i](s), s) = b[i]/(b[1]^2+b[2]^2+b[3]^2+b[4]^2+b[5]^2+b[6]^2+b[7]^2+b[8]^2)^.5, i = 1 .. n+1):

     

    In an old post, vv reported a bug in simpl/max, which has been "fixed" in Maple 2018. However, it seem that such repairs are not complete enough.
    For example, suppose it is required to find the (squared) distance between the origin and a point on x3 - x + y2 = ⅓ which is closest to the origin. In other words, one needs to minimize x²+y² among the points on this curve, i.e., 

    extrema(x^2 + y^2, {x^3 + y^2 - x = 1/3}, {x, y}, 's'); # in exact form

    Unfortunately, an identical error message appears again: 

    restart;

    extrema(x^2+y^2, {x^3+y^2-x = -2/(3*sqrt(3))}, {x, y})

    {4/3}

    (1)

    extrema(x^2+y^2, {x^3+y^2-x = 1/3}, {x, y})

    Error, (in simpl/max) complex argument to max/min: 1/36*((36+12*I*3^(1/2))^(2/3)+12)^2/(36+12*I*3^(1/2))^(2/3)

     

    `~`[`^`](extrema(sqrt(x^2+y^2), {x^3+y^2-x = 1/3}, {x, y}), 2)

    {4/3, 4/27}

    (2)

    extrema(x^2+1/3-x^3+x, {x^3+y^2-x = 1/3}, {x, y})

    {4/3, 4/27}

    (3)

    MTM[limit](extrema(x^2+y^2, {x^3+y^2-x = a}, {x, y}), 1/3)

    {4/3, 4/27}

    (4)

    Download tryHarder.mws

    How about changing the values of parameter ?

    for a from -3 by 3/27 to 3 do
        try
            extrema(x^2 + y^2, {x^3 + y^2 - x = a}, {x, y}); 
        catch:
            print(a); 
        end;
    od;
                                   -1
                                   --
                                   3 
    
                                   -2
                                   --
                                   9 
    
                                   -1
                                   --
                                   9 
    
                                   1
                                   -
                                   9
    
                                   2
                                   -
                                   9
    
                                   1
                                   -
                                   3
    

    By the way, like extrema, Student[MultivariateCalculus]:-LagrangeMultipliers also executes the Lagrange Multiplier method, but strangely, 

    Student[MultivariateCalculus][LagrangeMultipliers](y^2 + x^2, [x^3 + y^2 - x - 1/3], [x, y], output = plot):

    does not cause any errors.

    This is about functionality introduced in Maple 2022, which however is still not well known: Integral Vector Calculus and parametrization using symbolic (algebraic) vector notation. Four new commands were added to the Physics:-Vectors package, implementing the parametrization of curves, surfaces and volumes, as well as the computation of path, surface and volume vector integrals. Those are integrals where the integrand is a scalar or vector function. The computation is done from any description (algebraic, parametric, vectorial) of the region of integration - a path, surface or volume.
     
    There are three kinds of line or path integrals:

    NOTE Jan 1: Updated the worksheet linked below; it runs in Maple 2022.
    Download Integral_Vector_Calculus_and_Parametrization.mw

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

     

    I have been making animated 3d plots recently; the last time was perhaps three years ago, and I had some problems then.  If I recall correctly, I couldn't make an animated 3d plot that was plotted in non-Cartesian coordinates.

     

    I am very happy to report that this works very smoothly now in Maple 2022, and it's pretty fast, too.  I have a fairly complex function to plot, involving piecewise polynomials on a tensor product grid in the xi and eta variables (actually, I let plot3d pick out the grid; it seems happier to do so) and then plot them on an elliptical base, in coordinates x = d*cosh(xi)*cos(eta) and y=d*sinh(xi)*sin(eta)  (d is just a numerical constant, giving the location of the foci at (d,0) and (-d,0)), for 0 <= xi <= xi[0] (the outer elliptical boundary) and 0 <= eta <= 2Pi.  The straightforward command works, and building a sequence of plots and using plots[display] works.  I put option remember into my procedure w(xi,eta) and because the sample points are consistent for the time-dependent function exp(I*omega*t)*w(xi,eta) the xi-eta grid needs only to be done once and then one can compute (basically) as many frames as one wants in rapid succession.

     

    Works great.  Thanks, folks!

     

    for k to nplots do
        t := evalf(2*Pi*(k - 1)/nplots);
        plts[k] := plot3d([(xi, eta) -> focus*cosh(xi)*cos(eta), (xi, eta) -> focus*sinh(xi)*sin(eta), (xi, eta) -> Re(exp(omega*t*I)*w(xi, eta))], 0 .. xi[0], 0 .. 2*Pi, colour = ((xi, eta) -> Re(exp(omega*t*I)*w(xi, eta))), style = surfacecontour, lightmodel = "none");
    end do;
    plots[display](seq(plts[k], k = 1 .. nplots), insequence = true);
     

    With the winter solstice speeding towards us, we thought we’d create some winter themed documents. Now that they’re here, it’s time to show you all! You’ll see two new puzzle documents in this post, along with three informative documents, so keep reading.

    Let’s start with the tromino tree!

     

    First, what’s a tromino? A tromino is a shape made from three equal sized squares, connected to the next along one full edge. In this puzzle, your goal is to take the trominos, and try to fill the Christmas tree shape.

    There’s a smaller and larger tree shape, for different difficulties. Try and see how many ways you can fill the trees!

    Next, we’ll look at our merry modulo color by numbers.

    Table

Description automatically generated

    In this puzzle, your goal is to solve the modulo problems in each square, and then fill in the square with the color that corresponds to the answer. Have fun solving the puzzle and seeing what the image is in the end!

    Snowballs are a quintessential part of any winter season, and we’ve got two documents featuring them.

    A picture containing icon

Description automatically generated

    The first document uses a snowball rolling down a hill to illustrate a problem using differential equations. Disclaimer: The model is not intended to be realistic and is simplified for ease of illustration. This document features a unique visualization you shouldn’t miss!

    Our second document featuring snowballs talks about finding the area of a 2-dimensional snowman! Using the formula for the area of a circle and a scale factor, the document walks through finding the area in a clear manner, with a cute snowman illustration to match!

    Shape

Description automatically generated

    The final document in this mini-series looks at Koch snowflakes, a type of fractal. This document walks you through the steps to create an iteration of the Koch snowflake and contains an interactive diagram to check your drawings with!

    I hope you’ve enjoyed taking a look at our winter documents! Please let us know if there’s any other documents you’d like to see featured or created.

    This command should have been in Physics on day one. Being more familiar with functional differentiation, and Physics:-Fundiff was the first Physics command that ever existed, I postponed writing LagrangeEquations year after year. In general, however, functional differentiation is seen as a more advanced topic. So there is now a new command, Physics:-LagrangeEquations, taking advantage of functional differentiation on background, and distributed for everybody using Maple 2022.2 within the Maplesoft Physics Updates. This is the first version of its help page.


    Download LagrangeEquations.mw

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

    Welcome back to another Maple Learn blog post! Today we’re going to talk about the gift-wrapping algorithm, used to find the convex hull of a set of points. If you’re not sure what that means yet, don’t worry! We’re going to go through it with four Maple Learn documents; two which are background information on the topic, one that is a visualization for the gift-wrapping algorithm, and another that goes through the steps. Each will be under their own heading, so feel free to skip ahead to your skill level!

    Before we can get into the gift-wrapping algorithm we need to define a few terms. Let’s start by defining polygons and simple polygons.

    A picture containing text, person

Description automatically generated

    Polygon: A closed shape created by joining a series of line segments.

    Simple polygon: A polygon without holes and that does not intersect itself.

    Shape, polygon

Description automatically generated

    So, what are convex and concave polygons? Well, there are three criteria that define a convex polygon. A polygon that is not convex is called concave. The criteria are…

    1. Any line segment connecting any two points within the polygon stays within the polygon.
    2. Any line intersects a polygon’s boundary at most twice.
    3. All interior angles are less than 180 degrees or pi radians.

    A picture containing chart

Description automatically generated

    Because the criteria are equivalent, if any one is missing, the shape is concave. AKA, all three criteria must be present for a shape to be convex. Most “regular shapes”, such as trapezoids, are convex polygons!

    A shape that satisfies convex criteria but not the criteria for being a polygon is called a convex set.

    As mentioned at the start of this post, the gift-wrapping algorithm is used to find the convex hull of a set of points. Now that we know what convex polygons and convex sets are, we can define the convex hull!

    Convex hull: The convex set of a shape or several shapes that fully contains the object and has the smallest possible area.

    A picture containing text, stationary, envelope, businesscard

Description automatically generated

    Why was the convex polygon important? Well, the convex hull of a set of points is always a convex polygon. Some of the points in the set are the vertices of said polygon, and are called extreme points. You can find the convex hull of either concave or convex polygons.

    This document amazed me when I tried it for the first time. Here, you can generate a set of points with the “Generate Another” button, and then press the “Visualize” button. The document then calculates the perimeter of the convex hull of the set of points! The set can be further customized below the buttons, by changing the number of points. The other option below it allows you to slow down or speed up the visualization. Pretty cool, huh? It’s like it’s thinking!

    Try the document out a few times, or watch the gif below to get a quick idea of it.

    This final document walks you through the steps of how to use the gift wrapping algorithm. It is a simple loop of 4 steps, with one set-up step. Unlike the other documents in this post, I won’t be delving too far into the math behind the steps. I want to encourage you to check this one out yourself, as it’s really quite a fun problem to solve once you have some time!

    Chart, line chart

Description automatically generated

    I hope you check out the documents in this post. Please let us know below if there’s any other documents you’d like to see featured!

    A failing slinky is another intriguing physics phenome that can be easily reproduced with MapleSim.

    The bottom of a vertically suspended slinky does not move when the top is released until the slinky is fully collapsed.

     

     

    To model this realistically in MapleSim, it is necessary to

    • Establish a stretched equilibrium state at the start of the fall
    • Avoid penetration of windings when windings collapse (i.e. get into contact)

    The equilibrium state is achieved with the snapshot option. Penetration is avoided with the Elasto Gap component. Details can be found in the attached model.

    A good overview of “Slinky research” is given here. The paper provides a continuous description of the collapse process (using an inhomogenous wave equation combined with contact modeling!!!) and introduces a finite time for the collapse of all windings. Results for a slinky are presented that collapses after 0.27s. The attached model has sufficient fidelity to collapse at the same time.

    Real Slinkies also feature a torsional wave that precedes the compression wave and disturbs an ideal collapse. This can be seen on slowmo footage and advanced computer models. With a torsion spring constant at hand (are there formulas for coil springs?), it could also be modeled with MapleSim.

    Falling_slinky.msim

    Have you ever heard of the Maurer Rose?

    The Maurer Rose was demonstrated in 1987 by Peter Maurer and is created by connecting certain points on a rose curve. This creates petal-like patterns, caused by the oscillation of a sine curve.

     

    Chart, radar chart

Description automatically generated

    So, how are these created? A "rose curve" is created in polar coordinates with the equation sin(nt) for a (positive integer) value of n.  To create the Maurer Rose, straight line segments are drawn connecting points on the curve at incrementing angle values.  The size of this increment (called d in our examples) leads to different patterns of lines across the curve.

    This can be done in Maple Learn! One example of the Maurer Rose already exists, complete with a full interactive visualization and a more detailed overview of the Maurer Rose.

    Play around with it and look below at some of the different shapes that can be created using this document! The first is created with an n value of 31 and a d value of 65, with blue and red. The second uses an n value of 4 and a d value of 133, and purple and green.

    Chart, diagram, radar chart

Description automatically generated

    Are there any other concepts you’d like to see represented in Maple Learn’s document gallery? Please let us know in the comments below!

    Hello MaplePrimes community,

    We just created a Frequently Asked Question article that may address some Primes questions about updates to Physics in Maple 2022.2: 

    Why does Maple 2022.2 throw an error executing Physics:-Version(latest)?

    For searchability, the specific error in question is

    Error, (in Physics:-Version) unable to determine the Physics Updates version, could you please report the problem to support [at] maplesoft [dot] com

     

    • Maplesoft will work to improve package updating in future versions of Maple.
    • In Maple 2022.2, the workaround is to install and/or update the Maplesoft Physics Updates using the MapleCloud toolbar.

    An example of uniform motion along a generalized coordinate using the Draghilev method. (This post was inspired by school example in one of the forums.)
    The equations used in the program are very simple and, I think, do not require any special comments. DM is a procedure that implements the Draghilev method with "partial parameterization".

    DM_V.mw

    When K = 1, parameterization is carried out by changing the angle of rotation of the wheel. That is, uniform rolling is carried out.

    For K = 4, the coordinate corresponding to the position of the slider is parametrized.

     

    When K = 6, the slider moves with acceleration, according to a given equation. Hence, we have carried out the parameterization with respect to “time”.



    With the help of such techniques, we can obtain the calculation of the kinematics of both lever mechanisms and various types of manipulators.

    First 14 15 16 17 18 19 20 Last Page 16 of 304