自動解答系の作り方(18)

 今回は関数のグラフに限らず,一般の「曲線に囲まれた図形の面積」を定義し,それを求めるMathematicaのコマンドコードを公開します.

 R^2の部分集合からRへの関数族 {f_{k}}_{k∈{1,...,n}} に対して,∃k.( f_{k}(x,y) = 0 ∧ k∈{1,...,n} ) を満さない(つまり,どの曲線も通らない)点(x,y)全体をAとし,Aの連結成分のうち面積確定である(有界でなくてよい)ものの面積の総和を「{f_{k}(x,y) = 0}_{k∈{1,...,n}}に囲まれた図形の面積」と定めます.

 実装のアイデアは,各関数が定める正・負領域の交わり(空を含めて2^n個あり,それらの和集合がAとなる)を表す式をTuplesとInnerで用意,それらの個々の連結成分を表す式をReduce(CylindricalDecomposition+α)及び適当なRuleによる開集合の分解により生成し,それぞれのBooleをR^2でIntegrateして,Realsに属するものをSelect,Plusするというものです.

 以下,コードと使用例.

(* 分解規則 *)
rules=#//.{Or[P___,(A:Inequality[_,Less,_,Less,_]&&_),(B:Inequality[_,Less,_,Less,_]&&_),Q___]->Sequence[P||A,B||Q]}//.{Or[P___,(A:Inequality[_,Less,_,Less,_]&&_),(B:Inequality[_,Less,_,LessEqual,_]&&_),Q___]->Sequence[P||A,B||Q]}//.{Or[P___,(A:Inequality[_,LessEqual,_,Less,_]&&_),(B:Inequality[_,Less,_,Less,_]&&_),Q___]->Sequence[P||A,B||Q]}//.{Or[P___,(A:Inequality[_,LessEqual,_,Less,_]&&_),(B:Inequality[_,Less,_,LessEqual,_]&&_),Q___]->Sequence[P||A,B||Q]}//.{Or[(A:_<_&&_),(B:Inequality[_,Less,_,Less,_]&&_),B1___]->Sequence[A,B||B1]}//.{Or[(A:_<_&&_),(B:Inequality[_,Less,_,LessEqual,_]&&_),B1___]->Sequence[A,B||B1]}//.{Or[A1___,(A:Inequality[_,Less,_,Less,_]&&_),B:_>_&&_]->Sequence[A1||A,B]}//.{Or[A1___,(A:Inequality[_,LessEqual,_,Less,_]&&_),B:_>_&&_]->Sequence[A1||A,B]}//.{Or[(A:_<_&&_),B:_>_&&_]->Sequence[A,B]}&;

(* コマンド本体 *)
areaCAD0[fxys_List,cond_:True,x_Symbol:x,y_Symbol:y]:=Module[{pList=fxys/.L_==R_->L-R/.Sqrt[A_]->Piecewise[{{Sqrt[A],Reduce[0<=A]}}]},(Reduce[#1,{x,y},Reals]&)/@(Inner[#1*#2>0&,pList,#1,And]&)/@Tuples[{1,-1},Length[pList]]];areaCAD1[fxys_List,cond_:True,x_Symbol:x,y_Symbol:y]:=rules@areaCAD0[fxys,cond,x,y];areaCAD2[fxys_List,cond_:True,x_Symbol:x,y_Symbol:y]:=(Integrate[Boole[#1&&cond],{x,-Infinity,Infinity},{y,-Infinity,Infinity}]&)/@areaCAD1[fxys,cond,x,y];areaCAD[fxys_List,cond_:True,x_Symbol:x,y_Symbol:y]:=FullSimplify[Plus@@Select[areaCAD2[fxys,cond,x,y],Element[#,Reals]&]];areaCAD[fxys_List,x_Symbol:x,y_Symbol:y]:=areaCAD[fxys,True,x,y];

(* 曲線のリストの例のリスト *)
l={{y==-Sqrt[1-x^2],x==1,y==1+x},{y==0,y==x^2,y==(-1+x)^2},{y==0,y==(-1+x^2)^2},{y==0,y==x,y==2-x},{y==-1-2x,y==x^2,y==-1+2x},{y==(-1+x)^2x^3,y==0},{y==0,y==-1+x^2},{y==0,y==-x+x^3},{y^2-3x y+2x^2==0,x y==1},{y==x,y==x+x^2,y==(-1+x)^2+x},{y==x,y==-2+Abs[-3+(3x^2)/4]},{x^2+x y+y^2==1},{y==(x-2#)^2,y==-(x-#)^2}&/@{-1,0,2}//Flatten,{x^2y^4==1,x==0,x==1}};

In[2]:= areaCAD[#] & /@ l

Out[2]= {(4 + \[Pi])/2, 1/12, 16/15, 1, 2/3, 1/60, 4/3, 1/2, Log[2], 1/12, 208/27, (2 \[Pi])/Sqrt[3], 27/2, 4}

In[3]:= areaCAD1[{x^2 y == 1, y == 0, x^2 == 2}]

Out[3]= {x < -Sqrt[2] && y > 1/x^2,  x > Sqrt[2] && y > 1/x^2, -Sqrt[2] < x < 0 && y > 1/x^2,  0 < x < Sqrt[2] && y > 1/x^2, False, False,  x < -Sqrt[2] && 0 < y < 1/x^2,  x > Sqrt[2] &&   0 < y < 1/   x^2, (-Sqrt[2] < x < 0 && 0 < y < 1/x^2) || (x == 0 &&     y > 0) || (0 < x < Sqrt[2] && 0 < y < 1/x^2),  x < -Sqrt[2] && y < 0,  x > Sqrt[2] && y < 0, -Sqrt[2] < x < Sqrt[2] && y < 0}

In[4]:= areaCAD2[{x^2 y == 1, y == 0, x^2 == 2}]

Out[4]= {Infinity, Infinity, Infinity, Infinity, 0, 0, 1/Sqrt[2], 1/Sqrt[2], Integrate[1/x^2, {x, -Sqrt[2], Sqrt[2]}], Infinity, Infinity, Infinity}

In[5]:= areaCAD[{x^2 y == 1, y == 0, x^2 == 2}]

Out[5]= Sqrt[2]