2014計算機実験1_1

計算機実験1

自然情報学講座 瀬戸 繭美

seto@ics.nara-wu.ac.jp

数学モデリングのプロセス

問題点の抽出

定義・仮定

数式化

解析・数値計算

により解を得る

解釈

モデルの検証

モデルによる

説明・将来予測

解析・数値計算

により解を得る

数学モデリングのプロセス

問題点の抽出

定義・仮定

数式化

解釈

モデルの検証

モデルによる

説明・将来予測

画像引用_{: http://www.matumura-ringo.com/info.html} すべての物体は引き合う

M

m

r

F = G

mM

_r

₂ F: 万有引力 (kg m s-2_)! G: 万有引力定数 (m3 _s-2 _kg -1_)! M: 地球の質量 (kg)! m: 落下する物質の質量 (kg)! r: 2物質間の距離 (m)

万有引力の法則

法則と微分方程式

ニュートンの運動方程式

ma = F

物体の加速度 a は物体の質量 m に

反比例し、物体に働く力 F に比例する

m

d

2

x

dt

2

= F

2階常微分方程式

x(t): 位置

m

F

v(t) =

dx

dt

a(t) =

dv

dt

=

d

2

x

dt

2 Fの力で質量mの物体を押したら t秒後の速度vはどれくらい変わる?

速度

加速度

運動方程式を解くと何が解るか

?

x

x(0)

mg

ma(t) = mg

自由落下の運動方程式

v(t) = gt + v(0)

x(t) =

1

2

gt

2

+ v(0)t + x(0)

速度

位置

ある時点_{(初期条件)における} 速度と位置の情報が与えられれば ある時刻tにおける物体の速度と位置がわかる t = 0 において x = x(0), υ = υ(0)とし、   を t で積分すると

g: 重力加速度

d2x dt2 = g 時刻 t においてボールは どこ(x(t))にあるだろう?

常微分方程式

f (t,y,y ,y ,...,y

(n)

) = 0

f (t,y,

dy

dt

, ...,

d

n

y

dt

n

) = 0

n 階の微分を含むものは次のように表すことが可能

•

導関数の最高階数n をこの微分方程式の階数という

•

t を独立変数, yを従属変数と呼ぶこともある

•

2 つ以上の変数による微分方程式は偏微分方程式と呼ぶ

独立変数 t, 未知関数 y = y(t) とその導関数 yʹ′, yʹ′, ...

を含む方程式

初期値問題

問

_{: 方程式     , 初期条件    を満たす   を求めよ。}

dy

dt = f (t,y) y(t0) = y0 y = y(t)

!2 !1 0 1 2 !2 !1 0 1 2

f (t,y) = ky

k = 0.5 y

0

= 0.5

t

y

初期条件を満たす解曲線

dy

dt

=

ky

y(t) = y

₀

e

kt

一般的な解法

解析的に解けるのは線形

_{, 変数分離形, 同次形などの特}

殊な場合にのみ限られる、、

変数分離形

同次形

微分方程式を数値的に解き近似解を得る

= 微分方程式の数値解法

dy

dx =

f (x)g(y)

dy

dx =

f

y

x

オイラー法

漸化式を計算することで

微分方程式の近似解を得るアルゴリズム

tの刻み幅Δt = ti+1 - ti = h が十分小さいとき、

dx

dt =

f (t, x) = lim

t⇥0

x(t + t) x(t)

t

f (t

_i

,

x

_i

)

⇥

x

i+1

x

i

t

_i+1

t

_i

x

_i+1

_{⇥ x}

_i

_{+ h f (t}

_i

,

x

_i

)

初期値 x₀ を与えると時刻t_iにおけるx₁, x₂, x₃, ... が決定する 微分の定義から x(t) = xi, x(t + t) = xi+1

オイラー法のアルゴリズム

while (t <= t

n

){

t

i+1

= t

i

+ h,

x

i+1

= x

i

+ h f(t

i

, x

i

)

}

_t

₀

_t

₁

_t

₂

x

0

x

1

x

2

f(t

0

,x

0

)

f(t

1

,x

1

)

x

t

h

区間 [t

0

, t

n

] を等間隔 h で分割し、ステップ t

i

から

t

i+1

の

x の増加率を h f(t

i

, x

i

)で与える

解析解

数値解

dx

dt

t

0

t

1

t

2

x

t

x

0

x

1

x

2

t

0

t

1

t

2

x

0

x

1

x

2

f(t

0

,x

0

)

f(t

1

,x

1

)

x

t

誤差

真の値との誤差

h

オイラー法の精度

x(t)をt = ti の周りでテイラー展開すると

x(t

_i+1

) = x(t

_i

+ h) = x(t

_i

) + hx (t

_i

) +

h

2

2

x (t

i

) + ···

x(t

_i+1

) = x(t

_i

) + h f (t

_i

,

x

_i

)

オイラー法はテイラー展開の1次の項までしか考慮しない オイラー法

_誤差

区間[ t0, t0 + T ]を刻み幅hで分割した 時の終端 t0 + T での誤差の累積

h

2

2

x (t

i

)

T

h

=

h

2

x (t

i

)

刻み幅 h

を 1/10 にすれば累積誤差も約 1/10 になる

(ただし計算量は 10倍になる)

1回のステップ における誤差 テイラー展開による 多項式での近似 全ステップ数

プログラム例 (オイラー法, 1変数)

#define  DT  0.01              /*  時間刻み幅  */  

#define  STEP_MAX  1000  /*  実時間  DT*STEP_MAX  =  10.0  まで  */  

!

double  fn(double,  double);  /*  導関数  */  

void  euler(double,  double,  double*,  double);  /*  オイラー法  */  

! main(){      long  step;      double  t,  x,  x_next;   !    x=0.1;      /*  初期値設定  */   !

   for(i=0;  i<STEP_MAX;  i++){          t  =  i*DT;          euler(t,  x,  &x_next,  DT);          x  =  x_next;      }   }   !

void  euler(double  t,  double  x,  double  *x_out,  double  h){   ...  

}  

!

double  fn(double  t,  double  x){   ...   }

dx

dt

= f (t,x)

t = 0 より t = DT*STEP_MAX

までを計算

プログラム例 (オイラー法, 1変数)

void  euler(double  t,  double  x,  double  *x_out,  double  h){  

! double  tmp;   tmp  =  x  +  h*fn(t,x);   *x_out  =  tmp;   ! }   !

double  fn(double  t,  double  x){   ...   }

tmpを作成し, t

i

時のx(t

i

)からx(t

i+1

)

を計算した結果を一時的に保存,

x_outを介して返す

関数euler

t

i

, x(t

i

), h の引数を与えると t

i+1

= t

i

+ h 時の x(t

i+1

) を

引数を介して関数呼び出し元に返す

ti

x(ti

)

_アドレス

x(ti+1

)

の

h

Mathematicaによる数値解の視覚化 (1/2)

オイラー法により書き出したMathematica により読

み込み, 視覚化する

0.000 1.00   0.001 1.02   0.002 1.05   ...

t

i

x(t

i

)

ti と x(ti) の間に スペースを空ける $ mathematica SetDirectory[“~/keisanki-1-2010/”] 1) オイラー法により計算 した結果を右図の形式で 書き出す 2) GNOME端末からMathematicaを起動する 3) Mathematicaのシートでデータファイルのあるディレクトリ を【SetDirectory】コマンドで指定する

