LoginSignup
0
0
@mrrclb48z

アポロニウスの球？「2023岡山大学前期数学I・数学II・数学A・数学B第3問」をsympyとFreeCADでやってみたい。

Last updated at Posted at 2023-04-13

sympyは、2次式にならなくていいね。と思いながら、コーディングしていると、
だんだん雲行きが怪しくなってきました。アドバイスをお願いします。
①b=0? a^2=...?
②作図もしてみたい？

修正してみた。2023/05/10

オリジナル
T氏様

kamelink様

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

sympyで(オリジナルのやり方)

勉強中。

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

怪しいです。

from sympy import *
r    =symbols('r'    ,real=True,positive=True)
a,b,c=symbols('a b c',real=True)
cc   =symbols('cc'   ,real=True)
O=Point(0, 0,0)
A=Point(1, 1,0)
B=Point(1,-1,0)
P=Point(a, b,c)
#
r=1
ans=solve([
          Eq( A.distance(P),r*O.distance(P) ),
          Eq( B.distance(P),r*O.distance(P) )
          ],[a,b,c])[0]
OP=Point(ans[0],ans[1],c).distance(O)
c2=solve(Eq(OP, minimum(OP,c,Interval(-oo,oo)) ))[0]
print("#(1)",ans[0],ans[1],c2)
#########################################################
r=sqrt(3)/2
ans=solve([
           Eq( A.distance(P),r*O.distance(P) ),
           Eq( B.distance(P),r*O.distance(P) )
          ],[a,b,c])
amin=min(
         minimum(ans[0][0],c,Interval(-oo,oo)),
         minimum(ans[1][0],c,Interval(-oo,oo)) 
        )
amax=max(
         maximum(ans[0][0],c,Interval(-oo,oo)),
         maximum(ans[1][0],c,Interval(-oo,oo)) 
        )
print("#(2)",amin,",",amax)          
#########################################################
print("#   b =",ans[0][1],ans[1][1])          
c01=solve(Eq( A.distance(P),r*O.distance(P) ),c)          
print("#   a2=",c01[0].subs({b:0})**2)
print("#   a2=",c01[1].subs({b:0})**2)
mySub={b:ans[0][1],c:c01[0].subs({b:0})}
#########################################################
OPAP=P.dot(P-A).subs(mySub).expand()
print("#(3)",          
         maximum(OPAP,a,Interval(amin,amax)),",",
         minimum(OPAP,a,Interval(amin,amax))
     ) 
#(1) 1 0 0
#(2) 4 - 2*sqrt(2) , 2*sqrt(2) + 4
#   b = 0 0
#   a2= -a**2 + 8*a - 8
#   a2= -a**2 + 8*a - 8
#(3) 14*sqrt(2) + 20 , 20 - 14*sqrt(2)

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

FreeCADのマクロで作図(r=1)

import FreeCAD
import Part
import Draft
import Mesh
#########################################################################################################
# 計算
from sympy import *
r    =symbols('r'    ,real=True,positive=True)
a,b,c=symbols('a b c',real=True)
cc   =symbols('cc'   ,real=True)
O=Point(0, 0,0)
A=Point(1, 1,0)
B=Point(1,-1,0)
P=Point(a, b,c)
#
r=1
ans=solve([
          Eq( A.distance(P),r*O.distance(P) ),
          Eq( B.distance(P),r*O.distance(P) )
          ],[a,b,c])[0]
OP=Point(ans[0],ans[1],c).distance(O)
c2=solve(Eq(OP, minimum(OP,c,Interval(-oo,oo)) ))[0]
print("#(1)",ans[0],ans[1],c2)
#########################################################
PP=Point( ans[0],ans[1],c2 )
##########################################################################################################
### 3D作図
def my_Circle(x,y,z,a):
    zaxis = App.Vector(0, 0, 1)
    p5 = App.Vector(float(x),float(y),float(z))
    place5 = App.Placement(p5, App.Rotation(zaxis, 0))
    circle5 = Draft.make_circle(float(a), placement=place5)
def my_Sphere(myx,myy,myz,myr):
    sphere = doc.addObject("Part::Sphere", "mySphere2")
    sphere.Radius = float(myr)
    sphere.Angle1 = -90
    sphere.Angle2 =  90
    sphere.Angle3 = 360
    sphere.Placement = App.Placement(App.Vector(float(myx),float(myy),float(myz)), App.Rotation(0, 90, 0))
    FreeCADGui.getDocument('Unnamed').getObject('mySphere2').Transparency = 90
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.1 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
def myLine_H(*args):
    for i in range(1,len(args)):
        myLine(args[0],args[i])
    return
#
myTxtXYZ_S(PP,"P",A,"A",B,"B",O,"O")
myLine_H  (PP, A,B,O)
#
doc = App.activeDocument()
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")

0png.png

????? FreeCADのマクロで作図(r=sqrt(3)/2)

何がなんだか、
わかりませんでした。

import FreeCAD
import Part
import Draft
import Mesh
#########################################################################################################
# 計算
from sympy import *
r    =symbols('r'    ,real=True,positive=True)
a,b,c=symbols('a b c',real=True)
cc   =symbols('cc'   ,real=True)
O=Point(0, 0,0)
A=Point(1, 1,0)
B=Point(1,-1,0)
P=Point(a, b,c)
#
r=1
ans=solve([
          Eq( A.distance(P),r*O.distance(P) ),
          Eq( B.distance(P),r*O.distance(P) )
          ],[a,b,c])[0]
