[PR]当サイトはアフィリエイト広告を利用しています。
前回は味噌の発酵具合をロジスティックモデルで表してごちゃごちゃ計算をしました。
そして最後に数値積分するための径数μを微分方程式を手作業で解かずに求めることはできないだろうかということを考えて終わりました。
今回はμを微分方程式をきちんと解かずに求めてみようと思います。
ルンゲクッタ法は以下の式で表せます。
のとき刻み幅をhとして
今回は
です。
ここで考えてみたいのは、P’にμを適当に代入してもルンゲクッタ法でなんらかのtとPの対応を示した数表は得られるということです。
そしてμの値を色々と変化させれば、そのうちあるtに対して正しいPとなるμが求まるはずです。
この試行錯誤をどうやってコンピュータに求めさせるか。
今回使うのは最適化という手法です。
情報処理の試験で線形計画法というのをやりますが、それの仲間です。
詳しい解説は専門書をお読みください。
色々と出版されていますが、今回扱うExcelのソルバーのエボリューショナリーの参考書として以下を挙げておきます。
大谷紀子（2018）『進化計算アルゴリズム入門』オーム社.
今回の問題に対しては、
ルンゲクッタ法で求めた数値積分のt=270のときの数値解をP(270)として
目的関数は
ABS(P(270)-10)
制約条件は
0≦μ≦1
とします。
ABS()は()で囲まれた数値の絶対値を表します。
最適化のアルゴリズムは色々ありますが、いちいち作ると大変なのでエクセルのソルバーを使います。
手順
以下手順です。
まず以下のような表を作ります。
設定値
A2:0
B2:0.1
C2:=$H$1*(1-B2/11)*B2
D2: =$H$1*(1-(B2+C2/2)/11)*(B2+C2/2)
E3: =$H$1*(1-(B2+D2/2)/11)*(B2+D2/2)
F4: =$H$1*(1-(B2+E2)/11)*(B2+E2)
H1:0.1
A3: =A2+1
B3: =B2+(C2+2*D2+2*E2+F2)/6
次にA3とB3をドラッグしてからフィルハンドルをドラッグして272行目まで数式をコピーします。
A272: =A271+1
B272: =B271+(C271+2*D271+2*E271+F271)/6
この式が入っているか確認します。
次にC2からF2までをドラッグして、フィルハンドルをドラッグして272行目まで数式をコピーします。
フィルハンドルをダブルクリックすると早いです。
数式がコピーされているのを確認します。
そしてF274に以下の式を入力します。
F274: =ABS(B272-10)
数値積分の結果でt=270のときのPの値であるB272とその時の理論値の発酵度合いの値10の差をとってその絶対値を出力しています。
これが0になれば、積分の結果と理論値が等しいと考えられます。
さて次はソルバーの設定です。
リボンのデータから右端のソルバーをクリックします。
ダイアログを設定していきます。
目的セル：F274
目標値：最小値
変数セル：H1
制約条件：H1≦1,H1≧0
解決方法：エボリューショナリー
以上を設定したら解決ボタンを押します。
うまく計算できない場合何回か同じ設定で解決ボタンを押してみてください。
エボリューショナリーの初期点の生成がおそらくランダムなので答えが出る前に変な収束をする場合があるようです。
結果が出ると
F274: 7.21011E-06
H1: 0.0259034261347831
理論値のμは0.0259でしたから十分な値が求まりました。
以上「数値積分するための径数μを微分方程式を手作業で解かずに求める」でした。
