メインコンテンツに移動

常微分方程式 II

2階常微分方程式\[\frac{\text{d}^2 y}{\text{d}x^2}=-Cy\]の場合，解を一意に定めるには条件が2つ必要です.

初期値問題

\(x=0\)において, 初期条件\[ y(0)=1,\quad\frac{\text{d}y}{\text{d}x}=0\]が与えられている場合は, どうしましょうか. この場合, \[z=\frac{\text{d}y}{\text{d}x}\]とおいて, 連立常微分方程式\[\frac{\text{d}y}{\text{d}x}=z,\quad \frac{\text{d}z}{\text{d}x}=-Cy\]を初期値\[y(0)=1,\quad z(0)=0\]で解けばよろしい. 解を\(0\leq x\leq 2\pi\)で求めるとすると, プログラムはこんな感じですよね:

#include <cmath>  // 円周率πとか, 他にもいろいろ入ってて便利だし
...中略...
ODE(){ //めんどくさいし, コンストラクタでフォルダーを設定してしまおーっと
  boost::filesystem::path folder=getenv("HOME");
  boost::filesystem::current_path(folder/"tes_test");
};
...中略...
void solve2(){
  std::ofstream ofs("my.data");
  int N=100; //100等分してみる
  double C=1.;        //C=1ということにしよう
  double xmax=2*M_PI; //x=2piまで解く
  double dx=xmax/N; //ってこた区間幅はこうなる
  double y=1.,z=0.; //初期値は y=1, z=0
  ofs << 0.0 << " " << y << std::endl;
  for(int i=0;i<N;i++)
  {
    y=y+z*dx;
    z=z-C*y*dx;
    ofs << dx*(i+1) << " " << y << std::endl;
  }
}

結果はこんな感じか?

my

ですが数値計算というものは, 誤差がどれくらいであるのか？を知ることが必要です. 分割数を変えて結果を比較してみましょうね:

m3

これは10分割, 20分割, 40分割を比較したものです.

境界値問題

\(x=0, 1\)において, 境界条件\[ y(0)=1,\quad y(1)=0\]が与えられている場合は, どうしましょうか. この場合には, 2階微分に対する差分式\[\frac{\text{d}^2 y}{\text{d}x^2}\approx\frac{y_{i+1}-2y_i+y_{i-1}}{\Delta x^2}\]を連立させて解くことができます. すなわち \[ \begin{bmatrix}2-C\Delta^{2} & -1\\ -1 & 2-C\Delta^{2} & -1\\ & -1 & 2-C\Delta^{2} & -1\\ & & \ddots & \phantom{2-}\ddots\phantom{2-} & \ddots\\ & & & -1 & 2-C\Delta^{2} & -1\\ & & & & -1 & 2-C\Delta^{2} \end{bmatrix}\begin{bmatrix}y_{1}\\ y_{2}\\ y_{3}\\ \vdots\\ y_{N-2}\\ y_{N-1} \end{bmatrix}=\begin{bmatrix}y_{0}\\ 0\\ 0\\ \vdots\\ 0\\ y_{N} \end{bmatrix} \]を解けば良いわけですね!

というわけであとは行列をいかに解くか? という問題に帰着します.

三重対角行列スペシャル

上の行列を三重対角行列と呼びます. これには，格別に簡単な解法があります. まず, 行列の要素を配列データに格納しましょう. C++の配列は std::vector を用います. 

#include <vector>
...中略...
std::vector<double> A;

と定めておくと, Aに多数の double (浮動小数点数)の値を格納できます．

A.resize(100);

で, 100個の値を格納できるようになります. それぞれの値は, 添字 [i]で区別します. 具体的には

A[0], A[13]

は, 最初の値, 次は14番目の値です. 最初が0なので1ずれることに注意してください. 例えば

A[100]

は101番目の値です. 100個しかないのに101番目にアクセスすると, 運が良ければコンピュータがエラーして停止します. 大抵の場合, コンピュータはエラーせずに計算が終了し, 最終結果として誤った結果が得られます. で，それを論文に投稿してしまい，「こいつ，嘘つきで賞」を貰えます. こういうのを, バグといいます. 

嘘つくつもりはなかった, という言い訳が通るほど世間は甘くありません．間違える奴＝嘘つき→雑用係です．

初心者のうちは, これが恐ろしいのですが, A[i] の代わりに

A.at(i)

