./note007.html

ToC

6.1. ルンゲ・クッタ

という形の初期値問題を考える。

ルンゲクッタ法とは、非常に一般的には以下のような形に書ける方法である。

自然数

を段数（number of stages）という。

はパラメータであるが、

と

は普通

となるようにとる。これは、一般にそうでないような公式は不可能ではないがあまりいいことがない（精度がよくならない）からである。

と、こう、式に書いてしまうとすぐにはわからないが、例えばの場合に書き下してみると

と、まあ、こんな感じになる。こちらを良く見ればすぐにわかるように、

がすべて0ならば、から順に計算していくこと

ができる。つまり、「陽的」公式になっている。

が0のときは、各についての式にだけが入ってくる。これを半陰的 (semi-implicit) 公式という。この場合には、まずについての方程式をといて、次にについてのを解いて、、、と順番に計算出来る。
上のような制約が全くない時は、「陰的」公式ということになる。このときは、すべてのに対する（一般には非線形な）方程式を一度に解く必要がある。

今日はとりあえず陽的な公式だけを考える。形式的には、もっとも簡単な前進オイラー法もルンゲ・クッタ公式の一つということになるが、まあ、普通もっとも簡単なRK法というと、2次の公式

である。この方法がどのように働くかについては、図的な説明というものが可能である。元の点からまずオイラー法と同様に接線を引く。が、これを次の時刻まで延ばすのではなく、ステップの半分のところで止める。で、ここでもう一回微分方程式の右辺を評価する。ここでの導関数の値を使って、もとのところ

から直線を引くわけである。

Figure 4: 2次ルンゲクッタ法

普通は、RK法というといわゆる古典的 RK 公式

のことをさす。これは、いろいろ良い性質をもつ。例えば

が以外すべて0なので、右辺の計算が楽である。
次数が4次であり、4段陽的公式で到達可能な最高次数を達成している
係数が簡単な有理数なので、プログラムしやすい。また丸め誤差を小さくできる。

この公式が 4 次であることを示すのは、それほど簡単ではない。腕力に自信があるひとは挑戦してみて欲しい。

4次よりも高精度な公式というのも非常によく研究されていて無数にいろんなものが知られている。ただし、4次よりも高い次数では、次数よりも必要な段数が大きくなる。この2つの関係が一般にどうなるかは未解決の問題である。というようなこともあるので、よほど変わったことをしたいというのでない限り、専門家が作った公式を使うのが無難である。

Dormand と Prince は、段数／次数のさまざまな組合せについて、離散化誤差を非常に小さくした公式を求めている。これらは、

  http://www.unige.ch/math/folks/hairer/software.html

から入手できる (Fotran と C がある)。

普通の（というのがどういう意味かは今のところ明らかではないが）常微分方程式を手軽に解くにはこれらのRK公式を使えば十分である。

実用的には、こちらが欲しい計算精度をどうやって達成するかという問題がある。この問題を解決するためには、誤差を推定しながら積分の刻み幅を自動的に調節する必要がある。これについては次回にちょっと触れる。

6.2. ルンゲ・クッタ以外の方法

さて、「普通は」RKで十分というのは、もちろん時と場合によっては違う方法を使うべきであるという意味である。なぜそういう場合があるかというのは、以下の理由による。

RK法は（次数のいかんにかかわらず）ある問題を解くのにもっとも計算量が少ない方法ではないことが多い。
特に、問題の種類によっては、特別な RK 法を使うと非常にうまくいくことがある。
問題の種類によっては、普通のRK法ではどんなに計算時間をかけてもまともな答がでない。

まあ、このへんは詳しい話をやりだすときりがないので、簡単に。以下、上の３つを順番に扱う。

6.3. 線形多段階法

RK 法は、「一段階」法である。これはどういう意味かというと、微分方程式

の

での近似解

があったとして、次の

の値 <$x_{i+1} = x_i + \Delta x_i$> での近似解

を求めるのに、

と微分方程式そのものだけで十分であるということである。中間のややこしい値を計算する必要はあるが、それは RK 法のほうが勝手にやることで、使う方が入力を与える必要はない。

