【第7回】

<内容>

1.  ファイルに対する入出力

·       数値をファイルへ保存する方法と,ファイルから数値を取り込む方法について学ぶ.

·       C言語プログラムのファイル入出力方法によく似ているが,ベクトルを扱える点が異なる.

 

2.  常微分方程式の解法

·       常微分方程式の数値解法としてよく知られているオイラー法とルンゲクッタ法について学ぶ.

·       常微分方程式のソルバ(solver)であるlsodeを用いて解く方法について学ぶ.

 

 

<演習>

1.ファイルに対する入出力

 

11.ファイルへの出力

ベクトルxexp(x)からなる行列yをファイルexp.datへ出力する方法(教科書p.55

 

<解説>

(1)      01まで,0.1刻みの値を持つベクトルxを生成する.

(2)      ベクトルxと,そのxの値により計算されたexp(x)2つのベクトルを成分として持つ行列yを生成する.

(3)      ファイル名exp.datを書き込み用ファイル(オプション:’w’)として開き,これをfidという名前で認識させる.ファイルオープン命令がfopenである.

(4)      fidへ書式「%6.2f %12.4f」でyの値を出力する.\nUnix系)は改行を意味し,Windows系では\nと入力する.ファイル出力命令がfprintfである.

(5)      最後に,開いていたファイルfid,すなわち,書き込み用ファイルexp.datをファイルクローズ命令fcloseにより閉じて終了する.

 

 

練習7-1

 上記の操作(教科書p.55)を実行し,作成されたファイルexp.datの内容を確認せよ.

 

注)Octaveウィンドウ左のファイルブラウザにて該当ファイルを選択後,右クリックで表示されるサブウィンドウの「開く」を左クリックすれば,エディタに内容が表示される.

 

 

12.ファイルからの入力

 ファイルexp.datから値を呼び出して,変数aに入力する方法(教科書p.55

 

<解説>

(1)      ファイル名exp.datを読み出し用ファイル(オプション:’r’ は省略可能)として開き,これをfidという名前で認識させる.

(2)      fidから「%f %f」形式で読み込んだ値を,変数a[2, inf]形式で代入する.すなわち,変数aは行列となる.[2, inf]2inf(無限)列を意味し,変数a2行としてデータが無くなるまで読み込めということである.教科書では[2 inf]となっているが,通常は[2, inf]と書く.

(3)      行列aを転置することにより,ファイル内容と同じ形式とする.ファイルからの入力では,列単位で値を代入していくので,ファイル内容を転置した形式でしか読み込めない点に注意せよ.

(4)      最後に,開いていたファイルfid(すなわち,exp.dat)をfcloseで閉じる.

 

 

練習7-2

 上記の操作(教科書p.55)を実行し,作成された行列aの中身を確認してみよ(コマンドウィンドで「a」+「Enter」で表示できる).さらに,2行目をa = fscanf(fid, ’%f %f’, [11, 2]) とし,3行目のa=a' を削除した場合,aの中身がどのようになっているかを確認せよ.

 

<解説>

 exp.datのファイルを開くと112列に並んでいるように見えるが,ファイルのトラック上には以下のような流れ(ストリーム)で記憶されている.\nは改行を意味する特殊文字である.必ず,上流から下流へ向かってデータの処理が行われる.

 

 

―――――――――――――――――――――――――――――――――――――――――――――――――――――――

 

 

2.常微分方程式の解法

 常微分方程式の解法としてよく知られている,オイラー法ルンゲクッタ法のプログラム作成法を紹介する.また,Octaveで用意されているソルバ,lsodeの使用方法も説明する.

 

 

21.オイラー法

 連立一階常微分方程式

  

に対する解x(t)を,小さな正数Dtによりと近似し

  図面1.jpg

差分近似することにより求める方法である.オイラー法の近似精度は低く,Dtを十分小さくしないと精度良く解が得られないが,容易にプログラムできるという利点を持つ.

 

 

練習7-3

 例題2.19(教科書p.56を行え.さらに,入力をu(t) = tとしたときの応答も求めてみよ.

 

注)xx(:,i)=xは,ベクトルx21列)を行列xxに代入する際,列のみを指定して代入する操作を意味する.xxの行はxの行に対応づけられる.つまり,

  となる場合,

  というようになる.

 

 

22.ルンゲクッタ法

 4次ルンゲクッタ法はオイラー法より精度良く解を得ることができる.すなわち,時間刻みDtをオイラー法より大きくとることができ,プログラムの実行時間の短縮化を図ることができる.

 連立一階常微分方程式

  

に対する解x(t)

  図面1.jpg

により求める方法である.下図はルンゲクッタ法の概念図である.

runge-kutta

 

 

練習7-4

 例題2.20(教科書p.58を行え.さらに,入力をu(t) = tとしたときの応答も求めてみよ.

 

注)

¨    globalで宣言された変数は共通変数と呼ばれ,関数の引数として渡さなくても,関数内で使用することができる変数である.

¨    runge.mの中身は,教科書p.59function x = runge(fcn,xt,dt) 以下である.それまでの5行は関数rungeの説明である.

 

 

 2次系の状態変数x1x2をそれぞれx軸,y軸にプロットしてできた平面図を位相面図と呼び,2次系の挙動を調べるために用いられる.

 

 

練習7-5

 例題2.21(教科書p.60)を行え.

 

注)

¨     axisは座標軸の大きさを指定するものであり,これもplotの後に書く必要がある.

¨     メインプログラムの最後の行は,グラフ表示範囲の変更をデフォルトである’auto’(座標軸を描画するものに自動で合わせて表示するモード)に戻して終了するために書いてある.

 

 

23lsodeによる解法

 常微分方程式はソルバlsodeを用いることにより容易に解くことができる.ただし,lsodeOctaveのみで使用できるソルバでありMatlabでは使用できない.その代わり,Matlabでは,ルンゲクッタ法に基づくソルバ,ode23ode45などが用意されている.

 

 

練習7-6

 lsodeを用いて例題2.19(教科書p.56)を解いてみよ(教科書p.62参照).ただし,微分方程式を定義する関数f.m(教科書p.59)を用意しておく必要がある.

 

 

問題7-1

 問題1.3.(教科書p.62)を行え.

 

 

【レポート課題(7)(授業終了後にアップロード)

  RePORT