伊達 >> 生体情報システム >> 課題一覧 >> 課題2.5

課題2.5

目的

説明を読むよりは「例題」を実践し「学習」してみてください．

例

1. octave により微分方程式の解を数値実験により求める方法

   dx/dt = -3x， x(0)=1，という微分方程式の解を求める．
   # から始まる行や，各行の # から後ろはコメントです．

   lsode-sample.m

   function dx = f(x,t)
           dx(1) = -3.0*x
   end
   
   t=linspace(0,10,101);    # 横軸の作成． 0 から 10 を 101等分
   x=lsode("f",1,t);        # 1 は x の初期値 x(0)． lsode は関数．こう書けば解いてくれる．
   plot (t,x,";x(t);");
   pause
   

   % octave lsode-sample.m 
   

2. Fitzhugh- Nagumo モデルを記述

   nagumo.m

   #
   global theta gamma epsilon INPUT;
   theta= 0.01;
   gamma = 0.5;
   epsilon = 0.002;
   INPUT = 0.4;
   
   function dx = nagumo(x,t)
           global theta gamma epsilon INPUT;
           dx(1) = -x(1)*( x(1)- theta)*( x(1)-1.0 ) - x(2) + INPUT;
           dx(2) = epsilon*( x(1)  - gamma*x(2) );
   end
   
   t=linspace(0,2000,2000);  
   x=lsode("nagumo",[0.5; 0.2],t); # 0.5; 0.2 は x(1) x(2) の初期値
   __gnuplot_set__ key bottom
   plot( t,x(:,1), ";u;");
   hold on
   plot(t,x(:,2),";v;" );
   pause
   
   clearplot
   plot(x(:,1),x(:,2));   # もうひとつの図をみる．横軸 u 縦軸 v．
   pause
   
   # for creating PS figure
   #
   __gnuplot_set__ term postscript eps color "Helvetica" 20
   __gnuplot_set__ output "hoge.eps"
   replot
   

   % octave nagumo.m 
   

   
3. メモ
   プログラム上では u を x(1)， v を x(2) と記述している． v は u に比べ，非常にゆっくりと変化している場合を考える． epsilon の値が小さいのはそれを反映している． つまり，u が動いている間， v はある特定の値に固定されている と考えていい．