Quantifier Elimination with Mathematica #2
ー東大入試を ２行で  解くー

2018 年7月  Cabri研究会　生越　茂樹

QEのコマンド

§１．いろいろな方程式

g[x_]=x^2-(a+d)x+(a d-b c);
ForAll[{a,b,c,d,s,t},a>0&&b>0&&c>0&&d>0&&s>0&&t>0&&s (1-a)-t b>0&&-s c+t(1-d)>0]\[Implies]Exists[{x1,x2},g[x1]==g[x2]==0&&-1<x1<x2<1];  

Resolve[%,Reals]

(*(1)はExists{x1,x2,x3,p}とする*)
Clear[x1,x2,x3,M,p]
Exists[{x1,x2,x3,p},(x1|x2|x3|p|M)∈Reals,x1+x2+x3==0&&x1 x2+x2 x3+x3 x1==-3&&x1 x2 x3==p&&x1≤x2≤x3&&x1*x3==M]∨
Exists[{x1,x2,x3,p},(x1|p|M)∈Reals&&!x2∈Reals&&!x3∈Reals&&x1+x2+x3==0&&x1 x2+x2 x3+x3 x1==-3&&x1 x2 x3==p&&x1*x1==M];
Resolve[%,Complexes]

(*(2)はExists{x1,x2,x3}とする*)
Exists[{x1,x2,x3},(x1|x2|x3|p|M)∈Reals,x1+x2+x3==0&&x1 x2+x2 x3+x3 x1==-3&&x1 x2 x3==p&&x1≤x2≤x3&&x1*x3==M]∨
Exists[{x1,x2,x3},(x1|p|M)∈Reals&&!x2∈Reals&&!x3∈Reals&&x1+x2+x3==0&&x1 x2+x2 x3+x3 x1==-3&&x1 x2 x3==p&&x1*x1==M];
Resolve[%,Complexes]
RegionPlot[ImplicitRegion[%,{{p,-5,5},{M,-3,7}}],Axes->True]

§２．最大・最小値

Clear[x,y,t]
Exists[x,t>0&&((y-t/2≥(x-t/2)^2&&y-t/2≥(x+t/2)^2)||(y-t/Sqrt[2]≥x^2&&y≥(x+t/Sqrt[2])^2&&y≥(x-t/Sqrt[2])^2))]
reg=Resolve[%,Reals]
RegionPlot[reg,{t,0,5},{y,0,5},PlotPoints->50,Axes->True]

P (p, 0), A, B の中心を (a, b), (c, d)．半径を (Subscript[r, 1]), (Subscript[r, 2]) とおく． (p-c)(p-b)≤0　が 「一点のみ」共有の条件．

Clear[a,b,c,d,r1,r2,p]
Exists[{a,b,c,d,r1,r2,p},a^2+(p-b)^2==r1^2&&(p-c)^2+d^2==r2^2&&(p-c)(p-b)≤0&&a^2+b^2≤(1-r1)^2&&c^2+d^2≤(1-r2)^2&&0<r1<=1&&0<r2<=1&&-1<p<1&&k==r1+r2];

Resolve[%,Reals]

§３．領域

領域を表す式が求まらない，或いは複雑すぎる時でも ContourPlot(3D) , RegionPlot(3D) などを使って領域は図示できる．
また，面積・体積も RegionMeasure を使って 求められる事が多い．

ClearAll[x1,x2,y1,y2,d,a,x,y]
d[{x1_,y1_},{x2_,y2_}]=Abs[x1-x2]+Abs[y1-y2];
(*1*)
FullSimplify[d[{0,0},{x,y}]==d[{x,y},{1,1}]];
RegionPlot[%,{x,-2,3},{y,-2,3},Axes->True]

(*2*)
Exists[a,a>=0&&d[{0,0},{x,y}]==d[{x,y},{a,a^2+1}]];
Resolve[%,Reals]//FullSimplify
RegionPlot[%,{x,-1,1},{y,0,2}]

ClearAll["Global`*"]; 
a={-5,-5};b={5,-5};c={5,5};d={-5,5};p={-9,-3};q={9,3};
dP[x_,y_]=Which[(y≥2x+15||2y+x≤-15||x≤-5),ArcLength[Line[{p,{x,y}}]],y≥5||y≥1/2 Sqrt[x^2+10x+45],ArcLength[Line[{p,d,{x,y}}]],y<-5||y<1/2 Sqrt[x^2+10x+45]&&x≤5,ArcLength[Line[{p,a,{x,y}}]],y≥1/2Sqrt[x^2-10x+45],ArcLength[Line[{p,d,c,{x,y}}]],y<1/2Sqrt[x^2-10x+45],ArcLength[Line[{p,a,b,{x,y}}]]];
dQ[x_,y_]=dP[x,y]/.{x->-x,y->-y,p->q,q->p,a->c,c->a,b->d,d->b};
path=ContourPlot[dP[x,y]==dQ[x,y],{x,-12,12},{y,-12,12}];