data = ReadList[ “data_file”, {Real, Real}] 【ReadList[“file”, {type, type}]】

ファイル名”file”からtypeで指定した型のデータを用いてファイルの!

終わりに到達するまで読み込む (ここでは {x 座標(ti), y 座標(x(ti))}の対)。

ListPlot[ data, Joined->True, PlotRange->All ]

Mathematicaによる数値解の視覚化 (2/2)

4) データファイルを読み込む 読み込んだもの(リスト)に data という名前を付ける 5) データをグラフに描く 【ListPlot[list]】リスト形式になった値をxとy座標にプロットする データ各点を結ぶオプション グラフの範囲(縦軸)はすべて の範囲とするオプション PlotRange->{0, 10} とすると、! 縦軸の範囲が0から10となる

Mathematicaを用いてCプログラミングで求めた数値解と解析解を 比較してみよう。

fig1=ListPlot[ data, Joined->True, PlotRange->All ]

! ! ! !

fig2=Plot[Exp[t], {t, 0, 10}, PlotStyle ->Red]

! ! ! fig3=Show[fig1, fig2] 1. 数値解のプロットをfig1に代入 2. 解析解(ここでは et) を 区間 0 ≤ t ≤ 10  で赤線でプロットし、fig2に代入 3. fig2とfig2を重ねてプロットするし、  fig3に代入

