今後の予定 6/29 パターン形成第 11 回 7/6 データ解析第 12 回 7/13 群れ行動 ( 久保先生 ) 第 13 回 7/17 ( 金 ) 休講 7/20 まとめ第 14 回 7/27 休講?
数理生物学演習 第 11 回パターン形成
本日の目標 2 次元配列 分子の拡散 反応拡散モデル チューリングパタン
拡散方程式 拡散方程式 u t = D 2 u 拡散が生じる分子などの挙動を記述する. 集団遺伝学で出てくることも. 復習 空間微分演算子 ( ナブラ ) = x 1, x 2,, x n スカラー量 ( 例えば拡散性分子の濃度 ) の勾配 grad u = u = u x 1 e 1 + u x 2 e 2 + + u x n e n e i ( 単位ベクトル )
反応拡散モデル ギーラー マインハルト系 活性化因子 アクチベータ 抑制因子 インヒビター u t = D u 2 u + u2 v u v t = D v 2 v + u 2 v u 活性化因子 アクチベータ v 抑制因子 インヒビター 仮定 u と v は共に拡散する u は自己活性化する u は v の合成を促進する v は u の合成を抑制する
チューリングパタン ほぼ一様な構造のない状態から, 自発的に空間パタンがうまれる 近藤滋 Web ページより h0p://www.5s.osaka- u.ac.jp/labs/skondo/ M/Y/D/S 動物のイラスト集 より
有限差分法による離散化 時間方向への離散化は これまで通りオイラー法を使う 差分商により微分を近似することで, 微分方程式を離散化する方法 差分を刻み幅で割ったもの ある関数 u(x, y) を2 次元空間上で離散化する. " u% u i, j =u(x i, y j ), (x i, y j ) でのuのx 方向への偏微分を $ ',x 方向への刻み幅をΔxとすれば, # x & i 前進差分による近似 " $ # u x % ' & i = u i+1, j u i, j Δx 後退差分による近似 " $ # u x オイラー法と同じ % ' & i = u i, j u i 1, j Δx 中心差分による近似 " $ # u x % ' & i = u i+1, j u i 1, j 2Δx 2 階の中心差分による 2 階偏微分の近似 " 2 u % $ ' # x 2 & i = u i+1, j 2u i, j + u i 1, j Δx 2 y 方向へも同様に考えれば良し!
拡散方程式の離散化 有限差分法による離散化 時間方向の離散化 前進差分により近似 いつものオイラー法 y u i, j+1 t 空間方向の離散化 2 階の中心差分により近似 u i-1, j u i, j u i+1, j u t = D 2 u u i, j-1 u i, j+1 u i-1, j u i, j u i+1, j x 離散化 u i, j,t+δt = u i, j,t + Δt h 2 h: 空間の刻み幅 u i, j-1 ( ) D u i 1, j,t + u i+1, j,t + u i, j 1,t + u i, j+1,t 4u i, j,t x
情報処理センター講義室へ移動!
配列 同じ型をもつ変数のリスト たくさんの変数を個別に宣言するのは面倒! 型配列名 [ 配列サイズ ]; 型配列 1[ サイズ ], 配列 2[ サイズ ],, 配列 n[ サイズ ]; 11-1. 配列 #include <stdio.h> int main(void){ int i; int a[10]; 配列のなかのそれぞれの変数を 配列要素という a[0], a[1],, a[9] 配列 for(i=0;i<10;i++){ a[i]=i; i 番目の要素にiを代入 配列要素 for(i=0;i<10;i++){ printf("%d\n",a[i]); return 0; 各要素は添字によりアクセスする 特に注意!! 添字は 0 から始まり,( サイズ -1) で終わる int a[10]; で定義したならば, a[0] a[9] までの要素が存在する
2 次元配列 型配列名 [ 配列サイズ ][ 配列サイズ ]; 型配列 1[ サイズ ][ サイズ ], 配列 2[ サイズ ][ サイズ ],, 配列 n[ サイズ ][ サイズ ]; a[0][0], a[0][1],, a[0][n] a[1][0], a[1][1],, a[1][n] a[2][0], a[2][1],, a[2][n] a[m][0], a[m][1],, a[m][n] 添字 2, 1 の 配列要素 2 次元配列 11-2. 2 次元配列 #include <stdio.h> #include <stdlib.h> #include <time.h> int main(void){ int i,j; int a[10][10]; srand(time(null)); for(i=0;i<10;i++){ for(j=0;j<10;j++){ a[i][j]=rand()%10; for(i=0;i<10;i++){ for(j=0;j<10;j++){ printf("%d",a[i][j]); if(j!=9){ printf(", "); printf("\n"); 基本は 1 次元の場合と同じ. より多次元の配列も定義できます! return 0;
拡散方程式の離散化 有限差分法による離散化 時間方向の離散化 前進差分により近似 いつものオイラー法 y u i, j+1 t 空間方向の離散化 2 階の中心差分により近似 u i-1, j u i, j u i+1, j u t = D 2 u u i, j-1 u i, j+1 u i-1, j u i, j u i+1, j x 離散化 u i, j,t+δt = u i, j,t + Δt h 2 h: 空間の刻み幅 u i, j-1 ( ) D u i 1, j,t + u i+1, j,t + u i, j 1,t + u i, j+1,t 4u i, j,t x
周期境界条件 端 同士が張り合わされていると考える. プログラムを組むときも, この部分の処理は注意!
11-3. 2 次元の拡散方程式 #include <stdio.h> int main(void){ int i,j; int t; double dt=0.01; double u[50][50]; double utemp[50][50]; double D=0.2; FILE *fp; fp=fopen( 11-3.txt","w"); // 初期化 for(i=0;i<50;i++){ for(j=0;j<50;j++){ u[i][j]=0; u[24][24]=1; u[24][25]=1; u[25][24]=1; u[25][25]=1; // 初期値出力 for(i=0;i<50;i++){ for(j=0;j<50;j++){ fprintf(fp,"%f",u[i][j]); if(j!=49){ fprintf(fp, "); fprintf(fp,"\n"); for(t=1;t<1000;t++){ // 境界条件 //i=0,j=0 i=0; j=0; utemp[i][j]=u[i][j]+dt*(d*(u[49][j] +u[i+1][j]+u[i][49]+u[i][j+1]- 4*u[i] [j])); //i=0,j=49 i=0; j=49; utemp[i][j]=u[i][j]+dt*(d*(u[49][j] +u[i+1][j]+u[i][j- 1]+u[i][0]- 4*u[i][j])); //i=49,j=0 i=49; j=0; utemp[i][j]=u[i][j]+dt*(d*(u[i- 1][j] +u[0][j]+u[i][49]+u[i][j+1]- 4*u[i][j])); //i=49,j=49 i=49; j=49; utemp[i][j]=u[i][j]+dt*(d*(u[i- 1][j] +u[0][j]+u[i][j- 1]+u[i][0]- 4*u[i][j])); //i=0 i=0; for(j=1;j<49;j++){ utemp[i][j]=u[i][j]+dt*(d*(u[49][j] +u[i+1][j]+u[i][j- 1]+u[i][j+1]- 4*u[i] [j])); //i=49 i=49; for(j=1;j<49;j++){ utemp[i][j]=u[i][j]+dt*(d*(u[i- 1] [j]+u[0][j]+u[i][j- 1]+u[i][j+1]- 4*u[i] [j])); //j=0 j=0; for(i=1;i<49;i++){ utemp[i][j]=u[i][j]+dt*(d*(u[i- 1] [j]+u[i+1][j]+u[i][49]+u[i][j+1]- 4*u[i] [j])); //j=49 j=49; for(i=1;i<49;i++){ utemp[i][j]=u[i][j]+dt*(d*(u[i- 1] [j]+u[i+1][j]+u[i][j- 1]+u[i][0]- 4*u[i] [j])); // その他 for(i=1;i<49;i++){ for(j=1;j<49;j++){ utemp[i][j]=u[i][j]+dt*(d*(u[i- 1] [j]+u[i+1][j]+u[i][j- 1]+u[i][j+1]- 4*u[i] [j])); // 更新 for(i=0;i<50;i++){ for(j=0;j<50;j++){ u[i][j]=utemp[i][j]; // 出力 if(t%100==0){ for(i=0;i<50;i++){ for(j=0;j<50;j++){ fprintf(fp,"%f",u[i][j]); if(j!=49){ fprintf(fp, "); fprintf(fp,"\n"); fclose(fp); return 0;
エクセル 0. データを開いた後 全選択して 書式 セル 数値 から小数点を下 2 桁まで表示するようにする 表の左上の角を押して全選択をする 書式 セルをクリック
エクセル 数値を選んで 小数点 2 桁まで表示
エクセル 1. データを貼付けた後 列の幅を調整してセルを正方形にする 表の左上の角を押して全選択をする 列 A と B の間をダブルクリック
ズームアウトして全体像を見る 示する この数字を 35 に変更
条件付き書式で値が 1 のセルだけ強調表示する あとは OK を押す
こんな感じの拡散が見れます
反応拡散モデル 12-4. ギーラー マインハルト系についてプログラムを組み, さまざまな模様を描く. やるべきこと ギーラー マインハルト系の離散化 アクチベータ インヒビター u t = D u 2 u + u2 v u v t = D v 2 v + u 2 v 離散化 2 次元拡散方程式のプログラムを参考にプログラムを組む 基本的には拡散方程式を反応拡散方程式系に変更するだけ. ただし,2 種類の拡散性分子があることに注意! 初期値の設定 u=1.0, v=1.0 にわずかなノイズ (0.0 0.001 程度 ) を加える. パラメータを調整して, どんなパタンがあらわれるかを調べる 二つの拡散係数を変化させてみよう.?