裏 RjpWiki

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

算額(その1538)

2025年01月13日 | Julia

算額(その1538)

(6) 大阪府茨木市大字総持寺 総持寺 嘉永7年(1854)
近畿数学史学会:近畿の算額「数学の絵馬を訪ねて」,平成4年5月16日 初版第一刷,大阪教育図書株式会社,大阪市.

キーワード:円6個,外円
#Julia, #SymPy, #算額, #和算

問は「如図大円之内五箇円径問」とだけあり,答は「中円一尺二寸,脇円八寸,小円六寸」とある。

注:大円の直径は中円の直径の 2 倍ではあるが,問で大円の直径を述べて置かなければ問が成り立たない。また,算額の図がどうなっているかのわからないが横に 3 個並んでいる円は等円のようだし,脇円,小円なのだろうが,答えにある直径を加えると2尺2寸になり,大円の直径より小さくなる。とりあえず脇円,小円を区別して立式するが,脇円が中円,小円,大円に接しているということは譲れない。また,「只云大円径二尺四寸」があるものとする。

大円の半径と中心座標を R, (0, 0)
中円の半径と中心座標を r1, (0, r1)
脇円の半径と中心座標を r2, (R - r2, 0)
小円の半径と中心座標を r3, (0, 0)
とおき,以下の方程式を解く。

include("julia-source.txt")
# julia-source.txt ソース https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf

using SymPy
@syms R, r1, r2, r3
eq1 = r1 - R/2
eq2 = r3 - (2r1 - 2r2)
eq3 = r1^2 + (r3 + r3)^2 - (r1 + r2)^2
res = solve([eq1, eq2, eq3], (r1, r2, r3))[1]  # 1 of 2

    (R/2, R/3, R/3)

中円の半径は大円の半径の 1/2 倍,脇円は同じく 1/3 倍,小円も 1/3 倍ということで,脇円と小円は同じ大きさであるということになった。
大円の直径が 2 尺 4 寸のとき,中円の直径は 1 尺 2 寸,脇円,小円の直径は 8 寸である。

function draw(R, more=false)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    (r1, r2, r3) = R.*(1/2, 1/3, 1/3)
    plot()
    circle(0, 0, R)
    circle22(0, r1, r1, :green)
    circle2(R - r2, 0, r2, :blue)
    circle(0, 0, r3, :magenta)
    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, R, "R", :red, :center, :bottom, delta=delta/2)
        point(0, r1, "中円:r1,(0,r1)", :green, :center, delta=-delta/2)
        point(R - r2, 0, "脇円:r2,(R-r2,0)", :blue, :center, delta=-delta/2)
        point(0, 0, "小円:r3,(0,0)", :magenta, :center, delta=-delta/2)
    end
end;

draw(24/2, true)

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

算額(その1537)

2025年01月13日 | Julia

算額(その1537)

(17) 京都府京都市東山区清水 清水寺 明治25年(1892)
近畿数学史学会:近畿の算額「数学の絵馬を訪ねて」,平成4年5月16日 初版第一刷,大阪教育図書株式会社,大阪市.

キーワード:半円1個,扇,斜線2本
#Julia, #SymPy, #算額, #和算

1/3 円でできる扇に,半円 1 個と斜線 2 本を容れる。扇長(要から扇の先端までの長さ)が 20 寸,骨径(見えている骨の長さ)が 8寸のとき,斜線の長さ,半円の直径はいかほどか。

扇長を R, 骨長を r2
半円の半径と中心座標を r1, (0, R - r1)
とおき,以下の方程式を解いて r1 を求める。

include("julia-source.txt")
# julia-source.txt ソース https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf

using SymPy
@syms R, r1, r2
eq1 = r1/(√Sym(3)*R/2) - (R - r1 - r2)/(R/2 - r2);

ans_r1 = solve(eq1, r1)[1]
ans_r1 |> println

    3*R*(R - r2)/(sqrt(3)*R + 3*R - 2*sqrt(3)*r2)

半円の半径は 3R*(R - r2)/(√3R + 3R - 2√3r2) である。

2*ans_r1(R => 20, r2 => 8).evalf() |> println

    21.5155932850234

扇長が 20 寸,骨長が 8 寸のとき,半円の直径は 21.5155932850234 寸である。

斜線の端点は (0, r2), (√3R/2, R/2) なので,直接計算できる。

