0
0

三角形と円周角「2010 京都大学前期 理系甲,乙共通【3】」をChatGPTとWolframAlphaとsympyとFreeCADでやってみたい。

Last updated at Posted at 2023-09-27

・問題文は2次元ですが、3次元FreeCADのマクロで、XY平面上に円も、作図しました。

オリジナル 
shaitan 様より

望星塾 様より

上と同じです。大学入試数学問題集成 様>テキスト 理系甲,乙共通【3】

ChatGPT-3.5で(???できませんでした。???)

私の質問の仕方が悪いかもしれません。
入力文

xを正の実数とする.座標平面上の3点A(0,1),B(0,2),P(x,x)をとり,△APBを考える.
xの値が変化するとき,∠APBの最大値を求めよ.

以下の結果は抜粋です。

...
手計算では複雑な式となるため、数値計算を行うか、数値計算ソフトウェアを使用して最大値を
求めることがおすすめです。

WolframAlphaで

省略

sympyで(いつものスタイルで)

・angle_between関数とmaximum関数の組み合わせです。
・solveで終了しなかったので、solvesetにしました。
angle_between(l2)
https://docs.sympy.org/latest/modules/geometry/lines.html#sympy.geometry.line.LinearEntity.angle_between
maximum
https://docs.sympy.org/latest/modules/calculus/index.html#:~:text=*pi%5D)-,sympy,-.calculus.util.
solveset
https://docs.sympy.org/latest/modules/solvers/solveset.html

from sympy import *
var('x',real=True,positive=True)
A=Point(0,1)
B=Point(0,2)
P=Point(x,x)
ang_APB=Line(P,B).angle_between(Line(P,A))
max_APB=maximum(ang_APB,x,Interval(0,pi))
# print("#",max_APB,solve   ( Eq(ang_APB,max_APB),x                     )        )
print  ("#",max_APB,solveset( Eq(ang_APB,max_APB),x,Interval.Lopen(0,oo)).args[0])
# pi/4 {1}

上と同じです。mapにしただけです。 ← 資源の無駄使いするな。最初から、書けよ!!!です。

from sympy import *
var('x',real=True,positive=True)
A,B,P=map(Point,[(0,1),(0,2),(x,x)])
ang_APB=Line(P,B).angle_between(Line(P,A))
max_APB=maximum(ang_APB,x,Interval(0,pi))
print("#",max_APB,solveset( Eq(ang_APB,max_APB),x,Interval.Lopen(0,oo)).args[0])
# pi/4 1

sympyで(shaitan 様の方法を参考に)

・shaitan 様の解答に近づけました。円の作図のため、FreeCADのマクロへコピー利用しました。
projection(other)
https://docs.sympy.org/latest/modules/geometry/lines.html#sympy.geometry.line.LinearEntity.projection
distance(other)
https://docs.sympy.org/latest/modules/geometry/points.html#sympy.geometry.point.Point.distance

from sympy import *
var('x',real=True,positive=True)
A,B,P=map(Point,[(0,1),(0,2),(x,x)])
H=Line(Point(0,0),Point(x,x)).projection(B)
print("#",H)
print("#",Ge( maximum(B.distance(P),x,Interval(0,pi)) , B.distance(H)))
print("#",Line(H,A).angle_between(Line(H,B)) )
# Point2D(1, 1)
# True
# pi/4

sympyで(望星塾 様の方法を参考に)

・Eqを使わなくてもいいような気もしましたが、Eqだと、solveに放り込めるから、
 便利に思いました。
Derivatives
https://docs.sympy.org/latest/tutorials/intro-tutorial/calculus.html#derivatives