Mathematicaによる解析解との比較

1. 2. 3.

0 2 4 6 8 10t 0.5 1.0 1.5 2.0 2.5x 0 2 4 6 8 10t 0.5 1.0 1.5 2.0 2.5x

fig3出力の例

dx

dt =

x(1

x)

x(0) = 2

刻み幅 h が 0.5 の時 刻み幅 h が 0.05 の時

Mathematicaを用いて常微分方程式を数値的に解くことができる。 Cプログラミングで求めた結果と比較してみよう。 deq = { x’[t] == (1 - x[t])x[t], x[0] == 0.01} ! ! ! ! !

sol = NDSolve[ deq, x[t], {t, 0, 10}]

! ! !

Plot[ Evaluate[ x[t]/.sol ] , {t, 0, 10}]

! ! ! ! x[t]/.sol/.t->5 1. 右の初期値問題 を deq として定義 ★ Mathematica では代入が=, 等式が == で示される 2. 初期値問題 deq を 区間 0 ≤ t ≤ 10 で 数値的に解き、結果を sol と名付ける 3. 数値的に解いた結果を t の関数として描く 4. x(5) を数値的に求める ★ /. は右辺のルールを左辺に代入することを表す

Mathematicaによる数値解法

1. 2. 3. 4.

1. 初期値問題(x(0) = 1)を刻み幅 h = 0.1 及び h = 0.05 のオイラー法で解き、! 解析解(自分で求める)と比較し、下記の形式で出力し、 刻み幅 h を半分に!  すると解の精度はどの程度改善するか述べよ。! ! ! 2. 解析解と数値解をMathematicaを用いて同一グラフで視覚化せよ。 解析解 誤差 数値解

レポート課題1

下記(1)式、(2)式の常微分方程式について 1 - 3. の課題に取り組め (1) (2) 解析解 誤差 数値解 dx dt = x sin (t) 0 ⇥ t ⇥ 50 dx dt = x 0 t 1

昨年質問が多かったところ

Q1 データ書き出し行が多すぎてリダイレクト

によるファイル出力がうまくいきません。

A1 fprintf関数を使ってデータのファイル出力

をしましょう。

FILE  *fp  

fp  =  fopen(“data.txt”,  “w”);  

/*書き出し用ファイルを開く*/

 

! !

fprintf(fp,  “%f  %f  ¥n”,  t,  x);  

! !

fclose(fp);

 /*書き出し用ファイルを閉じる*/

…

…

(例)

_{char型で与えること}

書き込みモードで開く

Q2 Cプログラミングでsin(t)が使えません。!

A2 以下2点を確認してください。

•

#include  <math.h>  が宣言されているか  

•

コンパイル時に  -­‐lm  をつけているか

Q3 Mathematicaでsin(t)が使えません。!

A3 Sin[t]と入力してください。Sは大文字です。

1 階の連立常微分方程式

ベクトル表示の場合

オイラー法

プログラム例 (オイラー法, 多変数)

#define  DIM  2                    /*  次元  (変数の数)  */   #define  DT    0.01              /*  時間刻み幅  */  

#define  STEP_MAX  1000    /*  実時間  DT*STEP_MAX  =  10.0  まで  */  

!

void  derives(double,  double[],  double[],  int);  /*  導関数  */  

void  euler(double,  double[],  double[],  int,  double);  /*  オイラー法  */  

! main(){      long  step;      double  t,  x,  x_next;   !    x=0.1;      /*  初期値設定  */   !

   for(i=0;  i<STEP_MAX;  i++){          t  =  i*DT;          euler(t,  x,  &x_next,  DT);          x  =  x_next;      }   }   !