sqrt((√Sym(3)R/2)^2 + (R/2 - r2)^2)(R => 20, r2 => 8).evalf() |> println

    17.4355957741627

斜線の長さは 17.4355957741627 寸である。

function draw(R, r2, more=false)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    r1 = 3R*(R - r2)/(√3R + 3R - 2√3r2)
    plot([-√3R/2, 0, √3R/2], [R/2, r2, R/2], color=:green, lw=0.5)
    circle(0, 0, R, beginangle=30, endangle=150)
    plot!([-√3R/2, 0, √3R/2], [R/2, 0, R/2], color=:red, lw=0.5)
    circle(0, 0, r2, :magenta, beginangle=30, endangle=150)
    circle(0, R - r1, r1, :blue, beginangle=0, endangle=180)
    segment(-r1, R - r1, r1, R - r1, :blue)
    for θ = 30:15:150
        segment(0, 0, r2*cosd(θ), r2*sind(θ), :black, lw=2)
    end
    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, R, "R", :red, :center, :bottom, delta=delta)
        point(0, R - r1, "R-r1", :blue, :center, :bottom, delta=delta)
        point(0, r2, "r2", :magenta, :center, :bottom, delta=delta/2)
        circle(0, 0, r2/3, beginangle=0, endangle=30, lw=2)
        point(√3r2/6, 0, "30°", :red, :left, :bottom, delta=delta, deltax=4delta, mark=false)
    end
end;

draw(20, 8, true)

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

算額(その1536)

2025年01月13日 | Julia

算額(その1536)

(12) 京都市伏見区御香宮門前町 御香宮神社(ごこうのみやじんじゃ)文久3年(1863)
近畿数学史学会:近畿の算額「数学の絵馬を訪ねて」,平成4年5月16日 初版第一刷,大阪教育図書株式会社,大阪市.

キーワード:円4個,水平線,斜線2本
#Julia, #SymPy, #算額, #和算

3 本の直線で区画された領域に,大円,中円,小円,天円を容れる。大円,中円,天円の直径がそれぞれ 85 寸,25 寸,8 寸のとき,小円の直径はいかほどか。

中円の直径を d1
大円の直径を d3
天円の直径を d4
小円の直径を d2
として,算法助術の公式62を用いる。

using SymPy
@syms d1, d2, d3, d4
eq62 = -d1*d2*d3 + (d1 + d2)*d3*d4 + d1*d2*d4;
solve(eq62, d2)[1] |> println

    d1*d3*d4/(d1*d3 - d1*d4 - d3*d4)

小円の直径は d2 = d1*d3*d4/(d1*d3 - d1*d4 - d3*d4) である。

ans_d2 = solve(eq62(d1 => 25, d3 => 85, d4 => 8), d2)[1].evalf()
ans_d2 |> println

    13.6546184738956

d1, d3, d4 に具体値を代入して d2 = 13.6546184738956 である。

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

算額(その1535)

2025年01月13日 | Julia

算額(その1535)

滋賀県新旭町太田 大田神社 (1868)
(四)滋賀県 (8)大田神社
近畿数学史学会:近畿の算額「数学の絵馬を訪ねて」,平成4年5月16日 初版第一刷,大阪教育図書株式会社,大阪市.

キーワード:等差数列,部分和
#Julia, #SymPy, #算額, #和算

一辺の長さが a[i], i=1,2,...,n, a[k] > a[k+1] の等差数列に従う n 個の正方形がある。
それらの正方形の面積の和,一辺の長さの和および公差 K = a[k] - a[k+i] が与えられたとき,正方形の総数を求めよ。

算額(その1534)の類題である。こちらのほうが易しい。

まず最初に,面積の和と一辺の長さの和を一辺の長さの初項 a,公差 K, と項数 n で表す式を決定する。

1. n 個の正方形の面積の和

初項は a^2, 第 2 項までは a^2 + (a - K)^2, 第 3 項までは a^2 + (a - K)^2 + (a - 2K)^2 ,... なので,以下のようにプログラムすれば n ごとの式が得られる。

using SymPy
@syms s, a, K
s = 0
for n = 0:10
    s += (a - n*K)^2
    s |> expand |> println
