点と平面の距離の公式「2007 早稲田大学 人間科学部【8】」をChatGPTとWolframAlphaとsympyとFreeCADでやってみたい。

・私は、「一撃」とは、いきませんでした。
**・教えて下さい。２元(x,y,z)の平方完成は,ありますか?**
・「sympyで(解法6)」は、１変数z消去です。
**・申し訳ありません。「FreeCADの接平面の移動」の勉強中です。**


オリジナル　
(youtube 5:32) math karat 様より

https://youtu.be/FfVHFricBYA?si=iJu8tPbJ-YdeXRHL

(youtube 6:43) math karat 様より

https://youtu.be/8rRUgJnqxYA

上と同じです。大学入試数学問題集成 様＞テキスト 人間科学部【8】

https://mathexamtest.web.fc2.com/2007/200713591/2007135910400.html#top-0408

# ChatGPT-3.5で(???できませんでした。？？？)
以下の結果は抜粋です。

入力文 原稿
```
実数x,y,zの間にx+2y+3z=7という関係があるとき，x2+y2+z2は
x=y/u=z/vのとき最小値p/2をとる．
```
```
＞...この方程式を用いて、x^2 + y^2 + z^2 の最小値を求めるためにラグランジュ乗数法を適用します。...
...最終的に、u, v の値に関する具体的な情報が必要です。それらの値を提供いただければ、最小値を計算する手助けができます。...
```
入力文 修正1回目
```
実数x,y,zの間にx+2y+3z=7という関係があるとき，x2+y2+z2の
最小値を求める．
```
```
???＞...従って、制約条件x + 2y + 3z = 7のもとでx^2 + y^2 + z^2の最小値は6です。
```

入力文 修正2回目(OKです。)
```
平面の方程式x+2y+3z=7と原点と距離は
```
```
＞...この式を計算すると、原点から平面までの距離は 
7/√14です。
```

# WolframAlphaで

https://ja.wolframalpha.com/input?i=x%2B2y%2B3z%3D7%2Cx2%2By2%2Bz2%E3%81%AE%E6%9C%80%E5%B0%8F%E5%80%A4

最小値
(x, y, z) = (1/2, 1, 3/2) のとき, min{x^2 + y^2 + z^2|x + 2 y + 3 z = 7} = 7/2

# sympyで(解法1)(math karat 様の方法で) 平方完成の形
勉強中
# sympyで(解法2)(math karat 様の方法で) コーシー・シュワルツの不等式
・いい方法ありませんか？
factor() takes a polynomial and factors it into irreducible factors over the rational numbers. 
https://docs.sympy.org/latest/tutorials/intro-tutorial/simplification.html#factor
```python
from sympy import *
var('a,b,c',real=True)
var('x,y,z',real=True)
Sahen=(a*x+b*y+c*z)**2
Uhen =(a**2+b**2+c**2)*(x**2+y**2+z**2)
mySubs1={a:1,b:2,c:3}
mySubs2={x+2*y+3*z:7}
print("#",Sahen                               ,"≦",       Uhen                                )
print("#",Sahen.subs(mySubs1).subs(mySubs2)   ,"≦",factor(Uhen.subs(mySubs1).subs(mySubs2))   )
print("#",Sahen.subs(mySubs1).subs(mySubs2)/14,"≦",factor(Uhen.subs(mySubs1).subs(mySubs2))/14)
print("#","x/",1,"=y/",2,"=z/",3)
# (a*x + b*y + c*z)**2 ≦ (a**2 + b**2 + c**2)*(x**2 + y**2 + z**2)
# 49 ≦ 14*(x**2 + y**2 + z**2)
# 7/2 ≦ x**2 + y**2 + z**2
# x/ 1 =y/ 2 =z/ 3
```