OP=Point(ans[0],ans[1],c).distance(O)
c2=solve(Eq(OP, minimum(OP,c,Interval(-oo,oo)) ))[0]
print("#(1)",ans[0],ans[1],c2)
#########################################################
r=sqrt(3)/2
ans=solve([
           Eq( A.distance(P),r*O.distance(P) ),
           Eq( B.distance(P),r*O.distance(P) )
          ],[a,b,c])
amin=min(
         minimum(ans[0][0],c,Interval(-oo,oo)),
         minimum(ans[1][0],c,Interval(-oo,oo)) 
        )
amax=max(
         maximum(ans[0][0],c,Interval(-oo,oo)),
         maximum(ans[1][0],c,Interval(-oo,oo)) 
        )
print("#(2)",amin,",",amax)          
#########################################################
print("#   b =",ans[0][1],ans[1][1])          
c01=solve(Eq( A.distance(P),r*O.distance(P) ),c)          
print("#   a2=",c01[0].subs({b:0})**2)
print("#   a2=",c01[1].subs({b:0})**2)
#########################################################
#?????????????????????????????????????????????????
myEq=Eq((a-1)**2+c**2+1,(a**2+c**2)*Rational(3,4))
myEq1=myEq.subs({a:amin})
ans=solve(myEq1,c)
Pmin=Point(amin,0,ans[0])
myEq2=myEq.subs({a:(amin+amax)/2})
ans=solve(myEq2,c)
Pminmax=Point((amin+amax)/2,0,ans[0])
myEq3=myEq.subs({a:amax})
ans=solve(myEq3,c)
Pmax=Point(amax,0,ans[0])
#print(Pmin)
#print(Pminmax)
#print(Pmax)
##########################################################################################################
### 3D作図
def my_Circle(x,y,z,a):
    zaxis = App.Vector(0, 0, 1)
    p5 = App.Vector(float(x),float(y),float(z))
    place5 = App.Placement(p5, App.Rotation(zaxis, 0))
    circle5 = Draft.make_circle(float(a), placement=place5)
def my_Sphere(myx,myy,myz,myr):
    sphere = doc.addObject("Part::Sphere", "mySphere2")
    sphere.Radius = float(myr)
    sphere.Angle1 = -90
    sphere.Angle2 =  90
    sphere.Angle3 = 360
    sphere.Placement = App.Placement(App.Vector(float(myx),float(myy),float(myz)), App.Rotation(0, 90, 0))
    FreeCADGui.getDocument('Unnamed').getObject('mySphere2').Transparency = 90
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.1 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
def myLine_H(*args):
    for i in range(1,len(args)):
        myLine(args[0],args[i])
    return
#
myTxtXYZ_S(A,"A",B,"B",O,"O")
myLine_S  (A,B)
#
myTxtXYZ_S(Pmin,"Pmin",Pminmax,"Pminmax",Pmax,"Pmax")
myLine_H  (Pmin,    A,B,O)
myLine_H  (Pminmax, A,B,O)
myLine_H  (Pmax,    A,B,O)
#
doc = App.activeDocument()
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")

0png.png

参考

勉強中。
3Dプロット参考にしました。ありがとうございました。
2023岡山大学前期【３】解答＜岡山進研学院様

アポロニウスの球？

？勉強

修正してみた。2023/05/10

オリジナル(youtube)

from sympy import *
a,b,c,r=symbols('a b c r',real=True)
O=Point(0, 0,0)
A=Point(1, 1,0)
B=Point(1,-1,0)
P=Point(a ,b,c)
AP=Point(A).distance(P)
BP=Point(B).distance(P)
OP=Point(O).distance(P)
#
r_inp=1
ans=solve([Eq(AP,BP),Eq(BP,r*OP).subs({r:r_inp})],[a,b,c])
P_ans=P.subs({a:ans[0][0],b:ans[0][1]})
OP_ans=Point(O).distance(P_ans)
P_ans=P_ans.subs({c:solve(Eq(OP_ans, minimum(OP_ans,c,Interval(-oo,oo)) ))[0]})
print("#(1)",P_ans)
#
r_inp=Rational(1,2)*sqrt(3)
ans=solve([Eq(AP,BP),Eq(BP,r*OP).subs({r:r_inp})],[a,b,c])
mySub0={a:ans[0][0],b:ans[0][1]}
mySub1={a:ans[1][0],b:ans[1][1]}
P0=P.subs(mySub0)
P1=P.subs(mySub1)
# plot(P0.x,P1.x)
a_min=min( minimum(P0.x,c,Interval(-oo,oo)) , minimum(P1.x,c,Interval(-oo,oo))  )
a_max=max( maximum(P0.x,c,Interval(-oo,oo)) , maximum(P1.x,c,Interval(-oo,oo))  )
print("#(2)",a_min,",",a_max) 
dot_min=min( minimum((P0-O).dot(P0-A),c,Interval(-oo,oo)) , minimum((P1-O).dot(P1-A),c,Interval(-oo,oo))  ).simplify()
dot_max=max( maximum((P0-O).dot(P0-A),c,Interval(-oo,oo)) , maximum((P1-O).dot(P1-A),c,Interval(-oo,oo))  ).simplify()
print("#(3)",dot_min,",",dot_max) 
#(1) Point3D(1, 0, 0)
#(2) 4 - 2*sqrt(2) , 2*sqrt(2) + 4
#(3) 20 - 14*sqrt(2) , 14*sqrt(2) + 20

？？？？plotコマンドは、縦横比？異なります。

0png.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