end

    a^2
    K^2 - 2*K*a + 2*a^2
    5*K^2 - 6*K*a + 3*a^2
    14*K^2 - 12*K*a + 4*a^2
    30*K^2 - 20*K*a + 5*a^2
    55*K^2 - 30*K*a + 6*a^2
    91*K^2 - 42*K*a + 7*a^2
    140*K^2 - 56*K*a + 8*a^2
    204*K^2 - 72*K*a + 9*a^2
    285*K^2 - 90*K*a + 10*a^2
    385*K^2 - 110*K*a + 11*a^2

ついで,これらの係数の一般項を求める。

ai^2 の項の係数は (n + 1)
ai^1 の項の係数は -n*(n+1) * K
ai^0 の項の係数は n*(n+1)*(2*n+1)/6 * K^2
である。

まとめると,体積の和は以下のようになる。

@syms n::positive, K
S = n*(n+1)*(2*n+1)/6*K^2 - n*(n+1)*K*a + (n + 1)*a^2
S |> println

    K^2*n*(n + 1)*(2*n + 1)/6 - K*a*n*(n + 1) + a^2*(n + 1)

n = 10 のときの式は以下のようになり,前述の式のとおりになっていることが確認できる。

S(n => 10) |> println

    385*K^2 - 110*K*a + 11*a^2

2. n 個の正方形の一辺の長さの和

一辺の長さの和の一般項を求めるのも前節と同様にして以下のようになる。

using SymPy
@syms s, a, K
s = 0
for i = 0:10
    s += (a - i*K)
    s |> expand |> println
end

    a
    -K + 2*a
    -3*K + 3*a
    -6*K + 4*a
    -10*K + 5*a
    -15*K + 6*a
    -21*K + 7*a
    -28*K + 8*a
    -36*K + 9*a
    -45*K + 10*a
    -55*K + 11*a

ai^1 の項の係数は (n+1)
ai^0 の項の係数は n*(n+1)/2 * K

@syms n::positive, K
L = -n*(n+1)/2*K + (n + 1)*a
L |> println

    -K*n*(n + 1)/2 + a*(n + 1)

L(n => 10) |> println

    -55*K + 11*a

3. 算額問題を解く

1, 2 節で得られた一般項の式と具体的な面積の和と一辺の長さの和から正方形の個数を求める。

テストデータとして,初項 a = 20,公差 K = 3(実際は -3)のときの n = 4 個の正方形の一辺の長さが [20, 17, 14, 11] のとき,面積の和は 1006,一辺の長さの和は 62 である。

a = [20, 17, 14, 11]
sum(a.^2), sum(a)

    (1006, 62)

@syms n::positive, a::positive, K::positive, 末尾::positive
S = n*(n+1)*(2*n+1)/6*K^2 - n*(n+1)*K*a + (n + 1)*a^2
L = -n*(n+1)/2*K + (n + 1)*a
eq1 = S(K => 3) - 1006
eq2 = L(K => 3) - 62
solve([eq1, eq2], (n, a))[1]

    (3, 20)

n = 3 であるが,n は 0 から数えるので,正方形の個数は n + 1 = 4 個である。

もう一つのテストデータとして,初項 a = 100,公差 K = 5(実際は -5)のときの n = 16 個の正方形の一辺の長さは,
[100, 95, 90, 85, 80, 75, 70, 65, 60, 55, 50, 45, 40, 35, 30, 25]
である。
面積の和は 71000,一辺の長さの和は 1000 である。

a = collect(range(100, step=-5, length=16));
a |> println
(sum(a.^2), sum(a)) |> println

    [100, 95, 90, 85, 80, 75, 70, 65, 60, 55, 50, 45, 40, 35, 30, 25]
    (71000, 1000)

@syms n::positive, a::positive, K::positive, 末尾::positive
S = n*(n+1)*(2*n+1)/6*K^2 - n*(n+1)*K*a + (n + 1)*a^2
L = -n*(n+1)/2*K + (n + 1)*a
eq1 = S(K => 5) - 71000
eq2 = L(K => 5) - 1000
solve([eq1, eq2], (n, a))[1]

    (15, 100)

正方形の個数は n + 1 = 16 個である。

4. 一般解を求める

面積の和 K1,一辺の長さの和 K2,公差 K3 を記号定数として一般解を求める。

一度には解けないので,段階を追って解く。

