[Excel2013 で VBA] 微分方程式の数値的解法 (1) Excel のみの方法

前回、VBA から Excel のシートにアクセスする方法を学んだ。
この機能を用いると、様々な数学的な操作を Excel 上で行うことができるようになる。

ここでは「微分方程式の数値解法」を取り上げ、VBA (および Excel) 学習の総仕上げとする。

微分方程式の数値的に解くのは、「Excel のみを使う方法」と「VBA と Excel を組み合わせて使う方法」のどちらでもできる。
ここではその両方を取り上げ、その違いを体感してもらう。

準備

本日取り扱う微分方程式は以下のものである。



この微分方程式は、ばね系に周期的な外力が加わっている系の運動方程式



において、質量 m、ばね定数 k、外力の大きさ F、外力の振動数 ω に対し、 以下の制限を設けたものである。



この微分方程式 (運動方程式) は 2 階微分を含むが、これを数値的に解くためには 1 階微分のみの微分方程式に変換する必要がある。
そこで、以下の 2 変数 x1、 x2 を導入する。 (それぞれ、ばねの変位、速度を表す)



これらについて、以下の 1 階微分方程式が成り立つ (確認してみるとよい)。
(微分の階数が 2 階から 1 階に落ちている変わりに、変数が一つ増えている)



この形の微分方程式を一般化して以下のように書くことができる。



ここで、微分係数に対して以下の近似式を代入すると、



以下の式が得られる。これは、時刻 t における系の状態 (x1(t), x2(t)) から、時刻 t +Δt における系の状態
(x1(t+Δt), x2(t+Δt)) を求める式である。
この手法を Euler (オイラー) 法と呼ぶ。 この手続きを繰り返すことで、任意の時刻の系の状態を知ることができる。




また、このように微分方程式を数値的に解くことを数値積分する、などと言う。

なお、上記の議論は本日の内容の背景の解説としてつけたのものであり、本講義では覚える必要はない。



Euler 法で数値解を求めてみる

では、実際に Euler 法で数値解を求めてみよう。
冒頭で述べたように、ここでは Visual Basic によるプログラミングを使わず、 Excel の機能のみを使って実現してみよう。

Excel の操作法としては新しいことは学ばないので、これはExcel 利用法の復習と思って欲しい。
(詳細は Office 2007 の基礎 参照)

まず、準備として、第 1 列に以下のように「t」、「x1」、「x2」、「f1」、「f2」、「dt」と記述しよう。
この位置を合わせないと、以下の式も各自変更する必要があるので、図と同じ位置のセルを利用すること



次に、時間 t、位置 x1、速度 x2 に初期値として 0 を記述し、数値積分の時間刻み dtに 0.01 を記述しよう。



次に、f1(x1, x2, t) = x2 を実現するため、セル D2 に

=C2 

と記述しよう (上の枠つきの文字はコピー可能)。



次に、f2(x1, x2, t) = -4 x1 + 2 cos(2t) を実現するため、セル E2 に

=-4*B2+2*COS(2*A2) 

と記述しよう。



次に、時刻 t を dt だけ増やすため、セル A3 に

=A2+F$2 

と記述しよう。セル F2 を定数値として常に参照するため、絶対参照を用いていることに注意。



記述が終わると以下のように、セル A3 に 0.01 だけ進んだ時刻が表示される。
セル A3 が選択された状態で、セルの枠右下の■をマウスでつかみ、下方向へドラッグしてオートフィルしよう。



やや多めに、オートフィルで 1000 くらいコピーしておこう。



次に、x1(t+dt) = x1(t) + f1 (x1,x2,t)*dt を実現するため、
セル B3 に

=B2+D2*F$2 

と記述しよう。



同様に、x2(t+dt) = x2(t) + f2 (x1,x2,t)*dt を実現するため、
セル C3 に

=C2+E2*F$2 

と記述しよう。



セル D3 とセル E3 には、セル D2 とセル E2 に書いた式をコピーしよう。
2 つのセルをドラッグして選択し、右クリックでコピーし、



D3 セルを選択し、右クリックで貼り付ける。



すると、D3 と E3 に式がコピーされて値が表示される。



あとは、B3、C3、D3、E3 セルの式をオートフィルでコピーしていこう。
B3 から E3 までマウスでドラッグし、4 つのセルを選択する。
4 つのセルが選択された状態で、セルの枠右下の■をマウスでつかみ、下方向へドラッグしてオートフィルしよう。




A 列と同じ位置までコピーすればよいだろう。



さて、以上で Euler 法による微分方程式の数値解のデータが準備できた。
このグラフを描いて確認してみよう。

以下のように、A 列と B 列を共に選択する。
「A」を一回クリックした後「Ctrl キーを押しながら、『B』をクリック」または 「A、B をマウスでドラッグ」で選択可能である。



この状態で挿入タブから散布図を選択するのであった。

これにより、以下のようなグラフが描ける (これはまだラベルや軸名などを描いていない)。
これが微分方程式の数値解である。





数値解と解析解を比較してみる

以上で、Excel の機能のみを用いて微分方程式の数値解を Euler 法で求められることがわかった。

この数値解が、解析的な解とどれだけ異なっているかを調べてみよう。

この系の解析解は



で書けることが知られている (ただし、時刻 0 で位置と速度が 0 の場合)。
(なお、この解析解を得る方法を知るためには「微分方程式」を学ぶ必要がある)

数値解と解析解を同時にグラフに描いて、差を視覚化してみよう。

そのためには解析解のデータを準備する必要がある。

G 列に 以下のように解析解を用意しよう。書き込む式は

=0.5*A2*SIN(2*A2) 

である。



そしてオートフィルで必要数のデータを揃える。



このデータのグラフを描くには、まず、A、B、G 列を同時に選択する。
そのためには、「A」を一回クリックした後、「Ctrl キーを押しながら、『B』、『G』を一回ずつクリック」すると良い。



この状態で挿入タブから散布図を選択する。



すると下図のようなグラフが得られる。

重要なことは、Euler 法による数値解と解析解の間に大きな誤差があることである
これは Euler 法は (簡単ではあるが) 非常に性能の良くない数値積分法であるためである。

この誤差は次ページで学ぶ 4 次の Runge-Kutta 法を用いることでかなりの程度解消できる。





なお、グラフの凡例 (上の図で「Euler 法による数値解」や「解析解」と書かれた部分)
を変更するには、以下のように「元のデータ」から行う。



「データソースの選択」というダイアログが現れるので、 2たる系列のうち、編集したい方 (ここでは x1) をクリックして選択し、 「編集」ボタンを押す。



以下の「系列名」に表示したい内容を記述する。







←シートへの高速なアクセス微分方程式の数値的解法 (2) VBA + Excel→

Excel 2013 で学ぶ Visual Basic for Applications (VBA)に戻る