2008. 6. 23.
石立 喬

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

――― 定積分の計算と、正弦波形の平均値と実効値 ―――

 数値積分は、コンピュータによる数値計算で取り上げられるテーマの一つである。被積分関数が比較的単純な形で表現できる場合には解析的に解くことができるが、任意の波形に対しては、離散的なデータを用いた数値計算が必要になる。ここでは、直接的な実用性はあまり無いが、解説的な意味で定積分を取り上げた。

面積を求める定積分
 定積分の数値計算は離散的なデータを使用するので、誤差は免れない。区間a、b(aとbを含む)を等間隔に分割数nで割り、等間隔で割られた区間を多項式で近似して合計する。多項式近似にはラグランジュ補間が用いられる。実際には、その都度補間することなく、ニュートン・コーツの公式(Newton-Cotes Formula)を使用する。
 代表的なものは、台形法(台形則、Trapezoidal Rule、2点を使用)とシンプソン法(シンプソン則、Simpson’s Rule、3点を使用)である。これより精度の高い方法もあるが、複雑な割には効果が薄いので、あまり使われない。
 図1に、台形法、シンプソン法の他に、さらに単純な矩形法を併せて、原理を示す。グレイで塗りつぶした面積を計算し、関数の範囲全体に亘って合計する。

図1 定積分の求め方のいろいろ


◎ 矩形法
 矩形法は短冊法とも言われる。誤差が大きいので、定積分の数値計算には一般的でないが、比較のために紹介する。座標上の一点のみのデータと刻み幅hで矩形の面積を求め、それを累計する。図1(A)を参照して式で示すと、データ点f0からfnまでに対して、
   
となり、一番簡単である。図から分かるように、右上がりの曲線に対しては小さ目の、右下がりの曲線には大き目の結果が得られる。
◎ 台形法
 台形法は、座標上の隣接する二点のデータと刻み幅hで作る台形の面積を求め、それを累計する。これは、データの二点間を直線で結んで近似することに相当する。図2(B)を参照して式で示すと、データ点f0からfnまでに対して、
 
となり、矩形法に対して、f0 とfnを補正した形になっている。
◎ シンプソン法
 シンプソン法は、座標上の隣接する三点のデータから、それらを結ぶ二次曲線を補間で求め、それと刻み幅2hで作る面積を求め、それを累計する。これは、データの三点間を二次曲線で結んで近似することになる。刻み幅hを二つずつ一まとめにして計算するので、刻みの総個数を偶数に選ぶ必要がある(偶数個のため、分割数を2nで表してある。シンプソン法での2nは、台形法のnに相当する。)。図3(C)を参照して式で示すと、データ点f0からf2nまでに対して、
 
となり、やや複雑ではあるが、精度が上がる。図では、実際の曲線と近似した二次曲線の差異が明白に表現できていないが、誤差は存在する。

πの計算による各定積分方法の比較
前述した各定積分方法を比較する目的で、 を計算する。
  これは、半径1の円の面積の1/4に相当するので、これを4倍するとπが得られる。πを求める計算方法としては、効率が悪い、平方根ライブラリで精度が決まる、などの理由で使われることはない。

定積分を求めるプログラム
 非積分関数は、ユーザ関数functionとして用意した。これにより、別の関数に変えたときのプログラム変更を容易にした。
 台形法の場合の累積値sum1に、初期値として(f0 + fn)/2を入れた。この場合、f0=1.0、fn=0.0である。

#include "stdafx.h"

using namespace System;

double function(double x){

   return Math::Sqrt(1.0-x*x);

}

