solveで,Point3Dが使えました。
4点 A(1,1,3),B(-1,1,-1),C(1,-4,-2),D(-2,-1,1)
公式サイト
<問題>数学(Z)問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のマクロで
上のsympy計算に、FreeCAD作図を追加しました。
操作で、線2本を、赤色と緑色にしました。
import FreeCAD
import Part
import Draft
import Mesh
#########################################################################################################
# 計算
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]}
P=Point3D(P.x.subs(mySubs),P.y.subs(mySubs),P.z.subs(mySubs))
Q=Point3D(Q.x.subs(mySubs),Q.y.subs(mySubs),Q.z.subs(mySubs))
myFutosiki= 0 <= v.subs(mySubs) <= 1
print("#〔2〕線分上:",type(myFutosiki),myFutosiki)
print("# PQ :",P.distance(Q),"P:",P,"Q:",Q)
#
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]}
R=Point3D(R.x.subs(mySubs),R.y.subs(mySubs),R.z.subs(mySubs))
print("#〔3〕ABR :",ans[0],cosARB.subs({x:ans[0]}),
"R:",R)
#〔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
# #########################################################################################################
# 3D作図
def myXYZ2Txt(A):
return '(' + str(A.x) + ',' + str(A.y) + ',' + str(A.z) + ')'
def myTxt_XYZ(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.4mm'
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)
myTxt_XYZ(O,"O"+myXYZ2Txt(O))
myTxt_XYZ(A,"A"+myXYZ2Txt(A))
myTxt_XYZ(B,"B"+myXYZ2Txt(B))
myTxt_XYZ(C,"C"+myXYZ2Txt(C))
myTxt_XYZ(D,"D"+myXYZ2Txt(D))
myTxt_XYZ(P,"P"+myXYZ2Txt(P))
myTxt_XYZ(Q,"Q"+myXYZ2Txt(Q))
myTxt_XYZ(R,"R"+myXYZ2Txt(R))
myLine(A,B)
myLine(C,D)
myLine(P,Q)
myLine(B,R)
myLine(D,R)
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")