を用いて書くと, 範囲内であるかどうかをいちいちチェックしながら計算しますので，バグに遭遇するとコンピュータがエラーして停止します. コンピュータがエラーを出さなければ最終結果は正しいでしょう．ただし，毎回チェックしながら計算を行うので，計算速度はかなーり遅くなります. 

なお, 配列のサイズ(上の例では100)は

A.size()

で取得できますので, わざわざ別の変数に記録したりする必要はありません. では，コード例です：

...前略...
class BC
{
  public:
    class myMat
    {
      public:
        double A;
        double B;
        double C;
    };
    std::vector<double> myVEC;
    std::vector<myMat> myMAT;
    BC(int N)
    {
      ...中略...
      myMAT.resize(N+1);
      myVEC.resize(N+1);
      ...中略...
    };
    ...中略...
};

さてこのBCクラスでは，内部でmyMatクラスを定義しています．こういうのを内部クラスと言います. で，この部分では例の通りにdoubleクラスの配列を用意していますが, この部分ではmyMatの配列を用意しています. このように, std::vector は本当はクラスのオブジェクト(C++ではインスタンスと呼ぶことも多い)の配列を用意するものなのです. で，BCのコンストラクタのこのあたりでは，配列数を設定しています. その配列数は, BCのコンストラクタの引数Nになっています.

ところで, 変数の名前がmyMATとか微妙に長いです．これは，実戦のプログラム開発では，「検索・置換」が強力な武器になるためです．変数名がAとかCのように短いと，Aを含む文字列は無数にありますので，Aを検索しても引っかかりまくって使い物になりません．変数名が多少長いと，例えばmyMAT程度であると，検索によって，その変数がどこで使われているかを簡単に調べることができます．だからといって myFuckingAssholeVariableNeverUsedBefore とかに長い名前の変数は，入力が大変ですのでやめときましょう．

ええっとクラスの初期化において引数が必要であるということは？ BCクラスのインスタンスを作成するときには(要するにmain()では)次のように書かなければならない，ということを意味します:

def

さて配列の要素がクラスですので，そのメンバー変数であるA, B, Cにアクセスするには, 例えば

myMAT[12].A

というふうにアクセスします. そこで, 次のように定義しましょう:

for(int i=1;i<N;i++)
{
  if (i>1) myMAT[i].A=-1.;
  myMAT[i].B=2.-C*dx*dx;
  if (i<N-1) myMAT[i].C=-1.;
  if (i==1)
    myVEC[i]=y_min;
  else if (i==(N-1))
    myVEC[i]=y_max;
  else
    myVEC[i]=0.;
}

すると, 次のように定義できたことになります: \[ \begin{bmatrix}b_{1} & c_{1}\\ a_{2} & b_{2} & c_{2}\\ & a_{3} & b_{3} & c_{3}\\ & & \ddots & \ddots & \ddots\\ & & & a_{N-2} & b_{N-2} & c_{N-2}\\ & & & & a_{N-1} & b_{N-1} \end{bmatrix}\begin{bmatrix}y_{1}\\ y_{2}\\ y_{3}\\ \vdots\\ y_{N-2}\\ y_{N-1} \end{bmatrix}=\begin{bmatrix}v_{1}\\ v_{2}\\ v_{3}\\ \vdots\\ v_{N-2}\\ v_{N-1} \end{bmatrix} \]

1行目の定数倍を2行目に加えて, \(a_2\)の要素を0にできます.　もちろんそのとき\(b_2, v_2\)も変わってしまいますが, 計算機ですので計算するようにプログラムできます．これを繰り返すと，行列は次のようになります:\[ \begin{bmatrix}b_{1} & c_{1}\\ & b_{2} & c_{2}\\ & & b_{3} & c_{3}\\ & & & \ddots & \ddots\\ & & & & b_{N-2} & c_{N-2}\\ & & & & & b_{N-1} \end{bmatrix}\begin{bmatrix}y_{1}\\ y_{2}\\ y_{3}\\ \vdots\\ y_{N-2}\\ y_{N-1} \end{bmatrix}=\begin{bmatrix}v_{1}\\ v_{2}\\ v_{3}\\ \vdots\\ v_{N-2}\\ v_{N-1} \end{bmatrix} \]

もちろん各変数の値は，最初とは異なっています. 次に, \(N-1\)行目の定数倍を\(N-2\)行目に加えて, \(c_{N-2}\)の要素を0にできます. \[ \begin{bmatrix}b_{1}\\ & b_{2}\\ & & b_{3}\\ & & & \ddots\\ & & & & b_{N-2}\\ & & & & & b_{N-1} \end{bmatrix}\begin{bmatrix}y_{1}\\ y_{2}\\ y_{3}\\ \vdots\\ y_{N-2}\\ y_{N-1} \end{bmatrix}=\begin{bmatrix}v_{1}\\ v_{2}\\ v_{3}\\ \vdots\\ v_{N-2}\\ v_{N-1} \end{bmatrix} \] 

