今回はちょっとお遊び的かもしれません。まずは、以下のグラフを見てください。これは、ある式から作ったグラフです。一体これはどういう式から作られたものでしょうか。

「なんだ、ただの sin x のグラフじゃないか」と思うかもしれません。しかしこれは、dy/dx = cos x、y(0)=0 という微分方程式の初期値問題からコンピューターに数値データを作らせてそのデータから描いたグラフなのです。
　
　前回考えたような一階の微分方程式は全て解を出すために積分が必要でしたが、このようでは、手計算で微分方程式を解く問題においては、(現時点では) 積分によって原始関数を求められなければ即詰んでしまいます。積分できる関数、できない関数、どちらも無数にありますが、通常手計算では解くことができない方程式の方がずっと多く、そのような場合にはコンピューターを利用して近似的に解を推測するという方法が一般に用いられます。
　今回は、コンピューターによる微分方程式の図的解法の導入部分をかじってみようと思います。ここは学校の授業では触れたことすらなかったし、別に詳しくやらなくてもいい気がするのですが、先日購入した Shaum's Outlines の Differential Equations の演習を進めていたところ、Chapter 18 - Graphical and Numerical Methods Solving First-Order Differential Equations (一階微分方程式の図的、数値的解法) において、今回扱うような図的解法、数値的解法が紹介されていて、個人的にこれをコンピューターでやったら面白そうだと思った、言い換えれば、手で計算するのはすごい面倒だからコンピューターにやらせたら楽そうだと思ったので気分転換ということで。

方向場

　一階の微分方程式では、積分が出てきましたから、当然積分が出来なければ詰んでしまいます。そこで、積分が出来ないので手計算では明示的に解を出せなくても、グラフ的にはこういう動きを示す関数だ、というのを近似的に推測する手法も発展してきました。方向場 (Direction Fields) は、グラフのおおまかな流れを知る上では非常に有用な手段です。
　その手法はこうです。まず、一階微分方程式は導関数がただひとつしか出ないことを利用して、導関数 = x,y の関数という形に分離できます。つまり、

と表現できるということです。これは、x, y の組を適当に決めてやれば、その点における微分係数が分かることを意味しています。微分係数は接線の傾きですから、その点における関数の接線の傾きが分かるというわけですね。
　そして、関数の傾きというのは、その点における関数の直線近似でもあるわけですから、微分係数をもとにその点から短い線を延ばし、たくさんの点について同じように適用してやれば、関数のおおまかな流れが視覚的に分かるのではないか、ということになります。これが方向場の考え方です。
　ただ、あんまりデータを偏って採取してもしょうがないので、通常は、まんべんなく調べる、つまり等間隔で x, y の幅をとってその点について調べていきます。線の引き方は、全く同じ長さだけ、そして二つの点同士内向きにしても重ならない程度の短さで引くことにします。
　このように引けば、方向場が描けますが、あらゆる任意定数 C に対しての方向場であるため、書いてない空間含め無数のグラフが存在します。実用上はこのうち、初期条件を満たす 1 本を決めればいいだけなので、次に初期値問題からどのようにグラフを決定するかを考えます。

オイラー法

