Featured Post

21011

In the excellent book by W.G. Chinn, N.E. Steenrod "First Concepts of Topology" the theorem is proved which states that any bounded planar region can be cut into 4 regions of equal area by 2 perpendicular cuts (the pancake problem). This is an existence theorem which does not provide any way to find these cuts. In this post I made an attempt to find such cuts for any convex region on the plane bounded by a piecewise smooth self-non-intersecting curve.
The Into_4_Equal_Areas procedure returns a list of coordinates of 5 points: the first 4 points are the endpoints of the cutting segments, the fifth point is the intersection point of these segments. This procedure significantly uses my old procedure Area , which can be found in detail at the link  https://mapleprimes.com/posts/145922-Perimeter-Area-And-Visualization-Of-A-Plane-Figure-  . The formal argument of the Into_4_Equal_Areas procedure is a list  L specifying the boundary of the region to be cut. When specifying  L, the boundary can be passed clockwise or counterclockwise, but it is necessary that the parameter t (when specifying each link) should go in ascending order. This can always be achieved by replacing  t  with  -t  if necessary. The Pic procedure draws a picture of the source region and cutting segments. For ease of use, the code for the  Area  procedure is also provided. It is also worth noting that the procedure also works for "not too" non-convex regions (see examples below).

restart;
Area := proc(L) 
local i, var, e, e1, e2, P; 
for i to nops(L) do 
if type(L[i], listlist(algebraic)) then 
P[i] := (1/2)*add(L[i, j, 1]*L[i, j+1, 2]-L[i, j, 2]*L[i, j+1, 1], j = 1 .. nops(L[i])-1) else 
var := lhs(L[i, 2]); 
if type(L[i, 1], algebraic) then e := L[i, 1]; 
if nops(L[i]) = 3 then P[i] := (1/2)*(int(e^2, L[i, 2])) else 
if var = y then P[i] := (1/2)*simplify(int(e-var*(diff(e, var)), L[i, 2])) else 
P[i] := (1/2)*simplify(int(var*(diff(e, var))-e, L[i, 2])) end if end if else e1 := L[i, 1, 1]; e2 := L[i, 1, 2]; 
P[i] := (1/2)*simplify(int(e1*(diff(e2, var))-e2*(diff(e1, var)), L[i, 2])) end if end if end do; 
abs(add(P[i], i = 1 .. nops(L))); 
end proc:

Into_4_Equal_Areas:=proc(L::list)
local D, n, c, L1, L2, L3, f, L0, i, j, k, m, A, B, C, P, S, sol, Sol;
f:=(X,Y)->expand((y-X[2])*(Y[1]-X[1])-(x-X[1])*(Y[2]-X[2]));
L0:=map(p->`if`(type(p,listlist),[[p[1,1]+t*(p[2]-p[1])[1],p[1,2]+t*(p[2]-p[1])[2]],t=0..1],p), L);
S:=Area(L); c:=0;
n:=nops(L);
for i from 1 to n do
for j from i to n do
for k from j to n do
for m from k to n do
if not ((nops({i,j,k})=1 and type(L[i],listlist)) or (nops({j,k,m})=1 and type(L[j],listlist)))then
A:=convert(subs(t=t1,L0[i,1]),Vector): 
B:=convert(subs(t=t2,L0[j,1]),Vector):
C:=convert(subs(t=t3,L0[k,1]),Vector): 
D:=convert(subs(t=t4,L0[m,1]),Vector):
P:=eval([x,y], solve({f(A,C),f(B,D)},{x,y})):
L1:=`if`(j=i,[subsop([2,2]=t1..t2,L0[i]),[convert(B,list),P],[P,convert(A,list)]],`if`(j=i+1,[subsop([2,2]=t1..op([2,2,2],L0[i]),L0[i]),subsop([2,2]=op([2,2,1],L0[j])..t2,L0[j]),[convert(B,list),P],[P,convert(A,list)]], [subsop([2,2]=t1..op([2,2,2],L0[i]),L0[i]),L0[i..j][],subsop([2,2]=op([2,2,1],L0[j])..t2,L0[j]),[convert(B,list),P],[P,convert(A,list)]])):
L2:=`if`(k=j,[subsop([2,2]=t2..t3,L0[j]),[convert(C,list),P],[P,convert(B,list)]],`if`(k=j+1,[subsop([2,2]=t2..op([2,2,2],L0[j]),L0[j]),subsop([2,2]=op([2,2,1],L0[k])..t3,L0[k]),[convert(C,list),P],[P,convert(B,list)]], [subsop([2,2]=t2..op([2,2,2],L0[j]),L0[j]),L0[j..k][],subsop([2,2]=op([2,2,1],L0[k])..t3,L0[k]),[convert(C,list),P],[P,convert(B,list)]])):
L3:=`if`(m=k,[subsop([2,2]=t3..t4,L0[k]),[convert(D,list),P],[P,convert(C,list)]],`if`(m=k+1,[subsop([2,2]=t3..op([2,2,2],L0[k]),L0[k]),subsop([2,2]=op([2,2,1],L0[m])..t4,L0[m]),[convert(D,list),P],[P,convert(C,list)]], [subsop([2,2]=t3..op([2,2,2],L0[k]),L0[k]),L0[k..m][],subsop([2,2]=op([2,2,1],L0[m])..t4,L0[m]),[convert(D,list),P],[P,convert(C,list)]])):
sol:=fsolve({Area(L1)-S/4,Area(L2)-S/4,Area(L3)-S/4,(D-B).(C-A)},{t1=op([2,2],L0[i]),t2=op([2,2],L0[j]),t3=op([2,2],L0[k]),t4=op([2,2],L0[m])}) assuming real:
if type(sol,set(`=`)) then c:=c+1; Sol[c]:=convert~(eval([A,B,C,D,P],sol),list)
 fi; fi;
