2008. 4. 10.
2014. 2. 21. 訂正 delete → delegate、old_y 削除
石立 喬

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

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


 4次のRungeKutta法のやや複雑な場合として、二元連立微分方程式や、時間tを含む微分方程式を数値的に解く。例題は簡単なLCR共振回路の過渡現象とし、直流電圧を印加した場合と、交流電圧を印加した場合について紹介する。波形の表示が必要になるので、最小限の機能を持ったWindowsフォームアプリケーションを使用する。
 RungeKutta法については、すでに「Visual C++ 2008 Express Editionの易しい数値計算(6)」で簡単に述べてあり、「Visual C++ 2005 Express Editionの易しい使い方(15)」でも使われている。

二元連立微分方程式を解くための4次RungeKutta法
 二元連立微分方程式として、下記のような形を想定する。より一般化して、右辺に時間tを含めてある。
   
 上記の二元連立微分方程式の場合には、
   
を繰り返して計算することにより解が得られる。時間tが含まれていない場合は無視しておけばよい。

直流電圧を印加したLCR回路の過渡現象
 図1の回路で、電流iが流れていない状態でスイッチを入れた時の電流iを求める。

図1 LCR直列回路に直流源を印加する

図1の回路から、 の関係が成り立つので、
両辺を微分して、 を解けばよい。
これは二階微分方程式であるが、

として書き直すと、 となり、

の二元の連立一階微分方程式を解けば良いことが分かる。
 初期条件は、i=0、すなわち、x=0であり、最初はすべての電圧がコイルLの両端にかかると考えられるので、Ldi/dt=Eより、di/dt=E/L、すなわち、y=E/Lである。

抵抗Rの値による三通りの波形
の形の微分方程式は、解析的に解ける。
と書き換えて特性二次方程式を求めると、
 となり、この式の根の判別式
  の値により、微分方程式の解の形が決まる。
 LおよびCの値を一定とすると、Rの値によって、次の三通りの波形が得られる。
 D>0の場合 --- 過制動(Overdamping)となり、単調に減衰する指数関数になる。
 D=0の場合 --- 臨界制動(Critical Damping)となり、比較的急速に減衰する。
 D<0の場合 --- 不足制動(Underdamping)となり、包絡線が指数関数的に減衰する振動波形となる。

直流電圧を印加したLCR回路のプログラム
 「易しい数値計算(6)」と異なる点を中心に説明する。
1) 電源電圧および回路素子の定数は、表示したい波形によって最適のものを選ぶので、グローバル変数の宣言と同時に初期化していない。
2) delegateキーワードを使用してFunctionを定義する際、引数をt、x、yの三種類に一般化してある。ただし、直流電圧の印加の場合には、メソッド内部でtを無視する。
3) dx/dtを表す式のメソッドfunction_LCRx()と、dy/dtを表す式のメソッドfunction_LCRy()を用意する。
4) calculateRungeKutta()メソッドは、引数として、dx/dt用の関数func1、dy/dt用の関数func2、グラフの色、xの初期値、yの初期値を受け取り、RungeKutta法による数値計算を行い、グラフに表示する。グラフ描画に当たっては、一つ前の点を表すold_xとold_yを準備しておき、DrawLine()を使用する。
5) 「//-------- この部分だけを変える ----------」の中を書き換えれば、直流電圧印加以外にも一般化して使用できる。

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

static double H=1.0E-8;         //刻み幅.01μsec
static double TIME_MAX=1.0E-5;   //最大時間μsec
double E,R,C,L;

