数値計算@数値積分@その2@シンプソン法

参考文献及びWEB

参考文献

(1)「電気・電子工学のための数値計算法入門」橋本修著　総合電子出版社
(2)「ディジタル信号処理技術」玉井徳みち、長島厚、藤田泰弘、若井修造著 日経BP社
(3)「よくわかる有限要素法」福森栄次著 Ohmsha

参考WEB

http://szksrv.isc.chubu.ac.jp/java/physics/rlc/rlc0.html

http://www.asahi-net.or.jp/~jk2m-mrt/kiso_RLC.htm

http://ja.wikipedia.org/wiki/%E4%B8%89%E8%A7%92%E9%96%A2%E6%95%B0

http://www.cmplx.cse.nagoya-u.ac.jp/~furuhashi/education/CircuitMaker/chap1.pdf

http://www.akita-nct.ac.jp/~yamamoto/lecture/2003/5E/lecture_5E/diff_eq/node2.html

http://chemeng.on.coocan.jp/cemath/cemath08.html

http://homepage3.nifty.com/skomo/f6/hp6_3.htm

http://homepage1.nifty.com/gfk/rungekutta.htm

あわせて読む

(a) 微分の説明

微分方程式の解を求める

(b) 微分方程式の解を数値計算で求める

(c) 積分の説明

(d) 関数が微分の形 f'(x) の形で表されているときの積分値を求める

(一般的に数値積分をする意味は、関数自体が変化の式(例えば速度の関数)などで表されているとき
ある時刻から時刻までの距離を求めたいときなどで使用される。

もちろんその変化の式が、手計算によって計算可能ならば積分によって時間と距離の関係を導出し
その時間を入れることで直接距離を導いたほうが早い。

しかし、一般的に時間と距離の関係のような式自体が導出することが不可能な場合が世の中には多く
時間と速度の関係のような変化の式を見つけ出すことの方が容易なのだが、これを手計算で積分すること
は非常に難しいため数値積分が提案された)

台形法とシンプソン法の比較

ここでは積分を台形法とシンプソン法を用いてプログラムで計算してみます。
台形法は、1つの区間だけを考えます。
シンプソン法は、2つの区間を考えます。

台形法の理論

区間[a,b]で連続な関数 y = f(x) を等間隔にn分割して、その分割点を

x_0 , x_1 , x_2, \cdots , x_n (1-1)

としたとき、その関数曲線上の点

(x_0,f(x_0)), (x_1,f(x_1)), \cdots (x_n,f(x_n)) (1-2)

を通る曲線の方程式を近似式で表したもの

f(x) = \frac{ (x-x_1)(x-x_2)\cdots(x-x_n) }{ (x_0-x_1)(x_0-x_2)\cdots(x_0-x_n) } f(x_0) + \frac{ (x-x_0)(x-x_2)\cdots(x-x_n) }{ (x_1-x_0)(x_1-x_2)\cdots(x_1-x_n) } f(x_1) + \cdots + \frac{ (x-x_0)(x-x_1)\cdots(x-x_{n-1}) }{ (x_n-x_0)(x_n-x_1)\cdots(x_n-x_{n-1}) } f(x_n) (1-3)

となります。

ここで、2点 (x_i, y(x_i)) , (x_{i+1}, y(x_{i+1}))
を通る直線を考える。
y - f(x_{i+1}) = \frac{ f(x_{i+1}) - f(x_i)}{ x_{i+1} - x_{i} }(x - x_{i+1} )

y = \frac{ x - x_{i+1}}{ x_{i+1} - x_{i} }f(x_{i+1}) - \frac{ x - x_{i+1} }{ x_{i+1} - x_{i} }f(x_{i}) + \frac{ x_{i+1} - x_{i} }{ x_{i+1} - x_{i} }f(x_{i+1})

g(x) = \frac{ x - x_{i+1} }{ x_i - x_{i+1} }f(x_i) + \frac{ ( x - x_{i} )}{ ( x_{i+1} - x_{i} ) }f(x_{i+1}) (1-4)