using SymPy
@syms n::positive, a::positive, K::positive, K1::positive, K2::positive, K3::positive
S = n*(n+1)*(2*n+1)/6*K^2 - n*(n+1)*K*a + (n + 1)*a^2
L = -n*(n+1)/2*K + (n + 1)*a
eq1 = S(K => K3) - K1
eq2 = L(K => K3) - K2;
ans_a = solve(eq2, a)[1];

eq11 = eq1(a =>  ans_a) |> simplify
ans_n = solve(eq11, n)[3];  # 3 of 4
ans_n |> println

    (-6*K3^(5/2) + sqrt(3)*sqrt(K3^(13/3)*(144*K2^2 + K3^2)/(1944*K1^2 + 432*K2^2*K3^2 - K3^4 + sqrt(-K3^2*(144*K2^2 + K3^2)^3 + (1944*K1^2 + 432*K2^2*K3^2 - K3^4)^2))^(1/3) + K3^(11/3)*(1944*K1^2 + 432*K2^2*K3^2 - K3^4 + sqrt(-K3^2*(144*K2^2 + K3^2)^3 + (1944*K1^2 + 432*K2^2*K3^2 - K3^4)^2))^(1/3) + 2*K3^5) - sqrt(3)*sqrt(72*sqrt(3)*K1*K3^(22/3)/sqrt((2*K3^(26/3)*(1944*K1^2 + 432*K2^2*K3^2 - K3^4 + sqrt(-K3^2*(144*K2^2 + K3^2)^3 + (1944*K1^2 + 432*K2^2*K3^2 - K3^4)^2))^(1/3) + K3^(22/3)*(1944*K1^2 + 432*K2^2*K3^2 - K3^4 + sqrt(-K3^2*(144*K2^2 + K3^2)^3 + (1944*K1^2 + 432*K2^2*K3^2 - K3^4)^2))^(2/3) + K3^8*(144*K2^2 + K3^2))/(1944*K1^2 + 432*K2^2*K3^2 - K3^4 + sqrt(-K3^2*(144*K2^2 + K3^2)^3 + (1944*K1^2 + 432*K2^2*K3^2 - K3^4)^2))^(1/3)) - K3^(13/3)*(144*K2^2 + K3^2)/(1944*K1^2 + 432*K2^2*K3^2 - K3^4 + sqrt(-K3^2*(144*K2^2 + K3^2)^3 + (1944*K1^2 + 432*K2^2*K3^2 - K3^4)^2))^(1/3) - K3^(11/3)*(1944*K1^2 + 432*K2^2*K3^2 - K3^4 + sqrt(-K3^2*(144*K2^2 + K3^2)^3 + (1944*K1^2 + 432*K2^2*K3^2 - K3^4)^2))^(1/3) + 4*K3^5))/(6*K3^(5/2))

大変長い式であるが,ちゃんと数式解は求まる。

ans_n(K1 => 71000, K2 => 1000, K3 => 5).evalf() |> println

    15.0000000000000

 

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

算額(その1534)

2025年01月13日 | Julia

算額(その1534)

奈良県大和郡山市 庚申堂 明治13年(1880)

http://www.wasan.jp/nara/kosindo.html

(五)奈良県 (5)庚申堂
近畿数学史学会:近畿の算額「数学の絵馬を訪ねて」,平成4年5月16日 初版第一刷,大阪教育図書株式会社,大阪市.

キーワード:等差数列,部分和
#Julia, #SymPy, #算額, #和算

一辺の長さが a[i], i=1,2,...,n, a[k] > a[k+1] の等差数列に従う n 個の立方体がある。
それらの立体の体積の和,表面積の和および公差 K = a[k] - a[k+i] が与えられたとき,末尾の(最小の)立方体の一辺の長さを求めよ。

まず最初に,体積の和と表面積の和を一辺の長さの初項 a,公差 K, と項数 n で表す式を決定する。

1. n 個の立方体の体積の和

初項は a^3, 第 2 項までは a^3 + (a - K)^3, 第 3 項までは a^3 + (a - K)^3 + (a - 2K)^3 ,... なので,以下のようにプログラムすれば n ごとの式が得られる。

using SymPy
@syms s, a, K
s = 0
for n = 0:10
    s += (a - n*K)^3
    s |> expand |> println