これに対して、一段階ではない方法、つまり多段階法というのは、を計算するのに、以外の情報、例えばや、そこで計算した導関数の値、さらにもっと昔の情報を使うやりかたのことである。これは、プログラムとしてはもちろんRK法に比べれば面倒になる。昔の値をとっておかないといけないし、また、一番最初に計算を始める時にどうするかという問題もあるからである。

さらに、「一段階でない」というだけなので、可能な計算法にあまりに多様な可能性がある。一例として、以下のような方法が考えられる。

これは、陽的中点公式といわれるもので、式としては偏微分方程式の解法の時に空間微分の1次の項に使った中心差分と同じ形になっている。従って、理屈としては2次の精度をもった公式になっている。

この方法は、しかし、立派な名前までついているのにもかかわらず、実は使えない公式である。

どのような問題があるのかを示すために、以下の線形方程式

に陽的中点公式を適用したらなにが起きるかを考えてみる。刻みは固定の

とすると、離散化したものは

になる。これは線形差分方程式なので、前にもやったように固有値を調べればいい。今

ならどちらかは必ず絶対値が 1 より大きい。

ちょっと式を見ればわかるように、絶対値が 1 より大きい固有値は、つまりなら負である。したがって、必ず振動的に発散することになる。

なお、上のような、安定な線形微分方程式について振舞いを調べるというのが、常微分方程式の安定性解析の基本になる。非線形な方程式では違うとかいろんなことがあるかもしれないわけだが、まあ、少なくとも線形安定でないと話にならないし、それ以上のことは一般論としていうのは難しいからである。

6.3.1. アダムス法

さて、安定でちゃんと使える線形多段階法はじゃあどんなものかというわけだが、これも無限にいろんな作り方がある。そのへんの詳しい話はそういう本に譲ることにして、ここではもっとも広く使われているアダムス法について説明する。

原理は、いくつかのステップでの導関数（微分方程式の右辺）の値を憶えておいて、それを通る補間多項式を作り、それを積分して解を求めようというものである。

Figure 5:

図に概念を示す。ここでは、ラグランジュ補間（ニュートン補間）をして多項式を作る。で、その作った多項式を積分する。例えば、点からまで積分するのに、点からまでの関数値を使うとすれば、次の多項式で

簡単な例として、の場合を考えてみよう。この時、補間多項式は一次であって

となって、これを積分すれば、結局

となる。

一般に、アダムス法では任意段数の公式が構成でき、その次数は段数に等しいことがわかっている。これは、ルンゲ・クッタなどに比べればはるかによい性質をもっているということでもある。

6.3.2. 出発公式

アダムス法はいくらでも高次の公式が作れ、計算量もあまり多くないということがわかっているが、必ずしも広く使われているというわけでもない。その理由はいろいろあるが、一つは、

「どうやって計算を始めるべきかよくわからない」

ということである。つまり、初期値問題としてはもちろんにおけるしか知らないのに、多段階法ではその前の時刻での解が必要になるわけである。これに対する対応策はいくつかあるが、時間刻み一定の場合には、基本的にはルンゲ・クッタなどの別な方法で解を求めておくというやりかたが普通である。

というわけで結局プログラムを書く手間が2倍以上になるというのが、多段階法の実用上の問題である。

6.3.3. 陰的アダムス法

さて、前に述べた公式では、補間多項式を陽的に求めた。すなわち、時刻とそれより以前の値だけを使っていた。これに対し、陰的な補間多項式、つまりでの関数の値を使った公式というものも考えられる。刻みが定数であるとすれば、を決めれば前と同様に

の形に書けることになる。また手をぬいて

の場合を考えれば、これは単に台形公式

となる。

陰的公式の場合、例によってどうやって代数方程式を解くかが問題になる。通常の方法は、

初期値として同じ次数の陽的アダムス法の解を持ちいる。
そのあとは直接代入法で反復する。

というものである。ただし、特に線形多段階法の場合、反復を繰り返さないことが多い。反復回数を1回とか2回に固定してしまうのである。なぜそれでいいか、また、そもそもなぜ陰的公式を使うかというあたりはレポート課題ということで。

