(七-3) 奈良県大和郡山市小泉町 耳成山口神社 嘉永7年(1854)
近畿数学史学会：近畿の算額「数学の絵馬を訪ねて」，平成4年5月16日 初版第一刷，大阪教育図書株式会社，大阪市.
奈良県大和郡山市 庚申堂 明治13年(1880)
和算の館
http://www.wasan.jp/nara/kosindo.html
牧下英世：数学教育を通して取り組んだ総合的な学習とその実証的な研究-算額を―算額を用いた課題学習とそのフィールドワークの実践から―，2002筑波大学附属駒場論集第42集
https://tsukuba.repo.nii.ac.jp/record/6466/files/13.pdf
キーワード：累円
#Julia #SymPy #算額 #和算 #数学
大円と小円が交差してできる隙間に甲円から始まる累円を容れる。
大円，小円，甲円の直径をそれぞれ 7 尺 2 寸，6 尺 1 寸，1 尺 7 寸としたとき，乙円，丙円，丁円，戊円の直径を求めよ。
大円の半径と中心座標を \(R_1,\ (0,\ 0)\)
小円の半径と中心座標を \(R_2,\ (0,\ 2r_1 - R_1 + R_2)\)
甲円の半径と中心座標を \(r_1,\ (0,\ r_1 - R_1)\)
乙円の半径と中心座標を \(r_2,\ (x_2,\ y_2)\)
丙円の半径と中心座標を \(r_3,\ (x_3,\ y_3)\)
:
とおき，以下の連立方程式を解く。
include("julia-source.txt"); # julia-source.txt ソース
using SymPy
@syms R1, R2, r1, r2, x2, y2
(R1, R2, r1) = (72, 61, 17) .// 2
eq1 = x2^2 + y2^2 - (R1 -r2)^2
eq2 = x2^2 + (2r1 - R1 + R2 - y2)^2 - (R2 + r2)^2
eq3 = x2^2 + (y2 - r1 + R1)^2 - (r1 + r2)^2
#(r2, x2, y2) =
solve([eq1, eq2, eq3], (r2, x2, y2))[2]
(36465/4681, 204*sqrt(130845)/4681, -109509/4681)
漸化式的に累円のパラメータを数式で求めることが SymPy では困難なので，関数内で solve() を使う，半自動的なやり方を取る。
一つの累円のパラメータが得られたら，それを関数の引数に指定して，次の累円のパラメータを得る。
注意点として（Julia の SymPy 特有であるが），xx/yy の分数は xx//yy にすることと，sqrt(zz) は sqrt(Sym(zz)) にすること。そのようにしなくても良いが，処理時間がかかる。
function nextcircle(R1, R2, r1, r2, x2, y2)
@syms r3, x3, y3
eq4 = x3^2 + y3^2 - (R1 -r3)^2
eq5 = x3^2 + (2r1 - R1 + R2 - y3)^2 - (R2 + r3)^2
eq6 = (x3 - x2)^2 + (y3 - y2)^2 - (r2 + r3)^2;
solve([eq4, eq5, eq6], (r3, x3, y3))
end;
# r3
nextcircle(R1, R2, r1, 36465//4681, 204*sqrt(Sym(130845))/4681, -109509//4681)
2-element Vector{Tuple{Sym{PyCall.PyObject}, Sym{PyCall.PyObject}, Sym{PyCall.PyObject}}}:
(8690825/1404066, 50932*sqrt(130845)/702033, -19854559/1404066)
(17/2, 0, -55/2)
# r4
nextcircle(R1, R2, r1, 8690825//1404066, 50932*sqrt(Sym(130845))/702033, -19854559//1404066)
2-element Vector{Tuple{Sym{PyCall.PyObject}, Sym{PyCall.PyObject}, Sym{PyCall.PyObject}}}:
(18641819625/4108026529, 353486916*sqrt(130845)/4108026529, -18850643421/4108026529)
(36465/4681, 204*sqrt(130845)/4681, -109509/4681)
# r5
nextcircle(R1, R2, r1, 18641819625//4108026529, 353486916*sqrt(Sym(130845))/4108026529, -18850643421//4108026529)
2-element Vector{Tuple{Sym{PyCall.PyObject}, Sym{PyCall.PyObject}, Sym{PyCall.PyObject}}}:
(4442967010625/1379049831714, 62216188328*sqrt(130845)/689524915857, 4167487120889/1379049831714)
(8690825/1404066, 50932*sqrt(130845)/702033, -19854559/1404066)
# r6
nextcircle(R1, R2, r1, 4442967010625//1379049831714, 62216188328*sqrt(Sym(130845))/689524915857, 4167487120889//1379049831714)
2-element Vector{Tuple{Sym{PyCall.PyObject}, Sym{PyCall.PyObject}, Sym{PyCall.PyObject}}}:
(9530164237790625/4195069355668441, 378487752059964*sqrt(130845)/4195069355668441, 35723160673770891/4195069355668441)
(18641819625/4108026529, 353486916*sqrt(130845)/4108026529, -18850643421/4108026529)
# r7
nextcircle(R1, R2, r1, 9530164237790625//4195069355668441, 378487752059964*sqrt(Sym(130845))/4195069355668441, 35723160673770891//4195069355668441)
2-element Vector{Tuple{Sym{PyCall.PyObject}, Sym{PyCall.PyObject}, Sym{PyCall.PyObject}}}:
(2271355810006765625/1413029171738162306, 62689304552989212*sqrt(130845)/706514585869081153, 17460791512813260881/1413029171738162306)
(4442967010625/1379049831714, 62216188328*sqrt(130845)/689524915857, 4167487120889/1379049831714)
# r8
nextcircle(R1, R2, r1, 2271355810006765625//1413029171738162306, 62689304552989212*sqrt(Sym(130845))/706514585869081153, 17460791512813260881//1413029171738162306)
2-element Vector{Tuple{Sym{PyCall.PyObject}, Sym{PyCall.PyObject}, Sym{PyCall.PyObject}}}:
(4872058212464512265625/4255498594259851496689, 369959068275411941556*sqrt(130845)/4255498594259851496689, 63967589464505474522739/4255498594259851496689)
(9530164237790625/4195069355668441, 378487752059964*sqrt(130845)/4195069355668441, 35723160673770891/4195069355668441)
累円の直径と中心座標は以下のようになる。算額の答えとはずいぶん異なる。
乙円: 直径 = 15.58000, (15.7641, -23.3944)
丙円: 直径 = 12.37951, (26.2429, -14.1408)
丁円: 直径 = 9.07580, (31.1257, -4.58873)
戊円: 直径 = 6.44352, (32.6386, 3.02200)
己円: 直径 = 4.54351, (32.6356, 8.51551)
庚円: 直径 = 3.21487, (32.0960, 12.35700)
辛円: 直径 = 2.28977, (31.4472, 15.03170)
描画関数プログラムのソースを見る
function draw(more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
names = Char["甲乙丙丁戊己庚辛壬癸"...]
(R1, R2, r1) = (72, 61, 17) .// 2
plot(showaxis=false)
delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3 # size[2] * fontsize * 2
circlef(0, 0, R1, :lightgreen)
circlef(0, 2r1 - R1 + R2, R2, :cadetblue1)
circle(0, 2r1 - R1 + R2, R2, :lightslateblue)
circle(0, 0, R1, :mediumseagreen)
circlef(0, r1 - R1, r1, :gold)
number = 1
for (r, x, y) in [
(7.790002136295663, 15.764133064156171, -23.394360179448835)
(6.189755324892134, 26.242896581838256, -14.140759052637128)
(4.537901470061319, 31.12566720233952, -4.588734587745892)
(3.221759582902747, 32.63863611252851, 3.021998933649333)
(2.2717536779012537, 32.63557363626892, 8.515511340831878)
(1.607437309459633, 32.095998094235, 12.356992949646468)
(1.1448854005113114, 31.447186145793232, 15.031749640521548)]
number += 1
circlef(x, y, r, number)
circlef(-x, y, r, number)
point(x, y, names[number], :white, :center, :vcenter, mark=false)
@printf("%s円: 直径 = %8.5f, (%g, %g)\n", names[number], 2r, x, y)
end
if more
#hline!([0], color=:gray80, lw=0.5)
#vline!([0], color=:gray80, lw=0.5)
point(0, 0, "大円:R1,(0,0)", :red, :center, delta=-1)
point(0, 2r1 - R1 + R2, "小円:R2,(0,2r1-R1+R2)", :green, :center, delta=-1)
point(0, r1 - R1, " 甲円:r1\n(0,r1-R1)", :white, :center, delta=-1)
end
end;
以下のアイコンをクリックして応援してください
