裏 RjpWiki

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

マイナス保険証

2023年08月04日 | 雑感

https://www.47news.jp/9676569.html
> 【速報】マイナ保険証登録後でも解除できる方向検討
> 共同通信
> 速報
> 政治
> 2023年08月03日共同通信
>  

> 政府はマイナ保険証に関し、取得した後でも登録を解除して資格確認書を選択することを可能にする方向で検討に入った。関係者が3日、明らかにした。

 

何だこれ?
マイナ保険証の意味なくなったじゃん。

 

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

算額(その363)

2023年08月04日 | Julia

算額(その363)

山形県寒河江市大字柴橋字台下 伊豆神社
山形算額勝負-湯殿山神社を目指せ-

https://www.sci.yamagata-u.ac.jp/wasan/pdf/20180711SSEP.pdf

直線上に大円,中円,小円が互いに接している。大円,中円の直径をそれぞれ 9 寸,4 寸としたとき,小円の直径を求めよ。

大円の半径と中心座標を r1, (0, r1)
中円の半径と中心座標を r2, (x2, r2)
小円の半径と中心座標を r3, (x3, r3)
として以下の連立方程式を解く。

include("julia-source.txt");

using SymPy

@syms r1::positive, r2::positive, x2::positive, r3::positive, x3::positive;

eq1 = x2^2 + (r1 - r2)^2 - (r1 + r2)^2
eq2 = x3^2 + (r1 - r3)^2 - (r1 + r3)^2
eq3 = (x2 - x3)^2 + (r2 - r3)^2 - (r2 + r3)^2
res = solve([eq1, eq2, eq3], (x2, r3, x3))

   2-element Vector{Tuple{Sym, Sym, Sym}}:
    (2*sqrt(r1)*sqrt(r2), (-2*r1^(7/2)*r2^(3/2) + 4*r1^(5/2)*r2^(5/2) - 2*r1^(3/2)*r2^(7/2) + r1^4*r2 - r1^3*r2^2 - r1^2*r2^3 + r1*r2^4)/(r1^4 - 4*r1^3*r2 + 6*r1^2*r2^2 - 4*r1*r2^3 + r2^4), (2*r1^(5/2)*sqrt(r2) - 2*r1^(3/2)*r2^(3/2) - 2*r1^2*r2 + 2*r1*r2^2)/(r1^2 - 2*r1*r2 + r2^2))
    (2*sqrt(r1)*sqrt(r2), (2*r1^(7/2)*r2^(3/2) - 4*r1^(5/2)*r2^(5/2) + 2*r1^(3/2)*r2^(7/2) + r1^4*r2 - r1^3*r2^2 - r1^2*r2^3 + r1*r2^4)/(r1^4 - 4*r1^3*r2 + 6*r1^2*r2^2 - 4*r1*r2^3 + r2^4), (2*r1^(5/2)*sqrt(r2) - 2*r1^(3/2)*r2^(3/2) + 2*r1^2*r2 - 2*r1*r2^2)/(r1^2 - 2*r1*r2 + r2^2))

最初の組のものが適解である。

r1, r2 が 4.5, 2 のとき,x2 = 6, r3 = 0.72, x3 = 3.6 である。

(r1, r2) = (9, 4) .// 2
(2*sqrt(r1)*sqrt(r2), (-2*r1^(7/2)*r2^(3/2) + 4*r1^(5/2)*r2^(5/2) - 2*r1^(3/2)*r2^(7/2) + r1^4*r2 - r1^3*r2^2 - r1^2*r2^3 + r1*r2^4)/(r1^4 - 4*r1^3*r2 + 6*r1^2*r2^2 - 4*r1*r2^3 + r2^4), (2*r1^(5/2)*sqrt(r2) - 2*r1^(3/2)*r2^(3/2) - 2*r1^2*r2 + 2*r1*r2^2)/(r1^2 - 2*r1*r2 + r2^2))

   (6.0, 0.7199999999999942, 3.6)