Manipulate[ 
a={-5,-5};b={5,-5};c={5,5};d={-5,5};p={-9,-3};q={9,3};  x=r[[1]];y=r[[2]];
pts1=Which[(y≥2x+15||2y+x≤-15||x≤-5),{p,r},y≥5||y≥1/2 Sqrt[x^2+10x+45],{p,d,r},y<-5||y<1/2 Sqrt[x^2+10x+45]&&x≤5,{p,a,r},y≥1/2Sqrt[x^2-10x+45],{p,d,c,r},y<1/2Sqrt[x^2-10x+45],{p,a,b,r}];
pts2=Which[-y≥15-2 x||-x-2 y≤-15||-x≤-5,{q,r},-y≥5||-y≥1/2 (Sqrt[(((-x))^2)+10 (-x)+45]),{q,b,r},-y<-5||(-y<1/2 (Sqrt[(((-x))^2)+10 (-x)+45])&&-x≤5),{q,c,r},-y≥1/2 (Sqrt[(((-x))^2)-10 (-x)+45]),{q,b,a,r},-y<1/2 (Sqrt[(((-x))^2)-10 (-x)+45]),{q,c,d,r}];
delta=ArcLength[Line[pts1]]-ArcLength[Line[pts2]];
Show[
Graphics[{Point[r],Thick, Line[{a,d}],Line[{b,c}],Thin,Blue,Line[pts1],Red,Line[pts2],Black,
Point[p],Point[q],Text["P",p,{2,Center}],Text["Q",q,{-2,Center}],Text["R",r,{0,2}],Text["A",a,{Center,Top}],Text["B",b,{Center,Top}],Text["C",c,{Center,Bottom}],Text["D",d,{Center,Bottom}],Text[ Style["[青線長-赤線長]≒"<>ToString[Floor[100N[delta]]/100.],9,Bold],{7,12}]},PlotRange->{{-13,13},{-13,13}},Axes->True,ImageSize->240],path
],(*insert path*) {{r,{0,0}},Locator}]

１Tテキストに載っていた問題．Rが中点以外のときは複雑だが，QEなら簡単．

Clear[p,q,a,b]
Exists[{p,q},-1≤p≤1&&-1≤q≤1&&2p+q==3a&&2p^2+q^2==3b];
Resolve[%,Reals];
reg=FullSimplify[%]
RegionPlot[reg,{a,-1,1},{b,0,1},PlotPoints->50,Axes->True]

Pの軌跡（曲面)を求めてから 線分OPの通過領域を求める．さもないとReduceでは一晩以上かかる．  
また式は途中でできるだけ簡単にした方がよい．「(s^2) + (t^2) == 1 && (q^2) + (r^2) + (p^2) == 1」 のとき
 「(p^2) + (q - ((s))^2) + (r - ((t))^2) = 1 ⇔ 2 qs + 2 rt = 1」　このほうが実行時間が1/10 になる．

Clear[p,q,r,s,t,x,y,z,u]
Exists[{s,t},2q s+2r t==1&&s^2+t^2==1&&q^2+r^2+p^2==1];
surface=Resolve[%,Reals]

Exists[{p,q,r,u},x==u p&&y==u q&&z==u r&&0<u<1&&surface];(*{p,q,r}と{p,q,r,u}の間違いに注意*)
Resolve[%,Reals]
body=ImplicitRegion[%,{x,y,z}];

RegionPlot3D[body,PlotPoints->50, PlotStyle->Directive[Green,Opacity[0.5]],Mesh->None]

RegionMeasure で体積は求まるが時間は掛かる (約4分)

RegionMeasure[body]//Timing

[点射影]

a,b,cのままでは図示できないので，a = 2, b = 1, c = 3 として解答する. 
P(X,Y,c+1), R上の点をR’(x,y,c), 直線PR’とxy平面の交点をQ(u,v,0)とおくと,　(OverVector[QR'])= t (OverVector[QP])⇔(({{x-u}, {y-v
c}})) = t (({{X-u}, {Y-v
c+1}}))

Clear[a,b,c,u,v,x,y,X,Y];
a=2;b=1;c=3;

Exists[{x,y,X,Y,t},Abs[x]≤a&&Abs[y]≤b&&{x-u,y-v,c}==t{X-u,Y-v,c+1}&&X^2/a^2+Y^2/b^2==1];
pred=Resolve[%,Reals]//FullSimplify;

RegionPlot[pred,{u,-15,15},{v,-8,8},Axes->True,AspectRatio->Automatic]

面積は，文字のままで求まる

Clear[a,b,c,u,v,x,y,t,X,Y];
Exists[{x,y,X,Y,t},a>0&&b>0&&c>0&&Abs[x]≤a&&Abs[y]≤b&&{x-u,y-v,c}==t{X-u,Y-v,c+1}&&X^2/a^2+Y^2/b^2==1];
pred=Resolve[%,Reals];
{time,area}=RegionMeasure[ImplicitRegion[pred,{u,v}]]//Timing//FullSimplify