void  euler(double  t,  double  x,  double  *x_out,  double  h){   ...  

}  

!

double  fn(double  t,  double  x){   ...   }

t = 0 より t = DT*STEP_MAX

までを計算

d

_x

dt =

f (t, x)

プログラム例 (オイラー法, 多変数)

void  derives(double  t,  double  x[],  double  derivatives[],int  n){  

! ! ! ! ! ! ! }

関数derives

t

i

, x(t

i

), n の引数を与えると導関数を計算し、

derivatives[] 引数を介して計算結果を返す

ti

_x(t

_i

₎

dx/dt

or

f(t, x)

n

の導関数を計算し, derivatives[] に格納

プログラム例 (オイラー法, 多変数)

void  euler(double  t,  double  x[],  double  x_out[],  int  n,  double  h){   int i;!

double  x_tmp[DIM],  derivatives[DIM];  

!

derivs(t,  x,  derivatives,  n); /* 導関数の値を求める */!

!

for(i=0;  i<DIM;  i++) /* Euler 法 */! x_tmp[i]  =  x[i]  +  h*derivatives[i];!

!

for(i=0;  i<DIM;  i++)  

x_out[i]  =  x_tmp[i]; /* 引数を介して演算結果を返却 */! }

関数euler

t

i

, x(t

i

), n, h の引数を与えると t

i+1

= t

i

+ h 時の x(t

i+1

)

を引数を介して関数呼び出し元に返す

ti

_x(t

_i

₎

x(t

i+1

)

n

h

もう少しヒント

void  derives(double  t,  double  x[],  double  derivatives[],int  n){  

! ! ! ! ! ! }

の導関数を計算し,

一時的に

_{derivatives[] に格納}

double  fn(double  t,  double  x){  

! ! ! }

f (t, x)

の導関数を計算して返す

1変数のfn関数と多変数のderives関数の役割は同じ

double  fn1(double  t,  double  x1,  double  x2){   ! ! ! }   !

double  fn2(double  t,  double  x1,  double  x2){  

! ! ! }

f

1

(t, x

1

, x

2

)の導関数を計算して返す

...という説明でも???の場合は次のように2変数の

導関数を計算する場合を考えてみよう

f

2

(t, x

1

, x

2

)の導関数を計算して返す

もう少しヒント

1. 横軸をt, 縦軸をx1, x2とした図をMathematicaで作成せよ! ! 2. 横軸をx1, 縦軸x2とした相平面上の解軌道の図をMathematicaで作成せよ

レポート課題2

下記の連立常微分方程式のパラメータεを1.0, 0.5, 0.1と変化させたそ れぞれの場合についてオイラー法で数値的に解き、1-2.に取り組め

dx

₁

dt =

x

1

1

3

x

13

x

2

dx

₂

dt =

x

₁

ファン・デル・ポール振動子

Van der Pol oscillator equation

Mathematicaによる数値解の視覚化 (1/4)

オイラー法により書き出したMathematica により読

み込み, 視覚化する

_t

i

x

1

(t

i

)

スペースを空ける % /usr/local/bin/mathematica & SetDirectory[“~/keisanki-1-2010/”] 1) オイラー法により計算 した結果を右図の形式で 書き出す 2) ターミナルからMathematicaを起動する 3) Mathematicaのシートでデータファイルのあるディレクトリ を【SetDirectory】コマンドで指定する 0.000 1.00 2.00     0.001 1.02 1.97   0.002 1.05 1.94   ...

x

2

(t

i

)

data = ReadList[ “data_file”, {Real, Real, Real}]; データファイルより{ti, x1(ti), x2(ti)}の対を読み込み、読み込んだリストに! data という名前を付ける。

Mathematicaによる数値解の視覚化 (2/4)

4) データファイルを読み込む 5) グラフ描画用にデータを成形する - リストの転置 data は {{t₀, x₀, y₀}, {t₁, x₁, y₁}, {t₂, x₂, y₂}, ... } という形式になっているので! Transpose(転置)コマンドによりdataの成分を転置したものをdataTとする dataT = Transpose[data]; これにより dataT は {{t₀, t₁, t₂, ... }, {x₀, x₁, x₂, ... }, {y₀, y₁, y₂, ...}} となる

