LoginSignup
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.

@mrrclb48z

四面体「2023群馬大学前期【３】」をsympyとFreeCADでやってみたい。できませんでした。

Posted at

できませんでした。
FreeCADのマクロで作図しました。

オリジナル 公式ホームページより
医(医学科)学部 理工,情報学部【３】の類題
問題文と解答例

一般選抜（前期日程）、私費外国人留学生選抜＜群馬大学

上と同じです。大学入試数学問題集成＞【3】テキスト

sympyで(公式ページの解答例を参考に)

勉強のしなおしです。

sympyで(sympy的?安易なやり方)

できませんでした。座標計算です。
・Pointクラスの値の変更がわからなかったので、3回繰り返しました。
　(1),(2),(3)に分けました。
・「同一平面上は...」、点から平面への垂線(projection投影)です。
・絶対値のはずしかた？意味を教えて下さい。

from sympy import *
var('x y z')
ax,ay   =symbols("ax ay"   ,real=True)
cx,cy   =symbols("cx cy"   ,real=True)
dx,dy   =symbols("dx dy "  ,real=True)
dz      =symbols(      "dz",real=True,positive=True)
u,v,w   =symbols("u v w"   ,real=True)
s,t     =symbols("s t"     ,real=True,positive=True)
L       =symbols("L"       ,real=True,positive=True)
def myMenseki3D(O,A,B):
    P=A-O
    Q=B-O
    return Rational(1, 2) * sqrt(P.distance(Point(0, 0, 0)) ** 2 * Q.distance(Point(0, 0, 0)) ** 2 \
                                 - P.dot(Q) ** 2).simplify()
def mySphereFormula(myP,myR):
    return ((x-myP.x)**2+(y-myP.y)**2+(z-myP.z)**2-myR**2).expand()
def myTaiseki(PTO,PTA,PTB,PTC):
    return Matrix([[PTA.x-PTO.x, PTA.y-PTO.y, PTA.z-PTO.z],
                   [PTB.x-PTO.x, PTB.y-PTO.y, PTB.z-PTO.z],
                   [PTC.x-PTO.x, PTC.y-PTO.y, PTC.z-PTO.z]]
                  ).det()/6
def myVectorABC(OQ, OA,OB,OC):
    return solve([
                  Eq(OQ.x,u*OA.x+v*OB.x+w*OC.x),
                  Eq(OQ.y,u*OA.y+v*OB.y+w*OC.y),
                  Eq(OQ.z,u*OA.z+v*OB.z+w*OC.z)
                 ],
                 [u,v,w])
###########################################
O=Point(0 ,0 ,0 )
A=Point(ax,ay,0 )
C=Point(cx,cy,0 )
B=A+C
D=Point(dx,dy,dz)
X=B+Rational(2,3)*(D-B)
P=A+           s *(D-A)
Q=C+           t *(D-C)
Y=Plane(A,C,D).intersection(Line(O,X))[0]
print("#(1)",myVectorABC(Y, A,C,D))
########################################################
O=Point(0,0,0)
A=Point(L,0,0)
C=Point(0,L,0)
B=A+C
dum,Dxyz=solve([mySphereFormula(O,L),
                mySphereFormula(A,L),
                mySphereFormula(B,L)],
               [x,y,z]
              )
D=Point(Dxyz)
X=B+Rational(2,3)*(D-B)
P=A+           s *(D-A)
Q=C+           t *(D-C)
ft=solve(myTaiseki(O,X,P,Q),t)[0]
print("#(2)",solve(Eq(ft,0),s)[0],solve(Eq(ft,1),s)[0] )
########################################################
L=1
O=Point(0,0,0)
A=Point(L,0,0)
C=Point(0,L,0)
B=A+C
dum,Dxyz=solve([mySphereFormula(O,L),
                mySphereFormula(A,L),
                mySphereFormula(B,L)],
               [x,y,z]
              )
D=Point(Dxyz)
X=B+Rational(2,3)*(D-B)
P=A+           s *(D-A)
Q=C+           t *(D-C)
ft=solve(myTaiseki(O,X,P,Q),t)[0]
DPQ=myMenseki3D(D,P,Q).subs({t:ft})
print("#(3)XXXX",DPQ )
# 間違っています。絶対値をはずしました。
DPQabs=(s - 1)*((3*s - 2)/(4*s - 3) - 1)/2
ans=solve(diff(DPQabs,s))
print("#(3)",ans[0],DPQabs.subs({s:ans[0]}) )
plot(DPQ,(s,solve(Eq(ft,0),s)[0],solve(Eq(ft,1),s)[0]  ))
#(1) {u: 1/4, v: 1/4, w: 1/2}
#(2) 2/3 1
#(3)XXXX Abs((s - 1)*((3*s - 2)/(4*s - 3) - 1))/2
#(3) 1/2 1/8

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

1/3<=s<=1の△DPQの面積のグラフです。
plotコマンド。縦横比が異なります。

0png.png

