45 Julia で遊ぶ (2)

何か意味のあるプログラムを書かないと前進できない、と思ったので、 Julia で常微分方程式の初期値問題を解いてみることにした。

(埋め込み型公式のまともなプログラムとか、構造保存型数値解法とか、 本当は、常微分方程式で色々勉強しなければいけないことがたまっている。 ゼミ生誰かやらないかなあ、とか考えている。)

安直にボール投げ (このところ試しプログラムはこれを使っている)。 計算部分は、まあまあ納得出来た (http://nalab.mind.meiji.ac.jp/~mk/labo/text/julia-memo/node46.html)。

http://nalab.mind.meiji.ac.jp/~mk/misc/20191228/ballbound.jl

# ballbound.jl --- 投げたボールのバウンド, 空気抵抗ありのシミュレーション

# Runge-Kutta 法の1ステップ
function rungekutta(f,t,x,dt)
  k1=dt*f(t,x)
  k2=dt*f(t+dt/2, x+k1/2)
  k3=dt*f(t+dt/2, x+k2/2)
  k4=dt*f(t+dt, x+k3)
  x + (k1 + 2 * k2 + 2 * k3 + k4) / 6
end

function f(t,x)
  Gamma=1.0
  m=100.0
  g=9.8
  y=similar(x)
  y[1] = x[3]
  y[2] = x[4]
  y[3] = -Gamma/m*x[3]
  y[4] = -g-Gamma/m*x[4]
  y
end

function ballbound(n=1000) # 1000等分くらい
  t=0.0
  v0=50.0
  theta=50.0
  x=[0.0,0.0,v0*cos(theta*pi/180),v0*sin(theta*pi/180)]
  println("t x")
  println("$t $x")
  Tmax=20.0
  dt=Tmax/n
  println("Tmax=$Tmax, dt=$dt")
  s=@sprintf "%f %f %f %f %f\n" t x[1] x[2] x[3] x[4]
  print(s)
  of = open("ballbound.dat","w")
  print(of,s)
  for i=1:n
    # x=euler(f,t,x,dt)
    x=rungekutta(f,t,x,dt)
    if x[2]<0
      x[2] = - x[2]
      x[4] = - x[4]
    end
    t=i*dt
    s=@sprintf "%f %f %f %f %f\n" t x[1] x[2] x[3] x[4]
    print(s)
    print(of,s)
  end
  close(of)
  # 以下は工事中
  of = open("ballbound.gp","w")
  println(of,"plot \"ballbound.dat\" using 2:3 with lp")
  close(of)
  run(`gnuplot ballbound.gp`)
end

ターミナルで

curl -O http://nalab.mind.meiji.ac.jp/~mk/misc/20191228/ballbound.jl
julia

Julia の中で

using Printf
include("ballbound.jl")
ballbound(1000)

ボールの軌跡を描くためにとりあえず gnuplot, と考えたのだけど (学生は意外と gnuplot に慣れていることが分かった、 と言うのもある)、 Julia で pipe をどう扱えば良いのか分からなかったので (笑)、 原始的にまず全部データをファイルに書いて、 それから gnuplot を起動している。 安直 (締め切りに間に合わせるため)。

冬季休暇に入って数日、 やっと仕事が一段落したので、まともな可視化について考え始める。 最初は予定になかったけれど、(1) pipe で gnuplot をコントロールする、 というのをやって、それと (2) Plots を使って、 この2つくらいはとりあえず目処をつけよう。

(1) gnuplot に描かせる

ずっと以前 gnuplot を調べたとき、何とか見つけられたことは (「Cから gnuplot を呼び出す」)、 今では常識的なことになっているみたいだ。

暮れに出席した講習会 (http://www.isc.meiji.ac.jp/~akiyama_masakazu/#section16) のスライドから、

関数のグラフ描画 (C プログラムから gnuplot を呼び出す)

  FILE *gp;//For gnuplot
  gp = popen("gnuplot --persist","w");
  fprintf(gp, "set terminal x11\n");
  fprintf(gp, "plot sin(x)\n");
  fflish(gp);
  pclose(gp);

数値データのプロット (C プログラムから gnuplot を呼び出す)

  FILE *gp;//For gnuplot
  gp = popen("gnuplot --persist","w");
  fprintf(gp, "set terminal x11\n");
  ...
  fprintf(gp, "set yrange [ここから:ここまで]\n");
  fprintf(gp, "plot '-' with liens\n");

  ...
  fprintf("gp, "%f %f\n", なんとか, かんとか);
  ...

  fprintf("gp, "e\n");
  fflish(gp);
  pclose(gp);

Julia の version 0.6 くらいまでは、 readandwrite() とかいう関数があったらしいけれど、 今はなくなっていて、 どうすれば良いか説明した明快な文書が見つからなかった。

難しいときは、簡単なところから (Julia から bc を使う)

  p=open(pipeline(`bc`; stderr=Pipe()), "r+")

  nb=write(p.in,"2*3\n")
  s=readline(p.out)

原始的な電卓プログラム bc に、2*3 を入力して、 返事を読み取る。s には "6" が返る。

(bc を実際的な意味で使う人は今では少ないでしょうか。 bc には、良くこういう役が回ってきますね。 「たまに呼ばれると、いつもこんな役なんですよ」とか言いそう。)

では、gnuplot でやってみる。 関数のグラフを描く方は簡単で、例えば次のようにすれば良い。

関数のグラフ描画 (コピペすれば動くはず)

  p=open(pipeline(`gnuplot`; stderr=Pipe()), "r+")
  write(p.in,"plot sin(x)\n")

C に慣れていて、open(), read(), write() を使ったサンプルを見ると、 open() で返ってきた p に対して、 原始的なことしか出来ない気がするけれど、 println(p.in, "plot sin(x)") でもOKみたい (この辺は、一度ちゃんと調べないとダメな気がする)。

数値データのプロットについて、最初に成功したのは次のようなコード。

数値データのプロット (コピペすれば動くはず)

  using Printf
  p=open(pipeline(`gnuplot --persist`; stderr=Pipe()), "r+")
  write(p.in,"plot '-' with lp\n")
  n=100
  dx=2*pi/n
  for i=0:n
    x=i*dx
    @printf(p.in, "%f %f\n", x, sin(x)+sin(3*x)/5.0)
  end
  s=@sprintf "e\n";
  write(p.in,s);

最後は flush(p.in) とか、close(p) とかするのかな？

人間が読まないのならば、書式など不要だろうから (と考えるのだけれど、どうだろう？)、 次のようにしても良いかもしれない。 println() しか使っていない。 これなら using Printf も要らない。

  p=open(pipeline(`gnuplot --persist`; stderr=Pipe()), "r+")
  println(p.in, "plot '-' with lp")
  n=100
  dx=2*pi/n
  for i=0:n
    x=i*dx
    println(p.in,"$x $(sin(x)+sin(3*x)/5.0)")
  end
  println(p.in,"e")

ここまで進めて、以下が ballbound.jl の pipe バージョンである。

http://nalab.mind.meiji.ac.jp/~mk/misc/20191228/ballbound-v3.jl

# ballbound-v3.jl --- 投げたボールのバウンド, 空気抵抗ありのシミュレーション

# Runge-Kutta 法の1ステップ
function rungekutta(f,t,x,dt)
  k1=dt*f(t,x)
  k2=dt*f(t+dt/2, x+k1/2)
  k3=dt*f(t+dt/2, x+k2/2)
  k4=dt*f(t+dt, x+k3)
  x + (k1 + 2 * k2 + 2 * k3 + k4) / 6
end

function f(t,x)
  Gamma=1.0
  m=100.0
  g=9.8
  y=similar(x)
  y[1] = x[3]
  y[2] = x[4]
  y[3] = -Gamma/m*x[3]
  y[4] = -g-Gamma/m*x[4]
  y
end

function ballbound3(n=1000) # 1000等分くらい
  t=0.0
  v0=50.0
  theta=50.0
  x=[0.0,0.0,v0*cos(theta*pi/180),v0*sin(theta*pi/180)]
  println("t x")
  println("$t $x")
  Tmax=20.0
  dt=Tmax/n

  println("Tmax=$Tmax, dt=$dt")
  p=open(pipeline(`gnuplot --persist`; stderr=Pipe()), "r+")
  println(p.in, "plot '-' with lines title \"n=$n\"");

  println(p.in, "$(x[1]) $(x[2])")
  for i=1:n
    # x=euler(f,t,x,dt)
    x=rungekutta(f,t,x,dt)
    if x[2]<0
      x[2] = - x[2]
      x[4] = - x[4]
    end
    t=i*dt
    # @printf(p.in, "%f %f\n", x[1], x[2])
    println(p.in, "$(x[1]) $(x[2])")
  end
  println(p.in, "e")
  flush(p.in)
  println(p.in, "set term png");
  println(p.in, "set output \"ballbound3.png\"");
  println(p.in, "replot");
  println(p.in, "quit")
  flush(p.in)
  close(p)
end

Julia と gnuplot (と curl) が使えるならば、次のようにして試せる(はず)。

ターミナルで

curl -O http://nalab.mind.meiji.ac.jp/~mk/misc/20191228/ballbound-v3.jl
julia
include("ballbound-v3.jl")
ballbound3(1000)

図 2: Julia で常微分方程式を解いて、gnuplot を呼んで可視化

(2) Plots を使って

gnuplot は、(a) とりあえず描きたい場合、 (b) 他の手段で描けないときの最後の手段、のようなもので、 Plots を使って描けるようになりたい。

ballbound_v4.jl

# ballbound_v4.jl --- 投げたボールのバウンド, 空気抵抗ありのシミュレーション
#=
curl -O http://nalab.mind.meiji.ac.jp/~mk/misc/20200102/ballbound_v4.jl
julia
ENV["GKS_ENCODING"] = "utf-8"
using Plots
gr()
include("ballbound_v4.jl")
ballbound_v4()
=#

# Runge-Kutta 法の1ステップ
function rungekutta(f,t,x,dt)
  k1=dt*f(t,x)
  k2=dt*f(t+dt/2, x+k1/2)
  k3=dt*f(t+dt/2, x+k2/2)
  k4=dt*f(t+dt, x+k3)
  x + (k1 + 2 * k2 + 2 * k3 + k4) / 6
end

function f(t,x)
  Gamma=1.0
  m=100.0
  g=9.8
  y=similar(x)
  y[1] = x[3]
  y[2] = x[4]
  y[3] = -Gamma/m*x[3]
  y[4] = -g-Gamma/m*x[4]
  y
end

function ballbound_v4(n=1000) # 1000等分くらい
  t=0.0
  v0=50.0
  theta=50.0
  x=[0.0,0.0,v0*cos(theta*pi/180),v0*sin(theta*pi/180)]
  println("t x")
  println("$t $x")
  Tmax=20.0
  dt=Tmax/n

  println("Tmax=$Tmax, dt=$dt")
  xs=zeros(n+1)
  ys=zeros(n+1)

  #println("$(x[1]) $(x[2])")
  xs[1] = x[1]; ys[1] = x[2]
  for i=1:n
    # x=euler(f,t,x,dt)
    x=rungekutta(f,t,x,dt)
    if x[2]<0
      x[2] = - x[2]
      x[4] = - x[4]
    end
    t=i*dt
    # @printf("%f %f\n", x[1], x[2])
    xs[i+1]=x[1]; ys[i+1] = x[2]
  end
  p=plot(xs,ys,title="速度に比例する抵抗を受けるボール投げ",
         xaxis=("x",(0.0,xs[n+1])),yaxis=("z",(0.0,findmax(ys)[1])))
  display(p)
  println("図を保存する")
  savefig("ballbound.pdf")
  savefig("ballbound.png")
  display(p)
end

グラフを描くこと自体はとても簡単だったが、 日本語タイトルをつけようとして、ちょっと苦戦した。 Plots と GR バックエンドでは日本語表示は出来ない、という書き込みも見た。 ENV["GKS_ENCODING"] = "utf-8" とすれば良い、 という情報を得て、無事表示出来たのだけれど…

試してみよう

$ curl -O http://nalab.mind.meiji.ac.jp/~mk/misc/20200102/ballbound_v4.jl
$ julia
julia> ENV["GKS_ENCODING"] = "utf-8"
julia> using Plots
julia> gr()
julia> include("ballbound_v4.jl")
julia> ballbound_v4()

日本語文字列をタイトルにする: https://github.com/JuliaPlots/Plots.jl/issues/791

図は、タイトルも含めて画面にはきちんと描くことができるが、 保存するときに GKS: invalid bitmap size と表示されて、 タイトルが化けたり、空白になったりする。 以下の図は、savefig() で保存したものではない。

図 3: 空気抵抗を受ける場合のボール投げ

もう一息なのかな？