　初期値の組が一つ与えられることは、積分定数 C を決めることに等しく、上における方向場のなかから、条件に合致する C をもつ一つの曲線が描かれることとなります。したがって、初期値の組が分かっていれば、計算によって近似解を求めることが可能です。この方法が色々あるようですが、もっとも単純で分かりやすい方法としてオイラー法 (Eular's Method) というものがあります。
　このアイデアは非常に簡単です。まず、初期値の組、その点を出発します。方向場で見てきたように、

と変形し、f(x,y) を決定します。すると、初期値の点において、その点の微分係数が当然この式からわかります。接線は、その点からごく近い範囲ならば十分精度のよい直線近似となりえますので、その "ごく短い距離" だけ x を進めれば (これを Δx とする)、近似的に

と表現できることが分かります。Δx さえ小さければ、十分近い近似ですから、この調子でまた Δx 加え・・・とやっていけば、解 y のグラフの数値データが近似的に得られることがわかります。このときの Δx をステップ幅 (Step Size) といいます。個人的には分かりやすいので刻み幅と呼ぶことにします。

プログラミング環境とグラフ描画ソフトの準備

　上のような方法は、非常に有用ではありますが、実際これを高い精度で求めようとすると、幅は当然 1 や 2 といったおおざっぱすぎる値では使い物になりません。したがって、もっと何十、何百にも分割することが必要となります、しかし、今度は計算量が膨大すぎて鉛筆によってこれを実行することは不可能になってしまいます。というわけで、いまや誰もがパソコンを持つ時代ですから、上記の考え方をそのままパソコンにやらせてみます。これには、少々のプログラミングの知識が必要になるでしょうが、手法自体は単純なのでさほど難しくはならないでしょう。プログラミングのプの字も知らない私が書いているのだから大丈夫です。
　今回、以下のようにしてパソコンに図的解法を行わせます。まず、オイラー法のアルゴリズムを指定した関数について適用し、得られた数値データをファイルとして保存するようなプログラムを作成します。グラフを描くプログラムまで作るとすごく苦労するうえにここでやりたいこととは大幅に逸れてしまうため、ここではその出力されたデータに基づいて、グラフを描いてくれるソフトにデータを渡し、描画させます。
　グラフを描画してくれるソフトは多々ありますが、ここでは学校でもよく使われる GNUPLOT を用います。プログラムを作成するうえでは、これまた学校でよく使われる言語である C 言語を使って作成します。なんだか今回は学校のプログラミングの授業の題材としてありそうな感じがしますね。まあ、私のところはベクトルの内積と成分から角度を算出しろとか三角形の辺の長さが分かっているときの面積を求めろとか、1 年も授業したわりにすごい初歩的なのしか作らされませんでしたが・・。本当にプログラミング、またそれにより応用的な数学に興味を持たせたいなら、中学や高校でならったことの焼き直しよりは、こういうもうちょっと高度でかつコンピュータだからこそできるような面白い題材を取り上げるべきでしょうね。さて、C 言語のプログラムが作れるソフトなんていくらでもあると思いますが、無料で高機能となるとやっぱり Microsoft の Visual C++ 2010 Express しかないでしょう。と、いうわけで、今回はこの組み合わせで行きます。Windows 以外の場合は知りませんが、確か学校の授業では Emacs とかいうのを使わされてた記憶があります。
　GNUPLOT は、公式サイトからダウンロードしましょう。Download from SourgeForge というのを押すと、OS の選択に移動しますので、Windows binary とか書いてあるのを選択し、それでダウンロードできます。

これはインストーラー形式じゃないので、置きたいところにそのまま解凍しましょう。実行ファイルは、binary フォルダに入っている gnuplot.exe です。

試しに起動すると、次のような画面になります。

コマンドとかすごい面倒そうですが、軽くデータからプロットするだけならコマンド 1 個で済むのでそこまで苦痛ではありません。とりあえずここでは、まだデータすらできあがっていないので、そのまま終わります。× ボタンで終わるのもいいですがコマンドプロンプトの礼儀として、終了コマンドで終わってあげましょう (昔は × ボタンなどという便利なものはありませんでしたからね)。GNUPLOT は、exit コマンドを入力すると終了します。
　次に大仕事、Visual C++ 2010 Express の導入が残っています。公式サイトから直接ダウンロードできるので、ダウンロードしましょう。言語が違ってもやってることは一緒なので、Visual Basic か C# は扱えるが C 言語は駄目という場合は他でも多分よいでしょう (これらでもコンソールアプリが簡単に作れるかどうかは分かりませんが)。 別にダウンロード自体は大仕事ではないのですが、導入に非常に多くの時間を要します。インストールが終了するまで気長に待ちましょう。なお、こちら、一定期間後登録しなければ使えなくなります。しかし、製品登録も無料であるため Windows Live アカウントを取得して申請すれば簡単に登録キーが手に入るので、特に問題ありません。
　実際にプログラミングをする前に、Visual C++ の環境を整えます。Visual とついてるだけあってボタンとかいろんな部品を視覚的に取り付けてプログラミングできるのですが、データを出力するだけのプログラムにそんなものは必要有りませんから、昔ながらのコンソールアプリ (黒い画面に文字だけ) で作ります。Visual C++ を起動したら、"新しいプロジェクト" を選びます。

プロジェクト名を入力し、作成したら、次に出てくるウィザードで、空っぽのプロジェクトを作ります。

これをしないと、最初からなにやら色々書いてあるサンプルコードが付属してしまいます。最後に、完了を選べばおしまいですが、空っぽのプロジェクトなので、ソースコードを作らなければなりません。そこで、ソリューションエクスプローラーから、ソースコードのフォルダを右クリックして、"新しい項目を追加" を選びます。

ファイルを追加したら、試しに何かプログラムを作って実行しましょう。さすがにすごい基本的なことは、インターネット上にいくらでもあるであろう解説サイトに譲るとして (もしかしたらここで今後書くかもしれませんが・・)、とりあえず以下のようなものを書いて、正しく実行できるかどうか試します。

メニューバーのデバッグ → デバッグ開始から実行することができるのですが、この作業は今後何回も行うことになるので、一発でこれが出来るショートカットキー F5 は覚えておくべきです。なお、実行せずにファイルを作るだけでいいならソリューションのビルド (F7 キー) で大丈夫です。

このように、ちゃんとプログラムが作られていれば、何も問題がありません。getchar 関数は必要ないのですが、なぜ入れているかというと、実はこういったコンソールアプリケーションは MS-DOS プロンプトから実行されることが前提であるため、これ単体で起動してしまう (要するにデバッグで実行している今がまさにそう) と終了したとき勝手にウィンドウが閉じてしまいますので、それを防ぐためにわざとユーザーに何か入力させる関数を挟んで勝手にウィンドウを閉じるのを阻止しているだけです。なので、プログラムの内容的には、これ自体に特に重要な意味があったりするわけではありません。
　これで準備ができたので、どうやってオイラー法でデータを出力していくか、その方法を次で考えます。 　

プログラムの作成

　GNUPLOT の二次元グラフ描画に関しては、データファイルから点を打ってそれを線で結んでグラフにする、というようなことが可能であるため、そのデータファイルの書式にあわせてプログラムを作ればよいですね。GNUPLOT では、二次元の場合、x座標 y座標の組みを改行をはさんで並べていけばグラフを描くことができますので、それに沿って出力させます。要するに、書式は、

こんな感じになってればいいわけです。今回は、何でも解けるというような優秀なものを作るのは結構難しい、もとい、プログラミング超初心者の私には敷居が高すぎて無理 (おそらく、微分方程式の分類学を隅々まで知り尽くしているなら可能だが、ここでは入り口程度のことしか要求されていない) なので、陰関数は完全なかたちで解けなくてもいいという前提を作っておくことにします。さて、上のようなデータを出力するためには、以下の情報が事前に必要となります。

• x, y の初期値が 1 組
• x 座標の変域
• 刻み幅
• 関数 f(x,y)

　オイラー法は初期値から出発して、一定の刻み幅を決定したあと、左右隣はいくつだ、そしてまたその左右隣は・・・という手法なので、初期値が分かっていなければそもそも解くことができません。変域も必須です。どこからどこまで計算するのかを決定しなければ、コンピューターにはどうしようもありませんからね。ここでは、初期値の x 座標を出発して左右同じ距離だけ計算することにします。

イメージとしてはこんな感じですね。原点を起点としますから、どちらかの方向に行ったあと、また原点に戻ってもう片方を計算しないといけないのもわかります。なお、このとおりに計算するとひとつ不都合が起こってしまいますが、これについては後述とします。物理が題材の場合は、初期値はほとんど t=0 からでしかも t<0 は存在しないので、右一方だけでよいので、プログラムはたいへん簡単になります。この場合はこれから書くプログラムで左方向に関する記述を適宜切り捨ててください。
　刻み幅はたいへん重要です。これが唯一、精度に影響する材料となります。この幅が大きいと本来のグラフとは全く違った値が出てしまうことになりますが、理想の無限小に持って行くと今度は計算回数がコンピューターにとってすら膨大なものとなり、時間がかかってしまいますので、いいところで妥協しておきましょう。
　関数についてですが、関数も式のひとつの表現だと考えれば、#define で定義しておけばこの場合は十分ですね。つまり、いちいち関数を作らなくても、"#define function 式" で事足ります (function は任意の名前でよい)。しかし、#define で f(x,y) を表現する場合注意してほしいのは、この方法をとる場合、function は単純な置き換えですから、function のまわりにかっこをつけておくのが無難だということです。たとえば、function が -x だった場合は、5*function という表現は一見何も問題なさそうに見えますが、単純な置き換えなので 5*-function と、コンピューターにとっては意味の通らない内容になってしまいます。これを避けるため、そういう記述がしたい場合は、混乱が生じる可能性がある限り、5*(function) と記述しなければなりません。最初はお遊びで f(x,y) = x、つまり、y'=x なので二次関数を表す微分方程式で実験することにしましょう。function は x でいいですね。あと、使わないかもしれませんが微分方程式というと sin とか cos とかそんなのも良く出てくるのでねんのため math.h も #include しておきます。
　以下が、完成品です。さらっと書きましたが実は作るのに 3 時間かかりました。慣れないことはするものじゃないですね。

#include <stdio.h>
#include <math.h>

#define    function                2*x        // 関数f(x,y)にあたる表現をここで記述する
#define    step_size                0.01    // 刻み幅(小さすぎると処理時間が異常にかかるので注意)
#define    x_init_value            0.0        // 初期値のx座標x0
#define    y_init_value            0.0        // 初期値のy座標y0
#define    x_from_init_value        3.0        // x0から左右どこまでの範囲計算するか
#define WIDER_RANGE                        // 物理モデルの場合はt=0の左は考えないので、この記述は抜けばよい。

int main(){

    // data_dとdata_rの用途について
    // date_dでは、ファイル"data.dat.1"を作成し、マイナス側の計算情報を順に格納する。
    // 次に、data_rで扱う"data.dat"に、"data.dat.1"から行ごとに逆順で格納していく。
    // そうすれば、x座標が左の値-右の値へと綺麗に並ぶから、GNUPLOTで線をつないでプロットしても
    // ちゃんと線を引いてくれる。(GNUPLOTは上から順番に線でつなぐため順番が乱れているとおかしなことになる)
    FILE            *data_d, *data_r;

    data_d    =    fopen("data.dat.1", "w+"); // data.dat.1 は読み出しもあるから
    data_r    =    fopen("data.dat", "w");

    if(data_d == NULL || data_r == NULL){
        printf("ファイル作成に失敗");
        getchar();
        return -1;
    }

    double x,y;

#ifdef WIDER_RANGE
    x=x_init_value-step_size, y=y_init_value;y-=step_size*(function); // 初期値点の重複を避けるため
    while(x >= (x_init_value-x_from_init_value)){
        printf("%lf / %lf\n", x, y);
        fprintf(data_d, "%lf %lf\n", x, y);

        x-= step_size;
        y-= step_size*(function);
    }

    int pos;
    char buf[128];
    fseek(data_d, -3, SEEK_CUR); // Windowsでは、改行コードは2バイト扱いなので3バイト戻さないとおかしなことになる。
    while(true){
        if(fgetc(data_d) == '\n'){
            pos = ftell(data_d);
            fputs(fgets(buf, 128, data_d), data_r);
            fseek(data_d, pos-3, SEEK_SET);
            continue;
        }
        fseek(data_d, -2, SEEK_CUR);
        if(ftell(data_d)==0){
            fputs(fgets(buf, 128, data_d), data_r);
            break;
        }
    }
#endif
    x=x_init_value, y=y_init_value;
    while(x <= (x_init_value+x_from_init_value)){
        printf("%lf / %lf\n", x, y);
        fprintf(data_r, "%lf %lf\n", x, y);

        x+= step_size;
        y+= step_size*(function);
    }

    fclose(data_d);
    fclose(data_r);
    remove("data.dat.1");
    getchar();

    return 0;
}

ここでは、実行できさえすればいいのでエラー処理などそういうことはほぼやってありません。先ほど述べましたが、オイラー法は初期値座標を起点とするため、左右両方向に計算すると、左方向の場合、x が大きい座標から小さい座標へと推移していくことになります。ということは、順にファイルに書き出すと大きい → 小さい という順になります。逆方向の場合、また起点から起点からやり直さないといけないので、また起点に戻って今度は大きい方へとデータが書き出されていきます。こうなると、GNUPLOT は与えられたデータを順にプロットしていくため、一番小さい x 座標に達した時点で次にそこから起点に向かって線を引いてしまいます。これでは見た目的に悲惨なことになりますので、マイナス側のデータは逆順に並べ替えて出力する必要があります。
　今回は、データ量が膨大ではないので、ソースコード側の省エネのため上のコードのような方法をとりました。まず、本当に保存したいファイル、data.dat とは別のファイル、ここではdata.dat.1 にマイナス側へ順番にデータを書き出していきます。次に、data.dat.1 を最後から戻していって、改行コードを検出した時点で、そこから次の改行コードまでの行、つまり一行を読み込み、本来のファイル data.dat へ書き出していきます。これを繰り返せば、data.dat.1 のデータを data.dat へ逆順で格納することができます。プラス側はそのままでいいのでそのまま書き出しておしまいです。ファイル操作系の関数には "挿入" というのはなく、全部上書き扱いなため、面倒ですがこのような方法をとらなければうまくいきません。
　オイラー法自体についてはそこまで複雑じゃないので特に書くこともないと思います。ループ毎に刻み幅を x 座標に追加し、それに相当する増加量、刻み幅かける近似的なyの増加量(=微分係数)を y 座標に追加すればよいですね。
　物理モデルの場合は原点が起点である場合がほとんどなため、左側のプロットは必要ありません。このときは、#define WIDER_RANGE を取り除けば、コード内の #ifdef - #endif 間の部分、つまりマイナス側のデータ出力が丸ごと自動で取り除かれるため、この記述を取れば事足ります。

データ出力とグラフの確認

　上のコードは、dy/dx = 2x, y(0)=0 という微分方程式のオイラー法による近似解ですね。書くまでもないとおり、これは放物線です。本当に放物線になるのでしょうか。実際に実行してみましょう。

正常に終了すれば、このように大量のデータが表示されます。そして、

このように、data.dat が出力されます。場所は、デバッグで出力の場合はプロジェクトがあるフォルダの一番上の階層に出ますが、Debug や Release から直接実行する場合は実行ファイルと同じ場所に現れるため、その点に注意します。
　次にこれを gnuplot.exe と同じフォルダに移動します。そして、gnuplot.exe を起動します。起動したら、plot "data.dat" w l と入力します。

すると・・・

見事、放物線が現れました。
　次に、精度を確かめてみましょう。手計算で解けない場合有用だといっているのになんだか手計算で求めた結果で確認するというのも変な気分ですが、とりあえず x の二乗を同時にプロットします。これをするには、 , で区切ればいいだけです。つまり、「plot "data.dat" w l, x*x」 でいいですね。

このように、非常に高い精度で一致しているのが分かります。
　一方、刻み幅の重要性を見るために、step-size をわざと大きくしたものを用意してみましょう。刻み幅を 0.1 (最初の 10 倍) にした data.dat でもう一度同じことをしてみます。全く同じプロットがしたい場合は replot と入力すればいいので、data.dat を差し替えて replot します。

このように、刻み幅が大きいと誤差が相当大きくなってきてしまうことがわかります。なお、精度を上げるためには、刻み幅を小さくするという方法もありますが、別の計算方法をとるという方法もあります。このほか、修正オイラー法 (Modified Eular's Method) や、ルンゲ・クッタ法 (Runge-Kutta Method) といった数多くの方法があるため、より精度の高い方法についても今後試してみることにします。

いろいろな微分方程式を解かせよう

y'=cos(x), y(0)=0, 刻み幅 0.01, 初期値からの計算範囲が 3.14 (≒ π) のとき

これは、線形ではない上にベルヌーイの方程式等でもないので、今の段階では明示的に解くことができませんが、このようにコンピュータのおかげでだいたいの形が理解できます。

これは tan x であることは明白ですが、発散するような点を範囲に含んでしまうとうまく動かないため、そういう意味でもまだこのプログラムは未熟なことがわかります。ためしに π/2 を超える値を変域に設定してみるとよいでしょう。このような発散点を持つ関数は発散点よりも奥の部分は描くことができません。

一方、アークタンジェントは ∞ に行かないと発散しないため、このように高い精度で描くことが可能です。

すごくいびつな形ですね・・もはや元の関数がどういう式であるのか想像すらつきません。
　今回はここまでにしておきます。やっぱり数学も実験してみると楽しいですね。テイラー展開やフーリエ級数展開が実際の関数と一致するのを確かめるのも楽しかったけど、これはもっと楽しい。もし私が中学生の頃これを知っていたら間違いなくこれを自由研究の題材にしたでしょう・・・。まあ、自由研究なんて提出すらしなかったですけどね。