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

今日の演習

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

今日の目標

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


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

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

課題３で載せているグラフには、タイトルと軸、各曲線に名前がついてます．このように名前をつけるには、次のようにします.(データファイル名はそれぞれdata1、data2、data3、data4、とした.)

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)で"第三列から第二列を引いたもの"を描画することも可能です。
several


対数表示
指数的に大きくなる値は，そのまま描画するより
縦軸を対数座標にしたほうがわかりやすいことがある．
exp(x)
を縦軸のみ対数座標で描画するためには，
set logscale y
plot exp(x)
とする．

元のｙ座標に戻したい場合は，
unset logscale
で元に戻せる．

グラフの印刷の仕方

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

gnuplot 起動中に

set term postscript eps
set output "sin_graph.eps"
plot　sin(x)

set term x11

set output

 

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

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

作成した psファイルをファイルブラウザでダブルクリックすると、適当なpsファイルビューアが起動するので、それから印刷してください．（最寄りのプリンターから印刷されます.）ちなみに ps ファイルとは PostScriptファイルの略です．PostScriptとはプリンタの制御言語で，プリンタに絵を描かせる命令の集まりです．

 

色付きで保存したい場合：

最初の文を

set term postscript enhanced color

と入力

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

課題１

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

課題２

次の等式 π/4 = 1 - 1/3 + 1/5 - 1/7 + ... + (-1)n+1/(2n-1) + ... （ライプニッツ-グレゴリーの等式）を利用して、π の近似値を求めるプログラムを作成せよ。右辺の項の数とπの近似値 との関係を gnuplot で可視化せよ. （ n: 1 〜 1000 程度までで、おおよその傾向を見てください。）

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

課題３

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

X_{n+1} = 1 + 1/X_{n}

X_{0} = 1

（ i → ∞ で X{i} はいわゆる黄金比に収束します。収束の様子を確認してください。）

課題４ (Logistic map)

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

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

A < 1 なら X_i  → 0 （i → ∞）となるのは、すぐ分ると思います.

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},..,X{49} を描画せよ．Aの値により，一つの値に収束する場合，二つの値の間を振動する場合，など様々な振る舞いを示すことを確認せよ．


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

課題４の式で、A=2.0, 2.001,2.002,..,3.999,4.0の場合について
X_{1000},..,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.60001のとき）

となる数値データを出力し，

横軸が一列目の値，

縦軸が二列目と三列目の差の絶対値を出力せよ．

(gnuplot で絶対値を計算するには， abs()を使う）

また，縦軸を対数表示にすることにより，二つの値の離れ方を確認してみよ．

（興味のある人は，初期値鋭敏性というキーワードを調べてみてください）


課題おまけ

次の関数を描画してみてください.

plot [-7:7][-2:2] x-(x**3)/3/2/1+(x**5)/5/4/3/2/1-(x**7)/7/6/5/4/3/2/1+(x**9)/9/8/7/6/5/4/3/2/1-(x**11)/11/10/9/8/7/6/5/4/3/2/1+(x**13)/13/12/11/10/9/8/7/6/5/4/3/2/1, sin(x)

(長いですが、行を変えずに続けて打ち込んで下さい.) これはsin(x)のTaylor展開です. ちょっと長いですが、展開項数を変つつ、sin(x)とその多項式展開とを比べてみてください.

今日学んだ事

戻る