int main()
{

   Console::Title="定積分のいろいろ";

   int i,j;
   double h,sum1,sum2;
   double rectangle,trapezoidal,simpson;
   double a=0,b=1;
   array<int>^ n={100,200,500,1000,2000,5000,10000,20000,50000,100000};
   String^ string1;

   Console::WriteLine("   n        矩形法          台形法         シンプソン法    ");
   Console::WriteLine("------- ------------------ ------------------ ------------------");

   for(j=0;j<10;j++){
      h=(b-a)/n[j];
      //矩形法
      sum1=0.0;
      for(i=0;i<n[j];i++)
         sum1+=function(i*h);
      rectangle=h*sum1;
      //台形法
      sum1=0.5;
      for(i=1;i<n[j];i++)
         sum1+=function(i*h);
      trapezoidal=h*sum1;
      //シンプソン法
      sum1=0.0;
      for(i=1;i<=n[j];i+=2)
         sum1+=function(i*h);
      sum2=0.0;
      for(i=2;i<n[j];i+=2)
         sum2+=function(i*h);
      simpson=(1.0+4.0*sum1+2.0*sum2+0.0)*h/3.0;  //f0=1.0, f2n=0.0

      string1=String::Format("{0,6}   {1}   {2}   {3}",n[j],4.0*rectangle,4.0*trapezoidal,4.0*simpson);
      Console::WriteLine(string1);
   }

   Console::WriteLine();
   return 0;

}

得られた結果
 図2は得られた結果である。分割数nが同じ場合には(シンプソン法の場合は2nに相当する)、矩形法、台形法、シンプソン法の順に精度が高くなる。分割数n=100の場合、矩形法で0.6%、台形法で0.04%、シンプソン法で0.015%の精度があるので、用途によっては簡単な方法でも十分なことが分かる。


図2 得られた結果


電気回路に関係の深い定積分
 一般に、100Vの交流電圧と言えば、最大値が約141Vの正弦波の電圧を指す。このように、交流波形には、最大値(maximum value)、実効値(effective value、特に言及しなければこの値を使う)、平均値(average value)がある。
 平均値は、下式のように、正弦波の半周期、すなわち0からπまで、を積分し、半周期の期間で割って求める。全周期で積分すると0になるので、平均値と言えば半周期の平均を意味する。
  
 実効値は、その交流波形と電力的に等価な直流の値を求めるもので、下式のように、正弦波の二乗を半周期積分し、半周期の期間で割って、さらに平方根を求めたものである。二乗して、平均して平方根をとるので、r.m.s.(root mean square)とも呼ばれる。
  
 これらは解析的に解くことができ、最大値1の正弦波の平均値は 2/π≒0.636619772、実効値は 1/√2≒0.707106781 である。

正弦波形の平均値と実行値の計算
 念のために、プログラムを作成して結果を確認した。分割数nは200にとり(式の説明の上では、2nに相当)、シンプソン法で計算した。図3に示す通り、十分な精度が得られている。

int main()
{

   Console::Title="平均値と実効値";

   int i;
   double h,sum1,sum2;
   double simpson,average,effective;
   double a=0,b=Math::PI;
   int n=200;

   h=(b-a)/n;
   //平均値
   sum1=0.0;
   for(i=1;i<=n;i+=2)
      sum1+=Math::Sin(i*h);
   sum2=0.0;
   for(i=2;i<n;i+=2)
      sum2+=Math::Sin(i*h);
   simpson=(4.0*sum1+2.0*sum2)*h/3.0;  //f0=0.0, f2n=0.0
   average=simpson/Math::PI;
   //実効値
   sum1=0.0;
   for(i=1;i<=n;i+=2)
      sum1+=Math::Sin(i*h)*Math::Sin(i*h);
   sum2=0.0;
   for(i=2;i<n;i+=2)
      sum2+=Math::Sin(i*h)*Math::Sin(i*h);
   simpson=(4.0*sum1+2.0*sum2)*h/3.0;  //f0=0.0, f2n=0.0
   effective=Math::Sqrt(simpson/Math::PI);

   Console::WriteLine("平均値/最大値=" + average);
   Console::WriteLine("実効値/最大値=" + effective);

   Console::WriteLine();
   return 0;

}


図3 正弦波を積分した結果


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