2008. 7. 16.
石立 喬

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

――― 波形の数値積分と電気回路への応用 ―――

 電気回路で必要な積分演算の一つに、コンデンサに流れる電流i(t)と、その端子電圧v(t)の関係があり、下式で表される。
   
 ただし、Cはコンデンサの容量(キャパシタンス、単位はファラッド)である。
これを、離散的なデータを用い、コンピュータで時々刻々積分することによって、電気回路の波形を求めることができる。

波形を求める数値積分
 与えられた範囲の面積などを求める定積分ではなく、時々刻々変化する変数の積分値を、離散的に逐次出力するには、
    
に対して、hを刻み幅とし、刻みの個数nを、n=t/h として、離散的に下式で求めることができる。
    
 これを図で示すと、図1のように、時刻tiのときの関数値fiに刻み幅hを掛け、その値を累加算して、時刻tiの積分値とすることになる。面積を求める定積分の場合(「易しい数値計算(10)」参照)と異なり、形状は矩形法に似ているが、精度的には台形法に近い。


図1 不定積分で累加算する方法


 やや極端な例を図2によって説明する。(A)のような被積分波形に対して、時刻tiの関数値fiに刻み幅hを掛けると、(B)のグレイ部分の面積を求めていることになり、実際よりは面積Sだけ大きくなる。一方、無理して矩形法を適用したのが(C)で、時刻tiを中心に刻み幅hをとり、時刻ti-h/2における関数値fi-1/2と時刻ti+h/2における関数値fi+1/2で作る矩形の面積を求めると、実際よりも面積Sだけ小さくなる。(B)と(C)において、面積Sが等しいことは容易に分かることと思う。すなわち、誤差について考察すると、時刻tiの関数値fiに刻み幅hを掛けて矩形の面積を求めるのは、定積分における台形法と同程度の精度が得られることが分かる。
 見方を変えて考えると、波形の数値積分の場合は、面積を求める定積分に比べて、積分区間が長く、積分開始点と終了点の影響が無視でき、矩形法が台形法と同等になっているとも言える。これは、「易しい数値計算(10)」の、矩形法の式と台形法の式を比較すれば容易に理解できる。


図2 波形の数値積分は矩形法に相当するが、精度的には台形法と同等


数値積分の簡単な例
 電気回路で頻繁に使用する三角関数について、
    
の不定積分は良く知られているので、ここでは、sinθの波形の数値積分を試みる。Cは積分定数であって、任意の値を取ることができる。後述するプログラムでは、積分結果の初期値を-1とし、被積分関数のsinθ波形と積分結果の-cosθ波形のゼロレベルを揃えてある。

正弦波の数値積分プログラム(コンソールアプリケーション)
 結果の表示範囲を角度(ラジアン)0から2πまでとし、4096に分割して計算し、π/8刻みで表示している。出発点の0以外に16回表示するので、if(i % (n/16) == 0) の条件が成立する場合のみ表示することにした。刻み幅hは、関数値fiごとに掛けるのではなく、累加算した後に、表示の際に掛けている。そのため、初期値integralをあらかじめhで割っておく。

int main()
{
   Console::Title="三角関数の数値積分";

   double a=0,b=2*Math::PI;
   int n=4096;
   double h,integral,theta;
   String^ string1;

   h=(b-a)/n;

   Console::WriteLine("   角度(θ、xπrad)    被積分関数(sinθ)    積分結果(-cosθ)  ");
   Console::WriteLine(" ------------------ -------------------- --------------------");

   integral=-1.0/h;
   for(int i=0;i<=n;i++){
      integral+=Math::Sin(i*h);
      if(i % (n/16) == 0){
         theta=i*2/(double)n;
         string1=String::Format("     {0,4:F3}              {1,11:F8}              {2,11:F8} ",
                                              theta,Math::Sin(i*h),integral*h);
         Console::WriteLine(string1);
      }
   }

   Console::WriteLine();
   return 0;
}

正弦波の数値積分の実行結果(コンソールアプリケーション)
 分割数n=4096で一周期分を計算したが、非常に十分な精度とは言えない。分割数 n=16384 では、さらに改善されることを確認してある。


図3 正弦波の数値積分の結果


正弦波の数値積分プログラム(Windowsフォームアプリケーション)
 結果の表示範囲を角度0から4πまでの二周期とし、512分割で計算しているが、グラフ表示の精度としては十分である。ここでは、累加算結果を逐次出力するので、刻み幅hをその都度fiに掛けている。

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

   Graphics^ g=e->Graphics;

   int X0=10,Y0=100;
   int x,y1,y2,old_x,old_y1,old_y2;
   double a=0.0,b=4*Math::PI;
   int n=512;
   double h,integral;

   h=(b-a)/n;

   //グラフの背景を描く
   g->FillRectangle(Brushes::LightGray,X0,Y0-90,512,180);
   g->DrawLine(Pens::White,X0,Y0,X0+512,Y0);
   for(int i=1;i<=3;i++)
      g->DrawLine(Pens::White,X0+128*i,Y0-90,X0+128*i,Y0+90);

   integral=-1.0;
   for(int i=0;i<=n;i++){
      integral+=Math::Sin(i*h)*h;
      x=X0+i;
      y1=Y0-Math::Sin(i*h)*50;   //被積分波形
      y2=Y0-integral*50;       //積分波形
      if(i==0){
         old_x=x;
         old_y1=y1;
         old_y2=y2;
      }
      else{
         g->DrawLine(Pens::Red,old_x,old_y1,x,y1);    //被積分波形
         g->DrawLine(Pens::Blue,old_x,old_y2,x,y2);    //積分波形
         old_x=x;
         old_y1=y1;
         old_y2=y2;
      }
   }

}