これなら, すぐに\(y_i\)の値を求めることができますね! なお，最後のステップは, わざわざ行列を計算しなくても,

mySOL[N-1]=myVEC[N-1]/myMAT[N-1].B;
for(int i=N-2;i>0;i--) //i=N-2から始め，i=1まで上に進むループ
  mySOL[i]=(myVEC[i]-mySOL[i+1]*myMAT[i].C)/myMAT[i].B;

で十分ですよね．

で，これを計算すると, こんな感じになります(\(C=200\)としてみた):
my3

ここでは, \(N=10, 20, 40\)の結果を示しています. 初期値問題の時よりも, 解の収束が速い感じがしますね．

なぜでしょうか? 差分近似の精度をTaylor展開を用いて確かめてみましょう. つまり, \(y_{i+1}\)などを\(x_i\)の周りで展開し, どの程度の誤差を無視しているのか, 考えてみれば良いです.

一応，答えも出しておきましょう:

...中略...
for(int i=2;i<N;i++)
{
  double rate=myMAT[i].A/myMAT[i-1].B;
  myMAT[i].A+=-rate*myMAT[i-1].B;// 0になるだけなので無駄です
  myMAT[i].B+=-rate*myMAT[i-1].C;
  myVEC[i]+=-rate*myVEC[i-1];
}
mySOL[N-1]=myVEC[N-1]/myMAT[N-1].B;
for(int i=N-2;i>0;i--)
  mySOL[i]=(myVEC[i]-mySOL[i+1]*myMAT[i].C)/myMAT[i].B;
mySOL[0]=y_min;
mySOL[N]=y_max;
...中略...
//データを保存する部分
for(int i=0;i<=N;i++) ofs << dx*i << " " << mySOL[i] << std::endl;

途中の double rate のような，そこでしか使わない変数を中間変数と呼びます. C++では, 中間変数は，使いたくなったら，その場で定義して利用します．Fortranなどのように「プログラムの最初に宣言しておく」のは禁忌です．プログラムの可読性が下がるからです．C++コンパイラは，中間変数であると理解したものは，高速なメモリー領域に配置したり, 計算式に組み込んでしまったりして, 計算を高速にします．プログラムの最初に宣言すると，各所で利用されるグローバル変数であると理解し，低速なメモリー領域に配置する場合があります．

Jacobi法

連立方程式\[\text{A}\boldsymbol{x}=\boldsymbol{b}\]の解法には, 逐次近似法というものがあります. 行列\(\text{A}\)の対角成分を\(\text{D}\)とおき, \[\text{B}=\text{A}-\text{D}\]と定めると, 連立方程式は\[\text{D}\boldsymbol{x}=\boldsymbol{b}-\text{B}\boldsymbol{x}\]となります. Jacobi法では, 気分で定めた初期ベクトル\(\boldsymbol{x}^{(0)}\)から, 漸化式\[\text{D}\boldsymbol{x}^{(i)}=\boldsymbol{b}-\text{B}\boldsymbol{x}^{(i-1)}\]によってベクトル列\(\boldsymbol{x}^{(i)}\)を定めると, あらあら不思議，かなり多くの場合に解に収束するというものです.

Jacobi法の物理的イメージは熱方程式のセクションで説明することにして, 実際にやってみましょう．

上と同じように，3列だけ保存するのが正解ですが，行列演算の練習を兼ねて，三重対角行列であることを忘れて，密行列として計算プログラムを組んでみます．

例えば\(N=40000\)のように営業サイズの分割数になると, 非常に大量のメモリーを消費しますので，気を付けてください．具体的には，double 変数は1粒で8byte消費するので, 本当に必要なメモリーサイズは, 行列が3列, 解が1列, ベクトル1列で\(5N\)個あるので,  \[\frac{8\times5N}{10^9}\approx 0.0016 \text{GB}\]であるのに対し，密行列には\(N^2\)個の要素があるので, プログラムが消費するメモリーは\[\frac{8\times(N^2+2N)}{10^9}\approx 12 \text{GB}\]となり，大抵のデスクトップPCを暴走させるのに十分なサイズになってきます．つまり，うっかり密行列で解こうとした人には, 大変難しい問題に見えるわけです．正しい取扱いを行なった同業者には，なんの困難もない，ということになります．

