2012-04-13

歴史は「べき乗則」で動く　第５章

５章は砂山ゲームと地震のモデル化ゲームから地震の臨界状態への自己組織化を見ていきます。

まずは地震のモデル化について。
バリッジ＝ノポフの地震モデルPeriodicity, Chaos and Localization in a Burridge-Knopoff Model of an Earthquake with Dieterich-Ruina Frictionがべき乗則をうまく再現したモデルであった。
このモデルではブロックがそれぞれバネでつながっており、ブロックが周りのブロックに常に相互に影響している。
このゲームは実は砂山ゲームと数学的には共通した構造をもっている。
（とはいうものの、はじめ砂山ゲームの構造がいまいちよくわからなかった。
山田先生のブログにあったこちらの論文KAKEN — Research Projects | Physical elucidation of the behaviors of sand-pile avalanches like self-organized criticality and like characteristic earthquakes (KAKENHI-PROJECT-18340134)を書いた人の研究http://www.fgi.or.jp/research/research_yoshioka.pdfを見てわかったことだが、砂山内部に応力鎖が形成され、ちょうど地震モデルのようになるらしいです。）
しかしながら、このモデルはエネルギーが保存されており、地震でエネルギーが消費されていないといけない。

そこで、修正が入ったモデルがオラミ＝フェダー＝クリステンセンのモデルResearch: The Olami-Feder-Christensen Model of Earthquakesである。
このモデルの面白いところは、エネルギーロスの割合を変化させても常にべき乗則が成立することである。
加えて、グーテンベルク＝リヒターの法則が成立していた。
このモデルを用いてシュミレーションすると大森公式（余震）大森公式 - Wikipediaにほぼ一致する結果となった。
（つまりは規模のスケール不変性だけでなく、時間に対してもスケール不変性が成り立つようです。）

このモデルが果たして地震を正確にモデル化しているかだが、かなり簡略化しているため異論もあるらしい。
そのあたりはおいておくとして、べき乗則を作る臨界状態への自己組織化の一つのモデルなのは確かだ。

次の６章以降では臨界状態へといたる他の現象に移っていくみたいです。

2012-04-08

歴史は「べき乗則」で動く　第４章

歴史は「べき乗則」で動く

べき乗則についてさらに突っ込んでいきました。

べき乗則の性質スケール不変性は自然界・人間社会で多くみられる。
マンデルブローによる綿花価格変動などの発見。ブノワ・マンデルブロ - Wikipedia
安定分布安定分布 - Wikipediaになっており、スケール不変になっている。
また、エイリー・ゴールドバーガーによる心拍数の自己相似性の発見。

このような自己相似を数学の概念で「フラクタル」フラクタル - Wikipediaと呼ぶ。
文章中ではコッホ曲線コッホ曲線 - Wikipediaが例に挙げられていた。
（測度で勉強したカントールカントール集合 - Wikipediaもフラクタルだったとは意外なつながりがあった。）

続いて、ではフラクタルを編み出すものはなんなのか?に移っていく。
１つにはカオスがある。
本文中では加速器の電子の振る舞いを図にしたチリコフ=テイラーの標準マップを見て、カオスな領域と安定な領域がフラクタルになっていることを説明していた。
（『マクロ経済分析』岩波書店を読むと、数理経済学では複雑系をカオスから説明しようとしているのだと感じた。）

もう一つは「歴史」を見ていくしかない。
「歴史」の重要性を教えてくれるものに拡散律速凝集http://zeong.s.chiba-u.ac.jp/pdf/dla-.pdfというゲームがあり、「歴史」が違うために同じ結果は存在しない。
それでも、できた結晶にはスケール不変性が存在している。
これは生物学における偶然凍結説偶然凍結説とは何？ Weblio辞書にも見て取れる。

異なる歴史的過程が作られるゲームの中には現実世界にも存在するものがあるかもしれない。
べき乗則は、それを見つけるための一種の道しるべになっている。

５章は砂山ゲームと地震のモデル化ゲームから地震の臨界状態への自己組織化を見ていきます。

2012-03-29