FreeCADのマクロで作図

import FreeCAD
import Part
import Draft
import Mesh
#########################################################################################################
# 計算
from sympy import *
var('x y z')
ax,ay   =symbols("ax ay"   ,real=True)
cx,cy   =symbols("cx cy"   ,real=True)
dx,dy   =symbols("dx dy "  ,real=True)
dz      =symbols(      "dz",real=True,positive=True)
u,v,w   =symbols("u v w"   ,real=True)
s,t     =symbols("s t"     ,real=True,positive=True)
L       =symbols("L"       ,real=True,positive=True)
def myMenseki3D(O,A,B):
    P=A-O
    Q=B-O
    return Rational(1, 2) * sqrt(P.distance(Point(0, 0, 0)) ** 2 * Q.distance(Point(0, 0, 0)) ** 2 \
                                 - P.dot(Q) ** 2).simplify()
def mySphereFormula(myP,myR):
    return ((x-myP.x)**2+(y-myP.y)**2+(z-myP.z)**2-myR**2).expand()
def myTaiseki(PTO,PTA,PTB,PTC):
    return Matrix([[PTA.x-PTO.x, PTA.y-PTO.y, PTA.z-PTO.z],
                   [PTB.x-PTO.x, PTB.y-PTO.y, PTB.z-PTO.z],
                   [PTC.x-PTO.x, PTC.y-PTO.y, PTC.z-PTO.z]]
                  ).det()/6
def myVectorABC(OQ, OA,OB,OC):
    return solve([
                  Eq(OQ.x,u*OA.x+v*OB.x+w*OC.x),
                  Eq(OQ.y,u*OA.y+v*OB.y+w*OC.y),
                  Eq(OQ.z,u*OA.z+v*OB.z+w*OC.z)
                 ],
                 [u,v,w])
###########################################
O=Point(0 ,0 ,0 )
A=Point(ax,ay,0 )
C=Point(cx,cy,0 )
B=A+C
D=Point(dx,dy,dz)
X=B+Rational(2,3)*(D-B)
P=A+           s *(D-A)
Q=C+           t *(D-C)
Y=Plane(A,C,D).intersection(Line(O,X))[0]
print("#(1)",myVectorABC(Y, A,C,D))
########################################################
O=Point(0,0,0)
A=Point(L,0,0)
C=Point(0,L,0)
B=A+C
dum,Dxyz=solve([mySphereFormula(O,L),
                mySphereFormula(A,L),
                mySphereFormula(B,L)],
               [x,y,z]
              )
D=Point(Dxyz)
X=B+Rational(2,3)*(D-B)
P=A+           s *(D-A)
Q=C+           t *(D-C)
ft=solve(myTaiseki(O,X,P,Q),t)[0]
print("#(2)",solve(Eq(ft,0),s)[0],solve(Eq(ft,1),s)[0] )
########################################################
L=1
O=Point(0,0,0)
A=Point(L,0,0)
C=Point(0,L,0)
B=A+C
dum,Dxyz=solve([mySphereFormula(O,L),
                mySphereFormula(A,L),
                mySphereFormula(B,L)],
               [x,y,z]
              )
D=Point(Dxyz)
X=B+Rational(2,3)*(D-B)
P=A+           s *(D-A)
Q=C+           t *(D-C)
ft=solve(myTaiseki(O,X,P,Q),t)[0]
DPQ=myMenseki3D(D,P,Q).subs({t:ft})
print("#(3)XXXX",DPQ )
# wrong
DPQabs=(s - 1)*((3*s - 2)/(4*s - 3) - 1)/2
ans=solve(diff(DPQabs,s))
print("#(3)",ans[0],DPQabs.subs({s:ans[0]}) )
#plot(DPQ,(s,solve(Eq(ft,0),s)[0],solve(Eq(ft,1),s)[0]  ))
##########################################################################################################
### 作図用
# wrong
P=P.subs({s:1/2})
Q=Q.subs({t:1/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 = '0.05 mm'
    return
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(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)
    return
def myLine_S(*args):
    for i in range(1,len(args)):
        myLine(args[i-1],args[i])
    return
def myLine_C(*args):
    for i in range(1,len(args)):
        myLine(args[i-1],args[i])
    myLine(args[i],args[0])
    return
def myLine_H(*args):
    for i in range(1,len(args)):
        myLine(args[0],args[i])
    return
#
print("X",X)
print("P",P)
print("Q",Q)
#
myLine_C  (O,A,B,C)
myLine_H  (D, O,A,B,C)
myLine_C  (O,P,X,Q)
myLine_C  (D,P,Q)
myTxtXYZ_S(O,"O",A,"A",B,"B",C,"C",D,"D")
myTxtXYZ_S(X,"X",P,"P",Q,"Q")
#
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")

isometric方向です。
Originは、非表示です。マウス操作で、線分と文字を赤色にしました。

0png.png

上図の拡大図です。

1png.png

参考

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
Sign upLogin
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?