高松市松縄町 元祖わかめうどん 大島家
今月（今週?）は，わかめを練り込んだうどんと，伊予柑を練り込んだうどんの二種盛りが基本。それを他のメニューに使うということで，きつねうどんを頼んだ。
やや細麺で，こしがある。わかめも伊予柑も仄かな香りがある。
店先の瀬戸物に湛えられた水は野鳥のためだという...
算額（その1461）
七十五 群馬県吾妻町三島 馬頭観世音 嘉永4年(1851)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：円9個，正方形2個，斜線8本
#Julia, #SymPy, #算額, #和算
算題の文字は判読できないということであるが，「倒立した正方形の中に正方形を容れ，8 個の頂点を 8 本の斜線で結ぶ。区画された領域に，大円 1 個，中円 4 個，小円 4 個を容れる。」ということであろう。
一般性を損なわずに，倒立した正方形の対角線の長さを 2a
大円の半径と中心座標を r1, (0, 0)
中円の半径と中心座標を r2, (r1 + r2, 0)
小円の半径と中心座標を r3, (x3, x3)
とおき，以下の連立方程式を解く。
a は図全体の大きさであり，a が違う図形は全て相似である。
図の見た目は r1 の違いにより変化する。
include("julia-source.txt");
# # julia-source.txt ソース https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf
using SymPy
@syms a::positive, r1::positive, r2::positive, r3::positive, x3::positive
eq1 = dist2(r1, r1, a, 0, r1 + r2, 0, r2)
eq2 = dist2(r1, r1, a, 0, x3, x3, r3)
eq3 = dist2(0, a, a, 0, x3, x3, r3);
res1 = solve(eq1, r2)[2];
res2 = solve([eq2, eq3], (r3, x3))[4]; # 4 of 4
r2 = r1*(-r1 + sqrt(a^2 - 2*a*r1 + 2*r1^2))/(a - r1)
r3 = √2*a*sqrt(3a^2 - 4a*r1 - 2a*sqrt(2a^2 - 4a*r1 + 4r1^2) + 4r1^2)/(2a - 4r1)
x3 = a*(2a - 2r1 - sqrt(2a^2 - 4a*r1 + 4r1^2))/(2a - 4r1)
倒立正方形の対角線の長さが 2，大円（中心にある円）の半径が 0.35 のとき，全てのパラメータは以下のとおりである。
a = 1; r1 = 0.35; r2 = 0.209053; r3 = 0.103781; x3 = 0.426616
r1 は 0 ≦ r1 ≦ a/2 の範囲の値を取れる（実際にはそれ以外の値も取れる）。値によりイメージは異なる。
function draw(a, r1, more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
r2 = r1*(-r1 + sqrt(a^2 - 2a*r1 + 2r1^2))/(a - r1)
r3 = √2*a*sqrt(3a^2 - 4a*r1 - 2a*sqrt(2a^2 - 4a*r1 + 4r1^2) + 4r1^2)/(2a - 4r1)
x3 = a*(2a - 2r1 - sqrt(2a^2 - 4a*r1 + 4r1^2))/(2a - 4r1)
@printf("a = %g; r1 = %g; r2 = %g; r3 = %g; x3 = %g\n", a, r1, r2, r3, x3)
p = plot([a, 0, -a, 0, a], [0, a, 0, -a, 0], color=:magenta, lw=0.5, showaxis=false)
plot!(r1.*[1, 1, -1, -1, 1], r1.*[-1, 1, 1, -1, -1], color=:brown, lw=0.5)
plot!([a, r1, 0, -r1, -a, -r1, 0, r1, a], [0, r1, a, r1, 0, -r1, -a, -r1, 0], color=:orange, lw=0.5)
circle(0, 0, r1)
circle42(0, r1 + r2, r2, :blue)
circle4(x3, x3, r3, :tomato)
point(0.5, 0.7, @sprintf("r1 = %ga", r1/a), :black, :left, mark=false)
if more
delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3 # size[2] * fontsize * 2
hline!([0], color=:gray80, lw=0.5)
vline!([0], color=:gray80, lw=0.5)
point(a, 0, "a", :magenta, :left, :bottom, delta=delta/2)
point(r1, 0, " r1", :red, :left, :vcenter)
point(r1 + r2, 0, " r1+r2", :blue, :left, :vcenter)
point(x3, x3, "r3,(x3,x3)", :tomato, :left, delta=-delta/2, deltax=-3delta)
end
return p
end;
draw(1, 0.35, true)
算額（その1460）
七十五 群馬県吾妻町三島 馬頭観世音 嘉永4年(1851)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：円3個，四分円2個，正方形
#Julia, #SymPy, #算額, #和算
算題の文字は判読できないということであるが，「正方形の中に 90 度に開いた二面の扇を配し，隙間に3個の日の丸を容れる」という趣向の額面であろう。扇の骨の部分は装飾である。
正方形の一辺の長さを a
四分円の半径と中心座標を r1, (0, 0), (a, a)
日の丸円の半径と中心座標を r2, (a - r2, r2), (r2, a - r2)
とおき，以下の連立方程式を解く。
変数は a, r1, r2 の 3 個であるが，どれか 1 つを与えて他の 2 変数の値を求める。a, r1, r2 は比例関係にあるので，たとえば a = 1 として得られた解を元に，r1, r2 が与えられたときの残りの 2 変数の比を計算しておけばよい。
include("julia-source.txt");
# # julia-source.txt ソース https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf
using SymPy
@syms a::positive, r1::positive, r2::positive
eq1 = (a - r2)^2 + r2^2 - (r1 + r2)^2
eq2 = 2(a/2)^2 - (r1 - r2)^2;
1. a を与えて r1, r2 を求める場合
res1 = solve([eq1, eq2], (r1, r2))[1]; # 1 of 2
# r1
res1[1] |> simplify |> println
a*(-1 + sqrt(2*sqrt(2) + 4))/2
# r2
res1[2] |> simplify |> println
a*(-sqrt(2) - 1 + sqrt(2*sqrt(2) + 4))/2
2. r1 を与えて a, r2 を求める場合
res2 = solve([eq1, eq2], (a, r2))[1]; # 1 of 2
# a
res2[1] |> simplify |> println
2*r1*(-2*sqrt(2) + sqrt(20 - 14*sqrt(2)) + 3)
# r2
res2[2] |> simplify |> println
r1*(-3*sqrt(2) - 2*sqrt(10 - 7*sqrt(2)) + 5)
3. r2 を与えて a, r1 を求める場合
res3 = solve([eq1, eq2], (a, r1))[2]; # 2 of 2
# a
res3[1] |> simplify |> println
2*r2*(1 + sqrt(2) + sqrt(2*sqrt(2) + 4))
# r2
res3[2] |> simplify |> println
r2*(sqrt(2) + 3 + 2*sqrt(sqrt(2) + 2))
数値の場合は，横方向に 2 個の比を取ればどれも同じになる
a を与えて r1, r2 を求める場合
a = 1, r1 = 0.806563, r2 = 0.0994562
r1 を与えて a, r2 を求める場合
a = 1.23983, r1 = 1, r2 = 0.123309
r2 を与えて a, r1 を求める場合
a = 10.0547, r1 = 8.10973, r2 = 1
function draw(arg, type, more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
if type == 1
a = arg
r1 = a*(-1 + sqrt(2*sqrt(2) + 4))/2
r2 = a*(-sqrt(2) - 1 + sqrt(2*sqrt(2) + 4))/2
elseif type == 2
r1 = arg
a = 2*r1*(-2*sqrt(2) + sqrt(20 - 14*sqrt(2)) + 3)
r2 = r1*(-3*sqrt(2) - 2*sqrt(10 - 7*sqrt(2)) + 5)
elseif type == 3
r2 = arg
a = 2*r2*(1 + sqrt(2) + sqrt(2*sqrt(2) + 4))
r1 = r2*(sqrt(2) + 3 + 2*sqrt(sqrt(2) + 2))
else
println("draw の引数は, (arg, type, more) で，type は 1, 2, 3 のいずれか")
return
end
@printf("type = %d; a = %g; r1 = %g; r2 = %g\n", type, a, r1, r2)
plot([0, a, a, 0, 0], [0, 0, a, a, 0], color=:magenta, lw=0.5, showaxis=false)
circle(0, 0, r1, :blue, beginangle=0, endangle=90)
circle(a, a, r1, :blue, beginangle=180, endangle=270)
circle(0, 0, r1/2, :green, beginangle=0, endangle=90)
circle(a, a, r1/2, :green, beginangle=180, endangle=270)
for i = 1:3
θ = 90i/4
segment(0, 0, r1/2*cosd(θ), r1/2*sind(θ), :green)
segment(a, a, a + r1/2*cosd(θ + 180), a + r1/2*sind(θ + 180), :green)
end
circlef(a - r2, r2, r2)
circlef(r2, a - r2, r2)
circlef(a/2, a/2, r2)
if more
delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3 # size[2] * fontsize * 2
# hline!([0], color=:gray80, lw=0.5)
# vline!([0], color=:gray80, lw=0.5)
point(a, a, "(a,a)", :magenta, :right, :bottom, delta=delta/2)
point(0, 0, "(0,0)", :magenta, :left, delta=-delta/2)
point(r1, 0, "r1", :magenta, :center, delta=-delta/2)
end
end;
draw(1, 1, true)
算額（その1459）
藤田貞資(1789)：神壁算法巻上
http://www.wasan.jp/jinpeki/jinpekisanpo1.pdf
所懸于東都不忍辯才天堂者一事
天明8年戊申 四月
末術 關流神谷幸吉定令門人 東都本郷二丁目 佐久間直次郎善久
キーワード：三角形，矢，相似
#Julia, #SymPy, #算額, #和算
三角形内を 3 本の線分で区画して，元の三角形と相似な小さい三角形を作る。対応する頂点を結ぶ線分を，甲矢，乙矢，丙矢と呼ぶ。
元の三角形の 3 辺（大斜，中斜，小斜）が 6 寸，4 寸，3 寸で，甲矢，乙矢，丙矢の和が 5 寸のとき，甲矢の長さは何寸か。
図を構成する座標点を，(0, 0), (x1, 0), (x2, y2), (x3, y3), (x4, y4), (x5, y5) とおき，以下の連立方程式を解く。
解析解を求めるのは Julia の SymPy.solve では能力的に無理なので， nlsolve により数値解を求める。
甲矢は 2 寸である。ちなみに，乙矢は 2.25 寸，丙矢は 0.75 寸である。
全てのパラメータは以下のとおりである。
x1 = 6; x2 = 2.41667; y2 = 1.77756; x3 = 2.06342; y3 = 0.897116; x4 = 4.05; y4 = 0.44439; x5 = 2.99769; y5 = 1.30331
include("julia-source.txt");
# julia-source.txt ソース https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf
using SymPy
@syms x1, x2, y2, x3, y3, x4, y4, x5, y5,
大, 中, 小, 大2, 中2, 小2, 甲矢, 乙矢, 丙矢, 計
(大, 中, 小, 計) = (6, 4, 3, 5)
x1 = 大
大2 = sqrt((x4 - x3)^2 + (y3 - y4)^2)
中2 = sqrt((x4 - x5)^2 + (y5 - y4)^2)
小2 = sqrt((x5 - x3)^2 + (y5 - y3)^2)
甲矢 = sqrt((x1 - x4)^2 + y4^2)
乙矢 = sqrt(x3^2 + y3^2)
丙矢 = sqrt((x5 - x2)^2 + (y2 - y5)^2)
eq1 = 中^2 - ((x1 - x2)^2 + y2^2)
eq2 = 小^2 - (x2^2 + y2^2)
eq3 = 小2/小 - 中2/中
eq4 = 小2/小 - 大2/大
eq5 = 甲矢 + 乙矢 + 丙矢 - 計
eq6 = (甲矢 + 大2)^2 - ((x1 - x3)^2 + y3^2)
eq7 = (丙矢 + 中2)^2 - ((x4 - x2)^2 + (y2 - y4)^2)
eq8 = (乙矢 + 小2)^2 - (x5^2 + y5^2);
# 連立方程式の解は以下で求められるべきであるが，SymPy の能力的に不可能である
# res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8], (x2, y2, x3, y3, x4, y4, x5, y5))
# x2, y2 については，以下で求めることができる
res = solve([eq1, eq2], (x2, y2))[2]
(29/12, sqrt(455)/12)
# nlsolve により数値解を求める
using NLsolve
function nls(func, params...; ini = [0.0])
if typeof(ini) <: Number
r = nlsolve((vout, vin) -> vout[1] = func(vin[1], params..., [ini]), ftol=big"1e-40")
v = r.zero[1]
else
r = nlsolve((vout, vin)->vout .= func(vin, params...), ini, ftol=big"1e-40")
v = r.zero
end
return Float64.(v), r.f_converged
end;
function H(u)
(x2, y2, x3, y3, x4, y4, x5, y5) = u
return [
-y2^2 - (6 - x2)^2 + 16, # eq1
-x2^2 - y2^2 + 9, # eq2
sqrt((-x3 + x5)^2 + (-y3 + y5)^2)/3 - sqrt((x4 - x5)^2 + (-y4 + y5)^2)/4, # eq3
-sqrt((-x3 + x4)^2 + (y3 - y4)^2)/6 + sqrt((-x3 + x5)^2 + (-y3 + y5)^2)/3, # eq4
sqrt(x3^2 + y3^2) + sqrt(y4^2 + (6 - x4)^2) + sqrt((-x2 + x5)^2 + (y2 - y5)^2) - 5, # eq5
-y3^2 - (6 - x3)^2 + (sqrt(y4^2 + (6 - x4)^2) + sqrt((-x3 + x4)^2 + (y3 - y4)^2))^2, # eq6
-(-x2 + x4)^2 - (y2 - y4)^2 + (sqrt((-x2 + x5)^2 + (y2 - y5)^2) + sqrt((x4 - x5)^2 + (-y4 + y5)^2))^2, # eq7
-x5^2 - y5^2 + (sqrt(x3^2 + y3^2) + sqrt((-x3 + x5)^2 + (-y3 + y5)^2))^2, # eq8
]
end;
x1 = 6
iniv = BigFloat[30, 23, 21, 10, 55, 3, 37, 17] .* 6/74.6
res = nls(H, ini=iniv)
([2.4166666666666665, 1.7775607506417952, 2.0634154925612784, 0.8971156586851534, 4.04999554843813, 0.4443901876604488, 2.9976902975298994, 1.3033123554115043], true)
function draw(j, k, more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
x1 = 6
(x2, y2, x3, y3, x4, y4, x5, y5) = [30, 23, 21, 10, 55, 3, 37, 17] .* 6/74.6
(x2, y2, x3, y3, x4, y4, x5, y5) = res[1]
甲矢 = sqrt((x1 - x4)^2 + y4^2)
乙矢 = sqrt(x3^2 + y3^2)
丙矢 = sqrt((x5 - x2)^2 + (y2 - y5)^2)
@printf("甲矢 = %g; 乙矢 = %g; 丙矢 = %g\n", 甲矢, 乙矢, 丙矢)
@printf("x1 = %g; x2 = %g; y2 = %g; x3 = %g; y3 = %g; x4 = %g; y4 = %g; x5 = %g; y5 = %g\n", x1, x2, y2, x3, y3, x4, y4, x5, y5)
plot([0, x1, x2, 0], [0, 0, y2, 0], color=:blue, lw=0.5)
segment(0, 0, x5, y5)
segment(x3, y3, x1, 0)
segment(x2, y2, x4, y4)
segment(x1, 0, x4, y4, :green, lw=1.5)
segment(0, 0, x3, y3, :green, lw=1.5)
segment(x2, y2, x5, y5, :green, lw=1.5)
if more
delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3 # size[2] * fontsize * 2
hline!([0], color=:gray80, lw=0.5)
vline!([0], color=:gray80, lw=0.5)
point(0, 0, "0", :red, :center, delta=-3delta)
point(x1, 0, "x1", :red, :center, delta=-3delta)
point(x2, y2, "(x2,y2)", :red, :center, :bottom, delta=3delta)
point(x3, y3, "(x3,y3)", :red, :left, delta=-4delta, deltax=-8delta)
point(x4, y4, "(x4,y4)", :red, :left, delta=-3delta, deltax=-15delta)
point(x5, y5, "(x5,y5)", :red, :right, :vcenter, deltax=-5delta)
point(x2/2, y2/2, "小斜 ", :blue, :right, :vcenter, mark=false)
point((x1 + x2)/2, y2/2, " 中斜", :blue, :left, :vcenter, mark=false)
point(x1/2, 0, "大斜", :blue, :center, delta=-3delta, mark=false)
point((x4+x1)/2, y4/2, "甲矢", :green, :right, mark=false)
point(x3/2, y3/2, "乙矢", :green, :left, mark=false)
point((x5 + x2)/2, (y2 + y5)/2, "丙矢 ", :green, :right, :vcenter, mark=false)
ylims!(-13delta, y2 + 3delta)
end
end;
draw(125/2, 2000/2, true)
算額（その1458）
七〇 加須市大字外野 棘脱地蔵堂
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：円4個，外円，半円1個，楕円2個
#Julia, #SymPy, #算額, #和算
外円の中に外円と同じ直径の円弧を描き，甲円 1 個，乙円 2 個，等楕円 2 個を容れる。甲円の直径は外円の直径の半分であり，楕円と外円および楕円と弧は 1 点で接する。乙円の直径が 1 寸のとき，楕円の長径の最大値はいかほどか。
注：円弧は外円の中心を通る。また，外円（および円弧）は楕円の曲率円である。
外円の半径と中心座標を R, (0, 0)
弧の半径と中心座標を R, (0, -R)
甲円の半径と中心座標を r1, (0, R - r1); R = 2r1
乙円の半径と中心座標を r2, (x2, y2)
楕円の長径と短径を 2a, 2b; b = r1/2; 曲率円であるから a = sqrt(b*R)
とおき，以下の連立方程式を解く。
R = 2r1, b = r1/2 のとき， R = a^2/b ゆえ，a = sqrt(2r1*r1/2) = r1 である。
よって，「楕円の長径 = 甲円の直径」である。
include("julia-source.txt");
# # julia-source.txt ソース https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf
using SymPy
@syms R::positive, r1::positive, r2::positive,
x2::positive, y2::positive,
a::positive, b::positive
R = 2r1
b = r1/2
eq1 = x2^2 + (R - r1 - y2)^2 - (r1 + r2)^2
eq2 = x2^2 + (R + y2)^2 - (R + r2)^2
eq3 = x2^2 + y2^2 - (R - r2)^2
res = solve([eq1, eq2, eq3], (r1, x2, y2))[1]
(5*r2/3, 4*sqrt(3)*r2/3, r2/3)
r1 = 5r2/3である。
よって，楕円の長径は 2r1 = 10r2/3 で，r2 = 1/2 寸のとき 2r1 = 5/3 ≒ 1.67 寸。
以下は，図を描くためにパラメータを求める。
R = 2r1 = 10r2/3
b = 5r2/6 である。
蛇足
外円が楕円の曲率円であるから R = a^2/b より，
a = sqrt(b*R) = sqrt(5r2/6 * 10r2/3) = sqrt(25r2^2/9) = r2*5/3
楕円の長径 = 2a = 2r2*5/3 = 乙円の直径 * 5/3 ≒ 乙円の直径 * 1.67
乙円の直径が 1 寸のとき，楕円の長径は 5/3 寸 ≒ 1.67 寸である。甲円の直径と等しい。
function draw(r2, more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(r1, x2, y2) = (5*r2/3, 4*sqrt(3)*r2/3, r2/3)
R = 2r1
b = r1/2 # R/4
a = sqrt(b*R) # 曲率円: R = a^2/b
@printf("乙円の直径が %g のとき，楕円の長径は %g である。\n", 2r2, 2a)
@printf("r2 = %g; r1 = %g; x2 = %g; y2 = %g; R = %g; a = %g; b = %g\n",
r2, r1, x2, y2, R, a, b)
plot()
circle(0, 0, R, :magenta)
circle(0, -R, R, beginangle=30, endangle=150)
circle(0, R - r1, r1, :green)
circle2(x2, y2, r2)
ellipse(0, -R/4, a, R/4, color=:blue)
ellipse(0, -3R/4, a, R/4, color=:blue)
if more == true
delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3 # size[2] * fontsize * 2
hline!([0], color=:gray80, lw=0.5)
vline!([0], color=:gray80, lw=0.5)
point(0, R, "R", :green, :center, :bottom, delta=delta)
point(0, R - r1, "甲円:r1,(0,R-r1)", :green, :center, delta=-delta)
point(x2, y2, "乙円:r2,(x2,y2)", :red, :center, delta=-delta)
point(0, -R/4, "楕円:a,b,(0,-R/4)", :blue, :center, delta=-delta)
point(0, -3R/4, "楕円:a,b,(0,-3R/4)", :blue, :center, delta=-delta)
#plot!(xlims=(-0.2, 0.2), ylims=(-0.1, 0.05))
end
end;
draw(1/2, true)
PVアクセスランキング にほんブログ村