歴史は「べき乗則」で動く　第２＋３章

歴史は「べき乗則」で動く

今回から具体的な事象を扱っていきます。
まずは地震です。

地震には明確な前兆が（少なくとも発見されて）なく、
また巨大地震にも規則的な周期は（発見されて）ない。

そこで、世界中の地震の規模と回数の統計を取ると、
べき乗則冪乗則 - Wikipediaが成立していることがわかった。
べき乗則にはスケール不変性という性質が存在する。
つまり、「典型的な」値が存在せず、どんな2つの値を単純な規則が成立する。

とすると、大地震の周期や前兆を探すことは無意味になってしまう。
（仮に、地震の大きさに比例して大きく観測される前兆があれば別だと思いますが）

では、べき乗則はどのように考えればよいのか?　→5章へ続く
その前に、べき乗則についてもう少し詳しく　→4章へ続く

2012-03-22

第1章　なぜ世界は予期せぬ大激変に見舞われるか

歴史は「べき乗則」で動く

歴史、地震、火災、金融に見られる劇的な連鎖反応を挙げて、このような激変のあいだにある類似性の発見と、それが起こりうるネットワークの自発的構築の過程について考えていく本らしい。

気になった・よく知らない単語
臨界状態
カタストロフィー理論とは (Catastrophe Theory) カタストロフィーりろん： - IT用語辞典バイナリ
 カオス理論 - Wikipedia
複雑系 - Wikipedia
非平衡状態
歴史

カオスと複雑系の違いは歴史(history)にあり、
カオスは歴史に依存せず、初期条件と相互作用の完全な情報が知ることができれば予測可能になるが、
複雑系は非平衡状態のなかで組織化されたネットワーク（＝歴史）がいま（＝臨界状態）をつくっている。

2012-02-27

微分方程式で数学モデルを作ろう　2章

人工腎臓器、いわゆる透析装置のモデルについてです。
前回の薬のモデルは濃度によって排出量が決まりましたが、
今回は流速によって除去率が決まるモデルです。
老廃物のある血液濃度と透析液の濃度をを変数 $u,v$ で表し、 $x$ を透析開始点からの位置とすると、
$u,v$ は $u=u(x),v=v(x)$ で表されると仮定します。

ここでは以下のフィックの法則を用いて、膜を通過する物質の総量を考えます。

フィックの法則：単位時間に単位面積の幕を通過する物質の総量は、
　　　　　　　　　　膜のその位置における濃度差に比例する

位置 $x$ から $x+\delta x$ の微小区間について
単位時間当たり、血液から透析液まで移行する量は比例定数 $k$ を用いて、
$k\{u(x)-v(x)\}\delta x$
と与えられます。
さらに、(位置 $x$ に入ってきた血液の量)=(膜を通して移行した量)+( $x+\delta x$ で出ていく血液の量)が成立していますので、
微小区間内での血液の量を $Q_{b}$ とすると、
$Q_{b}u(x)=k\{u(x)-v(x)\}\delta x+Q_{b}u(x+\delta x)$
と書き換えられます。
$\delta x \rightarrow 0$ とすると、以下の微分方程式が成り立ちます。
$Q_{b}\frac{du}{dx}=-k(u-v)$
逆に、微小区間内での透析液の量を $Q_{d}$ とすると、
$Q_{d}v(x)=k\{u(x)-v(x)\}\delta x+Q_{d}v(x+\delta x)$
から、以下の微分方程式が成り立ちます。
$Q_{d}\frac{dv}{dx}=k(u-v)$
ここで、 $z=u-v , \alpha=\frac{k}{Q_{b}}-\frac{k}{Q_{d}}$ とすると、
$\frac{du}{dx}-\frac{dv}{dx}=-(\frac{k}{Q_{b}}-\frac{k}{Q_{d}})(u-v)$
より、
$\frac{dz}{dx}=-\alpha z$
となります。
よって、 $u,v$ は任意定数 $C_{1},C_{2}$ を用いて
$u=C_{1}+\frac{kC_{2}}{\alpha Q_{b}}e^{-\alpha x}$
$v=C_{1}+\frac{kC_{2}}{\alpha Q_{d}}e^{-\alpha x}$
と表されます。
ここで、初期条件について
$x=0$ で $u=u_{0}$
$x=L$ で $v=0$ 　( $L$ は透析の開始位置)
と定めると、
$u=u_{0}[\frac{(e^{-\alpha L}/Q_{d})-(e^{-\alpha x}/Q_{b})}{(e^{-\alpha L}/Q_{d})-(1/Q_{b})}]$
$v=\frac{u_{0}}{Q_{d}}[\frac{e^{-\alpha L}-e^{-\alpha x}}{(e^{-\alpha L}/Q_{d})-(1/Q_{d})}]$
となります。