# sympyで(解法3)(math karat 様の方法で) 原点Oと平面の距離の公式
＞原点Oと平面の距離、方向ベクトル(1,2,3)、垂直な時
Returns the n-th coefficient of f where N are the exponents of the generators in the term of interest.
https://docs.sympy.org/latest/modules/polys/reference.html#sympy.polys.polytools.Poly.nth
```python
from sympy import *
var('x ,y ,z'   ,real=True)
var('xh,yh,zh,k',real=True)
def myGentenToHeimen(f):
    O=Point(0,0,0)
    pol=poly(f,x,y,z)
    a=pol.nth(1,0,0)
    b=pol.nth(0,1,0)
    c=pol.nth(0,0,1)
    d=pol.nth(0,0,0)
    h=abs(a*O.x+b*O.y+c*O.z+d)/sqrt(a**2+b**2+c**2)
    # 
    k=1
    (xh,yh,zh)=k*(a,b,c)
    return h,xh,yh,zh
f  =x+2*y+3*z-7
ans=myGentenToHeimen(f)
print("#",ans[0]**2,ans[1],ans[2],ans[3] )
# 7/2 1 2 3
```
# sympyで(解法4)(math karat 様の方法で) 点と平面の距離関数(メソッド?)
(解法3)と同じです。
class sympy.geometry.plane.Plane(p1, a=None, b=None, **kwargs)
https://docs.sympy.org/latest/modules/geometry/plane.html
Distance between the plane and another geometric entity.
https://docs.sympy.org/latest/modules/geometry/plane.html#sympy.geometry.plane.Plane.distance
```python
from sympy import *
var('x,y,z',real=True)
def myXYZ2Plane(myeq):
    pol=poly(myeq.lhs-myeq.rhs,x,y,z)
    return Plane( Point( solve(myeq.subs({y:0,z:0}))[0] ,0,0),normal_vector=(
                  pol.nth(1,0,0),
                  pol.nth(0,1,0),
                  pol.nth(0,0,1)
                 )),pol.nth(1,0,0),pol.nth(0,1,0),pol.nth(0,0,1)
eq1=Eq(x+2*y+3*z,7)
eqPlane=myXYZ2Plane(eq1)
print("#",eqPlane[0].distance  (Point(0,0,0))**2 )
print("#",eqPlane[0].projection(Point(0,0,0))    )
print("#",eqPlane[1],
          eqPlane[2],
          eqPlane[3]
     )
# 7/2
# Point3D(1/2, 1, 3/2)
# 1 2 3
```

# sympyで(解法5)（math karat 様の方法で) 内積=0で
・(解法5)は、難しいです。
```python
from sympy import *
var('x ,y ,z'   ,real=True)
var('xp,yp,zp,k',real=True)
def myPoly(f):
    pol=poly(f,x,y,z)
    a=pol.nth(1,0,0)
    b=pol.nth(0,1,0)
    c=pol.nth(0,0,1)
    d=pol.nth(0,0,0)
    return a,b,c,d
f  =x+2*y+3*z-7
ans=myPoly(f)
O=Point(0 ,0 ,0 )
P=Point(xp,yp,zp)
# 
A=Point(solve(f.subs({y:0,z:0}),x)[0],0,0)
B=Point(0,solve(f.subs({x:0,z:0}),y)[0],0)
C=Point(0,0,solve(f.subs({x:0,y:0}),z)[0]) 
print("# A",A)
print("# B",B)
print("# C",C)
# 
eqOPAB=Eq((P-O).dot(B-A),0)
eqOPAC=Eq((P-O).dot(C-A),0)
print("#",eqOPAB)
print("#",eqOPAC)
# 
x_y=solve(eqOPAB,xp)[0]
x_z=solve(eqOPAC,xp)[0]
print("#",x_y)
print("#",x_z)
eqypzp=Eq(x_y,x_z)
print("# xp/1 ",eqypzp)
# A Point3D(7, 0, 0)
# B Point3D(0, 7/2, 0)
# C Point3D(0, 0, 7/3)
# Eq(-7*xp + 7*yp/2, 0)
# Eq(-7*xp + 7*zp/3, 0)
# yp/2
# zp/3
# xp/1  Eq(yp/2, zp/3)
```