（おまけ）これらも　３行で　十分

【楕円が円に含まれる条件】

(*1*)Clear[a,b]
ForAll[{x,y},x^2/a+(y-Sqrt[b])^2/b==1&&a>0&&b>0\[Implies]x^2+y^2≤1];
reg=Resolve[%,Reals]
RegionPlot[reg,{a,0,1},{b,0,1/2}]

(*2*)Exists[{a,b},Pi Sqrt[a b]==k&&reg&&a>0&&b>0];
Resolve[%,Reals]//FullSimplify

【三角形が正方形に含まれる条件】

手で解くときは正三角形は回転で考える事が多いが，QEの場合は長さで考える．

dis[x_]:=Total[x*x];(*Oとxの距離の2乗*)

(*1*)Exists[{c1,c2},dis[{a,0}-{c1,c2}]==dis[{0,b}-{c1,c2}]==dis[{a,0}-{0,b}]&&0<a≤1&&0<b≤1&&0≤c1≤1&&0≤c2≤1];
reg=Resolve[%,Reals]//FullSimplify
RegionPlot[%,{a,0,1},{b,0,1}]
(*2*)Exists[{a,b},Sqrt[3]/4(a^2+b^2)==k&&reg];
Resolve[%,Reals]

【3角関数を含んだ最小値】

三角関数も置き換えで多項式になる場合がある．そのときはQEで解ける．

(*1*)Clear[a,b,x]
TrigExpand[Cos[3θ]+a Cos[2θ]+b Cos[θ]]
f[x_]=%/.{Sin[θ]->Sqrt[1-x^2],Cos[θ]->x}//Expand;
g[x_]=Simplify[(f[x]-f[1])/(x-1)];
Print["f(θ)=", TraditionalForm[f[x]],"    ","g(θ)=", TraditionalForm[g[x]]]

(*2*)p=ForAll[x,(0<x<1\[Implies]g[x]≥0)]&&Exists[y,g[y]==0&&0<y<1];
(*p=Exists[y,g[y]==0&&0<y<1&&ForAll[x,0<x<1\[Implies]g[x]≥0]]もOK*)
r=Resolve[p,Reals]//FullSimplify
RegionPlot[ImplicitRegion[r,{a,b}]]
(*ContourPlot[(a^2)==4 (a+b),{a,-20,20},{b,-20,20},RegionFunction->Function[{a,b},((1+2 a+b>0&&-6<a≤-4)||(9+4 a+b>0&&-4<a<-2))]]*)
(*↑[別法]ContourPlotを使い RegionFunctionで (a,b)の範囲を制限する*)

積分計算してしまえば多項式の問題なのでQEで解ける．

Clear[x]
Exists[{α,β},0≤α≤β&&Integrate[(x-α)(x-β),{x,-1,1}]==1&&Integrate[(x-α)(x-β),{x,0,α}]==S];

Resolve[%,Reals]

虚数を扱うときは Resolve[, Complexes] と a ∈ Reals を組み合わせる．

Clear[x1,x2,p,q]
Exists[{x1,x2},(p|q)∈Reals&&x1+x2==-(p^2-q)/p&&x1*x2==-(2p-q-1)/p&&p≠0&&Re[x1]<0&&Re[x2]<0];
A=Resolve[%,Complexes]//FullSimplify;
Exists[x,-q x+q+1==0&&p==0&&x<0];
B=Resolve[%,Reals]//FullSimplify;
FullSimplify[A||B]
RegionPlot[A||B,{p,-2,2},{q,-2,4},PlotPoints->50,Axes->True]

【パラメータの多い通過領域】

残念ながら領域を表す式は綺麗にならないが，図示はでき，その上，面積も求まる．

Clear[p1,p2,r1,r2,x,y]
a=2(Sqrt[3]+1);
Exists[{p1,p2},p1^2+p2^2==4&&(p1-x)^2+(p2-y)^2==4]&&Exists[{r1,r2},(r1-a)^2+r2^2==2&&(x-r1)^2+(y-r2)^2==2] ;
Resolve[%,Reals];
reg=FullSimplify[%]

RegionPlot[reg,{x,2,5},{y,-2,2},Axes->True]

RegionMeasure[ImplicitRegion[reg,{x,y}]]//FullSimplify

[立体射影]

Clear[k,u,v,x,y,z];k=3;(*オリジナルは k　のまま*)
Exists[{x,y,z,t},x^2+y^2+z^2==1&&x≥1/k&&y≥1/k&&z≥1/k&&{x,y,z-1}==t{u,v,-1}&&0<t<1];
pred=Reduce[%,{u,v},Reals]//FullSimplify

reg=ImplicitRegion[pred,{u,v}];
RegionPlot[reg]