0
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 1 year has passed since last update.

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

Last updated at Posted at 2023-09-15

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

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

(youtube 6:43) math karat 様より

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

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で

最小値
(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

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

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

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)は、難しいです。

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

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

FreeCADのマクロで作図

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
正面からです。
1png.png

sympyの実行環境

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

参考

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

いつもと違うやつです。

0
0
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
0
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?