(途中)外サイクロイド<2/2>点Pが2周します。「2022 産業医科大学 医学部【3】」をWolframAlphaとsympyでやってみたい。

(2023-10-23タイトル変更済み)

前頁<1/2>より

sympyで(抜粋)"セ"から"ツ"までです。(作図 0≦t≦2.0に変更しました。)

　・ｎ本の折れ線を関数にしたいです。

・sympy.plotting.plot.plot_parametric(*args, show=True, **kwargs)
https://docs.sympy.org/latest/modules/plotting.html#sympy.plotting.plot.plot_parametric

# ver0.２
from sympy import *
var('f,a,b',real=True,nonnegative=True)
var('θ1,θ2',real=True,nonnegative=True)
var('ω1,ω2',real=True,nonnegative=True)
var('t,Φ'  ,real=True,nonnegative=True)
(Tp,a,b,f)   =(1,2,1,5)
(st0,st1,st2)=(0,1,2)
O=Point(0,0)
# ------------------------------------------------------------------
P=(Point(f,0)+Point(a*cos(θ1),a*sin(θ1)))
S=(P         +Point(b*cos(θ2),b*sin(θ2)))
print("#セ",P)
print("#ソ",S)
# ------------------------------------------------------------------
rep_θt={θ1:ω1*t,θ2:ω2*t+Φ}
print("#タ",solve(Eq(2*pi,ω1*Tp),Tp)[0])
# ------------------------------------------------------------------
(P,S)=(P.subs(rep_θt),S.subs(rep_θt))
PS   =O.distance(P)-O.distance(S)
Φ_sol=solveset( Eq(PS.subs({t:0}),1),Φ,Interval(0,2*pi)).args[0]
ω1_sol=(            2 *pi)/Tp
ω2_sol=(Rational(10,2)*pi)/Tp   #回転2.5倍(=5÷2)
print("#チ",ω1_sol)
print("#チ",ω2_sol)
print("#ツ", Φ_sol)
# 図1" ########################################################################
P=(Point(f,0)+Point(a*cos(θ1),a*sin(θ1))).subs(rep_θt).subs({ω1:ω1_sol,ω2:ω2_sol,Φ:Φ_sol})
S=(P         +Point(b*cos(θ2),b*sin(θ2))).subs(rep_θt).subs({ω1:ω1_sol,ω2:ω2_sol,Φ:Φ_sol})
ax1=plot_parametric(P.x,P.y,(t,st0,st1),aspect_ratio=(1.0,1.0),
                    # line_color="blue",
                    show=False
                  )
ax2=plot_parametric(S.x,S.y,(t,st0,st1),
                    # line_color="blue",
                    show=False
                  )
(st1,st2)=(1.0,2.0)
ax3=plot_parametric(S.x,S.y,(t,st1,st2),
                    # line_color="blue",
                    show=False
                  )
# ------------------------------------------------------------------
def myPlotLineOPS1(tOP,O,P,S):
    var('u'  ,real=True)
    OPu=O+u*Point(P.x,P.y).subs({t:tOP})
    PSu=  1*Point(P.x,P.y).subs({t:tOP})+u*Point(S.x-P.x,S.y-P.y).subs({t:tOP})
    ay1=plot_parametric(PSu.x,PSu.y,(u,0,1),line_color="black",
                        show=False
                       )
    return ay1
def myPlotLineOPS2(tOP,O,P,S):
    var('u'  ,real=True)
    OPu=O+u*Point(P.x,P.y).subs({t:tOP})
    PSu=  1*Point(P.x,P.y).subs({t:tOP})+u*Point(S.x-P.x,S.y-P.y).subs({t:tOP})
    ay2=plot_parametric(PSu.x,PSu.y,(u,0,1),line_color="black",
                        show=False
                       )
    return ay2
def myPlotLineOPS4(tOP,O,P,S):
    var('u'  ,real=True)
    OPu=O+u*Point(P.x,P.y).subs({t:tOP})
    PSu=  1*Point(P.x,P.y).subs({t:tOP})+u*Point(S.x-P.x,S.y-P.y).subs({t:tOP})
    ay4=plot_parametric(PSu.x,PSu.y,(u,0,1),line_color="black",
                        show=False
                       )
    return ay4
def myPlotLineOPS6(tOP,O,P,S):
    var('u'  ,real=True)
    OPu=O+u*Point(P.x,P.y).subs({t:tOP})
    PSu=  1*Point(P.x,P.y).subs({t:tOP})+u*Point(S.x-P.x,S.y-P.y).subs({t:tOP})
    ay6=plot_parametric(PSu.x,PSu.y,(u,0,1),line_color="black",
                        show=False
                       )
    return ay6
ax1.extend(ax2)
ax1.extend(ax3)
ax1.extend(myPlotLineOPS1(0.1,O,P,S))
ax1.extend(myPlotLineOPS2(0.2,O,P,S))
ax1.extend(myPlotLineOPS4(0.4,O,P,S))
ax1.extend(myPlotLineOPS6(0.6,O,P,S))
ax1.show()
# 図2" ########################################################################
Pt=(O.distance(P)).subs({ω1:ω1_sol})
St=(O.distance(S)).subs({ω1:ω1_sol,ω2:ω2_sol,Φ:Φ_sol})

tHani=(t,st0,st2)
pl_Pt=plot(Pt,tHani,show=False,aspect_ratio=(5.0,1.0),ylim=(2,9))
pl_St=plot(St,tHani,show=False)
pl_Pt.extend(pl_St)
pl_Pt.show()
#セ Point2D(2*cos(θ1) + 5, 2*sin(θ1))
#ソ Point2D(2*cos(θ1) + cos(θ2) + 5, 2*sin(θ1) + sin(θ2))
#タ 2*pi/ω1
#チ 2*pi
#チ 5*pi
#ツ pi

図1"　

・点Pが2周すると、t=0の状態に戻ります。
　青色 点Ｐ (0.0≦t≦2.0)
　緑色 点Ｓ (0.0≦t≦1.0) 　
　橙色 〃　 (1.0≦t≦2.0)

図2"

<1/2>のコメントと同じです。ありがとうございました。
　(0≦t≦2.0)

前頁<1/2>へ