Mathematica 演習
ver.2010.1.8
ここでは情報処理実習２（岩舘先生担当）の１つとして、新しい演算専用ソフトMathematicaを学習します。これまで勉強してきたソフトと共通する所もありますし、全く新しい使い方も数多くあります。ソフトの特徴を習得した後で、酵素反応などを表現した微分方程式を解く例を見ます。ここでは、プログラムを新しく作るよりも、計算結果を理解することを目標にしたいと思います。
STEP 1: 微分方程式について
微分・積分の定義やそれを使った方程式は、理解も大変ですが、紙の上の計算では、なかなか答えが出て来ないのが悩ましい所です。ここでは、微分方程式の本質的と思える箇所だけ、概要を理解するように努めて下さい。あとの作業は、すべてMathematicaが皆さんの代わりに解いてくれます。
まず、ある「状態」が変化する様を想像してください。状態がA→Bと変化するようなケースです。「酸素→二酸化炭素」でも良いし、「子供→大人」だったり、「後楽園駅にいる状態→池袋駅にいる状態」でも結構です。具体的なものをイメージすると考えやすいと思います。状態の変化は一瞬では起こりえないので、その速度を定義します。その速度とは、化学反応の速度、成長の速度、丸の内線の電車の速度だったりします。これをkとします。
ここで、ある単位時間（例えば１秒）の間に、A→Bへ変化するものの総数は、k×Aとなることに気づきます。酸素の分子がどれだけあるかで、酸素→二酸化炭素と変化する量が異なる（比例する）のと同じです。そもそも子供がいないと「子供→大人」の反応は起こりませんし、後楽園に人がいるので、その数に応じて「後楽園駅にいる状態→池袋駅にいる状態」の状態変化も発生します。
さて、変化した結果どうなるでしょうか？AからBが生まれます。二酸化炭素が生まれるし、大人が増え、池袋にいる人が増えます。増える速度を、１秒前と現時点での差、B(t)-B(t-1)として表現すると
k×A(t)=B(t)-B(t-1)
が成立します。本当は、わずか１秒間の短い時間でも、刻々と状態や速度が変化するので、上の式はそう正確ではありません。ある非常に短い時間Δtの間で起こる変化として正確に表現すると、上の式の左側は[B(t)-B(t-Δt)]/Δtとなります。Δtで割り算するのは、１秒間の平均速度に換算するためのものです。つまり
k×A(t)=[B(t)-B(t-Δt)]/Δt]
となり、右側は、Δtの短い時間の変化速度となります。ここでΔtをずっと短い時間へ、極限０に近づけると、微分と同じ定義となります。それが次の式です。
右の３つ、
はすべて同じことで、学問分野や習慣で表記が多少変わっているだけです。これで、状態の変化は、すべて微分方程式に置き換えて表現できることがわかりました。
さて、世の中はそう単純ではないことも多く、B→Cへとさらに別の状態に変化したり、B→Aで戻ったりもします。池袋で乗り換えたり、帰りの電車に乗ったりするのと同じです。これは、例えば、
と表記することができます。こういった複数の微分方程式の組み合わせで、多くの現象がうまく記述できます。より、正確には、A(t)+B(t)+C(t)は常に一定なので、表記は、下のようになります。
この４つの式の中の３つあれば、過不足なく現象を記述したことになります。増加・減少の関係と、kの前の＋や−の符号は間違わないように注意します。あとは、最初がどんな状態だったのか、すべてAか、すべて酸素か、すべて子供か、皆が後楽園にいるのか、といった初期条件を入れると良いことになります。もし、もっと複雑な現象、例えば、速度そのものが一定ではなく、他の条件で変化する場合などは、
などと、条件を追加すれば良いことになります。これは、例えば、後楽園から池袋に行く乗客の数の多いか少ないかで、電車の速度が変わるようなケースです。この計算が正確に解けるかどうかは別問題ですが、数式として正確に記述することは可能です。問題の正解は、Mathematicaが教えてくれます。さて、ここまでが概論です。まずは演算ソフトを使ってみます。
STEP 2: Mathematicaを起動して計算
下の手順に従って、Mathematicaを使った計算演習 1)〜4)を実施します。時間的な余裕のある場合だけ、先に進むようにしてください。プログラムを実行させ、バグなどがなく、計算結果が表示されたのを確認できたら、ファイル（.nbと拡張子が付いています）を指定のファイル名で保存します。保存したファイルの提出方法は、講義の時の指示に従ってください。
1) Mathematica(ver.7)を起動する。
1-a 情報処理実習室のコンピュータを各自の登録ID名で立ち上げます（立ち上がっていますね）。
1-b Mathematica 7.0を実行します。左上の「アプリケーション」→「数値計算/グラフ化ソフト」→「Mathematica」と進めます。この方法は、他のソフトと同じですが、不明の時は、手を挙げてTAの学生に聞いてください）。
1-c 「Untitled-1」や「名称未定義-1」と表示されたら、次の作業に移ります。
2) 簡単な計算の実行練習を行う。Mathematica独特のコマンド表記方法に慣れる。
下の参照ファイルにあるコマンド例をよく見て、間違いがないように注意しながら、キー入力します。スペースや大文字・小文字の違いも注意してください。スペースは×（かけ算）を意味します。(* ・・・・ *)はコメント行で、各コマンドの注意書きが記されています。これは入力する必要はありません。
１行入力したら、その度に、「Shift」+「Enter」を同時に押してコマンドの実行内容を確認します。その場で答えがでます。In[*] が入力した計算式、Out[*] が計算結果です。下の参照ファイルにはない計算ですが、例えばaa、a a、a+aなどの計算や、定数を少し変えた計算を平行してトライすると、式の意味がよく理解できるようになります。この作業はMathematica独特にコマンドになれるためのものです。
参照ファイル：Mathematica 基礎練習
リンク先にあるpdfファイル上のテキストは、テキストデータではなく画像ファイルです。そのままコピー＆ペーストできませんので、注意して下さい。練習なので、１つ１つ正確にキー入力します。計算をすべて終了して作成したファイルは、001.nbの名前で保存してWebClassから「ファイル１」として提出して下さい。ここまで、時間内に終わったら、次のステップに進みます。時間がなくなったら、途中までの実行結果でも構いませんので、１のファイルとして提出してください。最後は「ファイル」→「終了(Q)」で終了させます。
3) A→B方式の微分方程式を解く。
さて、ここまで来たらMathematicaに合わせて方程式を作るだけです。状態がA→Bと変化する単純なケースで解を求めてみましょう。微分方程式で表現すると、
となります。ここでは、Aをある決まった関数で与えて、その後のBの変化を調べるという課題を２つの条件で解くことにします。
２つの条件とは、それぞれ、A(t)が、以下の２つの関数で変化する条件です。
少し複雑ですが、それぞれ、
（最初の条件）
ClearAll["Global`*"];
k=10
T0=2
A0=1
AA=50
Plot[k*A0*(1-Exp[-AA*(t-T0)])*(1+Sign[t-T0]),{t,0,5}]
NDSolve[{B'[t]==k*A0*(1-Exp[-AA*(t-T0)])*(1+Sign[t-T0]),B[0]==0},B,{t,0,5}]
Plot[B[t]/.%,{t,0,5}]
（２つ目の条件）
ClearAll["Global`*"];
k=10
T0=2
A0=50
AA=100
Plot[A0*Sqrt[AA/Pi]*Exp[-AA (t-T0)^2],{t,0,5},PlotRange->{0,400}]
NDSolve[{B'[t]==k*A0*Sqrt[AA/Pi]*Exp[-AA (t-T0)^2],B[0]==0},B,{t,0,5}]
Plot[B[t]/.%,{t,0,5}]
をMathematicaの上で実行させると計算結果が得られます。これは、すべてを１つの段落内（右側に ] のマークがあります）にまとめて記入し、全体をいっぺんに実行させます。コピー＆ペーストで実行可能ですが、時々、文字化けでうまく作動しないケースもありますので、上の原文を確認して実行してください。はじめのk=10などの式は、速度などのパラメタの設定です。最後の３行で、計算結果を表示しますが、はじめにA(t)のグラフが、次にB(t)を示すグラフが表示されます。上のまま入力し、パラメタをいろいろ変えながら、実行してみてください。
＜問＞上にあげた２つの条件は何が異なることになりますか？Ａ→Ｂと言う化学反応を想定した時、この２つの条件の差は、具体的に、あるいは、実験の上では、どのような違いに相当するとかんがえられますか？ Mathematicaのファイル内に、(* ・・・・ *)のコメント行として答えを書き込み、計算の実行結果とともに002.nbの名前で保存します。WebClassから「ファイル２」として提出して下さい。ここまで、時間内に終わった学生は、次のステップに進みます。
4) ミカエリス・メンテン型酵素反応の応答を調べる。
さて、ここから、もう少し難しくなります。微分方程式を作ること、そしてパラメタを入れて計算する作業です。下の化学反応は、ミカエリス・メンテン型の酵素反応です。基質Ｓは酵素と結合して、ＳＥという中間体を作ります。その中から、新しく生成物Ｐができる反応です。この化学反応は、上の例にならって微分方程式で表現するとどうなるでしょうか？これはやがて２・３年生で習うかも知れませんが、頭の体操のつもりで、考えてみてください。これは特にレポートして提出する必要はありません。
この化学反応を計算するプログラムは、
条件Ａ：ここを右クリックして「名前を付けてリンク先を保存（K）...」で保存するnbファイル
式の解説：ここのjpegファイル
にあります。このいずれか一方をダウンロードして実行してみてください。上をクリックするだけでは実行できません。かならず一旦自分の管理下のフォルダ内に一旦ダウンロード・保存して、その後、Mathematicaのプログラムから「ファイル(F)」で開いて実行させます。下の問いに答えて、実行したMathematicaのファイルの中に(* ・・・・ *)のコメント行として答えを書き込み、計算の実行結果とともに003.nbの名前で保存します。WebClassから「ファイル３」として提出して下さい
＜問＞ プログラムの中のx(t)、y(t)、z(t)は、それぞれ、反応基質（Ｓ）、フリーの酵素（Ｅ）、生成物質（Ｐ）の濃度を示しています。速度定数k10を変えたとき、最終的な生成物の濃度の変化（時間経過）に、どのような差が出てくるかを計算で試しなさい。その結果を簡潔なコメント文としてファイル内に書き込みなさい。パラメタをあまり大きく変えるとグラフ表示ができなかったりするので、注意。
さらに詳しくやってみたい学生は下へ