常微分方程式\[\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.N=100;                          //   分割数を決定
   my_ode.solve();                        //   ODEクラスのsolve()を実行
   return 0;
}

これで, ODEクラスの仕様が決定したことになります. すなわち

ODE.N で，分割数を指定できる
ODE.solve() で, 解析を実行する
解析結果を，どのように処理するつもりなのかは，未定

クラスを書く時には，このように，「最初に仕様書を策定する」のが良いです．良い仕様書であれば，そのままプログラムになるからです．

書いて行く途中で，外部の情報が必要になったりしたら，それは「悪い仕様書」です．エンジンの開発部門が，数ヶ月おきにシャーシ開発部に「あ，ここにネジ穴追加」「あ，その重量配分やばいっす」とか伝達すると，そのうちシャーシ部門が反乱します．

ODEクラスの定義をODE_1.hpp に書くとしよう．前者の差分式は\[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 で良いんじゃないですかね．

では書いてみましょう．なお ODE は ordinary differential equation すなわち常微分方程式で, そこにsolve()関数と書いたところで, AIが勝負をかけてきました：

・・・いや，そう書きたかったんですけど状態ですねえ．TABキーを押すと，AIの言う通りになります.　

ま，違うところは, AIは\[\frac{\text{d}y}{\text{d}x}=f(x,y)\]を解こうとしていると読んだのですが，今は\[\frac{\text{d}y}{\text{d}x}=-y\]でいいから, そこは訂正しましょうね．

繰り返し：for文

上のサンプルで, 「ある作業を繰り返す」基本のやり方が書いてあります．

for(size_t i=0  ; i<N ; ++i){
    ....
}

これは，数えながら実行するfor文です．

この辺で，数を数える（非負）整数 i を定義します．

この辺は，繰り返しを続ける条件を書きます．繰り返し数iがNより小さければ続行

っこた，N=100 なら i=0, 1, 2, .., 99 まで実行するので100回ね！i=100では実行しないので注意

ここは，赤色の部分を1回終わったら何をするか，を定義．現在，1回終わったら繰り返し数を1増やす

++とはなにかというと, 1歩前に進む，という演算です．--だと一歩下がる. なのでC++はC言語より一歩前進したもの，ということらしい

繰り返し数の管理は当たり前っぽいですが自分で書かないとダメ！書き忘れると，無限に繰り返してしまうので注意

繰り返しは頻繁に使うので，いろんなバージョンがあります．

こいつらに，全部やれ
この条件が満たされる限り，際限なくやれ

とか，色々あります．数えながら行うのは基本ですね！

やってみる

では，これを呼び出す omain.cpp を作りましょう．数単語書いた段階でAIは先を読みました：

ん. まあ動くよね．

で，解けました:

可視化

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

ExcelとかNumbersとかの表計算ソフトに貼り付けて，グラフにする
- 残念でした. 流体力学界では図はjpeg, png などイメージファイルではなく，pdf や eps などベクトルファイルで提出することが常識です．で, Microsoft製品はベクトルファイルが極めて汚く, 役に立ちません.
- まあでも手っ取り早く描くのに慣れているから良い，という向きもあるでしょう．でも，それを論文で提出してはダメですよ．
Matlab あるいは Octave を利用する
- 京都大学はMATLABライセンスの正式ユーザーですので，ライセンス数無限で自由に利用できます．どんどんダウンロードして利用してください．
- 絶賛推奨しているにも関わらず，なかなか学生さんが使ってくれません．
  - とりあえず起動が重い．YouTubeを見ながら作業すると，再生が止まったりする
  - MATLABライセンスを利用するには研究室にVPN接続せねばならない. したがって，エロビデオを見ながら作業することができない．
    - 現在ではsshuttleの柔軟な設定で，エロビデオはご自宅の回線を使い，研究は大学の回線を使う，というふうに設定できます．頑張ってみてね！（失敗すると，「国民の皆様の税金を何に使っとるんだ?」というお怒りのメールを情報環境機構から送られるので十分に注意）
GNU-plotutilsを利用する
- 1980年代のソフトで，古い→低性能，という刷り込みによってユーザーが消えました．
  - Windowsでは動作しない，ってのがやばかったのかもしれない
- 実は高性能なので，流体力学研究室では現役です．実はOctaveの描画エンジンです．
Python/Matplotlibを利用する
- ユーザーが増えないですね．なぜでしょう？欧米人はこればっかりだよね，いま
Gnuplotを利用する
- それならMatlabやOctaveを利用した方が手っ取り早い気がする. 若者がこれを使うのをみたら，年寄りが喜ぶ以外の利点は無いような気がする.

データファイルの作成

多くの作図ソフトは１行に　X (空白) Y と書いてあるデータを理解します．なので，ODE_1.hpp の最後の方で，そんな出力になるように

std::cout << (i+1)*h << " " << y << std::endl;

にしとく必要がありますね. で，

　全てコピー，【新規テキストファイルの作成】で

ま，何かのファイル名で保存しましょう．左では「data.data」に保存してる

PlotUtilsで作図

ターミナルで，次のコマンドを入力し，適当な名前のPDFファイルに保存

% graph -Tps data.data  |  ps2pdf -dEPSCrop - new.pdf

% open new.pdf

GhostscriPDFをインストールしている場合，次のコマンドを入力し，適当な名前のEPSファイルに保存しても良い:

% graph -Tps data.data   > new.eps

% open new.eps

こう言うのが出る:

PDFファイルだから，論文にも使えるよ！

線の太さとか変えたい場合は，ここを参照．

MatPlotLibを使う

こっちに書いた

常微分方程式I

差分法

繰り返し：for文

やってみる

可視化

データファイルの作成

PlotUtilsで作図

MatPlotLibを使う