裏 RjpWiki

Julia ときどき R, Python によるコンピュータプログラム,コンピュータ・サイエンス,統計学

算額(その2100)

2024年09月19日 | Julia

算額(その2100)

百三 群馬県高崎市八幡町 八幡宮 安政7年(1860)
群馬県和算研究会:群馬の算額,上武印刷株式会社,高崎市,1987年3月31日.
キーワード:楕円3個,長方形

長方形の中に相似な楕円 3 個を容れる。長方形の短辺が与えられたとき,長辺はいかほどか。

長方形の長辺を c
楕円の長半径,短半径,中心座標を a, b, (b, a), (c - a, b)
とおき,以下の連立方程式を解く。
長方形の短辺は 2a = 4b である。

include("julia-source.txt");

using SymPy
@syms a::positive, b::positive, c::positive,
     x0::positive, y0::positive;
eq1 = (x0 - b)^2/b^2 + (y0 - a)^2/a^2 - 1
eq2 = (x0 -c + a)^2/a^2 + (y0 - b)^2/b^2 - 1
eq3 = a^2*(x0 - b)/(b^2*(y0 - a)) - b^2*(x0 - c + a)/(a^2*(y0 - b));

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)
   (c, x0, y0) = u
   return [
-1 + (-b + x0)^2/b^2 + (-a + y0)^2/a^2,  # eq1
-1 + (-b + y0)^2/b^2 + (a - c + x0)^2/a^2,  # eq2
a^2*(-b + x0)/(b^2*(-a + y0)) - b^2*(a - c + x0)/(a^2*(-b + y0)),  # eq3
   ]
end;
a = 1/2
b = a/2
iniv = BigFloat[2.6475964903602565, 0.8542484563784245, 0.5046106449806166] .*a
res = nls(H, ini=iniv)

   ([1.4708869390890313, 0.4745824757657914, 0.28033924721145365], true)

長方形の長辺は,楕円の長半径の 2.9417738781780627 倍である。
楕円の長半径が 1/2(すなわち,長方形の短辺が 1)のとき,長辺は 1.47089 である。

楕円の長半径,短半径は 0.5, 0.25,接点の座標は (0.474582, 0.280339) である。

function draw(短辺, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   a = 短辺/2
   b = a/2
   (c, x0, y0) = a .*[2.9417738781780627, 0.9491649515315828, 0.5606784944229073]
   @printf("長方形の短辺が %g のとき,長辺は %g である。\n楕円の長半径,短半径は %g, %g,接点の座標は (%g, %g) である。", 短辺, c, a, b, x0, y0)
   plot([0, c, c, 0, 0], [0, 0, 2a, 2a, 0], color=:green, lw=0.5)
   ellipse(b, a, b, a, color=:red)
   ellipse(c - a, b, a, b, color=:red)
   ellipse(c - a, 3b, a, b, color=:red)
   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(b, a, "(b,a)")
       point(c - a, b, "(c-a,b)")
       point(x0, y0, "(x0,y0)")
       point(c, 4b, "(c,2a)", :green, :right, :bottom, delta=delta/2)
   end

end;

draw(1, true)

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

算額(その2099)

2024年09月19日 | Julia

#和算 - ブログ村ハッシュタグ
#和算

算額(その2099)

八十七 群馬県碓氷郡松井田町峠 熊野神社 安政4年(1857)
群馬県和算研究会:群馬の算額,上武印刷株式会社,高崎市,1987年3月31日.
キーワード:直角三角形,鈎股弦

直角三角形(鈎股弦)において,「鈎*股」が 52440 歩,「股*弦」が 130065 歩のとき,鈎,股,弦を求めよ。

以下の連立方程式を解く。

using SymPy
@syms 鈎::positive, 股::positive, 弦::positive,
     鈎股相乗::positive, 股弦相乗::positive;
eq1 = 鈎*股 - 鈎股相乗
eq2 = 股*弦 - 股弦相乗
eq3 = 鈎^2 + 股^2 - 弦^2
eq3 = sqrt(鈎^2 + 股^2) - 弦
res = solve([eq1, eq2, eq3], (鈎, 股, 弦))[2]  # 2 of 2

   (鈎股相乗*(1/(股弦相乗^2 - 鈎股相乗^2))^(1/4), (1/(股弦相乗^2 - 鈎股相乗^2))^(-1/4), 股弦相乗*(1/(股弦相乗^2 - 鈎股相乗^2))^(1/4))

鈎 = res[1](鈎股相乗 => 52440, 股弦相乗 => 130065)
股 = res[2](鈎股相乗 => 52440, 股弦相乗 => 130065)
弦 = res[3](鈎股相乗 => 52440, 股弦相乗 => 130065)
(鈎, 股, 弦)

   (152, 345, 377)

「鈎*股」が 52440 歩,「股*弦」が 130065 歩のとき,鈎 = 152 寸,股 = 345 寸,弦 = 377 寸である。

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

算額(その2098)

2024年09月19日 | Julia

#和算 - ブログ村ハッシュタグ
#和算

算額(その2098)

七十六 群馬県桐生市天神町 天満宮 嘉永5年(1852)
群馬県和算研究会:群馬の算額,上武印刷株式会社,高崎市,1987年3月31日.
キーワード:年利

銀 76 匁を,複利で 7 年預け,元利合計が 357.62786865234375 匁であった。年利はいかほどか。

「術」は 400 文字以上にわたり述べられている。

「元本」,「年利」,「年数」,「元利合計」の関係は以下の式になる。

元利合計 = 元本*(1 + 年利)^年数

using SymPy
@syms 元本::positive, 年利::positive, 年数::positive,
     元利合計::positive;
eq = 元利合計 ⩵ 元本*(1 + 年利)^年数
eq |> println

   Eq(元利合計, 元本*(年利 + 1)^年数)

方程式を解いて,年利を求める。
年利 = (元利合計/元本)^(1/年数) - 1

res = solve(eq, 年利)[1]
res |> println

   元利合計^(1/年数)/元本^(1/年数) - 1

「元本」,「年数」,「元利合計」を式に代入し,「年利」を得る。
元本 76 匁を,複利で 7 年預け,元利合計が 357.62786865234375 匁のとき,年利は 0.25 = 2 割 5 分 である。

res(元本 => 75, 年数 => 7, 元利合計 => 357.62786865234375).evalf() |>  println

   0.250000000000000

 

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

PVアクセスランキング にほんブログ村

PVアクセスランキング にほんブログ村