なお、このやり方を、予測子・修正子法と呼ぶ。線形多段階法はほとんどこの形で使われるため、線形多段階法のことをさして予測子・修正子法と呼ぶ人もいる。

6.4. 構造体とクラス

ここからはちょっと C++ 言語の話題というか、今回以降の話で便利な機能について。

偏微分方程式とか、あるいは連立常微分方程式の数値解法をプログラムするのに、普通は配列を使う。しかし、これは割合に面倒だし、いろいろ妙な間違いをする可能性も高い。例えば、にを足すのは、数式としてはで済むのに、例えば C では

 for(i=0;i<n;i++) x[i] += dx[i];

Fortran なら

      do i = 1, n
          x(i) = x(i) + dx(i)
      enddo

という具合で、数式なら4文字で済むところをその何倍も書かないといけない。これは、 C にしても Fortran にしても、配列（あるいはベクトルとか行列）といった、数学では基本的な要素であるものに対する演算を直接に表現する記法を持っていないからである。

まあ、Fortran 95 とかそういったものを使うとそういう記法があるという話もないわけではないが、あまり普及していない。普及していない理由はいろいろあるが、その一つは C++ を使えば同様なことが実現できるからである。

C++ 言語自体は、配列に対する演算というものを用意しているわけではない。しかし、C++ では、プログラムの中で新しい「型」（実数型とか整数型というのと同じ意味での）を定義して、さらにそれに対する演算を定義することができる。これでベクトル型とかいったものを自分の使いやすいように定義すればいいことになる。

例えば、非常に基本的なベクトル型の定義は以下のようなものになる。ここでは加算と入出力くらいしか定義していないが、他の必要な演算も同様に定義できる。一応使いそうな演算を定義したものが

 http://grape.astron.s.u-tokyo.ac.jp/~makino/pcphysics/programs/vector.h

にあるので、実際にプログラムを作る時にはこれを使ってもよい。牧野が書いたプログラムは信用できないという向きは自分で書くこと。

 #ifndef  STARLAB_VECTOR_H
 #  define  STARLAB_VECTOR_H
 #define PR(x)  cerr << #x << " = " << x << " "
 #define PRC(x) cerr << #x << " = " << x << ",  "
 #define PRL(x) cerr << #x << " = " << x << "\n"
 
 #ifndef real
 # define real double
 #endif
 /*-----------------------------------------------------------------------------
  *  mvector  --  a class for 3-dimensional vectors
  *-----------------------------------------------------------------------------
  */
 class mvector
     {
     private:
         real element[VLEN];
 
     public:
 //Default: initialize to zero.
         mvector(real c = 0)
     {for (int i = 0; i<VLEN;i++)element[i] = c;}
 
 //  []: the return type is declared as a reference (&), so that it can be used
 //  on the left-hand side of an asignment, as well as on the right-hand side,
 //  i.e.  v[1] = 3.14  and  x = v[2]  are both allowed and work as expected.
 
         real & operator [] (int i)       {return element[i];}
 
 inline void print() {for(int i = 0; i<VLEN;i++)
      cout << element[i] << " ";
      cout << "\n";}
 
 //Unary -
 
         const mvector operator - ()
     {mvector v;for(int i = 0; i<VLEN;i++)v.element[i] = - element[i];
           return v;}
 
         friend const mvector operator + (const mvector &, const mvector & );
         friend const mvector operator - (const mvector &, const mvector & );
         friend mvector operator * (real, const mvector & );
         friend real operator * (const mvector &, const mvector & ); 
         friend mvector operator * (const mvector &, real);
         friend mvector operator / (const mvector &, real);
 
 //Mvector +=, -=, *=, /=
 
         mvector& operator += (const mvector& b)
     {for(int i = 0; i<VLEN;i++)element[i] += b.element[i];       
      return *this;}
 
 mvector& operator -= (const mvector& b)
     {for(int i = 0; i<VLEN;i++)element[i] -= b.element[i];
      return *this;}
 
 mvector& operator *= (const real b)
     {for(int i = 0; i<VLEN;i++)element[i] *= b;
      return *this;}
 
 mvector& operator /= (const real b)
     {for(int i = 0; i<VLEN;i++)element[i] /= b;
      return *this;}
 
 //      Input / Output
 
         friend ostream & operator << (ostream & , const mvector & );
 
 friend istream & operator >> (istream & , mvector & );
     };
 
 inline  ostream & operator << (ostream & s, const mvector & v)
 {
     for(int i = 0; i<VLEN;i++)s << v.element[i] << "  " ;
     return s;
 }
 inline  istream & operator >> (istream & s, mvector & v)
 {for(int i = 0; i<VLEN;i++)s >> v.element[i]; return s;}
 
 inline  const mvector operator + (const mvector &v1, const mvector & v2)
 {
     mvector v3;
     for(int i = 0; i<VLEN;i++)v3.element[i] = v1.element[i]+ v2.element[i];
     return v3;
 }
 inline  const mvector operator - (const mvector &v1, const mvector & v2)
 {
     mvector v3;
     for(int i = 0; i<VLEN;i++)v3.element[i] = v1.element[i]- v2.element[i];
     return v3;
 }
 
 inline  mvector operator * (real b, const mvector & v)
 {
     mvector v3;
     for(int i = 0; i<VLEN;i++)v3.element[i] = b* v.element[i];
     return v3;
 }
 inline  mvector operator * (const mvector & v, real b)
 {
     mvector v3;
     for(int i = 0; i<VLEN;i++)v3.element[i] = b* v.element[i];
     return v3;
 }
 
 inline  real operator * (const mvector & v1, const mvector & v2)
 {
     real x = 0;
     for(int i = 0; i<VLEN;i++) x+= v1.element[i]* v2.element[i];
     return x;
 }
 
 
 inline  mvector operator / (const mvector & v, real b)
 {
     mvector v3;
     for(int i = 0; i<VLEN;i++)v3.element[i] =  v.element[i]/b;
     return v3;
 }
 
 typedef const mvector (*vfunc_ptr)(const mvector &); 
 typedef const mvector (vfunc)(const mvector &); 
 
 #endif
 
 //=======================================================================//
 //  +---------------+        _\|/_        +------------------------------\\
 //  |  the end of:  |         /|\         |  inc/vector.h
 //  +---------------+                     +------------------------------//
 //========================= STARLAB =====================================\\

