計算数学演習第１０回

gnuplotを用いたデータのグラフ描画

今日の演習

今日は以下の内容で演習を行います。

  • Cプログラムとgnuplotを使った、グラフ、図形の描写 のつづき

今日の目標

gnuplot で計算結果を可視化する。


Gnuplotのあれこれ（軸、ラベル等）

タイトル(title)、軸の名前(xlabel , ylabel)、各データの名前

グラフにはタイトルや軸、各曲線に対して名前をつけることができます。
実際に名前をつけるには、次のようにします。

gnuplot を起動

set title "Sin waves"
set xlabel "x"
set ylabel "F(x)"

　ここでは１行目でグラフの名前を、２行目で横軸(x軸)の名前を、３行目で縦軸(y軸)の名前を設定しています。

　これで、図全体のタイトルと、各軸の名前の設定が終わりました。

plot [0:2*pi][-1.5:1.5] sin(x) title "wave 1", sin(2*x) title "wave 2", sin(3*x) title "wave 3"

　描写の際、最後に title “名前” と付けることで、各曲線に自由な名前が付けられます。
これはもちろんデータファイルを描画する際にも有効です。

実際に入力して結果を確かめてください。

関数の定義

また、gnuplotでは自分で関数を定義することも可能です。

f(x,a) = sin(a*x)
plot f(x,1), f(x,2)

という使い方も可能なので、各自調べてみてください。

再描画

replot

と入力すると、一番最後に入力した描画（今回の範囲ではplotの文）の命令文が実行されます。


データの値をそのまま出力するのではなく、さらに処理を行った上で出力する方法

“第二列の値に0.5を加えたもの”や”第二列の値を2倍したもの”を描画したいとき

plot [0:2*pi] "N=50.data" using 1:($2+0.5) w lp, "N=50.data" using 1:($2*2) w lp, sin(x)

と打ち込むと、下図のような結果が出ます。
これと同様の方法でusing 1:($3-$2)で”第三列から第二列を引いたもの”を描画することも可能です。

 

対数表示

指数的に大きくなる値は、そのまま描画するより
縦軸を対数座標にしたほうがわかりやすいことがあります。

exp(x)を縦軸のみ対数座標で描画するためには、

set logscale y
plot exp(x)

とする。
元のｙ座標に戻したい場合は、

unset logscale

で元に戻せる。

グラフのファイル出力の仕方

gnuplot で作成したグラフをファイルに出力する場合には、pdf ファイルに保存すると便利です。
その場合の手順は次のようになります。
（例えば sin(x) のグラフが描かれた pdf ファイル 「sin_graph.pdf」 を作成する場合.）

gnuplot 起動中に

set term pdf
set output "sin_graph.pdf"
plot sin(x)

とすれば「sin_graph.pdf」 というpdfファイルができあがります。
ここでは、一行目で出力形式の変更、２行目で出力先のファイル名の指定、をしています.
ファイル名や plot の部分を適宜変更して下さい.
（set output … と plot … の間で set title … 等 (上の「gnuplot 付け足し」参照.) をすると、タイトル、軸の名前付きの ps ファイルができます。）

この状態でさらにplotするとおかしなことになりますので、一旦出力を切るために

set out

とファイル名を空にして指定しておくと良いでしょう。

再度グラフを画面上に出力したい場合は、

set term qt

とすれば戻ります。(もし、gnuplot自体を終了していいのなら、set term qtの代わりにquitと入力してくれてかまいません)

余談ですが、以前は eps ファイルに出力することが多かったです。epsとは Encapsulated PostScriptファイルの略です。
PostScriptとはプリンタの制御言語で、プリンタに絵を描かせる命令の集まりです。
この場合は

set term post enh col eps

と出力形式を指定するとepsファイルができます（ファイルの拡張子もepsにすること）。
png形式の場合は、基本的に

set term png

で出力できますが、PCの環境に依存してサポートされていないことがあります（ファイルの拡張子もpngにすること）。

数列の収束の様子をみる.その１

課題１

次の数列X_{n}=(1+\frac{1}{n})^nを、n: 1 ～ 200 程度まで求めるプログラムを作成し、nとX_n の関係を gnuplot で可視化せよ.
（n\rightarrow \inftyでX_n \rightarrow e（ネイピア数）ですね。収束の様子を確認してください。）

課題２