from sympy import *
var('x'      ,real=True,positive=True)
var('θ,θA,θB',real=True,positive=True)
var('f'     ,real=True,positive=True)
A,B,P=map(Point,[(0,1),(0,2),(x,x)])
eq1=Eq(tan(θ) ,tan(θA-θB))                         ;print("#",eq1)
eq2=Eq(eq1.lhs,eq1.rhs.expand(trig=True))          ;print("#",eq2)
mySubs={tan(θA):(x-1)/x,tan(θB):(x-2)/x}       
eq3=eq2.subs(mySubs)                               ;print("#",eq3)
eq4=Eq(eq3.lhs,eq3.rhs.simplify())                 ;print("#",eq4)
eq5=Eq(f,(denom(eq4.rhs)/numer(eq4.rhs)).expand()) ;print("#",eq5)
eq6=eq3.subs({x:solve(diff(eq5.rhs))[0]})          ;print("#",eq6)
print("#",solve(eq6)[0])
# Eq(tan(θ), tan(θA - θB))
# Eq(tan(θ), tan(θA)/(tan(θA)*tan(θB) + 1) - tan(θB)/(tan(θA)*tan(θB) + 1))
# Eq(tan(θ), -(x - 2)/(x*(1 + (x - 2)*(x - 1)/x**2)) + (x - 1)/(x*(1 + (x - 2)*(x - 1)/x**2)))
# Eq(tan(θ), x/(x**2 + (x - 2)*(x - 1)))
# Eq(f, 2*x - 3 + 2/x)
# Eq(tan(θ), 1)
# pi/4

<解説>部分です。
factor
https://docs.sympy.org/latest/tutorials/intro-tutorial/simplification.html#factor

from sympy import *
var('x'           ,real=True,positive=True)
var('θ'            ,real=True,positive=True)
var('PA,PB,AB,ABP',real=True,positive=True)
eq1=Eq(cos(θ),(PA**2+PB**2-AB**2)/(2*PA*PB)         )    ;print("#",eq1)
eq2=Eq(sin(θ),2*ABP/(PA*PB)                         )    ;print("#",eq2)
eq3=Eq((eq2.lhs/eq1.lhs).simplify(),eq2.rhs/eq1.rhs )    ;print("#",eq3)
mySubs={PA**2:x**2+(x-1)**2,
        PB**2:x**2+(x-2)**2,
        AB**2:1,
        ABP  :x/2 
       }
eq4=eq3.subs(mySubs)                                     ;print("#",eq4)
eq5=Eq(eq4.lhs,factor(eq4.rhs)    )                      ;print("#",eq5)
# Eq(cos(θ), (-AB**2 + PA**2 + PB**2)/(2*PA*PB))
# Eq(sin(θ), 2*ABP/(PA*PB))
# Eq(tan(θ), 4*ABP/(-AB**2 + PA**2 + PB**2))
# Eq(tan(θ), 2*x/(2*x**2 + (x - 2)**2 + (x - 1)**2 - 1))
# Eq(tan(θ), x/(2*x**2 - 3*x + 2))

FreeCADのマクロで作図

・円を書きました。
・計算部分は、shaitan 様のソースコードをコピーです。