直線上にあり,互いに接している 2 円間の水平距離については
和算の心(その002)https://blog.goo.ne.jp/r-de-r/e/eeac1b6f4ed0f7c3a7dfdad299f066b3」
にも記したように,x2 = 2sqrt(r1*r2) = 6 で求めることができる。

2*sqrt(r1*r2)

   6.0

直線上にあり,互いに接している 3 円において,この場合の小円の半径については,デカルトの円定理により,1/(1/r1 + 1/r2 + 2sqrt(1/(r1*r2))) = 0.72 で求めることができる。

1/(1/r1 + 1/r2 + 2sqrt(1/(r1*r2)))

   0.72

なお,連立方程式の解で r3 は
(-2*r1^(7/2)*r2^(3/2) + 4*r1^(5/2)*r2^(5/2) - 2*r1^(3/2)*r2^(7/2) + r1^4*r2 - r1^3*r2^2 - r1^2*r2^3 + r1*r2^4)/(r1^4 - 4*r1^3*r2 + 6*r1^2*r2^2 - 4*r1*r2^3 + r2^4)
という長い式になっており,SymPy では簡約化されない。

大円と小円の中心の水平距離は x3 = 2sqrt(4.5 * 0.72) = 3.6 である。

2sqrt(4.5 * 0.72)

   3.6

r1, r2 を数値として与え連立方程式を解くと以下のようになる。

using SymPy

@syms r1::positive, r2::positive, x2::positive, r3::positive, x3::positive;
(r1, r2) = (9, 4) .// 2
eq1 = x2^2 + (r1 - r2)^2 - (r1 + r2)^2
eq2 = x3^2 + (r1 - r3)^2 - (r1 + r3)^2
eq3 = (x2 - x3)^2 + (r2 - r3)^2 - (r2 + r3)^2
res = solve([eq1, eq2, eq3], (x2, r3, x3))

   2-element Vector{Tuple{Sym, Sym, Sym}}:
    (6, 18/25, 18/5)
    (6, 18, 18)

(6, 18/25, 18/5)

   (6, 0.72, 3.6)

using Plots

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2) = (9, 4) .// 2
   (x2, r3, x3) = (6, 18/25, 18/5)
   @printf("x2 = %g;  r3 = %g;  x3 = %g\n", x2, r3, x3)
   plot()
   circle(0, r1, r1)
   circle(x2, r2, r2, :blue)
   circle(x3, r3, r3, :magenta)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       point(0, r1, " 大円:r1,(0,r1)", :red, :left, :vcenter)
       point(x2, r2, "中円:r2,(x2,r2)", :blue, :center, delta=-delta/2)
       point(x3, r3, "小円:r3,(x3,r3)", :magenta, :center, delta=-delta/2)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
       hline!([0], color=:black, lw=0.5)
       plot!(showaxis=false)
   end
end;

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

算額(その362)

2023年08月04日 | Julia

算額(その362)

山形県上山市金山峠 檜下不動堂
山形算額勝負-湯殿山神社を目指せ-

https://www.sci.yamagata-u.ac.jp/wasan/pdf/20180711SSEP.pdf

円内に同じ大きさの正方形が3個入っている。円の直径が 10 寸のとき,正方形の一変の長さを求めよ。

円の半径を r,正方形の一辺の長さを a とする。下にある正方形の下辺の y 座標を -b とする。
以下の連立方程式を解く。

include("julia-source.txt");

using SymPy

@syms r::positive, a::positive, b::positive;

eq1 = a^2 + b^2 - r^2
eq2 = (a/2)^2 + (2a - b)^2 - r^2;
res = solve([eq1, eq2], (a, b))

   1-element Vector{Tuple{Sym, Sym}}:
    (5*sqrt(17)*(-13*sqrt(17)*r/85 + r)*(13*sqrt(17)*r/85 + r)/(16*r), 13*sqrt(17)*r/85)

a は簡約化される。
a = 16*sqrt(17)*r/85

a = res[1][1] |> simplify
a |> println

   16*sqrt(17)*r/85

r = 10/2 のとき,a = 3.880570000581328 である。

r = 10/2
16*sqrt(17)*r/85

   3.880570000581328

