数値計算@数値積分@その1@台形法

参考文献及び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) の形で表されているときの積分値を求める

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

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

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

台形法

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

区間[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)

ここでいま、区間[a,b]内の一つの区間 $[x_i,x_{i+1}]$ を考えて、

その区間 (x_i,f(x_i)) と $(x_{i+1},f(x_{i+1}))$
の2点を結ぶ直線の方程式 g(x) を(1-3)式から求めると

$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)

これにより、区間 $[x_i,x_{i+1}]$ の積分値は、

$\int_{x_i}^{x_{i+1}}f(x)dx \approx \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+1}} _{x_i}$
$+ \frac{ f(x_{i+1})}{(x_{i+1} - x_i )}\Big[\frac{1}{2} x^2 - x_{i}x\Big]^{x_{i+1}} _{x_i}$

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

いま、
$x_{i+1} - x_{i} = h$ f(x_i) = y_i $f(x_{i+1}) = y_{i+1}$
とすると、

$\int_{x_i}^{x_{i+1}}f(x)dx \approx \frac{h}{2}(y_{i} + y_{i+1} )$ (1-5)

になるので、

これを一区間ではなく、区間[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-6)

となります。

台形法のコード

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

float f(float x); 

int _tmain(int argc, _TCHAR* argv[])
{
  errno_t err;
  FILE *fp_output;
  float a=1.0;
  float b=3.0;
  float ans=0;

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

  //printf("微小区間を代入してください\n");
  //printf("h << 0の値です\n");
  //scanf("%f",&h); //%fはfloatでないといけない
  //printf("h=%f\n",h);
  //rewind(stdin);

  //分割数を1から200まで分割する
  for(int n=1;n<=200;n++){
     float h = (b-a)/(float)n;
     ans=0;
      
     for(int i=1;i<=n;i++){
         ans=ans+f(a+(i-1)*h)+f(a+(i*h));
     }
     ans=h/2*ans;
     fprintf_s(fp_output,"%f %f\n",h,ans);
   }

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

float f(float x){
  float y;
  y=3.0*(x*x) + 10;
  return y;
}

Matlabのコード（結果描画）

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 3.3 数値積分
% written by embedded.samurai
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
echo off
clear all
close all

load 'data_sekibun.txt'

% comp1には行 comp2には列が入る
[comp1 comp2]=size(data_sekibun);

hi=data_sekibun(:,1);
iti=ones(size(hi));
sekibun_real=iti.*46.0;
sekibun_daikei=data_sekibun(:,2);
semilogx(hi,sekibun_real,'r-','linewidth',3);
hold on
semilogx(hi,sekibun_daikei,'m+-','linewidth',2);
hold on
grid on

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

h=legend('積分(理論)',...,
'積分(台形法)',...,
2);
set(h,'FontSize',18)

axis([0.01 2.0 45 50]);
h=gca
get(h)
set(h,'LineWidth',2,...,
'FontSize',18)

print -djpeg suuti_sekibun_daikei.jpg