透析装置の設計者は浄化値（クリアランス） $Cl$ を考えるそうです。
浄化値とは
$Cl=\frac{\int^{L}_{0}k\{u(x)-v(x)\}dx}{u_{0}}=\frac{Q_{b}}{u_{0}}\{u_{0}-u(L)\}$
で定義される数値のことで、具体的には
$Cl=Q_{b}[\frac{1-e^{-\alpha L}}{1-(Q_{b}/Q_{d})e^{-\alpha L}}]$
と表されます。

このモデルのパラメーターは $Q_{b}/Q_{d}$ と $kL/Q_{b}$ であるそうです。
なぜこの2つなのかと考えると、可変なのが $Q_{d}$ と $L$ であり、
これらを含み、モデルを考えるうえでわかりやすい独立変数がこの二つになるからだと思われます。

for i=1:1:800;
    k=1;
    L=200;
    Qd=400;
    a=k*(50/i-1/Qd);
    Cl=(i/50)*((1-exp(-a*L))/(1-(i/50*Qd)*exp(-a*L)));
    hold on;
    plot(i/50,Cl,'b-');
    axis([0 17 0 20]);
    title('浄化値の理論曲線')
    xlabel('Qb ml/分');
    ylabel('Cl ml/分');
    set(gca,'XTick',2:2:16);
    set(gca,'YTick',2:2:20);
    hold off;
end;
saveas(figure(1),'Dialysis','png');

グラフは本とはかなり異なった形になりました。
数式的には $Q_{b}/Q_{d}$ が小さい時は $Cl$ は $Q_{b}$ の一次近似になるので、
本のグラフとは $Q_{d}$ の値が異なるのかもしれません。

2012-02-15

微分方程式で数学モデルを作ろう　2章

今週のMIKUの内容とscilabの利用を少々。

条件
・濃度を $C(t)$ とする
・ $m \leq C(t)\leq M$ をみたす
・ $C(0)=m$ をみたす
・単位時間当たりの投与量を $Q(t)=($ 一定 $\)\geq 0$ とする
・ $C(t)$ と $Q(t)$ は $\frac{dC}{dt}=-kC+q$ をみたす(ただし、 $q(t)=\frac{Q(t)}{V}$ )
・ $min\{\int_{0}^{T}q(t)dt\}$ をみたす
・濃度が $C(t)=M$ となったときは $C(t)=m$ となるまで $Q(t)=0$ とする
をみたすような $q(t)$ を求めたいとします。

まず、 $km \leq q(t) \leq kM$ のとき、
$C(t)=(m-\frac{q}{k})e^{-kt}+\frac{q}{k}$
となり、
濃度　 $\frac{q}{k}$ 　にそのまま収束します。