using Plots

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r = 10/2
   (a, b) = (5*sqrt(17)*(-13*sqrt(17)*r/85 + r)*(13*sqrt(17)*r/85 + r)/(16*r), 13*sqrt(17)*r/85)
   @printf("a = %g;  b = %g\n", a, b)
   plot()
   circle(0, 0, r)
   rect(0, -b, a, a - b, :blue)
   rect(0, -b, -a, a - b, :blue)
   rect(-a/2, a - b, a/2, 2a - b, :blue)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       point(0, -b, " -b", :blue, delta=-delta/2)
       point(r, 0, " r", :red, :left, :bottom, delta=delta/2)
       point(a, 0, " a", :blue, :left, :bottom, delta=delta/2)
       point(a/2, 0, " a/2", :blue, :left, :bottom, delta=delta/2)
       point(0, a - b, " a-b", :blue, :left, :bottom, delta=delta/2)
       point(0, 2a - b, " 2a-b", :blue, :left, :bottom, delta=delta/2)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

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

算額(その361)

2023年08月04日 | Julia

算額(その361)

山形県山形市旅篭町 湯殿山神社
山形算額勝負-湯殿山神社を目指せ-

https://www.sci.yamagata-u.ac.jp/wasan/pdf/20180711SSEP.pdf

正方形の内部に 2 本の斜線と大円,中円,小円,甲円,乙円がある。中円の直径が与えられたときに甲円の直径を求めよ。


大円の半径と中心座標を r0, (0, r0)
中円の半径と中心座標を r1, (0, r1)
小円の半径と中心座標を r2, (x2, y2)
甲円の半径と中心座標を r3, (x3, y3) および (x5, y5)
乙円の半径と中心座標を r4, (r0 - r4, r4)
正方形の一辺の長さは 2r0 である。斜線と正方形の下辺の交点の x 座標を a とする。

r1 = 1 として,以下の連立方程式を立て,nlsolve() で数値解を求める。

include("julia-source.txt");

using SymPy

@syms r0::positive, a::positive, r1::positive,
     r2::positive, x2::positive, y2::positive,
     r3::positive, x3::positive, y3::positive,
     r4::positive, x5::positive, y5::positive;

r1 = 1
eq1 = x2^2 + (y2 - r0)^2 - (r0 - r2)^2
eq2 = x3^2 + (y3 - r0)^2 - (r0 - r3)^2
eq3 = (x2 - x3)^2 + (y2 - y3)^2 - (r2 + r3)^2
eq4 = 2(r0 - r4)^2 - (r0 +r4)^2
eq5 = r1^2*(4r0^2 + a^2) - a^2*(2r0 - r1)^2
eq6 = 2r0*(y2 - r0) - a*x2
eq7 = (r0 - 2r2)*(2r0 - r1) - r1*r0
eq8 = distance(0, 2r0, a, 0, r0 - r4, r4) - r4^2
eq9 = distance(0, 2r0, a, 0, x3, y3) - r3^2;
eq10 = x5^2 + (y5 - r0)^2 - (r0 -r3)^2
eq11 = distance(0, 2r0, a, 0, x5, y5) - r3^2;

# solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8, eq9])

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=1e-14)
       v = r.zero[1]
   else
       r = nlsolve((vout, vin)->vout .= func(vin, params...), ini, ftol=1e-14)
       v = r.zero
   end
   return v, r.f_converged
end;