private: delegate double Function(double,double,double);

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

   Graphics^ g=e->Graphics;

   int i;
   String^ string1;
   Function ^func1,^func2;

   //表示域をグレイで描く
   g->FillRectangle(Brushes::LightGray,X0,Y0,SIZE*10,SIZE*10);
   //目盛を入れる
   for(i=0;i<=10;i++){
      //縦軸
      g->DrawLine(Pens::Black,X0+i*SIZE,Y0,X0+i*SIZE,Y0+SIZE*10);
      string1=String::Format("{0}μs",i);
      g->DrawString(string1,Font,Brushes::Black,X0-20+i*SIZE,Y0+10+SIZE*10);
      //横軸
      g->DrawLine(Pens::Black,X0,Y0+i*SIZE,X0+SIZE*10,Y0+i*SIZE);
      string1=String::Format("{0,4:D}mA",(i-5)*20);
      g->DrawString(string1,Font,Brushes::Black,X0-10-SIZE,Y0-5+SIZE*5-(i-5)*SIZE);
   }

   //LCR回路に直流を加えたときの過渡現象
   func1=gcnew Function(this,&RungeKutta1::Form1::function_LCRx);
   func2=gcnew Function(this,&RungeKutta1::Form1::function_LCRy);

   //------------------------ この部分だけを変える ---------------------------

   C=2.0E-9,L=2.0E-5;   //共通の条件設定、C=0.002μF,L=20μH

   //過制動(Overdamping)
   E=50.0;   //E=50V
   R=500.0;  //R=500Ω
   calculateRungeKutta(func1,func2,Color::Blue,0.0,E/L);
   g->DrawString("Overdamping",Font,Brushes::Blue,X0+7*SIZE+10,Y0+8*SIZE+5);
   //臨界制動(Critical Damping)
   E=20.0;   //E=20V
   R=200.0;  //R=200Ω
   calculateRungeKutta(func1,func2,Color::Green,0.0,E/L);
   g->DrawString("Critical Damping",Font,Brushes::Green,X0+7*SIZE+10,Y0+8*SIZE+25);
   //不足制動(Underdamping)
   E=10.0;   //E=10V
   R=20.0;   //R=20Ω
   calculateRungeKutta(func1,func2,Color::Red,0.0,E/L);
   g->DrawString("Underdamping",Font,Brushes::Red,X0+7*SIZE+10,Y0+8*SIZE+45);

   //--------------------------------------------------------------------------

}

private: double function_LCRx(double t,double x,double y){

   return y;

}

private: double function_LCRy(double t,double x,double y){

   return -(x/C+R*y)/L;

}

private: void calculateRungeKutta(Function^ func1,Function^ func2,Color col,double xinit,double yinit){

   Graphics^ g=this->CreateGraphics();

   Pen^ pen1=gcnew Pen(col);

   double t,m1,m2,m3,m4,n1,n2,n3,n4;
   int i=0;
   double x=xinit;   //電流(i)の初期値
   double y=yinit;   //電流の微分(di/dt)の初期値
   double old_x;

   for(t=0.0;t<TIME_MAX;t+=H){
      m1=H*func1(t,x,y);
      n1=H*func2(t,x,y);
      m2=H*func1(t+H/2,x+m1/2,y+n1/2);
      n2=H*func2(t+H/2,x+m1/2,y+n1/2);
      m3=H*func1(t+H/2,x+m2/2,y+n2/2);
      n3=H*func2(t+H/2,x+m2/2,y+n2/2);
      m4=H*func1(t+H/2,x+m3/2,y+n3/2);
      n4=H*func2(t+H/2,x+m3/2,y+n3/2);
      x+=(m1+2*m2+2*m3+m4)/6;
      y+=(n1+2*n2+2*n3+n4)/6;
      if(i==0){
         old_x=x;
      }
      else{
         g->DrawLine(pen1,X0+(int)((i-1)*SIZE*0.01),Y0+SIZE*5
                             -(int)(old_x*SIZE*50),X0+(int)(i*SIZE*0.01),Y0+SIZE*5-(int)(x*SIZE*50));
         old_x=x;
      }
      i++;
   }

}

直流電圧を印加したLCR回路の場合の画面
 図2は得られた画面である。
 青色のグラフは、抵抗Rを500Ωと大きくして、過制動(Overdamping)にした場合で、波形を大きく表示するために、電圧Eを50Vに大きくとってある。
 緑色のグラフは、臨界制動(Critical Damping)の条件に合致した抵抗Rの値200Ωを用いた場合で、波形のサイズを適切にするために、電圧Eを20Vにしてある。最も短時間で収束する。
 赤色のグラフは、抵抗Rを20Ωと小さく選んで不足制動(Underdamping)にした場合で、きれいな減衰振動になっている。グラフから振動周期を読み取ると1.25μsecとなっていて、