この直線に対して積分すると
\int_{x_i}^{x_{i+1}}f(x)dx
\simeq \int_{x_i}^{x_{i+1}}g(x)dx
= \frac{f(x_i)}{x_i - x_{i+1}} \int_{x_i}^{x_{i+1}} ( x - x_{i+1} )dx + \frac{f(x_{i+1})}{x_{i+1} - x_{i}} \int_{x_i}^{x_{i+1}} ( x - x_{i} )dx

= \frac{f(x_i)}{x_i - x_{i+1}}\Big[ \frac{1}{2} x^2 - x_{i+1}x \Big]_{x_i}^{x_{i+1}} + \frac{f(x_{i+1})}{x_{i+1} - x_{i}}\Big[ \frac{1}{2} x^2 - x_{i}x \Big]_{x_i}^{x_{i+1}}

= \frac{f(x_i)}{x_i - x_{i+1}}(\frac{1}{2}x^2_{i+1} - x^2_{i+1} -\frac{1}{2}x^2_i + x_{i}x_{i+1}) + \frac{f(x_{i+1})}{x_{i+1} - x_{i}}(\frac{1}{2}x^2_{i+1} - x_{i}x_{i+1} -\frac{1}{2}x^2_i + x^2_{i})

= \frac{f(x_i)}{x_i - x_{i+1}}(-\frac{1}{2}x^2_{i+1} -\frac{1}{2}x^2_i + x_{i}x_{i+1}) + \frac{f(x_{i+1})}{x_{i+1} - x_{i}}(\frac{1}{2}x^2_{i+1} - x_{i}x_{i+1} +\frac{1}{2}x^2_i)

= -\frac{1}{2}\frac{f(x_i)}{x_i - x_{i+1}}(x^2_{i+1} +x^2_i - 2x_{i}x_{i+1}) + \frac{1}{2}\frac{f(x_{i+1})}{x_{i+1} - x_{i}}(x^2_{i+1} - 2x_{i}x_{i+1} +x^2_i)

= -\frac{1}{2}\frac{f(x_i)}{x_i - x_{i+1}}(x_{i} - x_{i+1})^2 + \frac{1}{2}\frac{f(x_{i+1})}{x_{i+1} - x_{i}}(x_{i+1} - x_{i})^2

= -\frac{1}{2}f(x_i)(x_i - x_{i+1}) + \frac{1}{2}f(x_{i+1})(x_{i+1} - x_{i})

= \frac{1}{2}(x_{i+1} - x_{i})\Big{f(x_i) + f(x_{i+1})\Big} (1-5)

ここで、
h = (x_{i+1} - x_{i})
f(x_i) = y_i
f(x_{i+1}) = y_{i+1}
とおくと

\int_{x_i}^{x_{i+1}}g(x)dx = \frac{1}{2}h\Big( y_i + y_{i+1} \Big) (1-6)

となる。

これを一区間ではなく、区間[a,b]にまで拡張します。
すると区間[a,b]を等間隔にn分割すると、

\LARGE h = \frac{b-a}{n}
ここで
a = x_0 < x_1 < x_2 < \cdots < x_{n-1} < x_n = b

区間[x_0,x_{1}],[x_1,x_{2}],[x_2,x_{3}],\cdots,[x_{n-1},x_{n}]
の各微小台形の面積の総和は、

\LARGE \int_{a}^{b} f(x) dx \approx \frac{h}{2}\Big{\sum_{i=1}^{n}(y_{i-1} + y_i) \Big}
= \LARGE \frac{h}{2}\Big{(y_{0} + y_{1}) + (y_{1} + y_{2}) + (y_{2} + y_{3}) + \cdots + (y_{n-1} + y_{n})\Big}
= \LARGE \frac{h}{2}\Big{(y_{0} + y_{n}) + 2(y_{1} + y_{2} + y_{3} + \cdots + y_{n-1} )\Big} (1-7)

となります。

シンプソン法

3点 (x_{i-1}, y(x_{i-1})), (x_i, y(x_i)) , (x_{i+1}, y(x_{i+1}))
を通る直線を考える。

