Maple Questions and Posts

These are Posts and Questions associated with the product, Maple

Hello,

I use Maple 2022 on a MacBook Pro. In most of the plots I create I use symbol=solidcircle and symbolsize=12.

Is there a way to make those (and perhaps oher) settings user default so I don't need them in every plot command?

Thanks.

Jose


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


 

 

After studying the plottools:-transform command, I intend to visualize the following regions with constrained parameters in 
 

(plottools[transform](proc (u, v) options operator, arrow; [u^3-v^2, u^2-v^3] end proc))(plots[inequal](`or`(u^2+4*v^2 <= 4, `and`(u^2+v^2 < 4, 4*v >= (u+2)^2+2*v^2)), nolines))

 

(plottools[transform](proc (s, t) options operator, arrow; [s^2*sqrt(t)*cos(t), s^2*sin(t)] end proc))(plots[inequal](`and`(`and`(s >= 1, 5*s <= 5+t), t < 5), s = 1 .. 2, t = 0 .. 5))

 

 

But Mma gives 

The first instance (with default settings) is the same, but as for the second instance, which graph is correct? 

restart;
with(plottools):
with(plots):
transform((u, v) -> [u^3 - v^2, u^2 - v^3])(inequal(Or(u^2 + 4*v^2 <= 4, And(u^2 + v^2 < 4, (u + 2)^2 + 2*(v - 1)^2 <= 2)), nolines));
transform((s, t) -> [s^2*sqrt(t)*cos(t), s^2*sin(t)])(inequal(`and`(1 <= s, 5*s <= 5 + t, t < 5), s = 1 .. 2, t = 0 .. 5));


Download TransformedRegion.mws

I take the liberty to rephrase my previous question as I believe the title was not very clear and so maybe some power users did not look at it. I am making the transition from Mathcad towards Maple and get stuck solving the equation in the attached worksheet. In mathcad I would solve it like this:

How can I achieve results in Maple? I know it is a very powerful program but for me the learning curve is at this moment quite steep. Any help would be very much appreciated.

Multiple_input.mw

plot-problem.mw

I have done something but what?

I want to use maple notation for input and  output.

I have done something to mess this up.

How do I get rid of the typesetting messages?

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

L:=[["O",3.85090000,0.45160000,0.00120000],
["O",−2.59990000,1.40410000,−0.00180000],
["N",−1.57050000,−0.71710000,0.00010000],
["C",−0.20660000,−0.42310000,−0.00020000],
["C",0.22050000,0.90470000,0.00040000],
["C",0.72980000,−1.45700000,−0.00070000],
["C",1.58410000,1.19860000,0.00020000],
["C",2.09330000,−1.16290000,−0.00070000],
["C",2.52040000,0.16480000,−0.00030000],
["C",−2.64850000,0.17820000,0.00090000],
["C",−3.97350000,−0.54200000,0.00100000],
["H",−0.44360000,1.75770000,0.00120000],
["H",0.41130000,−2.49630000,−0.00100000],
["H",−1.80100000,−1.70860000,0.00010000],
["H",1.90530000,2.23700000,0.00090000],
["H",2.81800000,−1.97260000,−0.00080000],
["H",−4.06550000,−1.14630000,−0.90580000],
["H",−4.79040000,0.18440000,0.02880000],
["H",−4.04450000,−1.18860000,0.88020000],
["H",3.96500000,1.41760000,0.00170000]];

S:={{1,9},{1,20},{2,10},{3,4},{3,10},{3,14},{4,5},{4,6},
{5,7},{5,12},{6,8},{6,13},{7,9},{7,15},{8,9},{8,16},{10,11},
{11,17},{11,18},{11,19}}

Now I have list L above a Set S below 

Now we delete from L all having "H" in first position and also  the subsets of the set S which as all contain those  list position number where L had "H"

Kind help with function which can

Takes input of the list and sets returns the new list and set where the above are done 

First of all I would like to wish all of you a happy, prosperous but especially healthy 2023! I have again a beginner question. Why is test2 not working in the attached document?

Thank you so much for your assistance!

QuestionFDS.mw

I remember a section “Tell us what we can do better” at the bottom of online help pages. I used this whenever I came across a potential error worth investigating.

Has this section disappeared (I hope not) or do I have a browser issue?

