計算数学演習第１０回 – Masashi Shiraishi Homepage

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 で作成したグラフを印刷したい場合には，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+\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$ の近似値を求めるプログラムを作成せよ。
右辺の項の数と $\pi$ の近似値との関係を gnuplot で可視化せよ.
（ n: 1 ～ 1000 程度までで、おおよその傾向を見てください。）

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

課題３

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

$X_{n+1} = 1 + \frac{1}{X_n}$

$X_0=1$

（ $i \rightarrow \infty$ で $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 \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()を使う）

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

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

課題おまけ

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

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)とその多項式展開とを比べてみてください.

今日学んだ事

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