g(x) = \frac{ (x - x_{i})(x - x_{i+1}) }{ (x_{i-1} - x_{i})(x_{i-1} - x_{i+1}) }f(x_{i-1}) + \frac{ (x - x_{i-1})(x - x_{i+1}) }{ (x_{i} - x_{i-1})(x_{i} - x_{i+1}) }f(x_{i}) + \frac{ (x - x_{i-1})(x - x_{i}) }{ (x_{i+1} - x_{i-1})(x_{i+1} - x_{i}) }f(x_{i+1})

\int_{x_{i-1}}^{x_{i+1}}f(x)dx
\simeq \int_{x_{i-1}}^{x_{i+1}}g(x)dx
= \int_{x_{i-1}}^{x_{i+1}} \Big{ \frac{ (x - x_{i})(x - x_{i+1}) }{ (x_{i-1} - x_{i})(x_{i-1} - x_{i+1}) }f(x_{i-1}) + \frac{ (x - x_{i-1})(x - x_{i+1}) }{ (x_{i} - x_{i-1})(x_{i} - x_{i+1}) }f(x_{i}) + \frac{ (x - x_{i-1})(x - x_{i}) }{ (x_{i+1} - x_{i-1})(x_{i+1} - x_{i}) }f(x_{i+1}) \Big} dx

この直線の方程式は、以下となる。

\int_{x_{i-1}}^{x_{i+1}}g(x)dx
= \frac{h}{3}\Big( y_{i-1} + 4y_{i} +y_{i+1}\Big) (2-1)

ただし、
x_{i+1} - x_i = x_i - x_{i-1} = h
f(x_{i-1}) = y_{i-1}
f(x_{i}) = y_{i}
f(x_{i+1}) = y_{i+1}

これを一区間ではなく、区間[a,b]にまで拡張します。
台形法では、区間の距離をhとしましたが、
シンプソン法の場合は、区間の距離は2hとなります。すると区間[a,b]をn回の積分で行うとすると、

\LARGE 2h = \frac{b-a}{n}
となるので、
\LARGE h = \frac{b-a}{2n}

また、
a = x_0 < x_1 < x_2 < \cdots < x_{2n-1} < x_{2n} = b
つまり、
n回の計算においてxの数はx_{2n}個となります。

区間[x_0,x_{1},x{2}],[x_{2},x_{3},x_{4}],[x_{4},x_{5},x_{6}],\cdots,[x_{2n-2},x_{2n-1},x_{2n}]
の各微小台形の面積の総和は、

\LARGE \int_{a}^{b} f(x)
= \LARGE \int_{x_0}^{x_2} f(x) dx \approx \frac{h}{3}\Big{(y_{0} + 4y_{1} + y_{2})\Big}
+ \LARGE \int_{x_2}^{x_4} f(x) dx \approx \frac{h}{3}\Big{(y_{2} + 4y_{3} + y_{4})\Big}
+ \LARGE \int_{x_4}^{x_6} f(x) dx \approx \frac{h}{3}\Big{(y_{4} + 4y_{5} + y_{6})\Big}

\LARGE :

\LARGE :

+ \LARGE \int_{x_{2n-2}}^{x_{2n}} f(x) dx \approx \frac{h}{3}\Big{(y_{2n-2} + 4y_{2n-1} + y_{2n})\Big}

= \LARGE \frac{h}{3}\Big{(y_{0} + 4y_{1} + y_{2}) + (y_{2} + 4y_{3} + y_{4}) + (y_{4} + 4y_{5} + y_{6}) + \cdots + (y_{2n-2} + 4y_{2n-1} + y_{2n})\Big}

= \LARGE \frac{h}{3}\Big{ y_{0} + 2(y_{2} + y_{4} + y_{6} + \cdots + y_{2n-2} ) + \LARGE 4(y_{1} + y_{3} + y_{5} + \cdots + y_{2n-1}) + y_{2n} \Big}

となります。

[例題] 関数 y = \frac{4}{1+x^2} を区間 [0,1] においてシンプソン法と台形法を用いて解け。

ソースコード

// suuti_sekibun.cpp : コンソール アプリケーションのエントリ ポイントを定義します。
//