import FreeCAD
import Part
import DraftTools
import Draft
import Mesh
############################################################################
# 計算
from sympy import *
var('x',real=True,positive=True)
A,B,P=map(Point,[(0,1),(0,2),(x,x)])
H=Line(Point(0,0),Point(x,x)).projection(B)
print("#",H)
print("#",Ge( maximum(B.distance(P),x,Interval(0,pi)) , B.distance(H)))
print("#",Line(H,A).angle_between(Line(H,B)) )
# Point2D(1, 1)
# True
# pi/4
############################################################################
# 作図用
O=Point(0,0)
L=Point(3,3)
C=(B+H)/2
ci=Circle(C,C.distance(H))
############################################################################
# 3D作図 z=0 XY平面に作図しました。
############################################################################
############################################################################
# 円の作図 FrecCADのdocより
# https://wiki.freecad.org/Macro_Circle
def Freecad3D_circle(x=0.0,y=0.0,z=0.0,radius=0.0,diameter=0.0,circumference=0.0,area=0.0,startangle=0.0,endangle=0.0,arc=0.0,anglecenter=0.0,cord=0.0,arrow=0.0,center=0,placemObject=""):
    from math import sqrt, pi
    if placemObject == "":
        pl = FreeCAD.Placement()
        pl.Rotation = FreeCADGui.ActiveDocument.ActiveView.getCameraOrientation()   
        pl.Base = FreeCAD.Vector(x,y,z)
    else:
        pl = FreeCAD.Placement()
        pl = placemObject                                  # placement imposted
    if diameter != 0:                                      # with diameter
        radius = diameter / 2.0
    elif circumference != 0:                               # with circumference
        radius = (circumference / pi) / 2.0
    elif area != 0:                                        # with area
        radius =  sqrt((area / pi))
    elif (cord != 0) and (arrow != 0):                     # with cord and arrow
        radius = ((arrow**2) + (cord**2) / 4.0) / (arrow*2) 
    elif (arc != 0) and (anglecenter != 0):                # with arc and anglecenter central in degrees
        radius = ((360/anglecenter)*arc) / pi/2.0
        if endangle != 0:
            startangle  = endangle - anglecenter
        endangle   = anglecenter + startangle
        startangle  = endangle - anglecenter
    if radius != 0:
        try:
            Draft.makeCircle(radius,placement=pl,face=False,startangle=startangle,endangle=endangle,support=None)
            if center != 0:
                x=pl.Base.x
                y=pl.Base.y
                z=pl.Base.z
                Draft.makePoint(x,y,z)
        except Exception:
            App.Console.PrintError("Unexpected keyword argument" + "\n")
        App.ActiveDocument.recompute()
    else:
        App.Console.PrintMessage("Unexpected keyword argument" + "\n")
        App.Console.PrintMessage("circle(x,y,z,radius,diameter,circumference,area,startangle,endangle,[arc,anglecenter],[cord,arrow],center,placemObject)" + "\n")
        App.Console.PrintMessage("circle(radius=10.0,placemObject=App.Placement(App.Vector(11,20,30), App.Rotation(30,40,0), App.Vector(0,0,0)))" + "\n")
    return
def myCircle_2D(myCi):
    x=myCi.center.x
    y=myCi.center.y
    r=myCi.radius
    Freecad3D_circle(
           x=float(x),y=float(y),z=0.0,
           radius=float(abs(r)),
           center=1,
           placemObject=App.Placement(App.Vector(float(x),float(y),0),
           App.Rotation(0,0,0),App.Vector(0,0,0)))
    return
############################################################################
def myXYZ2Txt_2D(A):
    return '(' + str(A.x) + ',' + str(A.y) +  ')'
def myTxtXYZ_2D(A,myWedgei):
    P5x=float(A.x)
    P5y=float(A.y)
    P5z=0.0
    p5 = FreeCAD.Vector(P5x, P5y, P5z)
    myText = Draft.makeText(myWedgei, p5)
    myText.Label = myWedgei
    FreeCADGui.ActiveDocument.ActiveObject.FontSize = '0.15 mm'
    return
def myTxtXYZ_S_2D(*xy_tx):
    for i in range(1,int(len(xy_tx)/2)+1):
        myTxtXYZ_2D(xy_tx[2*i-2],xy_tx[2*i-1]+myXYZ2Txt_2D(xy_tx[2*i-2]) )
    return
def myLine_2D(A,B):
    Ax,Ay,Az=float(A.x),float(A.y),0.0
    Bx,By,Bz=float(B.x),float(B.y),0.0
    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_2D(*args):
    for i in range(1,len(args)):
        myLine_2D(args[i-1],args[i])
    return
def myLine_C_2D(*args):
    for i in range(1,len(args)):
        myLine_2D(args[i-1],args[i])
    myLine_2D(args[i],args[0])
    return
#######################################################################################
myCircle_2D  ( ci )
myLine_S_2D  ( O,L )
myLine_C_2D  ( A,B,H )
myTxtXYZ_S_2D( O,"O",A,"A",B,"B",H,"H",C,"C"  )
#######################################################################################
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方向?です。
0png.png

・isometric方向?です。拡大図です。
1png.png

・上空から?です。Originは、CAD操作で非表示です。
2png.png

sympyの実行環境

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

参考

以下、いつもの?おすすめです。

私がやってみた感想

>∠APBが最大となるのは,△APBの外接円の直径dが最小となるとき...(shaitan 様より)
>図を書いてみると,4APBが45,のとき...ABPが直角二等辺三角形のとき,...(望星塾 様より)
・ABより、円周角の角の点が?向こう側(直線側)へ行くと細くなるから、垂線の足で角度が最大
です。
意味が通じます。でしょうか?

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