次に、　 $kM < q(t)$ 　のとき、
$0 \leq t \leq T_{1}$ 　で、
経路： $C(t)=(m-\frac{q}{k})e^{-kt}+\frac{q}{k}$ 　を通り、
$T_{1} \leq t \leq T_{1}+T_{2}$ 　で、
経路： $C(t)=(M-\frac{q}{k})e^{-k(t-T_{1})}+\frac{q}{k}$ 　を通り、
これを繰り返すことになります。
ちなみに、 $T_{1}=\frac{1}{k}log\{\frac{q-km}{q-kM}\}$ 、 $T_{2}=\frac{1}{k}log\{\frac{M}{m}\}$ です。
このとき、投与量の時間平均は $q(t)\times\frac{T_{1}}{T_{1}+T_{2}}$ 　となり、
具体的には　 $q\times\frac{log\{\frac{q-km}{q-kM}\}}{log\frac{M}{m}+log\{\frac{q-km}{q-kM}\}}$ 　となります。
$q\longrightarrow\infty$ のとき、
$q(t)\times\frac{T_{1}}{T_{1}+T_{2}}\longrightarrow k\times\frac{M-m}{log\frac{M}{m}}$ 　となりますが、
これは瞬間的に濃度 $m$ から $M$ になるまで投与した場合の投与量の時間平均 $\frac{(M-m)}{T_{2}}$ と同じになります。

