2008. 4. 4.
2014. 2.21. 訂正 delete → delegate
石立 喬

Visual C++ 2008 Express Edition の易しい数値計算(6)

――― RungeKutta法で一元微分方程式を解き、回路の過渡現象を計算する ―――

 微分方程式の代表的な解法である4次のRungeKutta法を用いて、一元微分方程式を数値的に解く。例題は簡単な電気回路の過渡現象とし、波形の表示が必要になるので、最小限の機能を持ったWindowsフォームアプリケーションを使用する。
 プログラムの解説と実証に重点を置き、実行画面の見栄えや使い勝手は重視しない。Windowsアプリケーションの詳細は、すでに「Visual C++ 2005 Express Editionの易しい使い方」で述べてある。

Windowsフォームアプリケーションの始め方
 すでに、「Visual C++ 2005 Express Editionの易しい使い方(2)」で紹介してあるが、あらためて簡単に説明する。
1)統合開発環境のメニューから、「ファイル」→「新規作成」→「プロジェクト」と順に選択する。
2)「新しいプロジェクト」ウインドウが開くので、
  プロジェクトの種類 --------------------------- CLR
  Visual Studioにインストールされたテンプレート --- Windowsフォームアプリケーション
  プロジェクト名 ------------------------------ RungeKutta(入力する、一例)
  ソリューションのディレクトリを作成 ------------- チェックが付いていることを確認
  ソリューション名 ---------------------------- RungeKutta(自動的に入る、一例)
 を入力または確認して、「OK」をクリックする。図1は、「OK」をクリックする直前の「新しいプロジェクト」ウインドウを示す。


図1 「新しいプロジェクト」ウインドウ


3)しばらくして現れた「Form1.h[デザイン]」の画面で、Windowsフォームの内部を右クリックし、「プロパティ」を選択し、「プロパティ」ウインドウを表示させる。
 「配置」欄の「Size」の右に「300,300」とあるのを、とりあえず「500,500」に書き換える(後でコードエディタから修正できる)。
 「表示」欄の「BackColor」をクリックし、右側の選択メニューから「Window」を選択する。
 「表示」欄の「Text」の右に「text1」とあるのを、「RungeKutta」に書き換える。
 図2は、書き換えた「プロパティ」ウインドウを示す。


図2 書き換えた「プロパティ」ウインドウの一部


4)「プロパティ」ウインドウの上部にある「イベント」(稲妻の絵がついている)アイコンをクリックし、「表示」欄の「Paint」をダブルクリックすると、右に「Form1_Paint」が表示され、「コードエディタ」ウインドウが開いて、Form1_Paint()のスケルトンができる。簡単なプログラムであれば、Form1_Paint()に記述するだけで完成する。
5)ビルドには、統合開発環境のメニューから「ビルド」→「RungeKutta(一例)のビルド」をクリックし、実行には、「デバッグ」→「デバッグなしで実行」をクリックする。

微分方程式を解くための4次RungeKutta法
 微分方程式をコンピュータで数値的に解く方法は、古くから知られている。
   
のような形式で記述できる、最も簡単な一元微分方程式があるとき、
   
を繰り返して求める方法を、Euler法と呼び、微分方程式の最も簡単な数値計算法である。ただし、これは、刻み幅hを十分小さくしないと誤差が大きく、hを小さくすると計算時間が長くなる欠点がある。
   
のように、中間点の傾斜を用いたのが改良Euler法で、hの幅をやや大きく選んでも、ある程度の精度が得られる。
 Euler法をさらに進歩させたのが4次Runge-Kutta法(ドイツの数学者、C.RungeとM.W.Kuttaが1901年に発表)で、実用性に優れているので、これを使用する。
 4次RungeKutta法は、
   
で表されるような一元の微分方程式に対して、
   
を繰り返しで、解を求める。K1からk4までを用いるので、4次と呼ばれる。
 より複雑な、二元以上の微分方程式や、右辺の関数に時間tを含むものに対してもRungeKutta法を使用できるが、これらについては、稿を改めて「Visual C++ 2008 Express Editionの易しい数値計算(7)」で紹介する。

一元微分方程式で解ける簡単な回路の過渡現象
◎CR回路
 図3の回路で、コンデンサCが放電している状態(端子間の電圧が0)でスイッチを入れた時の電流i(t)を求める。

図3 CR回路

 図3の回路から、 の関係が成り立つので、その両辺をtで微分すると、
なる微分方程式が得られる。これを変形して として、RungeKutta法で解くことができる。

 t=0においては、コンデンサCの両端の電圧がなく、すべての電圧が抵抗Rにかかるので、i=E/Rが、iの初期条件である。
 これを解析的に解く(解法は省略する)と、 i = [定数]・e-1/(C・R)・tが得られ、t=0でi=E/Rの条件を入れると、i = E/R・e-1/(C・R)・t となる。t=CRでは、i = E/R・1/e となって、初期値の約0.368倍になることが分かっているので、この値を時定数(τ=CR)と呼び、波形の概要を把握するのに使用する。