end

    a^3
    -K^3 + 3*K^2*a - 3*K*a^2 + 2*a^3
    -9*K^3 + 15*K^2*a - 9*K*a^2 + 3*a^3
    -36*K^3 + 42*K^2*a - 18*K*a^2 + 4*a^3
    -100*K^3 + 90*K^2*a - 30*K*a^2 + 5*a^3
    -225*K^3 + 165*K^2*a - 45*K*a^2 + 6*a^3
    -441*K^3 + 273*K^2*a - 63*K*a^2 + 7*a^3
    -784*K^3 + 420*K^2*a - 84*K*a^2 + 8*a^3
    -1296*K^3 + 612*K^2*a - 108*K*a^2 + 9*a^3
    -2025*K^3 + 855*K^2*a - 135*K*a^2 + 10*a^3
    -3025*K^3 + 1155*K^2*a - 165*K*a^2 + 11*a^3

ついで,これらの係数の一般項を求める。

ai^3 の項の係数は (n + 1) * K^0
ai^2 の項の係数は -3*n*(n + 1)/2 * K^1
ai^1 の項の係数は n*(n+1)*(2*n+1)/2 * K^2
ai^0 の項の係数は -(n*(n+1)/2)^2 * K^3
である。

まとめると,体積の和は以下のようになる。

@syms n::positive, K
V = -(n*(n+1)/2)^2*K^3 + n*(n+1)*(2*n+1)/2*K^2*a - 3*n*(n + 1)/2*K*a^2 + (n + 1)*a^3
V |> println

    -K^3*n^2*(n + 1)^2/4 + K^2*a*n*(n + 1)*(2*n + 1)/2 - 3*K*a^2*n*(n + 1)/2 + a^3*(n + 1)

n = 10 のときの式は以下のようになり,前述の式のとおりになっていることが確認できる。

V(n => 10) |> println

    -3025*K^3 + 1155*K^2*a - 165*K*a^2 + 11*a^3

2. n 個の立方体の表面積の和

体積の和の一般項を求めるのも前節と同様にして以下のようになる。

using SymPy
@syms s, a, K
s = 0
for i = 0:10
    s += 6(a - i*K)^2
    s |> expand |> println
end

    6*a^2
    6*K^2 - 12*K*a + 12*a^2
    30*K^2 - 36*K*a + 18*a^2
    84*K^2 - 72*K*a + 24*a^2
    180*K^2 - 120*K*a + 30*a^2
    330*K^2 - 180*K*a + 36*a^2
    546*K^2 - 252*K*a + 42*a^2
    840*K^2 - 336*K*a + 48*a^2
    1224*K^2 - 432*K*a + 54*a^2
    1710*K^2 - 540*K*a + 60*a^2
    2310*K^2 - 660*K*a + 66*a^2

ai^2 の項の係数は 6(n + 1)
ai^1 の項の係数は -6n*(n+1) * K
ai^0 の項の係数は n*(n+1)*(2n+1) * K^2

@syms n::positive, K
S = n*(n+1)*(2n+1)*K^2 - 6n*(n+1)*K*a + 6(n + 1)*a^2
S |> println

    K^2*n*(n + 1)*(2*n + 1) - 6*K*a*n*(n + 1) + a^2*(6*n + 6)

S(n => 10) |> println

    2310*K^2 - 660*K*a + 66*a^2

3. 算額問題を解く

1, 2 節で得られた一般項の式と具体的な体積と表面積の値から立方体の個数を求める。

テストデータとして,初項 a = 20,公差 K = 3(実際は -3)のときの n = 4 個の立方体の一辺の長さが [20, 17, 14, 11] のとき,体積の和は 16988,表面積の和は 6036 である。

a = [20, 17, 14, 11]
sum(a.^3), sum(6a.^2)

    (16988, 6036)

@syms n::positive, a::positive, K::positive, 末尾::positive
V = -(n*(n+1)/2)^2*K^3 + n*(n+1)*(2*n+1)/2*K^2*a - 3*n*(n + 1)/2*K*a^2 + (n + 1)*a^3
S = n*(n+1)*(2n+1)*K^2 - 6n*(n+1)*K*a + 6(n + 1)*a^2;
eq1 = V(K => 3) - 16988
eq2 = S(K => 3) - 6036
eq3 = (a - n*K)(K => 3) - 末尾
solve([eq1, eq2, eq3], (末尾, n, a))[1]

    (11, 3, 20)

一番小さい立方体の一辺の長さは 11 である。
n = 3 であるが,n は 0 から数えるので,立方体の個数は n + 1 であることに注意。
一番大きい立方体の一辺の長さ(初項 a) は 20 である。