実際にこれを使うには、

 const int VLEN = 2;
 #include  "vector.h"

 double k;
 mvector dxdt(mvector & x)
 {
     mvector d;
     d[0] = x[1];
     d[1] = -k*x[0];
     return d;
 }

    .....
    double h;
    mvector kx1;
    kx1 = dxdt(x)*h;
    ......

というような具合に、ベクトルの大きさを指定する（これはここでは固定である）。もちろん、可変にするとかいろんなことができるが、繁雑になるのでここでは固定サイズにする。

ここで、注意してほしいのは、 + という演算子や [] という、これも「演算子」が、新しく定義されていることである。これらは、もちろん

 + なら実数や整数に対してすでに定義されているが、ここではベクトル

同士の演算に対して新しい意味を持つように拡張されていることになる。 [] も同様で、配列の他にここで定義したベクトル型について新しい意味を持つようになったわけである（といっても、こちらは普通の配列に対してとまったく同じように働くが）。

これは、偉そうにいうと C++ の演算子多重定義 (overload)という強力な機能である。まあ、落ち着いて考えてみると、この機能の実現はそんなに難しいわけではなくて、例えば + という演算子が出てきたところで、それが適用されているデータ型を見て、そのデータ型に対して定義されている演算を呼ぶようにするというだけのことである。これは、コンパイラにそういう機能を付け加えるだけで実現できる。

C++ の場合には、このような、ある意味での機能拡張は、「クラス」というものを使って実現されている。これはもともとは「オブジェクト指向」とかそういった難しい機能を実現するためのものであるが、ここではとりあえずあんまりそういうことは考えないでベクトル型を使うことにする。

なお、名前が素直に vector ではなく mvector (mathematical vector のつもり)になっているのは、 C++ の新しい標準では vector という語に別の意味を与えているからである。

上のコードの意味を簡単に説明しておく。

 class mvector は mvector という名前のクラスの定義を始めますという

宣言である。その後に \{\} で囲って中身を書く。

