next up previous contents
Next: H..3 matplotlib で日本語タイトルを用いる Up: H. 本文中の図の描き方 Previous: H..1 図1 の描き方


H..2 図2の描き方

gnuplot を使って描いた。

  1. 次のようなファイルを作る。

    euler.data

    # FIRSTN,LASTN,r=#
10 1.245394e-01
20 6.498412e-02
40 3.321799e-02
80 1.679689e-02
160 8.446252e-03
320 4.235185e-03
640 2.120621e-03
1280 1.061069e-03
    runge-kutta.data

    # FIRSTN,LASTN,r=#
  10 2.084324e-06
  20 1.358027e-07
  40 8.666190e-09
  80 5.473058e-10
 160 3.438538e-11
 320 2.153833e-12
 640 1.358913e-13
1280 1.154632e-14

    euler.dat を作るには、 https://m-katsurada.sakura.ne.jp/labo/text/num-ode/node27.htmlにある reidai5-2.c が使える。 
    % cc reidai5-2.c
% ./a.out > euler.data
10 1280 2
%

    参考までそれを Python, Julia に書き換えたものを載せておく。
    euler-error.py

    # euler-error.py ---dx/dt=x (0<t<1), x(0)=1 を Euler 法で解く
import math

# 微分方程式 x'=f(t,x) の右辺の関数 f の定義
def f(t,x):
    return x

a=0.0; b=1.0
e=math.exp(1.0)
x0=1.0

str=input("N0,N1,r=")
s=str.split()
N0=int(s[0])
N1=int(s[1])
r=int(s[2])
print('#')

n=N0
while n <= N1:
  h=(b-a)/n
  # 開始時刻と初期値のセット
  t=a
  x=x0
  for i in range(0,n):
    x += h * f(t,x)
    t += h
  print('%4d %e' % (n,abs(e-x))) # C言語の printf("%4d %e\n", n, fabs(e-x))相当
  n *= r
    euler-error.jl

    # euler-error.jl --- dx/dt=x (0<t<1), x(0)=1 を Euler 法で解く
#  euler-error.c の真似

using Printf

function f(t,x)
  x
end

function compute_euler_error(N0, N1, r)
  a=0.0; b=1.0;
  e=exp(1.0);
  x0=1.0;
  n=N0;
  while (n <= N1)
    h=(b-a)/n;
    t=a;
    x=x0;
    for i=0:n
      x += h * f(t, x);
      t += h;
    end
    @printf("%4d %e\n", n, abs(e-x))
    n *= r;
  end
end

# julia euler-error.jl と実行した場合に実行する。
if abspath(PROGRAM_FILE) == @__FILE__
  s=split(Base.prompt("# FIRSTN,LASTN,r="))
  N0=parse(Int64,s[1]);
  N1=parse(Int64,s[2]);
  r=parse(Int64,s[3]);
  println("#")
  compute_euler_error(N0,N1,r)
end
    euler-error.py を実行する
    % python euler-error.py > euler.data
10 1280 2
% cat euler.data
# FIRSTN,LASTN,r=: #
  10 1.348349e-01
  20 6.768076e-02
  40 3.390861e-02
  80 1.697167e-02
 160 8.490220e-03
 320 4.246211e-03
 640 2.123381e-03
1280 1.061760e-03
%
    euler-error.jl を実行する
    % julia euler-error.jl > euler.data

    runge-kutta.dat を作るプログラムは公開していないが、 reidai5-2.c を書き換えれば良い。Runge-Kutta法のコードは、 https://m-katsurada.sakura.ne.jp/labo/text/num-ode/node28.htmlが参考になる (分かってしまえば簡単, ほぼ該当部分のコピペで済む)。

  2. gnuplot を起動して次のようにすれば図が(画面に)描ける。
    gnuplot で描く
    set logscale xy
plot "euler.data" with lp
replot "runge-kutta.data" with lp
    (lp の l は L の小文字。 lp は linespoints の略。)
    こうして表示されたグラフを保存するには、続いて
    gnuplot で描く (続き)
    set term pdf
set output "compare_data.pdf"
replot
    とすればよい。 ここでは PDF 形式で保存したが、 set term png とか set term postscript eps color とか、 形式は色々なものが選択可能である。

補足1     ここでは標準出力のリダイレクション (コマンド > ファイル名) を使って、 データ・ファイル euler.data を作成したが、 Cのプログラムの中で fopen(), fprintf(), flose() などの関数を使って直接ファイルを作成することも出来る (その方がずっと融通が効く)。 「常微分方程式の初期値問題を解くプログラムの書き方     2.3.1 Euler法のCプログラム例」 の問2を参考にすると良い(解答プログラム付き)。

補足2     Cのプログラムから直接 gnuplot を起動して、 データを渡してグラフを描くことも出来る。


桂田 祐史