Mathematica 教學 Non-uniform circle packing

circlepacking[ntemp_Integer]:=
Module[{k=ntemp,c,dist,cons,radi,vars,obj,cord,mypoint,testpoint,x,y,r,myrange},

(*定義圓心*)
c[i_]:={x[i],y[i]}/.s_[u_]:>ToExpression[ToString[s]<>ToString[u]];

(*定義任兩圓限制式*)
dist[i_Integer,j_Integer]:={Sqrt[(c[i]-c[j]).(c[i]-c[j])]>=1/Sqrt[i]+1/Sqrt[j]};
radi[i_Integer]:={Sqrt[c[i].c[i]]+1/Sqrt[i]
(*定義最佳解座標*)
cord[n_Integer]:={x[#],y[#],1/Sqrt[#]}&/@Range[n]/.s_[u_]:>ToExpression[ToString[s]<>ToString[u]];

(*定義問題的限制式*)
cons[n_]:=Flatten[{Apply[dist,Subsets[Range[n],{2}],1],radi[#]&/@Range[n]}];

(*定義變數*)
vars[n_]:=Flatten@{r,{x[#],y[#]}&/@Range[n]}/.s_[u_]:>ToExpression[ToString[s]<>ToString[u]];

(*變數最大變動區間*)
myrange=Total[1/Sqrt[#]&/@Range[k]]//N;

(*最佳化主程式*)
obj[n_Integer]:=NMinimize[{r,cons[n]},vars[n],
Method->{"DifferentialEvolution",
"InitialPoints"->RandomReal[{-myrange,myrange},{20n+20,2n+1}],
"PenaltyFunction"->Function[{d,i},2^i*d^2]}];

mypoint[n_]:=({#[[1]],cord[n]/.#[[2]]}&/@{obj[n]})[[1]];

(*求解k的圓的最佳解*)
testpoint=mypoint[k];

(*畫出圖形*)
Print@Show[Graphics[Circle[{0,0},testpoint[[1]]]],
Graphics[{RGBColor[Random[],Random[],Random[]],
Disk[{#[[1]],#[[2]]},#[[3]]],Black,PointSize[0.0125],
Point[{#[[1]],#[[2]]}]}]&/@testpoint[[2]]];
testpoint
]

circlepacking[20]

Tags: Mathematica 教學, Teaching

• 編輯文章
• 此文章網址
• 轉寄此文章
• 目前有0 意見