常微分方程式\[\frac{\text{d}y}{\text{d}x}=-y\]を, 初期値\[y(x=0)=1\]を用いて,\[0<x<1\]の区間で解いてみよう.

差分法

差分近似とは, 微分を有限区間$\Delta x$の間の関数値の差で近似するものです:\[\frac{\text{d}y}{\text{d}x}\approx\frac{y(x+\Delta x)-y(x)}{\Delta x}\quad\text{あるいは}\quad\frac{\text{d}y}{\text{d}x}=\frac{y(x+\Delta x)-y(x)}{\Delta x}+O(\Delta x^n)\] これを利用するため, まず方程式を解きたい領域$0<x<1$を$N$等分しましょう:\[\Delta x=\frac{1}{N},\quad x_i=i\Delta x (0\leq i\leq N)\]この「離散点」$x_i$における関数値を\[y_i=y(x_i)\]と定めます. 方程式の差分近似は\[\frac{y_{i+1}-y_i}{\Delta x}=-y_i\]あるいは\[\frac{y_{i}-y_{i-1}}{\Delta x}=-y_i\]となります.

プログラムとして, ODEクラスで上の解析を行うことにします. その場合，main()は次のようになるでしょう:

#include <iostream>
#include "my_first_solver.hpp"      //  自分の解析を書き込むクラス定義を読み込む
int main(int argc, const char * argv[]) {
   ODE my_ode;                            //   ODEクラスのオブジェクト my_ode を作成
   my_ode.solve();                         //   ODEクラスのsolve()を実行
   return 0;
}

my_first_solver.hpp には, ODEクラスの定義を書くことになります. 前者の差分式は\[y_{i+1}=y_i-y_i\Delta x\]と書けるので, 次の作戦になるでしょう:

まず y0=1 としておく
i=1 では, y1=y0-y0*dx となる
i=2 では, y2=y1-y1*dx となる
面倒くさいので, y=y-y*dx で良いんじゃないですかね．

なお，プログラムを書くときに「めんどくせえ」と思う心は極めて重要です．なぜならば，一番面倒くさいと思うはずなのは寡黙なコンピューターの方であり，あなたがめんどくせえと思わないと，高速な計算にならないからです．というわけで

...省略...
class ODE
{
public:
   ODE(){};
   void solve(){
      std::cout << "ほんじゃま, 解いてみっか？" << std::endl;
      int N=100; //100等分してみる
      double xmax=1.0; //x=1まで解く
      double dx=xmax/N; //ってこた区間幅はこうなる
      double y=1.; //初期値は y=1
      std::cout << 0.0 << " " << y << std::endl;
      for(int i=0;i<N;i++)
      {
          y=y-y*dx;
          std::cout << dx*(i+1) << " " << y << std::endl;
      }
      std::cout << "とけたぞコラ" << std::endl;
    }
};
...省略...

これを再生してみると, 結果が得られるでしょう:

解析解は$y=e^{-x}$ですから, 厳密な値は$y=0.367879\cdots$であります. 結構誤差があることが分かりますね!

可視化

結果をグラフに描くことは大切です．いろいろな方法があります．

ExcelとかNumbersとかの表計算ソフトに貼り付けて，グラフにする
- 残念でした. 流体力学界では図はjpeg, png などイメージファイルではなく，pdf や eps などベクトルファイルで提出することが常識です．で, Microsoft製品はベクトルファイルが極めて汚く, 役に立ちません.
- まあでも手っ取り早く描くのに慣れているから良い，という向きもあるでしょう．でも，それを論文で提出してはダメですよ．
Matlab あるいは Octave を利用する
- 京都大学~~流体力学研究室~~はMATLABライセンスの正式ユーザーですので，ライセンス数無限で自由に利用できます．どんどんダウンロードして利用してください．
- 絶賛推奨しているにも関わらず，なかなか学生さんが使ってくれません．
  - とりあえず起動が重い．YouTubeを見ながら作業すると，再生が止まったりする
  - MATLABライセンスを利用するには研究室にVPN接続せねばならない. したがって，エロビデオを見ながら作業することができない．
GNU-plotutilsを利用する
- 1980年代のソフトで，古い→低性能，という刷り込みによってユーザーが消えました．
  - Windowsでは動作しない，ってのがやばかったのかもしれない
- 実は高性能なので，流体力学研究室では現役です．実はOctaveの描画エンジンです．
Python/Matplotlibを利用する
- ユーザーが増えないですね．なぜでしょう？
Gnuplotを利用する
- それならMatlabやOctaveを利用した方が手っ取り早い気がする.