もう一つテストデータとして,初項 a = 100,公差 K = 5(実際は -5)のときの n = 16 個の立方体の一辺の長さは,
[100, 95, 90, 85, 80, 75, 70, 65, 60, 55, 50, 45, 40, 35, 30, 25]
である。
体積の和は 5500000,表面積の和は 426000 である。

a = collect(range(100, step=-5, length=16));
a |> println
(sum(a.^3), sum(6a.^2)) |> println

    [100, 95, 90, 85, 80, 75, 70, 65, 60, 55, 50, 45, 40, 35, 30, 25]
    (5500000, 426000)

@syms n::positive, a::positive, K::positive, 末尾::positive
V = -(n*(n+1)/2)^2*K^3 + n*(n+1)*(2*n+1)/2*K^2*a - 3*n*(n + 1)/2*K*a^2 + (n + 1)*a^3
S = n*(n+1)*(2n+1)*K^2 - 6n*(n+1)*K*a + 6(n + 1)*a^2;
eq1 = V(K => 5) - 5500000
eq2 = S(K => 5) - 426000
eq3 = (a - n*K)(K => 5) - 末尾
solve([eq1, eq2, eq3], (末尾, n, a))[1]

    (25, 15, 100)

最小の立方体の一辺の長さは 25,最大の立方体の一辺の長さは 100, 立方体の個数は n + 1 = 16 である。

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

算額(その1533)

2025年01月13日 | Julia

算額(その1533)

奈良県奈良市 円満寺 天保 15 年
http://www.wasan.jp/nara/enman.html

奈良県奈良市下山町 円満寺 天保15年(1844)
近畿数学史学会:近畿の算額「数学の絵馬を訪ねて」,平成4年5月16日 初版第一刷,大阪教育図書株式会社,大阪市.

キーワード:等差数列
#Julia, #SymPy, #算額, #和算

簡単なものも,絨毯爆撃していく。

甲,乙,丙の正方形において,甲,乙の一辺の二乗和が 724,乙,丙の一辺の二乗和が 580 である。また,それぞれの辺の長さは 2 寸ずつの差がある。

甲,乙,丙の一辺の長さを a, b, c; a > b > c
甲,乙の一辺の二乗和を K1
乙,丙の一辺の二乗和を K2
辺の長さの差 K3
とおき,以下の連立方程式を解く。

using SymPy
@syms a::positive, b::positive, c::positive,
      K1::positive, K2::positive, K3::positive
eq1 = a^2 + (a - K3)^2 - K1
eq2 = b^2 + (b - K3)^2 - K2
res2 = solve([eq1, eq2], (a, b))[4];  # 3 of 4

# 甲
res2[1] |> println
res2[1](K1 => 724, K2 => 580, K3 => 2) |> println

    K3/2 + sqrt(2*K1 - K3^2)/2
    20

# 乙 = 甲 - K3
res2[2] |> println
res2[2](K1 => 724, K2 => 580, K3 => 2) |> println

    K3/2 + sqrt(2*K2 - K3^2)/2
    18

# 丙 = 甲 - 2K3
res2[2] - K3 |> println
(res2[2] - K3)(K1 => 724, K2 => 580, K3 => 2).evalf() |> println

    -K3/2 + sqrt(2*K2 - K3^2)/2
    16.0000000000000

甲,乙,丙の一辺の長さはそれぞれ,20 寸,18 寸,16 寸である。

辺の長さの差 K3 を使わないで連立方程式を立ててもよい。
当然,同じ数値解が得られるが,数式解の見た目は複雑になる。

using SymPy
@syms a::positive, b::positive, c::positive,
      K1::positive, K2::positive, K3::positive
eq1 = a^2 + b^2 - K1
eq2 = b^2 + c^2 - K2
eq3 = (a - b) - (b - c)
res = solve([eq1, eq2, eq3], (a, b, c))[1];  # 1 of 2

# 甲
res[1] |> println
res[1](K1 => 724, K2 => 580).evalf() |> println

    sqrt(2)*(-K1 + K2 - 2*sqrt(-K1^2 + 6*K1*K2 - K2^2))*sqrt(-K1 + 7*K2 - sqrt(-K1^2 + 6*K1*K2 - K2^2))/(4*(K1 - 5*K2))
    20.0000000000000

 

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

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

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