その間の振る舞いなのですが、
まず、 $x=\frac{m}{M}$ 　 $q(t)=(l+1)kM$ とおきます。（ $x\in(0,1)$ , $l>0$ ）
すると、 $q(t)\times\frac{T_{1}}{T_{1}+T_{2}}=kM\times\frac{log(1+\frac{1-x}{l})^{l+1}}{log\{(1+\frac{1-x}{l})\times\frac{1}{x}\}}$ となります。
$(1+\frac{1-x}{l})^{l+1}$ 　、　 $\{(1+\frac{1-x}{l})\times\frac{1}{x}\}$ とも1より大きくなります。
そこで、 $(1+\frac{1-x}{l})^{l+1}\div\{(1+\frac{1-x}{l})\times\frac{1}{x}\}=x(1+\frac{1-x}{l})^{l}$ をとり、
$l>0$ をみたす $x$ の関数とみると、 $x=0$ で0、 $x=1$ で1となります。
また、微分すると $(1+\frac{1-x}{l})^{l-1}\{1-x+\frac{1-x}{l}\}>0$ となり単調増加なので
$x\in(0,1)$ , $l>0$ のもとで、 $1<(1+\frac{1-x}{l})^{l+1}<\{(1+\frac{1-x}{l})\times\frac{1}{x}\}$
よって、[tex:kM\times\frac{log(1+\frac{1-x}{l})^{l+1}}{log\{(1+\frac{1-x}{l})\times\frac{1}{x}\}}

2012-02-07

微分方程式で数学モデルを作ろう　2章

２章の「薬の吸収」と「水の過熱と冷却」についてmatlabで書いてみました。
まずは「薬の吸収」からです。

function dYdt=medicine_model(t,Y)
k=0.5
dYdt=[-k*Y];

%一定の濃度幅を維持し続ける場合
k=0.5;
y0=100;
T=0.5;
Ys=y0/(1-exp(-k*T));
Ym=Ys*exp(-k*T);
for i=[0:1:29];
    [t,Y]=ode45(@medicine_model,[i*T (i+1)*T],Ys);
hold on;
plot(t,Y,'b');
title('薬剤投与モデル');
xlabel('時間（day）');
ylabel('投与量(mg)');
axis([0,15,0,500]);
set(gca,'XTick',[1:1:15]);
plot([i*T i*T], [Ym Ys], 'r');
%一定の投与量を維持し続ける場合
Yi=y0*sum(exp([0:1:i]*(-k*T)));
[t,Y]=ode45(@medicine_model,[i*T (i+1)*T],Yi);
plot(t,Y,'g');
plot([i*T i*T], [(Yi-y0) Yi], 'r');
hold off;
end;
saveas(figure(1),'medicine_model','png');

次に「水の過熱と冷却」についてです。
液体の入った容器をヒーターを使って温めるケースを想定しています。
液体の比熱、質量をcw,mw,容器のをcc,mc、さらにC(t):時刻tにおける周囲との温度差とすると、
容器が時刻t=0から得た総熱量h(t)について
$h(t)=(cw*mw+cc*mc)*C(t)$ 　　　　　・・・・・・(1)
と表せます。
さらに、j:伝熱係数、a:容器の表面積とすると、
$\frac{dh(t)}{dt}=-jaC(t)$ 　　　　　・・・・・・(2)
と表せます。
(1)の微分形と(2)から
$\frac{dC}{dt}=[\frac{-ja}{cw*mw+cc*mc}]C$ 　　　　　・・・・・・(3)
と導けます。
これから、
$C(t)=C(0)*e^{$-k_{1}*t$}$ $k_{1}=[\frac{-ja}{cw*mw+cc*mc}]$
となります。
本では初期値45から何もせず冷まして、最後に一気にヒーターで加熱した場合と、42まで下がったら45まで加熱する場合を考え、
どちらのほうが熱量が少なくてすむかを計算しています。
結果は一気に加熱した場合のほうが安くすむそうです。

function dCdt=water_model(t,C);
a=100;j=20;
cw=10;mw=900;
cc=5;mc=200;
k1=j*a/(cw*mw+cc*mc);
dCdt=-k1*C;


a=100;j=20;
cw=10;mw=900;
cc=5;mc=200;
k1=j*a/(cw*mw+cc*mc);
[t,C]=ode45(@water_model,[0,9],45);
plot(t,C,'b');
title('水の加熱と冷却');
xlabel('時刻');
ylabel('周囲より高い温度');
axis([0,15,0,50]);
%最後に温める
hold on;
Cm1=45*exp(-9*k1);
x=9:14;
y=12*(x-9)+Cm1;
plot(x,y,'r');
%一定の温度帯を保つ
T1=5*log(15/14);
T2=3/12;
for i=[0:1:20];
    [t,C]=ode45(@water_model,[i*(T1+T2) i*(T1+T2)+T1],45);
    plot(t,C,'g');
    plot([i*(T1+T2)+T1 (i+1)*(T1+T2)],[42 45],'r');
end;
hold off;
saveas(figure(1),'water_model3','png');

ただ、本のモデルだと加熱中の熱の放出を無視しているので、それは具体的にコスト計算してるのにどうなのかと。
加熱中はq:ヒーターの単位時間当たりの加熱量が加わるので、
(2)式が
$\frac{dh(t)}{dt}=-jaC(t)+q$
と変わるので、(3)式は
$\frac{dC}{dt}=[\frac{-ja}{cw*mw+cc*mc}]C+[\frac{q}{cw*mw+cc*mc}]$
となるので、
$C(t)=C(0)*e^{$-k_{1}*t$}+\frac{k_{2}}{k_{1}}$ $k_{2}=[\frac{q}{cw*mw+cc*mc}]$
ということで、それを組み込んでもう一度matlabで書いてみました。

function dCdt=water_model(t,C);
a=100;j=20;
cw=10;mw=900;
cc=5;mc=200;
k1=j*a/(cw*mw+cc*mc);
dCdt=-k1*C;

function dCdt=water_model2(t,C);
a=100;j=20;
cw=10;mw=900;
cc=5;mc=200;
k1=j*a/(cw*mw+cc*mc);
q=120000;
k2=q/(cw*mw+cc*mc);
dCdt=-k1*C+k2;


a=100;j=20;
cw=10;mw=900;
cc=5;mc=200;
k1=j*a/(cw*mw+cc*mc);
[t,C]=ode45(@water_model,[0,9],45);
plot(t,C,'b');
title('水の加熱と冷却');
xlabel('時刻');
ylabel('周囲より高い温度');
axis([0,15,0,50]);
%最後に温める
hold on;
q=120000;
k2=q/(cw*mw+cc*mc);
Cm1=45*exp(-k1*9);
[t,C]=ode45(@water_model2,[9,15],Cm1);
plot(t,C,'r');
%一定の温度帯を保つ
T1=5*log(15/14);
T2=5*log(6/5);
for i=[0:1:10];
    [t,C]=ode45(@water_model,[i*(T1+T2) i*(T1+T2)+T1],45);
    plot(t,C,'g');
   [t,C]=ode45(@water_model2,[i*(T1+T2)+T1 (i+1)*(T1+T2)],42);
    plot(t,C,'r');
end;
hold off;
saveas(figure(1),'water_model2','png');

２つの図を比べると、結構違いがあったようです。
２つの具体的な違いの計算と、腎臓のモデルは後日あげます。