od: od: od: od:
convert(Sol,list);
end proc:

Pic:=proc(L,Sol)
local P1, P2, T;
uses plots, plottools;
P1:=seq(`if`(type(s,listlist),line(s[],color=blue, thickness=2),plot([s[1][],s[2]],color=blue, thickness=2)),s=L):
P2:=line(Sol[1],Sol[3],color=red, thickness=2), line(Sol[2],Sol[4],color=red):
T:=textplot([[Sol[1][],"A"],[Sol[2][],"B"],[Sol[3][],"C"],[Sol[4][],"D"],[Sol[5][],"P"]], font=[times,18], align=[left,above]);
display(P1,P2,T, scaling=constrained, size=[800,500], axes=none);
end proc: 


Examples of use:

L:=[[[0,0],[1,4]],[[1,4],[6,7]],[[6,7],[12,0]],[[12,0],[0,0]]]:
Sol:=Into_4_Equal_Areas(L);
Pic(L, Sol[1]);

# Check (areas of all 4 regions)
Area([[L[1,1],Sol[1,4],Sol[1,5],Sol[1,1],L[1,1]]]),
Area([[Sol[1,4],Sol[1,5],Sol[1,3],L[4,1],Sol[1,4]]]),
Area([[Sol[1,5],Sol[1,2],L[3,1],Sol[1,3],Sol[1,5]]]),
Area([[Sol[1,5],Sol[1,2],L[1,2],Sol[1,1],Sol[1,5]]]); 

        


 

L:=[[[1+cos(-t),1+sin(-t)],t=-3*Pi/2..-Pi],[[0,1],[-1,0]],[[cos(t),sin(t)],t=Pi..2*Pi]]:
Sol:=Into_4_Equal_Areas(L);
Pic(L,Sol[1]);

    

# The boundary is the Archimedes spiral and the arc of a circle

L:=[[[t*cos(t),t*sin(t)],t=0..2*Pi],[[Pi+5*cos(-t),sqrt(25-Pi^2)+5*sin(-t)],t=arccos(Pi/5)..Pi-arccos(Pi/5)]]:
Sol:=evalf(Into_4_Equal_Areas(L));
Pic(L,Sol[1]);

     

More examples can be found in the attached file.

4_Equal_Area.mw

Featured Post

I am very happy to announce the first public release of a project which I have been working on for the last couple of years.

NODEMaple consists of a set of Maple workbooks and a library for structural design based on the Eurocode.

Currently the main development of the workbooks is focused on "Eurocode 5: Design of timber structures" with the Norwegian Annex.

This software has been made public in the hope of that it might be useful for other structural designers, professionals as well as students. Everyone interested is very Welcome to contribute to this project. The code is published under the GPLv3 license.

For more information see https://github.com/Anthrazit68/NODEMaple.



A task for geometry fans

Maple asked by Kitonum 21011 January 11

Odegenerator plots

Maple asked by janhardo 430 Today