正弦波の数値積分の実行結果(Windowsフォームアプリケーション)
 図4で、赤が被積分関数のsinθ波形で、青は積分結果の-cosθ波形である。


図4 正弦波の数値積分の結果


半波整流電流をCR並列回路に流し込んだ場合の電圧波形
 交流電源から直流を取り出すには、交流電源を整流して(半波整流または全波整流)、コンデンサで平滑化する。実際には、半導体素子による帰還回路などを用いて、より高能率で安定な直流電圧を得ることができる。
 ここでは、極めて簡単な例として、半波整流電流源をCR回路に印加する場合の出力電圧を計算する。
 図5は、ここで計算する回路を示したもので、半波整流電流源の電流をiS、平滑コンデンサCに流れる電流をiC、負荷抵抗Rに流れる電流をiRとすると、
   
から、時間とともに変化する平滑コンデンサCの電圧vCを逐次求めることができる。


図5 電気回路の簡単な例として選んだCR回路


半波整流波形の数値積分のプログラム(Windowsフォームアプリケーション)
 最初は初期値としてvC=0を与えるので、抵抗に流れる電流iR=vC/Rはゼロとなり、定電流源iSのすべてがコンデンサに流れ込み、iC=iSとなる。このコンデンサ電流iCを積分して新しいvCを求め、iRを求めて、iC=iS-iR の数値積分を繰り返す。
 グラフ上の全表示時間8msecに対して、8サイクルの半波整流波形を発生させるために、全時間範囲を16で割り、商が偶数の場合は正弦波を出力し、奇数の場合には出力しないようにした。

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

   Graphics^ g=e->Graphics;

   int X0=10,Y0=100,Y1=200;
   int x,y1,y2,old_x,old_y1,old_y2;

   double C=1.0E-6;        //容量 C=1μF
   double R=1.5E+3;        //抵抗 R=1.5kΩ、他に 200Ω、2kΩの場合も試した
   double Im=1.0E-3;        //半波整流電流波形の最大値 Im=1mA
   double frequency=1.0E+3;   //周波数 frequency=1kHz
   double omega=2*Math::PI*frequency;   //角周波数 ω=2πf
   double a=0.0,b=8.0E-3;     //グラフ表示範囲 0〜8msec
   int n=512;
   double h;
   double is;    //電流源電流
   double ir;    //抵抗Rを流れる電流
   double vc;    //コンデンサCの端子電圧

   h=(b-a)/n;

   //グラフの背景を描く
   g->FillRectangle(Brushes::LightGray,X0,Y0-90,512,280);
   g->DrawLine(Pens::White,X0,Y0,X0+512,Y0);
   g->DrawLine(Pens::White,X0,Y1,X0+512,Y1);
   for(int i=1;i<=7;i++)
      g->DrawLine(Pens::White,X0+64*i,Y0-90,X0+64*i,Y1+90);

   vc=0.0;
   for(int i=0;i<=n;i++){

      //定電流波形(半波整流波形)の生成
      if((i/32) % 2 ==0)
         is=Im*Math::Sin(omega*i*h);
      else
         is=0.0;

      //抵抗Rを流れる電流(これが定電流波形から差し引かれる)
      ir=vc/R;

      vc+=1.0/C*(is-ir)*h;  //コンデンサCの端子電圧vcを積分で求める

      x=X0+i;          //時間軸
      y1=Y0-is*50.0E+3;   //定電流波形(半波整流波形)の表示
      y2=Y1-vc*100;     //コンデンサCの端子電圧(積分波形)の表示
      if(i==0){
         old_x=x;
         old_y1=y1;
         old_y2=y2;
      }
      else{
         g->DrawLine(Pens::Red,old_x,old_y1,x,y1);   //定電流波形
         g->DrawLine(Pens::Blue,old_x,old_y2,x,y2);   //コンデンサCの端子電圧
         old_x=x;
         old_y1=y1;
         old_y2=y2;
      }
   }

}

半波整流波形の数値積分の実行結果(Windowsフォームアプリケーション)
 図6で、上の赤いグラフは、半波整流電流源波形を示し、下の青い線は、コンデンサCの端子電圧を示す。負荷抵抗R=1.5kΩの場合には、8周期、すなわち8msecまでの間に、ほぼ安定的な出力が得られる。リップル(直流に載った交流成分)が多く、実用的ではないが、計算例として示した。
 図7は、負荷抵抗Rを減らして200Ωにした場合で、半周期の間にコンデンサCに充電された電荷は、次の半周期の間に、ほぼ完全に放電している。
 図8は、負荷抵抗を増やした場合で、負荷抵抗の増大につれて、出力電圧が増加する様子が分かる。負荷抵抗をさらに増やして行けば(軽い負荷しか使わない)、出力電圧値はさらに高まり、相対的にリップルの少ない直流電圧が得られる。


図6 R=1.5kΩで得られた結果



図7 R=200Ωの場合



図8 R=2kΩの場合



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