| [ QuizWit ] in KIDS 글 쓴 이(By): parsec ( 먼 소 류 ) 날 짜 (Date): 2002년 6월 14일 금요일 오후 07시 12분 43초 제 목(Title): 양치기 문제 -- Numerical Approach 양치기 문제를 풀기 위한 Mathematica Program입니다. 우선 필요한 함수와 연산을 정의합니다. ---------------------------------------------------------------------- part3[l_,i_]:=Take[RotateLeft[l,i-1],3] (* 삼각형 추출 *) part4[l_,i_]:=Take[RotateLeft[l,i-1],4] (* 사각형 추출 *) subs3[l_List,p_List,i_]:=Module[{q}, (* 삼각형 치환 *) q=RotateLeft[l,i-1]; q=Drop[q,3]; q=Insert[q,p,1]; q=Delete[q,{1,0}]; q=RotateRight[q,i-1]; Return[q]]/;Dimensions[p]=={3,2} subs4[l_List,p_List,i_]:=Module[{q}, (* 사각형 치환 *) q=RotateLeft[l,i-1]; q=Drop[q,4]; q=Insert[q,p,1]; q=Delete[q,{1,0}]; q=RotateRight[q,i-1]; Return[q]]/;Dimensions[p]=={4,2} rot[{x_,y_},th_]:={{Cos[th],-Sin[th]},{Sin[th],Cos[th]}}.{x,y} (* 벡터 회전 *) dist[v1_,v2_]:=Sqrt[(v2-v1).(v2-v1)] (* 두 점 사이의 거리 *) angle[{x_,y_}]:=Im[N[Log[x+I y]]] (* 벡터의 회전각 *) baseangle[p_]:=angle[p[[3]]-p[[1]]] /;Dimensions[p]=={3,2} (* 밑변의 기운 각 -- 삼각형 *) baseangle[p_]:=angle[p[[4]]-p[[1]]] /;Dimensions[p]=={4,2} (* 밑변의 기운 각 -- 사각형 *) basecenter[p_]:=(p[[3]]+p[[1]])/2 /; Dimensions[p]=={3,2} (* 밑변의 중심 -- 삼각형 *) (* 주어진 삼각형을 변의 길이의 합이 같은 이등변 삼각형으로 변환 -- 타원의 공식을 이용*) equilaterize3[p_]:= Module[{q,v0,th,a,b,c},q=p;v0=basecenter[q];th=baseangle[q]; q=rot[#,-th]& /@ (Map[(-v0+#)&,q]); a=(dist[q[[2]],q[[1]]]+dist[q[[2]],q[[3]]])/2; c=dist[q[[3]],q[[1]]]/2; b=Sqrt[a^2-c^2]; q=ReplacePart[q,{0,b},2]; q=(v0+#)& /@ (rot[#,th]& /@ q); Return[q] ] (* 2차원 벡터의 외적 *) cross2D[v1_,v2_]:=Apply[Cross,Map[Append[#,0]&, {v1,v2}]][[3]] (* 다각형의 넓이 *) area[l_List]:= Abs[Apply[Plus,Apply[cross2D,Transpose[{l,RotateLeft[l,1]}],{1}]]]/2 (* 다각형의 둘레 길이 *) peri[l_List]:=Apply[Plus,Apply[dist,Transpose[{l,RotateLeft[l,1]}],{1}]] (* Poly line의 길이의 합 *) extent[l_List]:= Apply[Plus,Apply[dist,Drop[Transpose[{l,RotateLeft[l,1]}],-1],{1}]] (* 다각형에서, 주어진 인덱스 사이의 점들을 움직여 이웃한 두 변의 길이를 길이의 합이 변하지 않고, 길이가 같도록 변환하는 작업을 반복*) Equilaterize[l_,n_,a0_,a1_]:= Module[{actv,l0,len,i,p}, l0=l; len=Length[l0]; actv=False; For[i=1 ,i<n * len,i++, idx=Mod[i,len]+1; If[idx==a0,actv=True]; If[actv==False,Continue[]]; p=part3[l0,i]; p=equilaterize3[p]; l0=subs3[l0,p,i]; If[idx==a1,actv=False]; ]; l0 ] (* 사각형을, 각 변의 길이를 고정시킨 채 원에 내접하는 사각형으로 변환. -- 내접 사각형의 대각선 공식을 이용 *) cyclify[p_List]:= Module[{p4,v0,th0,q2,a,b,c,d,psol,qsol,pval,qval,th1,th2,tx1,tx2}, p4=p; v0=p4[[1]]; th0=baseangle[p4]; q4=rot[#,-th0]& /@ ((-v0+#)& /@ p4); a=dist[q4[[4]],q4[[1]]]; b=dist[q4[[2]],q4[[1]]]; c=dist[q4[[3]],q4[[2]]]; d=dist[q4[[4]],q4[[3]]]; pval=Sqrt[(a b+c d)(a c+b d)/(a d+b c)]; qval=Sqrt[(a c+b d) (a d+b c)/(a b+c d)]; psol=FindRoot[Sqrt[(a+d Cos[th2])^2+(d Sin[th2])^2]==pval,{th2,1.55}]; qsol=FindRoot[Sqrt[(a-b Cos[th1])^2+(b Sin[th1])^2]==qval,{th1,1.55}]; tx1=th1/.qsol; tx2=th2/.psol; If[Sin[tx1]<0,tx1*=-1;tx2*=-1]; q4=ReplacePart[ q4,b{Cos[tx1],Sin[tx1]},2]; q4=ReplacePart[q4, {a,0}+d{Cos[tx2],Sin[tx2]},3]; p4=(v0+#)& /@ (rot[#,th0]& /@ q4); Return[p4] ]/;Dimensions[p]=={4,2} (* 위 연산을 이용하여 다각형 상의 이웃한 3개의 변들을 내접 사각형으로 변환하는 작업을 반복함.*) MaximizeArea[p_,n_,i0_,i1_]:= Module[{i,len,p2,p4,idx1,idx2,actv}, actv=False; p2=p; len=Length[p2]; For[i=1,i<n len,i++, idx1=Mod[i,len]+1; idx2=Mod[i+1,len]+1; If[idx1==i0,actv=True]; If[actv==False,Continue[]]; p4=part4[p2,i]; p4=cyclify[p4]; p2=subs4[p2,p4,i]; If[idx2==i1,actv=False]; ]; Return[p2] ] ClosedLine[l_List]:=Line[Append[l,l[[1]]]] -------------------- 이상은 함수 정의 ----------------- 위에서 정의한 연산들을 이용해서 원래의 양치기 문제를 풀면 poles=3 wirelen=100 (* 먼저 첫번째 양치기의 울타리 *) (* 둘레 길이의 합이 wirelen(=100)이 되도록 리스트를 초기화 *) sh1=Table[{(i*wirelen/2)/(poles-2+1),0},{i,0,poles-1}]; 등변 다각형으로 만들면: sh1=Equilaterize[sh1,100,1,poles] (* 결과입니다 *) {{40.1286,-35.1611},{20.7529,-8.03745},{53.9305,-4.81947}} 면적: area[sh1] ----결과------- 481.125 둘레: peri[sh1] ----결과--- 100. (* 일반화할 경우에 대해서 등변다각형이라 해도 최대 면적을 갖지 않을 수 있으므로 *) If[poles>=4,sh1=MaximizeArea[sh1,100,1,poles]] (* 그림으로 보려면 ... *) Show[Graphics[{RGBColor[1,0,0],ClosedLine[sh1]}],AspectRatio->Automatic, Frame->True] --> http://community.webshots.com/storage/1/v1/8/81/21/40688121jpYQZQ_ph.jpg (* 두번째 양치기의 울타리 *) (* 첫번째 양치기의 말뚝 중에서 사용할 말뚝은 먼저 첫번째와, 다음 pick1번째 *) pick1=2 (* 두번째 양치기의 울타리를 초기화한다. 첫번째 양치기의 울타리에서 이용할 말뚝 두개의 방향으로 울타리 길이가 100이 되도록 말뚝을 배열한다. *) v0=sh1[[pick1]]-sh1[[1]] -->{-19.3757,27.1237} k0=Sqrt[v0.v0] --> 33.3333 k1=(wirelen+k0)/2/poles --> 22.2222 v1=k1 v0/k0 --> {-12.9171,18.0825} sh2=(sh1[[1]]+#)& /@ Table[i v1,{i,1,poles}]; (* 울타리를 정의하게 될 말뚝을 모두 포함시킨다: 면적 계산시 이용 *) sh2=RotateRight[Delete[Append[sh2,Reverse[Take[sh1,pick1]]],{-1,0}],1] --> {{40.1286,-35.1611},{27.2115,-17.0787},{14.2943,1.00378},{1.37723,19.0862},{ 20.7529,-8.03745}} (* 초기화한 모습 *) Show[Graphics[{RGBColor[1,0,0],ClosedLine[sh1]}], Graphics[{RGBColor[0,0,1],Line[Take[sh2,{1,poles+2}]]}], AspectRatio->Automatic,Frame->True] http://community.webshots.com/storage/1/v2/8/81/25/40688125sOkxji_ph.jpg (* 등변화... *) sh2=Map[Re[#]&,Equilaterize[sh2,200,2,poles+1],{2}] Show[Graphics[{RGBColor[1,0,0],ClosedLine[sh1]}], Graphics[{RGBColor[0,0,1],Line[Take[sh2,{1,poles+2}]]}], AspectRatio->Automatic,Frame->True] http://community.webshots.com/storage/1/v1/8/81/22/40688122TrQWQM_ph.jpg (* 면적 증가 변환...*) sh2=MaximizeArea[sh2,100,2,poles+1] Show[Graphics[{RGBColor[1,0,0],ClosedLine[sh1]}], Graphics[{RGBColor[0,0,1],Line[Take[sh2,{1,poles+2}]]}], AspectRatio->Automatic,Frame->True] http://community.webshots.com/storage/1/v1/8/81/26/40688126RdHKqH_ph.jpg area[sh2] --> 1212.32 extent[Take[sh2,{1,poles+2}]] --> 100. (* 세번째 양치기의 울타리 *) (* 먼저 빌려 쓸 말뚝을 고른다. 첫번째 양치기의 말뚝 중 pick2번째, 두번째 양치기의 말뚝 중 pick3번째의 것을 사용하기로 한다. *) pick2=poles; pick3=2; (* 나머지 과정은 두번쩨 양치기의 경우와 비슷하다. 면적계산을 위해 울타리를 정의하게 될 말뚝들을 모두 리스트에 포함시키는 것을 잊지 말자...*) v2=sh2[[pick3]]-sh1[[pick2]] k2=Sqrt[v2.v2] k3=(wirelen+k2)/2/poles v3=k3 v2/k2 sh3=(sh1[[pick2]]+#)& /@ Table[i v3,{i,1,poles}]; sh3=Delete[Append[sh3,Reverse[Take[sh2,{2,pick3}]]],{-1,0}]; sh3=Delete[ Append[sh3,Reverse[Take[RotateLeft[sh1,pick2-1],poles-pick2+2]]],{-1,0}]; sh3=RotateRight[sh3,1] (* 초기화한 모습 *) Show[Graphics[{RGBColor[1,0,0],ClosedLine[sh1]}], Graphics[{RGBColor[0,0,1],Line[Take[sh2,{1,poles+2}]]}], Graphics[{RGBColor[0,1,0],Line[Take[sh3,{1,poles+2}]]}], AspectRatio->Automatic,Frame->True] http://community.webshots.com/storage/1/v2/8/81/32/40688132AAImoJ_ph.jpg sh3=Map[Re[#]&,Equilaterize[sh3,100,2,poles+1],{2}] (* 등변화한 모습 *) Show[Graphics[{RGBColor[1,0,0],ClosedLine[sh1]}], Graphics[{RGBColor[0,0,1],Line[Take[sh2,{1,poles+2}]]}], Graphics[{RGBColor[0,1,0],Line[Take[sh3,{1,poles+2}]]}], AspectRatio->Automatic,Frame->True] http://community.webshots.com/storage/1/v2/8/81/31/40688131wFYPyc_ph.jpg (* 면적 증가 연산...*) sh3=MaximizeArea[sh3,100,2,poles+1] Show[Graphics[{RGBColor[1,0,0],ClosedLine[sh1]}], Graphics[{RGBColor[0,0,1],Line[Take[sh2,{1,poles+2}]]}], Graphics[{RGBColor[0,1,0],Line[Take[sh3,{1,poles+2}]]}], AspectRatio->Automatic,Frame->True] http://community.webshots.com/storage/1/v1/8/81/34/40688134QgObRP_ph.jpg area[sh3] --> 1341.42 (* 말뚝의 갯수와 다른 울타리의 말뚝 중 이용할 말뚝의 인덱스를 변화시키고, 등변화 및 면적 증가 연산의 반복횟수를 적당히 늘려가며 시도하면 다른 일반적인 경우에 대해서도 근사적인 해를 구할 수 있읍니다. 다음 페이지를 보시면 말뚝 3개인 경우와 10개인 경우에 대한 그림을 볼 수 있습니다.*) http://community.webshots.com/album/40688066sbKcFo 참고로 말뚝이 10개인 경우: area[sh1] --> 769.428 area[sh2] --> 1050.63 (pick1-->5) area[sh3] --> 1322.67 (pick2-->9, pick3-->4) ◇ ~~~_ _ ∴ ~|~| | _/__, SEP. 11. 2001 _ ∴∴ _ ~ | | ＼ ` Armorica under a tat ,-| `,-,_| |__ | | | A ______|_|__|_|___|__|| | |__|_|_____________________________________ |