数値計算@数値積分@その3@ガウスルジャンドル法

参考文献及び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) 微分方程式の解を数値計算で求める

積分@その1@積分とはなにか

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

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

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

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

はじめに

図1-1 積分する関数

関数 f'(x) の x_1 から x_2 までの積分は、

$Area = \int^{x_3}_{x_1}f'(x)dx$ (1-1)

で表されます。
ガウスルジャンドル法とは、左図の x_1 から x_2 の座標を
右図のように-1から1までの $\xi$ 座標に座標変換して積分する手法である。

座標変換

図1-2 座標変換

x座標軸と $\xi$ 座標軸の変換式は、

$x(\xi) = h\xi + x_2$ (1-2)
$dx = hd\xi$ (1-3)

となり、

図1-1 積分する関数

(1-1)式は、

$Area = \int^{x_3}_{x_1}f'(x)dx = \int^{+1}_{-1}f'(x(\xi))hd\xi$ (1-4)

となる。

(1-4)式は、 w_i (ウェイティングポイント)を使用して

$Area = \int^{+1}_{-1}f'(x(\xi))hd\xi \simeq h\sum^2_{i=1} f(\xi_i)w_i$ (1-5)

と表される。

サンプリングポイントが2つの場合、
サンプリングポイント $\xi_i$ とウェイティングポイント w_i は、

表1-1 2点ガウス・ルジャンドルの係数

i	\xi_i	w_i
1	-0.577350	1.000000
2	0.577350	1.000000

となる。
2点のときのサンプリングポイント、ウェイティングポイントの決め方は別途書籍を参考に
してください。
ここでは、どのような曲線でも座標変換することで-1から1までの積分に落とし込むことができて
そしてその積分の計算方法は、表1-1のサンプリングポイントでの曲線の値とウェイティングポイントを
かけて、足し合わせることで簡単に積分を計算することができるという点である。

さきほどの台形法、シンプソン法に比べて2点を計算するだけで積分が求められる点が非常に魅力的である。
精度を上げるためにはサンプリングポイントを増やしていけばよいが、その場合のサンプリングポイントに
おけるウェイティングポイントは決まっているので別途他の参考書籍を参考のこと。

手計算によるガウス　ルジャンドル法

手計算で行ってみる。

まず領域を

0から1までを領域①
1からpi/2までを領域②

とし、それぞれを積分して足し合わせることで積分を計算する。

■領域1の計算

領域1の全体長は1なので

h=1/2=0.5
x1=0
x2=0.5
x3=1

となる。

(a) ガウスルジャンドル法を適用して $\xi_1$ , $\xi_2$ を計算する。

$x(\xi_1) = h\xi_1 + x_2$
x(-0.577350) = 0.5*(-0.577350) + 0.5
x(-0.577350) = 0.211325

$x(\xi_2) = h\xi_2 + x_2$
x(0.577350) = 0.5*(0.577350) + 0.5
x(0.577350) = 0.7887

(b) (a)を元にsin(x)に代入して $f(x(\xi))$ を求める。

$f(x(\xi_1))=f(0.211325)=\sin(0.211325)=0.2098$

$f(x(\xi_2))=f(0.7887)=\sin(0.7887)=0.7094$

$f(x(\xi_1))w_1 + f(x(\xi_2))w_2 = 0.2098+0.7094=0.9192$
h*0.9192= 0.5*0.9192 = 0.4596

■領域2の計算

領域1は、

h= (pi/2 - 1)/2  = 0.5708/2 = 0.2854
x1=1
x2=1.2854
x3=1.5708

となる。

(a) ガウスルジャンドル法を適用して $\xi_1$ , $\xi_2$ を計算する。

$x(\xi_1) = h\xi_1 + x_2$
x(-0.577350) = 0.2854*(-0.577350) + 1.2854
x(-0.577350) = 1.1206

$x(\xi_2) = h\xi_2 + x_2$
x(0.577350) = 0.2854*(0.577350) + 1.2854
x(0.577350) = 1.4502

(b) (a)を元にsin(x)に代入して $f(x(\xi))$ を求める。

$f(x(\xi_1))=f(1.1206)=\sin(1.1206)=0.9004$

$f(x(\xi_2))=f(1.4502)=\sin(1.4502)=0.9927$

$f(x(\xi_1))w_1 + f(x(\xi_2))w_2 = 0.9004+0.9927=1.8931$
h*1.8931= 0.2854*1.8931 = 0.5403

■領域1と2を足す。

0.4596+0.5403=0.9999

Matlab コード

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ガウス・ルジャンドル法
% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
echo off
clear all
close all

y=[];
x=[]; 
for i=0:0.1:pi
  x=[x,i];
  y=[y,sin(i)];
end

area1x=[];
area1y=[];
for i=0:0.01:1
  area1x=[area1x,i];
  area1y=[area1y,sin(i)];
end

area2x=[];
area2y=[];
for i=1:0.01:pi/2
  area2x=[area2x,i];
  area2y=[area2y,sin(i)];
end

plot(x,y,'r-','linewidth',2);
hold on
stem(area1x,area1y,'bo','linewidth',2);
hold on
stem(area2x,area2y,'go','linewidth',2);
grid on

% xラベルとその文字の大きさ、線の太さの設定
str={'0','1','p/2','p'}
set(gca,'FontName','symbol','xtick',[0,1,pi/2,pi],'xticklabel',str)
set(gca,'LineWidth',2,'FontSize',18) 

% x-y範囲
axis([-Inf Inf -Inf Inf]); 

% xラベル、yラベル
xlabel('x','Fontsize',20,'FontName','Times');
ylabel('sin(x)','Fontsize',20,'FontName','Times')
h=legend('sin(x)',...,
1);
set(h,'FontSize',15,'FontName','Times')

print -djpeg sin_sekibun_gause_rugendol_byouga.jpg

[例題] 関数 $y = \frac{4}{1+x^2}$ を区間においてGauss-Legendre法を用いて解け。

3.141588 3.141593 3.147541

ソースコード

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

#include "stdafx.h"
#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 gauss_legendre_method(double (*function)(double x),double a,double b);
double func(double x); 
 

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

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


  double a=0.;
  double b=1.0;

  int n=200;
  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);
  ans_gauss_legendre = gauss_legendre_method(&func,a,b);
  fprintf_s(fp_output,"%f %f %f\n",ans_daikei,ans_sympson,ans_gauss_legendre);
  printf("%f %f %f\n",ans_daikei,ans_sympson,ans_gauss_legendre);

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


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

  double w1=1.;
  double w2=1.;
  double xi1=-0.577350;
  double xi2= 0.577350;
  double h=(b-a)/2.;

  double x1=a;
  double x2=a+h;
  double x3=b;

  double x_xi1=h*xi1 + x2;
  double x_xi2=h*xi2 + x2;

  double ans1=0.;
  double ans2=0.;
  double ans=0.;

  ans1=(*function)(x_xi1);
  ans2=(*function)(x_xi2);

  ans = h*(ans1*w1 + ans2*w2);

  return ans;
}


//積分の範囲[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;
  double 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;
}