Mathematicaによる数値解の視覚化 (3/4)

4) グラフ描画用にデータを成形する - 描画用リストの作成

dataX = Transpose[ {dataT[[1]], dataT[[2]] } ]; dataT の第 1 成分と第 2 成分を組み合わせて転地し、 新たに( t, x1(t) ) のデータ dataX とする

dataY = Transpose[ {dataT[[1]], dataT[[3]] } ];

同様にして ( t, y(t) ) のデータ dataY 、( x(t), y(t) ) のデータ dataXY をつくる

Mathematicaによる数値解の視覚化 (4/4)

5) データをグラフに描く

時系列データをグラフに描くコマンド ListPlot を用いて、横軸を時 刻 t 、縦軸をそれぞれ x(t), y(t) としたグラフを描き、名前を付ける

gx = ListPlot[ dataX, Joined->True, PlotRange->All] gy = ListPlot[ dataY, Joined->True, PlotRange->All] 2 つのグラフを Show コマンドで重ね合わせて表示

Show[gx, gy]

横軸を x(t), 縦軸を y(t) としたグラフ（相平面上の解の軌道）を描く ListPlot[ dataXY, Joined->True, PlotRange->All]

deq2 = { x1’[t] == –x2[t] – x1[t]^2, x2’[t] == 2 x1[t] – x2[t], x1[0] == x2[0] == 1} ! ! !

sol2 = NDSolve[ deq2, {x1[t], x2[t]}, {t, 0, 20}]

! ! !

Plot[ Evaluate[ x1[t]/.sol2 ], {t, 0, 20}] Plot[ Evaluate[ x2[t]/.sol2 ], {t, 0, 20}]

! ! !

ParametricPlot[ Evaluate[ {x1[t], x2[t]}/.sol2 ], {t, 0, 20}] 1. 右の初期値問題 を deq として定義 2. 初期値問題 deq2 を 区間 0 ≤ t ≤ 20 で 数値的に解き、結果を sol2と名付ける 3. 数値的に解いた結果をtの関数として描く 4. 相平面上の解軌道 ( x₁(t), x₂(t) ) を描く

Mathematicaによる数値解法

1. 2. 3. 4.

グラフを LaTeX に貼り込む方法

EPS: Encapsulated PostScript

LaTeX ドキュメントのプリアンブルに右文を追加 \usepackage{graphicx}

EPS グラフ “ファイル名.eps” を取り込む \includegraphics{graph.eps}

\begin{cetner}

\includegraphics[width = 5cm]{graph.eps} \end{center}

(例) 図を中央そろえ、横幅 5cm に縮小拡大して表示する場合 保存したいグラフを選択して、Mathematica のメニューから、!

Edit ---> Save Selection As ---> EPS

1) Mathematicaで作成した図をEPS形式で保存

レポート作成のTips

•

図や表に番号を振ろう

•

図や表にキャプション(説明)をつけよう

ポイント

その図を初めて見る人を想定して説明する

図1 dx/dt = x(1-x)の解析解(赤線)並びにオ イラー法時間刻み幅 h = 0.5 で計算した数 値解(黒線)。初期値はx(0) = 2。 0 2 4 6 8 10t 0.5 1.0 1.5 2.0 2.5x

•

理科系の作文技術

_{/ 木下是雄 (2002) 中央公論社 (中央新書)}

•

レポートの組み立て方

/ 木下是雄 (1994) ちくま学芸文庫

•

これからレポート・卒論を書く若者のために

_{/ 酒井聡樹}

(2007) 共立出版

今後のレポート課題、卒業論文作成にも役立つ

ので図書館等で探して読んでみよう

レポート作成の指南書

等々

これは私の学生の時の必携の書

参考文献: !

1) Burghes, D. and Borrie, M. (1990)「微分方程式

で数学モデルを作ろう」, 株式会社 日本評論社!

2) 高須 夫悟 (2009) H21年度 計算機実験1 講義レ

ジュメ!

3) 原 進 (2000) はやわかりMathematica 第2版,

株式会社 共立出版