ここでは, 時間の節約のためGNU-plotutilsを利用します．これ (あるいは, このページ末尾から plotutils-2.6-新しげ.pkg) をダウンロードしてインストールしてください. また，論文の作成ページの

GhostScript たぶんこっちは, 以下で必須だな. ここの添付のやつ入れておけばOK
ImageMagick

をインストールしておいた方が良いです. どっちにしろ卒業するまでにはインストールする必要がありますので．

データファイルの作成

グラフを描くには, データをファイルに保存する必要があります．

コピペ：XCodeの出力画面で【CMD】+Aすると全選択ができます．【CMD】+Cでコピーして，エディターで【CMD】+Vすると，データファイルを作成できます.　100万行程度までは，これで全然大丈夫です．

ですが，毎回コピペでは作業が遅いです．ここでは，プログラム中でファイルに出力することを覚えましょう．

ファイルをフォルダーに入れる，という概念が希薄な初心者は，ここを読みましょう．

コンピュータには，「あんたのファイルは基本ここに置く」という「ホームフォルダー」ってのがあります．そこにファイルを置いたほうが便利です．ホームフォルダー？という人はここを読みましょう

ファイルを作成する例文を示しましょう．以下では，$HOME/tes_test というフォルダーに，my.data というファイルを作成して出力しています．

...前略...
#include <boost/filesystem.hpp>
...中略...
void solve(){
    //$HOME/tes_testにファイルを出す(フォルダーはあらかじめ作っとけ)
    boost::filesystem::path folder=getenv("HOME"); // あなたの基準フォルダーHOMEを取得
    boost::filesystem::current_path(folder/"tes_test");
    //ファイルを開く
    std::ofstream ough("my.data");
    std::cout << "ほんじゃま, 解いてみっか？" << std::endl;
    ...略...
    double y=1.; //初期値は y=1
    ough << 0.0 << " " << y << std::endl;
    for(int i=0;i<N;i++)
    {
        y=y-y*dx;
        ough << dx*(i+1) << " " << y << std::endl;
    }
    ough.close();  //プログラムが終了すると自動で行われるので省略可能
    std::cout << "とけたぞコラ" << std::endl;
    ...後略...

ではターミナルを開いて,

cd データファイルを置いたフォルダーを実行しておきます．

head というのを使うと, 長大なデータファイルの先頭を表示することができます(tailだとケツを表示します)．このように，データファイルにはX Yのペアを羅列するようにしておいてください.

EPSファイルを作成

EPSってのは古式ゆかしいファイル形式で，今でもこの形式で論文投稿を求められることが多いものです．上のデータは, 次のコマンドでEPSファイルに変換できます:

graph -Tps my.data > my.eps

作図できたか確認するには, 次のどちらかを行います.

open my.eps を実行
open . でfinderを開き，my.epsをクリックします:

MacではEPSを開くとなんでかPDFに変換されます．ので，【File】【PDFで保存】を行うとPDFファイルができます．

PDFファイルを作成

いや面倒だから最初っからPDFを作成して欲しい場合，次のコマンドを実行します:

graph -Tps my.data | eps2pdf - > my.pdf   (graph -Tps my.data | eps2pdf -> my.pdf でもOK)

なにがそんなに面倒なのかというと, 何かミスって再作成した場合，

EPSであると，もう一度クリックしないと再表示されない
PDFであると，再作成を検出して再表示される

という違いがあって，EPSだとクリックの手数が多いんすよね・・・オイラ生まれて一度も間違ったことをしたことがない，って向きならEPSでどうぞ

pdfファイルを表示したければ

open my.pdf

クリックしないと精神の安定が保てない場合は

open .

でFinderを開いて, my.pdf をクリックしてください.

JPGファイルを作成

おれはPowerPointを使うという場合, pdfファイルやepsファイルでは役に立ちません. そういう場合には，

pdfやepsファイルをクリックし，【ファイル】【書き出し】でJPEGファイルにする
convert my.pdf my.jpg コマンドで変換する　（少々オプションを工夫しても, 上の方が画質が高い）

ことになるでしょう．

作図の雰囲気を変えたい

それはここを見ながら変えてください. これくらいならすぐです:

常微分方程式 I

差分法

可視化

データファイルの作成

EPSファイルを作成

PDFファイルを作成

JPGファイルを作成

作図の雰囲気を変えたい