# sympyで(解法6)zを消去後、偏微分の極値で。
・zをx,yで表して、与式に代入して与式のzを消去した。
　2変数関数に変形した。偏微分を使って、2変数関数の最小値を計算した。。
・sympyの標準関数に、2変数関数の最小値は見当たりませんでした。
　mathemticaにあったような気がします。
・ユーザー定義関数にしたほうがいいですネ。
sympy.plotting.plot.plot3d(*args, show=True, **kwargs)
Plots a 3D surface plot.
https://docs.sympy.org/latest/modules/plotting.html#sympy.plotting.plot.plot3d
```python
from sympy import *
var('x,y,z',real=True)
eq1 =Eq(x+2*y+3*z,7)
h1  =x**2+y**2+z**2
#
z_xy=solve(eq1,z)[0]
h2  =h1.subs({z:z_xy})
ans =solve( [ diff(h2,x),diff(h2,y) ],[x,y])
print("#",h2.subs( {x:ans[x],y:ans[y]} ))
# 7/2
# 以下省略

#　作図
from sympy.plotting import plot3d
print ("#",h2)
plot3d(    h2, (x, -5, 5), (y, -5, 5))
# x**2 + y**2 + (-x/3 - 2*y/3 + 7/3)**2
```
x＊＊2 + y＊＊2 + (-x/3 - 2＊y/3 + 7/3)＊＊2
![0png.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/109698/ade49ac1-52f7-aeae-a66e-b88d5945d08e.png)

# FreeCADのマクロで作図