restart;
with(geometry);
with(plots);
_EnvHorizontalName = 'x';
_EnvVerticalName = 'y';
Vdot := proc(U, V) local i; add(U[i]*V[i], i = 1 .. 2); end proc;
dist := proc(M, N) sqrt(Vdot(M - N, M - N)); end proc;
EQ := proc(M, N) local eq, a, b, c; eq := simplify(expand((y - M[2])/(x - M[1]) - (N[2] - M[2])/(N[1] - M[1]))*(x - P1[1])*(P2[1] - P1[1])); a := coeff(eq, x); b := coeff(eq, y); c := tcoeff(eq, [x, y]); RETURN(-a*x/c - b*y/c - 1); end proc;
R := 5;
ang := [2/3*Pi, -3*Pi*1/4, -Pi*1/6];
seq(point(`||`(P, i), [R*cos(ang[i]), R*sin(ang[i])]), i = 1 .. 3);
seq(dsegment(`||`(seg, i), [`||`(P, i), `||`(P, irem(i, 3) + 1)]), i = 1 .. 3);
circle(cir, [point(OO, [0, 0]), R]);
sol := solve(subs(x = 2, Equation(cir, [x, y])), y);
point(A, [2, sol[1]]);
triangle(Tri, [P1, P2, P3]);
incircle(inc, Tri, 'centername' = oo);
circle(Cr, [A, oo]);
sol := solve({Equation(Cr, [x, y]), Equation(inc, [x, y])}, {x, y});
point(H1, [subs(sol, x), subs(sol, y)]);
line(L, [A, oo]);
reflection(H2, H1, L);
line(L1, [A, H1]);
line(L2, [A, H2]);
Equation(cir, [x, y]);
Equation(L1, [x, y]);
sol := solve({Equation(L1, [x, y]), Equation(cir, [x, y])}, {x, y});
evalf(%);
point(M1, [subs(sol, x), subs(sol, y)]);
sol2 := solve({Equation(L2, [x, y]), Equation(cir, [x, y])}, {x, y});
evalf(%);
point(M2, [subs(sol2, x), subs(sol2, y)]);
triangle(TR, [M1, M2, A]);
display([draw([P1(symbol = solidcircle, symbolsize = 8, color = blue), P2(symbol = solidcircle, symbolsize = 8, color = blue), P3(symbol = solidcircle, symbolsize = 8, color = blue), A(symbol = solidcircle, symbolsize = 8, color = black), H1(symbol = solidcircle, symbolsize = 8, color = black), H2(symbol = solidcircle, symbolsize = 8, color = black), L1(color = black), L2(color = black), seg1(color = magenta), seg2(color = magenta), seg3(color = magenta), Cr(color = black), cir(color = magenta), inc(color = blue)]), textplot([seq([coordinates(`||`(P, i))[], convert(`||`(P, i), string)], i = 1 .. 3)], 'align' = {'above', 'left'})], view = [-6 .. 10, -15 .. 6], scaling = constrained, size = [800, 800], axes = none);
It seems that there is conusion between M1 and M2. How to write letters A, M1 ,M2, H1, H2 ? Thank you.

The range is wrong. For details, see below, please.
 

restart;

assume(x, RealRange(0, 1))

plot([sqrt(x*(2 - x)/3), 1 - sqrt((1 - x^2)/3)], legend = InertForm:-Display~([sqrt(x*(2 - x) %/ 3), 1 - sqrt((1 %- x^2) %/ 3)], 'inert' = false));

 

smartplot([sqrt(x*(2 - x)/3), 1 - sqrt((1 - x^2)/3)]);

 

smartplot([''piecewise'(And(0 <= x, x <= 1), sqrt(x*(2 - x)/3))', ''piecewise'(And(0 <= x, x <= 1), 1 - sqrt((1 - x^2)/3))']);

 

smartplot(['piecewise(And(0 <= x, x <= 1), sqrt(x*(2 - x)/3))', 'piecewise(And(0 <= x, x <= 1), 1 - sqrt((1 - x^2)/3))']);

 

smartplot([''piecewise''(And(0 <= x, x <= 1), sqrt(x*(2 - x)/3)), ''piecewise''(And(0 <= x, x <= 1), 1 - sqrt((1 - x^2)/3))]);

 

x := 'x'NULL


 

Download SmartPlots.mw

The help page claims that smartplot(..) will call 2-D plot procedures ultimately, but why is the smartplots command incompatible with the use of assume?

Good day everyone,

I have been having problems with a system of PDE solution using `numeric`.

It's giving me the error code "Error, (in pdsolve/numeric/par_hyp) input system is too far from a 'standard' form (see ?pdsolve,numeric for more detail)" and I have checked to the best of my ability for the error but could not see anything.

The code is attached below.

Please, anyone with useful information should help. Thanks

Numeric.mw

How do I find the solutions "left" and "right" with only answers in the range 0 to +1?

The domain of eq1 and eq2 is 0 <=beta <= 1.

So I'm only interested in left[1] and right[2]:

Here is the code:

eq[1]:=(1/6)*beta^3+(1/2)*xi^2*beta-(1/2)*beta^2-(1/2)*xi^2+(1/3)*beta = 0;

eq[2]:=(1/6)*beta*(beta^2+3*xi^2-6*xi+2) = 0;

left:=solve(eq[1],xi);

right:=solve(eq[2],xi);

left:=plot(left[1],beta=0..1,thickness=2,color=blue);

right:=plot(right[2],beta=0..1,thickness=2,color=blue);

How do I plot a volume of revolution? I can plot other volumes using Student[Calculus1]]:-VolumeOfRevolution, but not this one. I get a blank plot. I do get the correct volume from output=value..

How do I plot this using plot3d?

restart;
a := 0; b := 1;
f := (x) -> x^2+2;
g := (x) -> 1/2*x+1;
V := int(f(x)^2 - g(x)^2,x=a..b)*Pi;
Student[Calculus1]:-VolumeOfRevolution(f(x),g(x),x=a..b,output=value);
Student[Calculus1]:-VolumeOfRevolution(f(x),g(x),x=a..b,output=plot);
 

I cannot plot the **erf** function, please see below, what's wrong ?

 

First 7 8 9 10 11 12 13 Last Page 9 of 1984