中身には private: と書いてから書くものと public: と書いてから書くものがあって、 private: のほうは、「普通には外から見えない」ものを定義するのに使う。典型的なのはこの例のように、あるクラスの「オブジェクト」が持つデータを宣言することである。データの宣言は普通の変数宣言と同じ。

 public:  のほうでは、主に「メンバー関数」というものを宣言するこ

とになる。ここで、自分の名前と同じ名前の関数は特別な意味を持っていて、 (constructor)、このクラスの変数がプログラムの中で使われる時に呼ばれるものである。上の例では、デフォルトでは vector a; みたいに書くと

 element の各要素に 0 が代入される。0以外にしたければ  vector

a = vector(1); とか書く。

その次の double \& operator [] (int i)\{ ...\} というのははっきりいってなんだかわからないが、 C++ ではいくつかの演算記号や上の配列の記号等に、クラス毎に別の意味を与えることができる。これが「演算子の多重定義」や「オーバーロード」と呼ばれる機能である。

というとなんだか恐ろしげに聞こえるが、考えてみると ``+'' だって実数と整数では全く違うことをするわけで、別のデータ型を定義したらそれの ``+'' はまた違うことをして欲しいというのは当然であろう、というより、ここで我々がしたいのはそもそもそういうことであった。上の例では、結局プログラムの中で

 mvector a; ... a[i] ... という感じに書くと、  a[i] が  a
 の要素である  element の  i 番目の要素、つまり

element[i] と全く同じことになる。

なお、 C 言語では構造体 struct というものがあり、例えば

 struct mvector{
     double element[NDIM];
 };

といったふうに使う。ここで struct mvector a; のように変数を宣言すると、 a.element[i] といった形で ``.'' を使って要素にアクセスすることになる。これはちょっと見苦しい。また、C には演算子のオーバーロードの機能はないので、折角構造体を作っても、 ``+'' で足すというわけにはいかない。これはある意味みかけだけの問題ではあるが、見てわかりやすいということはプログラムを書いたりデバッグしたりする時には割合大事なことである。

その後の friend const mvector operator + ... が、実際に mvector 同士の加算を宣言する。但し、ここでは形だけを宣言して、実体は下のほうに書く。但し、その次の += のように全部書いてもかまわない。

 mvector a, c; ... a += c; と書くと、`` +='' の中で いきなり

element と書かれているものは a の要素、 b.element になっているのは c の要素になって、結局 a の各要素に c の各要素が足されるという、普通に期待したい動作が実現されることになる。なお、このような、宣言してないクラス変数を使う関数を「メンバー関数」という。

`` +'' も事情は同様だが、ここでは 2 つ引数を取る普通の関数として定義されている。 friend を付けると、 private: と宣言したものもこの関数の中では使えるようになるので、ここでは friend を付けている。

`` <<'', `` >>'' はメンバ関数として実現されている。

6.5. 練習

6.5.1. 練習1

一次元調和振動子

について、初期条件

からの、

までの数値解を、古典的ルンゲ・クッタ公式を使って求めよ。刻み幅をいろいろ変えて、精度が刻みにどのように依存するか、また、どこまで高い精度が実現可能か調べ、なぜそうなるかを考察せよ。

なお、牧野が書いたプログラムがここの hermonic5a.C にあるので、今回はこれをコピーしてコンパイル・実行するだけでも一応いいことにする。但し、これは古い C++ の規格にそって書かれているので色々変更が必要である。

6.5.2. 練習2

刻み幅が等しい場合について、以下の公式を導け

4次の陽的アダムス公式
4次の陰的アダムス公式
4次の陽的シュテルマー公式

但し、シュテルマー公式とは、2階の微分方程式用の公式である。時間の2階微分に、偏微分方程式の時に使ったのと同じ3点を使う中心差分を使う。

6.5.3. 練習3

上で導いた4次の陽的・陰的アダムス公式を使って、1と同じ単振動の方程式について同様な解析を行なえ。出発公式にはルンゲ・クッタを使うか、あるいは厳密解をつかってもよい。

6.5.4. 練習4

陰的線形多段階法について、直接代入法が収束する速さが何で決まるかを考察せよ。

Previous ToC Next