今後の予定 6/29 パターン形成第 11 回 7/6 データ解析第 12 回 7/13 群れ行動 ( 久保先生 ) 第 13 回 7/17 ( 金 ) 休講 7/20 まとめ第 14 回 7/27 休講?

(1)

今後の予定

• _{6/29
パターン形成第11回}

• _{7/6
データ解析第12回}

• _{7/13
群れ行動（久保先生）第13回}

• 7/17 (金）休講

• _{7/20
まとめ第14回}

• _{7/27
休講？}

(2)

数理生物学演習

(3)

本日の目標

•  ２次元配列

•  分子の拡散

•  反応拡散モデル

•  チューリングパタン

(4)

拡散方程式

∂u

∂t

= D∇

2 u

拡散方程式

拡散が生じる分子などの挙動を記述する．

集団遺伝学で出てくることも．

復習

∇ = ∂

∂x

₁

,

∂

∂x

₂

,

,

∂

∂x

_n

⎛

⎝⎜

⎞

⎠⎟

空間微分演算子　∇(ナブラ）

スカラー量（例えば拡散性分子の濃度）の勾配

grad u

= ∇u = ∂

u

∂x

₁

e

1 + ∂

u

∂x

₂

e

2 ++ ∂

u

∂x

_n

e

n

e

_i

(単位ベクトル）

(5)

反応拡散モデル

∂u

∂t

= D

u

∇

2 u

+

u

2 v

− u

∂v

∂t

= D

v

∇

2 v

+ u

2 − v

ギーラー−マインハルト系

活性化因子

アクチベータ

抑制因子

インヒビター

v

u

抑制因子

インヒビター

活性化因子

アクチベータ

仮定

・　

u

と

v

は共に拡散する

・　

u

は自己活性化する

・　

u

は

v

の合成を促進する

・　

v

は

u

の合成を抑制する

1. u が自己活性化により増える

2.同じ場所で v が u により活性化される

3. v が周囲に拡散して u の合成を抑えこむ

4. v の少ない場所までu がわずかに拡散する

(6)

チューリングパタン

ほぼ一様な構造のない状態から，自発的に空間パタンがうまれる

「M/Y/D/S 動物のイラスト集」より

近藤滋Webページより

(7)

有限差分法による離散化

前進差分による近似

後退差分による近似

中心差分による近似

2階の中心差分による2階偏微分の近似

差分商により微分を近似することで，微分方程式を離散化する方法

ある関数

u(x, y)を2次元空間上で離散化する．

u

_{i, j}

=u(x

_i

, y

_j

), (x

_i

, y

_j

)でのuのx方向への偏微分を　　，x方向への刻み幅をΔxとすれば，

差分を刻み幅で割ったもの

∂u

∂x

"

#

$

%

&

'

i

=

u

i+1, j

− u

i, j

Δx

∂u

∂x

"

#

$

%

&

'

i

∂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 u

∂x

2 "

#

$

%

&

'

i

=

u

i+1, j

− 2u

i, j

+ u

i−1, j

Δx

2 y方向へも同様に考えれば良し！

時間方向への離散化は

これまで通りオイラー法を使う

オイラー法と同じ

(8)

拡散方程式の離散化

_{有限差分法による離散化}

時間方向の離散化

空間方向の離散化

前進差分により近似

いつものオイラー法

2階の中心差分により近似

u

_{i, j}

u

_{i+1, j}

u

_{i-1, j}

u

_{i, j-1}

u

_{i, j+1}

u

_{i, j}

u

_{i+1, j}

u

_{i-1, j}

u

_{i, j-1}

u

_{i, j+1}

x

y

x

t

u

_{i, j,t}

_+Δt

= u

_{i, j,t}

+ Δ

t

h

2 ⎡⎣

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

)

⎤⎦

∂u

∂t

= D∇

2 u

離

散

化

h:空間の刻み幅

(9)

(10)

配列

型　配列名[配列サイズ];

型　配列1[サイズ], 配列2[サイズ],…, 配列n[サイズ];

たくさんの変数を個別に宣言するのは面倒！

同じ型をもつ変数のリスト

11-‐1. 配列

#include <stdio.h>

int main(void){

int i;

int a[10];

for(i=0;i<10;i++){

a[i]=i;

}

for(i=0;i<10;i++){

printf("%d\n",a[i]);

}

return 0;

}

i番目の要素にiを代入

•  配列のなかのそれぞれの変数を

配列要素という

a[0], a[1], …, a[9]

特に注意！！

•  添字は

0から始まり

，

(サイズ-1)で終わる

int a[10];　で定義したならば，

a[0]〜a[9]までの要素が存在する

•  各要素は添字によりアクセスする

配列要素

配列

(11)

２次元配列

型　配列名[配列サイズ][配列サイズ];

型　配列1[サイズ][サイズ], 配列2[サイズ][サイズ],

…, 配列n[サイズ][サイズ];

2次元配列

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の

配列要素

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");

}

return 0;

}

基本は1次元の場合と同じ．

より多次元の配列も定義できます！

(12)

拡散方程式の離散化

_{有限差分法による離散化}

時間方向の離散化

空間方向の離散化

前進差分により近似

いつものオイラー法

2階の中心差分により近似

u

_{i, j}

u

_{i+1, j}

u

_{i-1, j}

u

_{i, j-1}

u

_{i, j+1}

u

_{i, j}

u

_{i+1, j}

u

_{i-1, j}

u

_{i, j-1}

u

_{i, j+1}

x

y

x

t

u

_{i, j,t}

_+Δt

= u

_{i, j,t}

+ Δ

t

h

2 ⎡⎣

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

)

⎤⎦

∂u

∂t

= D∇

2 u

離

散

化

h:空間の刻み幅

(13)

周期境界条件

“端”同士が張り合わされていると考える．

プログラムを組むときも，

この部分の処理は注意！

(14)

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; }

(15)

エクセル

0．データを開いた後、

全選択して

”

書式

”→”

セ

ル

→

数値

”

から小数点を

下２桁

まで表示するよ

うにする。

表の左上の角を押して

全選択をする

書式

_{→セルをクリック}

(16)

エクセル

数値を選んで

(17)

エクセル

１ ．

データを貼付けた後

、

列の幅

を調整してセ

ルを正方形にする

表の左上の角を押して

全選択をする

列

_{AとBの間を}

ダブルクリック

(18)

ズームアウト

して全体像を見る

•  示する

(19)

条件付き書式

で値が「1」のセルだけ強調表示する

(20)

(21)

反応拡散モデル

ギーラー−マインハルト系についてプログラムを組み，

さまざまな模様を描く．

12-‐4.

やるべきこと

・ギーラー−マインハルト系の離散化

・パラメータを調整して，どんなパタンがあらわれるかを調べる

∂u

∂t

= D

u

∇

2

u +

u

2

v

− u

∂v

∂t

= D

v

∇

2

v + u

2

− v

アクチベータ

インヒビター

離散化

?

・２次元拡散方程式のプログラムを参考にプログラムを組む

基本的には拡散方程式を反応拡散方程式系に変更するだけ．

今後の予定 6/29 パターン形成第 11 回 7/6 データ解析第 12 回 7/13 群れ行動 ( 久保先生 ) 第 13 回 7/17 ( 金 ) 休講 7/20 まとめ第 14 回 7/27 休講?