◎LR回路
 図4の回路で、コイルLに電流が流れてなく、電流変化もない状態(端子間の電圧が0)でスイッチを入れた時の電流i(t)を求める。

図4 LR回路

 図4の回路で、 の関係より、これを変形して として、RungeKutta法で解く。

 t=0において急激に電流を流そうとした時は、コイルLのインピーダンスを無限大とみなすことができ、電流が流れないので、i=0が、iの初期条件である。
 これを解析的に解く(解法は省略する)と、i = [定数]・(1 - e-R/L・t)が得られ、t=∞でi=E/Rの条件を入れると、i = E/R・(1 - e-R/L・t) となる。t=L/Rでは、i = E/R・(1 - 1/e) となって、最終到達値の約0.632倍になることが分かっているので、この値を時定数(τ=L/R)と呼び、波形の概要を把握するのに使用する。

delegateキーワードを使う
 仮引数の並びと型、戻り値の型が同一のメソッドが複数個あるとき、それらを切り替えて使用できれば便利である。これにはdelegateキーワードが使用できる。delegateとは、「代表する」の意味で、実際のメソッドの代わりに別の名称で呼び出すことができる。
 下記は、delegate宣言の一例で、ここで使用するものである。

   private: delegate double Function(double);

 これは、仮引数が一個でdouble型であり、戻り値がdoubleであるメソッドに、Functionと言う型の名称をつけることを意味する。Functionは、参照クラスと同様に、^(カレット)を付けて使用する。
 ここで使用した、Function型のオブジェクトを生成する(インスタンス化する)例を下記に示す。

   Function^ function1=gcnew Function(this,&RungeKutta::Form1::function_CR);

 これは、名前空間RungeKuttaの中にあるクラスForm1のメンバーメソッドfunction_CRを、Function型のfunction1として登録することを意味する。thisは、使用するメソッドが現在属しているオブジェクト、すなわちForm1自身を指し、&RungeKutta::Form1::function_CR は、使用するメソッドへのポインタである。C言語などで用いられてきた、関数へのポインタに比べて、オブジェクトを指定しているので、安全性が高い。

画面上に点を打ってグラフを描く方法について
 Java API では、
   g.drawLine(x,y,x,y); または、g.drawRect(x,y,0,0);
によって、1ピクセルの点を描画することができ、これによってグラフを描くことができる。ただし、あらかじめ、
   g.setColor(Color.red);
などで、色を指定しておく必要がある。
 一方、Visual C++ で使用する .NET Framework の場合には、
   g->DrawLine(pen,x,y,x,y); や g->DrawRectangle(pen,x,y,0,0);
では、何も描いてくれない(.NET Framework 2.0から.NET Framework 3.5になっても改善されていない)。
   g->DrawLine(pen,x,y,x+1,y+1); や g->DrawRectangle(pen,x,y,1,1);
は、一個の点ではなく、二個ないしは四個のピクセルが描画されるので太くなる。
   g->DrawEllipse(pen,x,y,1,1);
は、DrawRectangeよりは、すこし小さい点が描かれるが、場合により形状が変わる欠点がある。

 これらの問題を解決する方法には、下記の2種類があり、すでに「Visual C++ 2005 Express Editionの易しい使い方」などで、多用しているので、参照されたい。
1) 点を打ってグラフを作成することを避け、二点間を線で結ぶ。新しく点を打ちたい場所を点(x,y)とし、一つ前に打った点の場所を点(old_x,old_y)として、
   g->DrawLine(old_x,old_y,x,y);
  を使用する。ただし、変数old_xとold_yを用意し、最初に打つ点だけを特別扱いにしてold_xとold_yに格納するなどの面倒がある。
2) Imageクラスのオブジェクトに1ピクセルずつ順次点を打ってから、それを一括して描画する。すなわち、
   Bitmap^ bmap=gcnew Bitmap(width,height);   //オブジェクトbmapを生成する
   bmap->SetPixel(x,y,color);              //bmapに1ピクセルだけ点を打つ
   ---------------------
   ----- 順次点を打つ ----
   ---------------------
   g->DrawImage(bmap,x,y);              //bmapを画面に描画する
を利用する。ただし、描画の過程を見ることができず、一定期間待たされてから、瞬間的に画像が現れる欠点がある。

プログラム
1)画面上の表示位置、表示寸法、計算する電気回路の定数、RungeKutta法による計算のための刻み幅、最大時間は、Form1_Paint()以外のメソッドからも使われる定数であるので、プログラムの上方にある、「///必要なデザイナ変数です。」の
   System::ComponentModel::Container ^components;
 の下に、グローバル変数として下記を記述する。

   static int X0=50,Y0=30;
   static int SIZE=40;

   static double E=10.0,R=100.0,C=1.5E-8,L=2.2E-4;     //E=10V,R=100Ω,C=0.015μF,L=0.22mH
   static double H=1.0E-8;                    //刻み幅.01μs
   static double TIME_MAX=1.0E-5;              //最大時間μs

