(空間のベクトル)「2021東京農工大学数学(Z)問1」をsympyでやってみた。Point3Dでsolve


solveで,Point3Dが使えました。

<問題>数学(Z)問1
4点 A(1,1,3),B(-1,1,-1),C(1,-4,-2),D(-2,-1,1)

<解答例>

Pycharmで

答えは、あっていると思います。

from sympy import *
var('s t u')
def myInnerProduct(myPA,myPB):
    myA=myPoint3DtoMatrix(myPA)
    myB=myPoint3DtoMatrix(myPB)
    mySum=0
    for i in range(0, shape(myA)[1]):
        mySum = mySum+myA[i]*myB[i]
    return mySum
def myPoint3DtoMatrix(myPoint3D):
    myx = Matrix(1, 3, range(3))
    myx[0]=myPoint3D.x
    myx[1]=myPoint3D.y
    myx[2]=myPoint3D.z
    return myx
myZahyo=((0,0,0),(1,1,3),(-1,1,-1),(1,-4,-2),(-2,-1,1))
O=Point3D(myZahyo[0])
A=Point3D(myZahyo[1])
B=Point3D(myZahyo[2])
C=Point3D(myZahyo[3])
D=Point3D(myZahyo[4])
ans=solve(O-D-(s*(A-D)+t*(B-D)+u*(C-D)),[s,t,u])
print("#〔1〕s,t,u:",ans[s],ans[t],ans[u])
#
var('v w')
P=A+v*(B-A)
Q=C+w*(D-C)
ans=solve([myInnerProduct(B-A,Q-P),myInnerProduct(Q-P,D-C)],[v,w])
mySubs={v:ans[v],w:ans[w]}
myFutosiki= 0 <= v.subs(mySubs) <= 1
print("#〔2〕線分上:",type(myFutosiki),myFutosiki)
print("#    PQ   :",P.distance(Q).subs(mySubs),
               "P:",P.x.subs(mySubs),P.y.subs(mySubs),P.z.subs(mySubs),
               "Q:",Q.x.subs(mySubs),Q.y.subs(mySubs),Q.z.subs(mySubs))
#
var('x')
R=C+x*(D-C)
cosARB=myInnerProduct(A-B,R-B)/(Point(B).distance(A)*Point(B).distance(R))
ans=solve(diff(cosARB))
mySubs={x:ans[0]}
print("#〔3〕ABR  :",ans[0],cosARB.subs({x:ans[0]}),
               "R:",R.x.subs(mySubs),R.y.subs(mySubs),R.z.subs(mySubs))
#〔1〕s,t,u: 3/10 1/2 1/5
#〔2〕線分上: <class 'sympy.logic.boolalg.BooleanTrue'> True
#    PQ   : 5*sqrt(14)/7 P: -3/7 1 1/7 Q: -13/7 -8/7 6/7
#〔3〕ABR  : 5/4 sqrt(39)/13 R: -11/4 -1/4 7/4

参考

wolframalphaでグラフを作図

x = 5/4 のとき, max{(3 sqrt(5) x)/(5 sqrt((1 - 3 x)^2 + (5 - 3 x)^2 + (3 x - 2)^2))} = sqrt(3/13)

FreeCAD

検討中