Last modified: Thu Apr 26 21:47:52 JST 2007

伊達 >> 卒業研究 >> 課題一覧 >> 課題2

課題2

目的

細胞が他の細胞から入力を受け取り, 興奮してパルスをだす様を表現した微分方程式を考え, その数理モデルの挙動を計算機シミュレーションにより確認する. 各パラメータに依存して挙動がどう変化するか観察する. 非線形振動子の基礎. もう少し進めば,リズム現象をも理解できるようになる. octave と C言語で gsl (The GNU Scientific Library)を使った手法を学ぶ.

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

  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 はある特定の値に固定されている と考えていい. 大学院生は,より深く考えよう.

nagumo.m lsode-f1.m