理論的に求めた、


 と良く一致している。




図2 得られた画面


交流電圧を印加したLCR回路の過渡現象
 図3の回路で、電流iが流れていない状態でスイッチを入れた時の電流iを求める。

図3 LCR直列回路に交流電圧源を印加する

図3の回路から、 の関係が成り立つので、
両辺を微分して、 を解けばよい。
直流電圧の印加の場合と同様に変形すると、

の二元連立微分方程式を解けば良いことが分かる。
 初期条件は、i=0、すなわち、x=0であり、di/dtも0であるので、y=0である。

共振周波数(Resonant Frequency)と共振角周波数(Resonant Angular Frequency)
LCR直列回路には、最も電流が流れやすい周波数がある。これを共振周波数または自然周波数と言う。普通に使われている周波数fの単位はヘルツ(Hz)であるが、これは、一秒間に何回正負の方向が変わるかを表す。しかし、電気回路の計算では、角周波数(Angular Frequency)の方が便利であり、角周波数ωの単位はラジアン/秒である。周波数fと角周波数ωとの間には、ω=2πfの関係がある。特に、共振周波数をf0、共振角周波数をω0と呼ぶ。
LCR直列回路の共振角周波数は、 の関係があるので、
  これを確認するために、ここでは、与えられた回路素子の定数から共振角周波数を計算し、交流電源として使用した。共振角周波数の前後の交流電源を用いた場合も電流波形を計算し、比較した。

交流電圧を印加したLCR回路のプログラム
 直流電圧を印加した場合のプログラムと異なる点のみを記述する。

double OMEGA;   //角周波数ω

   //------------------------ この部分だけを変える ---------------------------

   //LCR回路に正弦波を加えたときの過渡現象
   func1=gcnew Function(this,&RungeKutta1::Form1::function_LCRx);
   func2=gcnew Function(this,&RungeKutta1::Form1::function_LCRSINy);
   E=1.0,R=20.0,C=1.5E-9,L=1.5E-5;   //共通の設定条件、E=1V,R=20Ω,C=0.0015μF,L=15μH
   //共振角周波数
   OMEGA=Math::Sqrt(1/(L*C)-R*R/(4*L*L));
   calculateRungeKutta(func1,func2,Color::Red,0.0,0.0);
   g->DrawString("Resonant Frequency",Font,Brushes::Red,X0+6*SIZE+10,Y0+8*SIZE+5);
   //共振角周波数の1.5倍
   OMEGA=1.5*Math::Sqrt(1/(L*C)-R*R/(4*L*L));
   calculateRungeKutta(func1,func2,Color::Green,0.0,0.0);
   g->DrawString("Resonant Frequency * 3/2",Font,Brushes::Green,X0+6*SIZE+10,Y0+8*SIZE+25);
   //共振角周波数の2/3倍
   OMEGA=0.667*Math::Sqrt(1/(L*C)-R*R/(4*L*L));
   calculateRungeKutta(func1,func2,Color::Blue,0.0,0.0);
   g->DrawString("Resonant Frequency * 2/3",Font,Brushes::Blue,X0+6*SIZE+10,Y0+8*SIZE+45);

   //--------------------------------------------------------------------------


private: double function_LCRSINy(double t,double x,double y){

   return -(x/C+R*y-OMEGA*E*Math::Cos(OMEGA*t))/L;

}

交流電圧を印加したLCR回路の画面
 図4は得られた画面である。
 赤色のグラフは、交流電源として、共振周波数ω0を使用した場合の波形で、共振現象によって、時間とともに振幅が増大し、一定振幅に達している。
 緑色のグラフは、共振周波数の1.5倍の交流電源を用いたもので、共振現象が起きない。青色のグラフは、共振周波数の2/3倍の交流電源を用いたもので、これも共振現象が起きていない。



図4 得られた画面


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