2)グローバル変数のさらに下に、下記のようにdelegateキーワードを使用して、Functionを定義する。

   private: delegate double Function(double);

3)最終的に、画面の出力結果を見ながら、Form1のクライアント領域(枠を含まない、ユーザが使える内部領域)を最適化するために、下記のようにカッコ内を書き換える(最初は余裕を見て大きくとってあったが、できたグラフに合わせて狭めた)。これは、「プロパティ」ウインドウにも反映し、「プロパティ」ウインドウ上で書き換えたのと同じ効果がある。「このメソッドの内容をコードエディタで変更しないでください。」の注意があるが、大体は大丈夫のようである。

   this->ClientSize = System::Drawing::Size(480, 270);

4)CR回路とLR回路の二種類の回路に対してRungeKutta法を使用するので、微分方程式を表現するfunction_CRとfunction_LRを、それぞれ順次にFunction型のfunction1に登録して、calculateRungeKutta( )を呼び出す。
5)calculateRungeKuttaメソッドは、Function型の引数、すなわち、微分方程式を表現するメソッド、グラフィック表示に用いる色、および電流の初期値を受け取って、画面上にグラフを描画する。
 4次のRungeKutta法は、すでに説明したものを、そのままコード化したものである。
 Form1_Paint()以外では、Graphicsオブジェクトを取得するのに、
    Graphics^ g=this->CreateGraphics();
 が必要である。

private: System::Void Form1_Paint(System::Object^ sender, System::Windows::Forms::PaintEventArgs^ e) {

   Graphics^ g=e->Graphics;

   int i;
   double time_constant;
   String^ string1;
   Function^ function1;

   //表示域をグレイで描く
   g->FillRectangle(Brushes::LightGray,X0,Y0,SIZE*10,SIZE*5);

   //目盛を入れる
   for(i=0;i<=10;i++){
      g->DrawLine(Pens::Black,X0+i*SIZE,Y0,X0+i*SIZE,Y0+SIZE*5);
      string1=String::Format("{0}μs",i);
      g->DrawString(string1,Font,Brushes::Black,X0-10+i*SIZE,Y0+10+SIZE*5);
   }
   for(i=0;i<=5;i++){
      g->DrawLine(Pens::Black,X0,Y0+i*SIZE,X0+SIZE*10,Y0+i*SIZE);
      string1=String::Format("{0,3:D}mA",i*20);
      g->DrawString(string1,Font,Brushes::Black,X0-5-SIZE,Y0-5+SIZE*5-i*SIZE);
   }

   //CR回路の過渡現象をRungeKutta法で解く
   function1=gcnew Function(this,&RungeKutta::Form1::function_CR);
   calculateRungeKutta(function1,Color::Blue,E/R);

   //CR回路の時定数を表示する
   time_constant=R*C*1.0E+6;
   g->DrawLine(Pens::Blue,X0+time_constant*SIZE,Y0+SIZE*5,X0+time_constant*SIZE,Y0);

   //LR回路の過渡現象をRungeKutta法で解く
   function1=gcnew Function(this,&RungeKutta::Form1::function_LR);
   calculateRungeKutta(function1,Color::Red,0.0);

   //LR回路の時定数を表示する
   time_constant=L/R*1.0E+6;
   g->DrawLine(Pens::Red,X0+time_constant*SIZE,Y0+SIZE*5,X0+time_constant*SIZE,Y0);

}

private: double function_CR(double i){

   return -i/(C*R);

}

private: double function_LR(double i){

   return (E-R*i)/L;

}

private: void calculateRungeKutta(Function^ func,Color col,double current){

   Graphics^ g=this->CreateGraphics();

   Pen^ pen1=gcnew Pen(col);

   int i=0;
   double t,k1,k2,k3,k4;

   for(t=0.0;t<TIME_MAX;t+=H){
      k1=H*func(current);
      k2=H*func(current+k1/2);
      k3=H*func(current+k2/2);
      k4=H*func(current+k3/2);
      current+=(k1+2*k2+2*k3+k4)/6;
      g->DrawRectangle(pen1,X0+i*SIZE*0.01,Y0+SIZE*5-current*SIZE*50,1,1);
      i++;
   }

}

得られた画面
 図5に得られた結果を示す。青はCR回路の電流波形、赤はLR回路の電流波形である。時定数が異なるので、上下対称にはなっていない。点を打つのに2x2の4ピクセルを使用しているので、曲線が太く描かれている。青と赤の縦の線は、それぞれの波形の時定数を解析的に求めた結果を示す。



図5 得られた画面



「Visual C++ の勉強部屋」(目次)へ