次の等式（Leibniz-Gregoryの等式）

    \[\frac{\pi}{4} = 1 - \frac{1}{3} +\frac{1}{5}-\frac{1}{7}+ \cdots +\frac{(-1)^n}{2n+1}+ \cdots\]

を利用して、\piの近似値を求めるプログラムを作成せよ。
右辺の項の数と\frac{\pi}{4}の近似値 との関係を gnuplot で可視化せよ.
（ n: 0 ～ 1000 程度までで、おおよその傾向を見てください。）

数列の収束の様子をみる.その２（漸化式）

課題３

次の漸化式で表される数列X_nを求めるプログラムを作成し、n (0 ~ 100)とX_nの関係を gnuplot で可視化せよ。

    \[X_{n+1} = 1 + \frac{1}{X_n}\]

    \[X_0=1\]

（ n \rightarrow \inftyでX_nはいわゆる黄金比に収束します。収束の様子を確認してください。）

課題４　Logistic map

次の漸化式で表される数列 (0<X_i<1) を求めるプログラムを作成し、iとX_iの関係を gnuplot で可視化せよ.

    \[X_{i+1} = A*X_{i}*(1.0-X_{i})\]

A < 1 ならX_i \rightarrow 0(i \rightarrow \infty)となるのは、すぐ分ると思います.

A = 0.5, 2.0, 3.2, 3.5, 3.9 等（A : 0 ～ 4）それぞれの場合について、初期値X_0を0.3, 0.6, 0.61 等（0 ～ 1）と変化させ、実行、描画して見比べてください.　X_{0}, \dots, X_{49}を描画せよ。Aの値により、一つの値に収束する場合、二つの値の間を振動する場合、など様々な振る舞いを示すことを確認せよ。

課題５　Logistic mapの分岐ダイアグラム

課題４の式で、A=2.0, 2.001,2.002,..,3.999,4.0の場合についてX_{1000}, \dots, X_{1050}を計算し、
一列目をAの値、二列目をそのAにおけるX_nとして、以下のような形式で出力するプログラムを作成せよ。
ただし、どのAにおいてもX_0=0.1であるものとする、
出力例：
2.000 X_{1000}
..
2.000 X_{1050}
2.001 X_{1000}
…
2.001 X_{1050}
2.002 X_{1000}
…
4.000 X_{1050}
（合計 2001*51行のデータ）
このデータをもとにgnuplotでグラフを作成してみること。
散布図( データファイルをdiagram.txtとした場合、plot “diagram.txt” w d　と書くと綺麗に見える。）

課題６　n周期解

課題４の漸化式で与えられるx_iで8周期解、3周期解、６周期解に収束するものは存在するだろうか。
存在するなら、それぞれの場合に対応するaの値を見つけ、横軸i　縦軸x_iのグラフを描画してみよ。
（初期値x_0は0.1とする。）
ただし、ここでいうn周期解とはx_i=x_{i+n}を満たす解である、
（a=3.2のときはiが大きいときこの性質を持つ解に収束するので、２周期解に収束するといえる。）

課題７　カオス入門（発展問題）

課題４の式で、A=3.9として、
一列目　i
二列目　X_i（X_0=0.6のとき)
三列目　X_i（X_0=0.600001のとき）
となる数値データを出力し、横軸が一列目の値、縦軸が二列目と三列目の差の絶対値を出力せよ。
(gnuplot で絶対値を計算するには、 abs()を使う）

また、縦軸を対数表示にすることにより、二つの値の離れ方を確認してみよ。
（興味のある人は、初期値鋭敏性というキーワードを調べてみてください）

Gnuplotのおまけ

描画設定ファイルの作成

色々設定した内容は、save "plot_setting.gp"とすれば保存でき（ファイル名はなんでも良い）、load "plot_setting.gp"とすれば呼び出すことができる。同じ設定で複数の画像ファイルを作成する際などに便利（ファイルの中に何が書かれているのか、テキストエディタで確認しよう。）

for文、if文の定義

実はgnuplotでもfor文やif文の定義ができる。

例：

f(x,a)=x>0?sin(a*x):0
p for [i=1:5] f(x,i) t sprintf("a=%d",i)

1行目が三項演算子（if文に相当）を使った関数の定義
2行目がfor文によるグラフの描画（titleにも数値を入れることができる）

他にも様々なことができるので、各自調べてみてください。

今日学んだ事

  • gnuplot とCプログラムでグラフ、図形を描画しました。

Department of Mathematical and Life Sciences, Graduate School of Integrated Sciences for Life, Hiroshima University