消費メモリーサイズが計算機のメモリーサイズを超えると，計算速度は恐ろしく低下します．10秒で済むものが1週間かかったりします．コンピュータは「メモリーが足りませんよ」と連絡してくることはありません．あなたが，確認する必要があります．ここでは，消費メモリーサイズを確認する方法も覚えましょう．

線形代数ライブラリを使う

密行列を取り扱うのは, Eigenライブラリーを用いるのが簡単です（Eigenは疎行列も扱えます．念のため）．

Eigenはとにかく計算が段違いに早いです．まずはここに従ってインストールを済ませましょう.

\(N\times N\)の行列を準備し，係数行列にするには次の通り：

#include <Eigen/Dense>
class BIG_MEM
{
private:
  Eigen::MatrixXd myMAT;
  Eigen::VectorXd myVEC;
  ...中略...
public:
  BIG_MEM(int N)
  {
    ...中略...
    myMAT.resize(N-1,N-1);
    myMAT.setZero();
    myVEC.resize(N-1);
    myVEC.setZero();
    for(int i=0;i<N-1;i++)
    {
      if (i>0) myMAT(i,i-1)=-1.;
      myMAT(i,i)=2.-200*dx*dx;
      if (i<N-2) myMAT(i,i+1)=-1.;
      if (i==0) myVEC(i)=y_min;
      else if (i==(N-2)) myVEC[i]=y_max;
    }
    ...後略...

private: というのが初めて出てきています．これは，クラス外部からアクセスを禁止する変数です．「禁止しなくってもよいのではないか？」といえばそうなのですが，禁止しておかなければ，よくわかっていない人が色々いじくり回して，動かないよ？と文句を言ってくる可能性が出てきます．よって，他人とプログラムを共有する場合，外から見える必然性がない変数はprivate:にしておくほうが良いです．赤色部分がそれぞれ行列とベクトルの宣言です．Xdは浮動小数であることを意味しています．次元は青色のところで設定しています．setZero()をしないと，要素が0になりませんので注意してください. 使い方には特に注意はないです．数学とは異なり, 添字が0からですので注意. 添字は[i]と書いても(i)と書いても同じですが, (i)は範囲チェックを行いますので, 計算は安全ですが遅いです．

Eigen行列には色々なオペレータが定義されているので，便利に使ってください．例えば inverse() は逆行列を与えます．ので，

void INV()
{
  auto SOL=myMAT.inverse()*myVEC;
  std::ofstream ofs("my.data");
  ofs << dx*0 << " " << y_min << std::endl;
  for(int i=0;i<N-1;i++)
    ofs << dx*(i+1) << " " << SOL[i] << std::endl;
  ofs << dx*N << " " << y_max << std::endl;
}

この auto というのは，オブジェクト指向につきもののやつで，右辺を見たら計算結果がベクトルであることはクラスの定義から明らかなんですね．だから「コンパイラーさん，勝手に推定して適当な型宣言してくれろ」ということができるわけです．これで, \[ \text{A}^{-1}\boldsymbol{b}\]を計算することができます．ってゆうか，それ答えなんで, ちゃんと出ます:

re

・・・まあしかしなんというか，\(N\)が1000とか10000になると, 逆行列の計算はまぢ遅いので，あまり実用しないほうが良いです(逆に言えば，それ以下なら，これが一番早い)．実際，N=5000で走らせてみると，

mem

このように, 400MBものメモリーを利用し，あなたのノートPCの電池を結構食いつぶしていることがわかります．消費メモリーを調べる別な方法は，アクティビティモニタです：

acta

あー，それから，今はJacobi法を書くのでした．こういうのは，まず作戦が必要です．

  1. 行列を設定する Ready() 関数
  2. \(\boldsymbol{x}_{i-1}\)から\(\boldsymbol{x}_{i}\)を求め, \[||\boldsymbol{x}_{i-1}-\boldsymbol{x}_{i}||\]を出力する関数Fire()
  3. 最後に結果を保存するCease()関数

と言った感じですね．すると，main()は

BIG_MEM myObj(number_N);
double err_max=1e-7;
myObj.Ready();
while (myObj.Fire()>err_max);
myObj.Cease();

こんな感じになりますね．ずれがerr_max以下になるまで頑張る，ってわけです．Ready()はBIG_MEMクラスのコンストラクタに書いても良いし, Cease()はデストラクタに書いてもOK. あとはパーツを書けば良い．

プログラムを書くとき，このように，大枠から定めていくのはコツです．パーツから積み上げると，だいたい発散して出来上がりません．

まず行列設定です:

Eigen::MatrixXd MAT;
Eigen::VectorXd VEC;
Eigen::VectorXd SOL;
std::ofstream ofs;    //最初にはファイル名を設定していないことに注意
void Ready()
{
  MAT.resize(N-1,N-1); MAT.setZero();
  VEC.resize(N-1);     VEC.setZero();
  SOL.resize(N-1);     SOL.setZero();
  double dia=2.-C*dx*dx;
  for(int i=0;i<N-1;i++)
  {
    if (i>0) MAT(i,i-1)=-1./dia;
    if (i<N-2) MAT(i,i+1)=-1./dia;
    if (i==0) VEC(i)=y_min/dia;
    else if (i==(N-2)) VEC[i]=y_max/dia;
  }
  ofs.open("my.data");　　//ここで my.data を開く．どんどん追加できるようにする
}

対角要素\(\text{D}\)で初めから割って行列を作成しています．すると\(\text{D}=\text{E}\)だから計算が楽になるわけですね．

しかしながら dia の計算を見ると,\[C\Delta x^2 = 2\]で，なんともやばそうであることがわかります．想像できるように, Jacobi法は必ずしも収束するとは限らないわけです.

関数Fire()は，そうすると

double Fire()
{
  auto NEW=VEC-MAT*SOL;
  double err=(NEW-SOL).squaredNorm();
  SOL=NEW;
  std::cout << "ERR= " << err << std::endl;
  return err;
}

ここは，関数がdouble数値を返すと定義しています．返却する値は，ここの err です. このあたりはちょいと難解かもしれませんね．NEWとSOLがベクトルですので，この引き算の結果はベクトルです．すると, その結果には，ベクトルクラスに定義されている関数を適用することができる！だから, squaredNorm() つまりえっと成分の2乗の和が err に設定されます．

今回のmain()プログラムは, 5000回ごとにデータを書き出せるように, 以下のものを用いました:

BIG_MEM myObj(number_N);
double err_max=1e-7;
myObj.Ready();
for(int i=1;i<50000;i++) //myObj.Fire()内の i 変数とはダブらないので気にしない
{
   myObj.Fire();
   if (i%2000==0) myObj.Write();
}

myObj.Write()は, iが2000で割り切れる場合に呼び出されるのですが, ofsにデータを追加する関数です. このために, ファイル ofs は初期にはファイル名が無設定であり，Ready()でファイル名を設定, その後ファイルにどんどん追加できる，という作戦を取っています．

void Write()
{
  ofs << dx*0 << " " << y_min << std::endl;
  for(int i=0;i<N-1;i++) ofs << dx*(i+1) << " " << SOL[i] << std::endl;
  ofs << dx*N << " " << y_max << std::endl << std::endl;
}

ここで std::endl が2連です．この場合，データを書いて改行して改行するので，空行が1つ開きます．それによって，異なるデータを区別するわけです．実は, graph プログラムは, 空行があると新しい別の曲線である，と理解します．よって，次のような結果が得られます．なお\(C=2\)のケースです．それ以上の\(C\)はちょっと難しいか.

my

5万回も繰り返したのに，なんやゆっくり変化を続けています．これでも，最後には \[||\boldsymbol{x}_{i-1}-\boldsymbol{x}_{i}|<10^{-10}|\]なんですけどね・・・

共役勾配法

逐次近似は収束するまで回数が掛かる場合があります. 場合によっては, \(N\times N\)行列なら最大\(N\)回で必ず収束するという共役勾配法もあります．ま，でも\(N\)が馬鹿でかいと, なかなか\(N\)回とは行かないんですけど, 行列を解くのが仕事になったら, 考慮しても面白いでしょう.

Gaussの消去法

三重対角行列のように, 一部分にしか非ゼロ要素がない場合には利用しないのですが, 中学生と同じように

  • 連立方程式を解くには, 一個づつ変数を消去していけば良い

を実行しても良いです. これをGaussの消去法と言います. 自分で書いても良いですが，もちろん，Eigenにはプロが書いたものが用意されていますので，それより良い性能を出すのは無理でしょう．

Eigenを使って，どのように解くのかは，「Eigen 連立方程式」とかでGoogleに聞いたら良いです．ここにも少しだけあります．でも，日本語で聞いてもぱっとしないことがありますので，英語で聞きましょう．日本人1億人よりも，世界人70億に聞いたほうが，情報量が70倍です．