```python
import FreeCAD
import Part
import Draft
import Mesh
#########################################################################################################
# 計算
from sympy import *
var('x,y,z',real=True)
def myXYZ2Plane(myeq):
    pol=poly(myeq.lhs-myeq.rhs,x,y,z)
    return Plane( Point( solve(myeq.subs({y:0,z:0}))[0] ,0,0),normal_vector=(
                  pol.nth(1,0,0),
                  pol.nth(0,1,0),
                  pol.nth(0,0,1)
                 )),pol.nth(1,0,0),pol.nth(0,1,0),pol.nth(0,0,1)
eq1=Eq(x+2*y+3*z,7)
eqPlane=myXYZ2Plane(eq1)
print("#",eqPlane[0].distance  (Point(0,0,0))**2 )
print("#",eqPlane[0].projection(Point(0,0,0))    )
print("#",eqPlane[1],
          eqPlane[2],
          eqPlane[3]
     )
# 7/2
# Point3D(1/2, 1, 3/2)
# 1 2 3
##########################################################################################################
# 3D
r=sqrt(7/2)
O=Point3D(0  ,0,  0)
H=Point3D(1/2,1,3/2)
##########################################################################################################
### 3D作図
def myXYZ2Txt(A):
    return '(' + str(A.x) + ',' + str(A.y) + ',' + str(A.z) + ')'
def myTxtXYZ(A,myWedgei):
    P5x=float(A.x)
    P5y=float(A.y)
    P5z=float(A.z)
    p5 = FreeCAD.Vector(P5x, P5y, P5z)
    myText = Draft.makeText(myWedgei, p5)
    myText.Label = myWedgei
    # FreeCADGui.ActiveDocument.ActiveObject.FontSize = '40.0 mm'
    # FreeCADGui.ActiveDocument.ActiveObject.FontSize = '3.0 mm'
    FreeCADGui.ActiveDocument.ActiveObject.FontSize = '0.5 mm'
def myLine(A,B):
    Ax,Ay,Az=float(A.x),float(A.y),float(A.z)
    Bx,By,Bz=float(B.x),float(B.y),float(B.z)
    pl = FreeCAD.Placement()
    pl.Rotation.Q = (0.4247081540122249, 0.17592004639554645, 0.33985110062924484, 0.8204732460821097)
    pl.Base = FreeCAD.Vector(-3.9166066876399563, -2.1670824762243774, 1.7495260956243028)
    points = [FreeCAD.Vector(Ax,Ay,Az), FreeCAD.Vector(Bx,By,Bz)]
    line = Draft.make_wire(points, placement=pl, closed=False, face=True, support=None)
    Draft.autogroup(line)
def myLine_S(*args):
    for i in range(1,len(args)):
        myLine(args[i-1],args[i])
    return 0
def myTxtXYZ_S(*xy_tx):
    for i in range(1,int(len(xy_tx)/2)+1):
        myTxtXYZ(xy_tx[2*i-2],xy_tx[2*i-1]+myXYZ2Txt(xy_tx[2*i-2]) ) 
    return
def myLine_C(*args):
    for i in range(1,len(args)):
        myLine(args[i-1],args[i])
    myLine(args[i],args[0])
    return
#################################################################################
var('myX,myY,myZ')
def myDegree(myX,myY,myZ):
    r_x=0
    r_y=0
    r_z=0
    if myX==0:
       r_z=0    
    else:
       r_z=atan( myZ/myX ) *180/pi    
    if myY==0:
       r_z=90    
    else:
       r_z=atan( myZ/myY ) *180/pi    

    r_x=float(r_x)
    r_y=float(r_y)
    r_z=float(r_z)
    return r_x,r_y,r_z 
#
def my_Sphere_name(myCenter,myr,myname):
    sphere = doc.addObject("Part::Sphere", myname)
    sphere.Radius = float(myr)
    sphere.Angle1 = -90
    sphere.Angle2 =  90
    sphere.Angle3 = 360
    sphere.Placement = App.Placement(App.Vector(float(myCenter.x),float(myCenter.y),float(myCenter.z)), App.Rotation(0, 90, 0))
    FreeCADGui.getDocument('Unnamed').getObject(myname).Transparency = 80
doc = App.activeDocument()
plane = doc.addObject("Part::Plane", "myPlane1")
(plx,ply,plz)             =(H.x,H.y,H.z)
(plane_Length,plane_Width)=(10,10)
plane.Length = plane_Length
plane.Width =  plane_Width
plane.Placement = App.Placement(App.Vector(plx,ply,plz), App.Rotation( 0, 0, 0 ))
#
plane = doc.addObject("Part::Plane", "myPlane2")
(plx,ply,plz)             =(H.x,H.y,H.z)
(plane_Length,plane_Width)=(10,10)
plane.Length = plane_Length
plane.Width =  plane_Width
plane.Placement = App.Placement(App.Vector(plx,ply,plz), App.Rotation( 0, 0, myDegree(H.x,H.y,H.z)[2]+90 ))
#
my_Sphere_name(O,r,"mySphere_AB")
myTxtXYZ_S( O,"O",H,"H" )
myLine_S  ( O,H   )
#################################################################################
doc = App.activeDocument()
App.ActiveDocument.addObject("App::Origin", "Origin")
# App.ActiveDocument.getObject('Origin').Visibility = True
App.ActiveDocument.recompute()
Gui.activeDocument().activeView().viewAxonometric()
Gui.SendMsgToActiveView("ViewFit")
```
**申し訳ありません。「FreeCADの接平面の移動」の勉強中です。**
接点が平面の角の点です。
datum planeでできませんでした。 planeでしました。
planeの隅の4点の座標計算を探しています。

isometric方向?です。
CADの操作で、transparencyを80にしました。
![0png.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/109698/6138b301-3949-4c18-60f0-c126e1358037.png)
正面からです。
![1png.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/109698/29a73a77-24e9-9f56-22dc-9a14d05e231b.png)

# sympyの実行環境
①私の環境は,pycharmです。
②よく聞くのは、Jupyterです。
③web上で、上記のソースを「SymPy Live shell」に、コピー貼り付けでもできました。
黒背景の右上に、マウスを移動すると、コピーマークが発生します。
???　タブレット環境で、コピー貼り付けが実行できませんでした。???

https://live.sympy.org/

# 参考

以下、いつもの?おすすめです。

https://qiita.com/mrrclb48z/items/75904cdc6e84ecef4d03#%E5%8F%82%E8%80%83wolframalpha-%E9%96%8B%E7%99%BA%E7%89%A9%E8%AA%9E

いつもと違うやつです。

https://qiita.com/mrrclb48z/items/38e849f41f46bf0c8f9c