計算物理学２

第 4 回演習

ver. 2018/5/10
以下の例題を Fortran もしくは C 言語を用いてプログラム化せよ。
ファイル入出力（１５） 配布資料のページに載せてある matrixa.dat と matrixb.dat から３ｘ３行列をそれぞれ amat, bmat として読み込み、それらの行列の積をファイル matrixc.dat に出力
・matrixa.dat と matrixb.dat をダウンロードして計算物理学 2 のディレクトリに移動
・テキストファイルのファイルの中身を確認するには cat コマンド$ cat matrixa.dat
1.000 2.000 3.000-1.000 1.000 2.000
2.000 1.000 -2.000$ cat matrixb.dat-1.9000

2.00002.00005.0000 -4.9000 -5.0000-7.0000

7.00007.1000とするとファイルの中身が表示される。ファイルサイズが大きすぎて一画面に収まらない
ような場合は less コマンドを使う。
$ less matrixa.dat
十字キーで順に表示できる。less を終了する場合は q とタイプ。
・matrixc.dat ファイルが作成されたかどうかはプログラムを実行したあとに ls コマンドで確認、cat コマンドでファイルの中身を確認する。
・emacs で開いても構わないが、ファイルの編集をせず表示するだけの場合は cat や lessを使うことが多い。
グラフ作成（１６）x が 0 から 10 まで 0.1 刻みで x, exp(-x/2), (sin 3x)*exp(-x/2), (cos 3x)*exp(x/2) を ex16data.dat ファイルに出力。
一番上の行にx が 0.0 のときの x, exp(-x/2), (sin 3x)*exp(-x/2), (cos 3x)*exp(-x/2)
x が 0.1 のときの x, exp(-x/2), (sin 3x)*exp(-x/2), (cos 3x)*exp(-x/2)
x が 0.2 のときの x, exp(-x/2), (sin 3x)*exp(-x/2), (cos 3x)*exp(-x/2)
1

x が 0.3 のときの x, exp(-x/2), (sin 3x)*exp(-x/2), (cos 3x)*exp(-x/2)...
x が 10.0 のときの x, exp(-x/2), (sin 3x)*exp(-x/2), (cos 3x)*exp(-x/2)というように一行に４つの数字を、x が 0 から 10 まで 0.1 刻みで(およそ 100 行程度)出力。cat や less コマンドで作成された ex16data.dat ファイルの中身を確認したあと、
gnuplot を 用 いて こ のフ ァイ ル を読 み込 ん で exp(-x/2), (sin 3x)*exp(-x/2), (cos3x)*exp(-x/2)の曲線を一枚のグラフに表示。
（１７）ロジスティック写像はxi+1 = axi (1<latexit sha1_base64="hfBBL9oLewMZtT3JwRlhodzYcAU=">AAACnHichVFfKz1RFF3G/+uHiwdKaXIjktseKVJKeFBS/t1LoWlmHJzMnZlm5t5w8wV8AQ+ekAfh2Qfw4gt48BHkkfLi4bfv3Ckh7OnMXmedvfY5q216tgxCoscKpbKquqa2rj7R8K+xqTnZ0poN3LxviYzl2q6/ahqBsKUjMqEMbbHq+cLImbZYMXenSucrBeEH0nWWw31PbOSMbUduScsImdKT7Xt6UQ5oh+q4aqh7ulT7tEFO/XoyRWmKQv0OtBikEMe8m7zFOjbhwkIeOQg4CBnbMBDwtwYNBI+5DRSZ8xnJ6FzgEAnW5rlKcIXB7C7/t3m3FrMO70s9g0ht8S02L5+VKnrogS7phe7pip7o/cdexahH6S37nM2yVnh681HH0tufqhznEDsfql8UJlf/7inEFkYjL5K9eRFTcmmV+xcOjl+WxhZ7ir10Rs/s75Qe6Y4dOoVX62JBLJ4gwQPSvo7jO8gOpTVKawvDqYnJeFR16EQ3+ngeI5jADOaRie49xzVulC5lWplV5sqlSkWsacOnULL/AXn5l2U=</latexit>

xi )で与えられる。初期値 x1 は 0 から 1 の間の適当な値とする。
パラメータ a の値を 1 から 4 までの間を 0.001 刻みで変え、各 a の値において、x500 から x600 までの値をaｘ500

aｘ501

...
aｘ600

a+0.001 x500
a+0.001 x501...
a+0.001 x600
a+0.002 x500...
として logistic.dat ファイルに出力。gnuplot で横軸 a、縦軸は xi としてプロット(分岐図)
ファイルサイズが大きくなるため、描画に時間がかかることがある。プロットするときはデータを線で結ばず サイズの小さい丸(pointtype 7 pointsize 0.1)を用いた
plot "logistic.dat" using 1:2 with points pointtype 7 pointsize 0.1を使う。
a が 3.5 付近にフラクタル構造が見えるがこの領域を拡大してみる。点数が足らない場合
は領域を区切って刻みを小さくしてプロットする(例えば 3.568 から 3.571 までを-710 刻みで見る)カオスはどのあたりから始まっているか。

2（１8）sin x の x=0 でのテーラー展開を 1 次から 11 次まで行い、組み込み関数で計算した SIN(x)とテーラー展開したものをグラフに表示。x がどのあたりまで近似が良いかを示す。ただし sin x の x=0 でのテーラー展開は
x3x5

sin x = x<latexit sha1_base64="PMk5EGs89LgehzNVhg7kswxpuhs=">AAACzHichVFNSxxBEH1O4tfG6Ca5BLxMshgCwaVGIxsCAYkHPYlfq4Jrlpm2VxtnZ4aZ3kUzzNVD/oAHTwl4kODZH+Ald/Eg6A8IOSrkkkNqZ4dIItFqurv69XvV/SgncFWkic46jHv3O7u6e3pzD/oe9g/kHz1ejPxGKGRZ+K4fLjt2JF3lybJW2pXLQSjtuuPKJWdzonW/1JRhpHxvQW8HcrVur3uqpoStGarmJyuR8swt8x1Pc9is1EJbxFsfRpN49FlivvoDjCXxGAPXjFISl9oMsebrqJovUJHSMG8mVpYUkMWMnz9CBWvwIdBAHRIeNOcubEQ8VmCBEDC2ipixkDOV3kskyLG2wSzJDJvRTV7X+bSSoR6fWzWjVC34FZdnyEoTQ3RKB3RJ3+grfadf/60VpzVaf9nm3WlrZVAd+PR0/uedqjrvGhvXqlsUDrNv96RRw5vUi2JvQYq0XIp2/ebH3cv5t3ND8Qv6Qj/Y32c6o2N26DWvxP6snNtDjhtk/duOm8niSNGiojX7ujD+PmtVDwbxHC+5HyWMYwozKPO7hzjBOS6MaUMbsZG0qUZHpnmCv8LY+Q1lPKm7</latexit>3!

+で与えられる。

35!

x7+ ···
7!