#include "stdafx.h"
#include<iostream>
#include<fstream>
#include<io.h>

double sympson_method(double (*function)(double x),double a,double b,int num);
double daikei_method(double (*function)(double x),double a,double b,int num);
double func(double x); 

int _tmain(int argc, _TCHAR* argv[])
{
  errno_t err; 
  FILE *fp_output;
  double ans_sympson=0.;
  double ans_daikei=0.;

  if( err = fopen_s(&fp_output,"simpsonAndDaikei_sekibun.txt","w")  != 0){ exit(2);}

  double a=0.;
  double b=1.0;

  //分割数を1から200まで分割する
  for(int n=1;n<=200;n++){
     double h2 = (b-a)/(double)n;
     double h=h2/2.;

     ans_daikei  = daikei_method(&func,a,b,n);
     ans_sympson = sympson_method(&func,a,b,n);

     fprintf_s(fp_output,"%f %f %f\n",h,ans_daikei,ans_sympson);
  }

  if(fp_output != NULL) fclose(fp_output);
  //getchar();
  return 0;
} 

//積分の範囲[a,b]とその範囲の分割数num
double sympson_method(double (*function)(double x),double a,double b,int num){


  //積分の範囲とその範囲を分割する数numから求められる離散間隔
  double h=0.;
  double h2=0.;

  //シンプソン法で使われる変数
  double ans=0;
  double term1=0.;
  double term2=0.;
  double term3=0.;
  double term4=0.;

  h2 = (b-a)/(double)num;
  h=h2/2.;
  ans=0.;
  term1=0.;
  term2=0.;
  term3=0.;
  term4=0.;

  ///// ここからsympson法 ///////////////
  //第1項の計算
  term1=(*function)(a);

  //第2項の計算 y2+y4+y6+.....+y_{2n-2}
  for(int i=2;i<=(2*num-2);i=i+2){
     term2=term2+(*function)(a+i*h);
  }

  //第3項の計算 y1+y3+y5+....+y_{2n-1}
  for(int i=1;i<=(2*num-1);i=i+2){
     term3=term3+(*function)(a+i*h);
  }

  //第4項の計算
  term4=(*function)(b);

  ans=h/3.*(term1+2*term2+4*term3+term4); 

  return ans;
} 

//積分の範囲[a,b]とその範囲の分割数num
double daikei_method(double (*function)(double x),double a,double b,int num){ 

  double h = (b-a)/(float)num;
  float ans=0;

  for(int i=1;i<=num;i++){
     ans=ans+(*function)(a+(i-1)*h)+(*function)(a+(i*h)) ;
  }

  return ans=h/2.*ans;
}

double func(double x){
  double y;
  y=4.0/(1.0+x*x);
  return y;
}

Matlabコード

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% シンプソン法
% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
echo off
clear all
close all

load 'simpsonAndDaikei_sekibun.txt'

hi=simpsonAndDaikei_sekibun(:,1);

iti=ones(size(hi));
sekibun_real=iti.*3.14159;
daikei=simpsonAndDaikei_sekibun(:,2);
simpson=simpsonAndDaikei_sekibun(:,3);

semilogx(hi,sekibun_real,'r-','linewidth',2);
hold on
semilogx(hi,daikei,'bo-','linewidth',2);
hold on
semilogx(hi,simpson,'go-','linewidth',2);
hold on
grid on

%x軸を逆にする
set(gca,'XDir','reverse');

set(gca,'XTick',[0.01 0.1 0.5])
set(gca,'XTickLabel',{'0.01','0.1','0.5'})

xlabel('微小区間(h)','Fontsize',18)
ylabel('積分値','Fontsize',18)

h=legend('積分(理論)',...,
'積分(台形法)',...,
'積分(シンプソン法)',...,
4);
set(h,'FontSize',15)

axis([0.01 0.5 3.100 3.145]);
h=gca
get(h)
set(h,'LineWidth',2,...,
'FontSize',18)

print -djpeg suuti_sekibun_simpson_daikei.jpg
Generated by ひとりWiki version 1.15.1