function H(u)
   (r0, a, r2, x2, y2, r3, x3, y3, r4, x5, y5) = u
   return [
       x2^2 + (-r0 + y2)^2 - (r0 - r2)^2,  # eq1
       x3^2 + (-r0 + y3)^2 - (r0 - r3)^2,  # eq2
       -(r2 + r3)^2 + (x2 - x3)^2 + (y2 - y3)^2,  # eq3
       2*(r0 - r4)^2 - (r0 + r4)^2,  # eq4
       -a^2*(2*r0 - 1)^2 + a^2 + 4*r0^2,  # eq5
       -a*x2 + 2*r0*(-r0 + y2),  # eq6
       -r0 + (r0 - 2*r2)*(2*r0 - 1),  # eq7
       -r4^2 + (-2*r0*(a^2 - a*r0 + a*r4 + 2*r0*r4)/(a^2 + 4*r0^2) + r4)^2 + (-a*(a*r0 - a*r4 + 4*r0^2 - 2*r0*r4)/(a^2 + 4*r0^2) + r0 - r4)^2,  # eq8
       -r3^2 + (-a*(a*x3 + 4*r0^2 - 2*r0*y3)/(a^2 + 4*r0^2) + x3)^2 + (-2*r0*(a^2 - a*x3 + 2*r0*y3)/(a^2 + 4*r0^2) + y3)^2,  # eq9
       x5^2 + (-r0 + y5)^2 - (r0 - r3)^2,  # eq10
       -r3^2 + (-a*(a*x5 + 4*r0^2 - 2*r0*y5)/(a^2 + 4*r0^2) + x5)^2 + (-2*r0*(a^2 - a*x5 + 2*r0*y5)/(a^2 + 4*r0^2) + y5)^2,  # eq11
   ]
end;

iniv = [big"2.1", 1.41, 0.67, 1.26, 2.44, 0.44, 1.41, 1.34, 0.34, 0.68, 3.4]
res = nls(H, ini=iniv);
println(res);

  (BigFloat[2.000000000000000000000000000320970390259289336928171646885307055660835194112333, 1.414213562373095048801688724010844988825749298548807658294693530490981696954068, 0.6666666666666666666666666668515396947153689160851700548263921689551120494197228, 1.257078722109417821157056643922483726053671707035922941395627750409082770279752, 2.444444444444444444444444444692495891824273215136865120549810095106586884337016, 0.4444444444444444444444444445433658437202953697363922783119495000794485144891305, 1.410452971059059754400848547640885881999130012137408290274003768454414456979157, 1.343969891811035677909661427194763924578403051134536294523138065736756827990007, 0.3431457505076198047932451032161729608070508078480639448376432554356357666557888, 0.6846782324566366141942458588087396030640151534840055426317378671228636978868817, 3.396770848929705062831079314009928699095283423608233263143025331890301809943354], true)

   r0 = 2;  a = 1.41421;  r2 = 0.666667;  x2 = 1.25708;  y2 = 2.44444;  r3 = 0.444444;  x3 = 1.41045;  y3 = 1.34397;  r4 = 0.343146;  x5 = 0.684678;  y5 = 3.39677

中円の直径が1のとき,甲円の直径は 0.4444(4/9) である。

using Plots

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r1 = 1
   (r0, a, r2, x2, y2, r3, x3, y3, r4, x5, y5) = res[1]
   @printf("r0 = %g;  a = %g;  r2 = %g;  x2 = %g;  y2 = %g;  r3 = %g;  x3 = %g;  y3 = %g;  r4 = %g;  x5 = %g;  y5 = %g\n", r0, a, r2, x2, y2, r3, x3, y3, r4, x5, y5)
   #@printf("甲円の直径が %g 寸のとき,乙円の直径は %g 寸\n", 2r1, 2r2)
   plot([r0, r0, -r0, -r0, r0], [0, 2r0, 2r0, 0, 0], color=:black, lw=0.5)
   circle(0, r0, r0)
   circle(0, r1, r1, :brown)
   circle(x2, y2, r2, :blue)
   circle(x3, y3, r3, :magenta)
   circle(x5, y5, r3, :magenta)
   circle(r0 - r4, r4, r4, :green)
   plot!([a, 0, -a], [0, 2r0, 0], color=:black, lw=0.5)
   if more
       point(0, r0, " r0", :red, :left, :bottom)
       point(0, r1, " r1", :brown)
       point(x2, y2, "(x2,y2)")
       point(x3, y3, "(x3,y3)", :magenta)
       point(x5, y5, "(x5,y5)", :magenta)
       point(r0 - r4, r4, "(r0-r4,r4)")
       point(a, 0, "a ", :black, :right, :bottom)
       point(r0, 0, " r0", :black, :left, :bottom)
       point(0, 2r0, " 2r0", :black, :left, :bottom)
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

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

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

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