微分@その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) の形で表されているときの積分値を求める

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

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

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

連続と不連続

微分や積分をコンピュータで解く場合、連続(アナログ)な現象を不連続(ディジタル)な現象として解く必要がある。

連続モデル　　　　不連続モデル
   微分　　　　　　　 差分
   積分　             和分

数値解法におけるある点xでの微分を求める方法は1つしかない。
微分方程式の解を求めることと、たんに微分を求めることとは区別すること。
ここでは微分の方法だけを示し、微分方程式の解を求めることは
オイラー、テイラー、ルンゲクッタ法で説明する。
ここでは数値微分における差分の方法をプログラム交えて説明する。

差分の定義

ある連続関数 y=f(x) , $a\lt x \lt b$ において、区間[a,b]をn等分し、その間隔を $h=\frac{b-a}{n}$ とする。
このように定義すると、第1差分と第2差分は次のように表される。

○第1差分
$\Delta y_0 = y_1 - y_0$
$\Delta y_1 = y_2 - y_1$

$\Delta y_{n-1} = y_n - y_{n-1}$

○第2差分
$\Delta^2 y_0 = \Delta y_1 - \Delta y_0$
$\Delta^2 y_1 = \Delta y_2 - \Delta y_1$

$\Delta^2 y_{n-2} = \Delta y_{n-1} - \Delta y_{n-2}$

この考え方を一般化すると、第n差分は以下のようになる。
$\Delta^n y_j = \Delta^{n-1} y_{j+1} - \Delta^{n-1}y_j$

ここで、
$y_{j+m} = K^m y_j$
となるオペレータKを導入すると、

$\Delta y_j = y_{j+1} - y_j = Ky_j - y_j = (K-1)y_j$
$\Delta^2 y_j = \Delta y_{j+1} - \Delta y_j$
$= (y_{j+2} - y_{j+1}) - (y_{j+1} - y_{j})$
$= y_{j+2} - 2 y_{j+1} + y_{j}$
$= (K^2 - 2K + 1)y_{j} = (K-1)^2 y_{j}$

$\Delta^n y_{j} = (K-1)^n y_{j}$
$=\[K^n - \Bigl(\begin{array}{GC+23}n\\1\end{array}\Bigr)K^{n-1}+\Bigl(\begin{array}{GC+23}n\\2\end{array}\Bigr)K^{n-2}+ \cdots +(-1)^n \] y_j$
$=y_{j+n}-\Bigl(\begin{array}{GC+23}n\\1\end{array}\Bigr)y_{j+n-1}+\Bigl(\begin{array}{GC+23}n\\2\end{array}\Bigr)y_{j+n-2}+ \cdots +(-1)^{n}y_{j}$

数値微分

微分オペレータ D=d/dx を導入して差分を表すと、

$Dy=y^{(1)}$
$D^{2}y=D \cdot Dy=Dy^{(1)}=y^{(2)}$
$\vdots$
$\vdots$
$D^{n}y=y^{(n)}$ (式1-1)

となる。ここで、 $h \ll 1$ に対してテーラー展開すると

$y(x+h) = y(x) + hy^{(1)}+\frac{h^2}{2!}y^{(2)}(x) + \cdots$
$= (1+hD+\frac{(hD)^2}{2!}+ \cdots)y(x)$
$= e^{hD} \cdot y(x)$ (マクローリン展開)

となる。ここで

$\Delta{y} = y(x+h) - y(x) = (e^{hD} -1)y(x)$

ゆえに

$\Delta = e^{hD} - 1$ 　(式1-2)

が導出される。

さらに、(式1-2)は、

$e^{hD} = 1+\Delta$ (式1-2)

だから、この式の両辺の $log_{e}$ を取ると、

$hD = ln(1+\Delta)$ (式1-3)

この式をテーラー展開すると、

$= \Delta - \frac{1}{2}\Delta^{2} + \frac{1}{3}\Delta^{3} - \cdots$
$(\frac{dy}{dx})_j = D^{n} y_j$
$= \frac{1}{h}(\Delta - \Delta^{2} + \frac{1}{3}\Delta^{3} - \cdots)y_j$
$= \frac{1}{h}\Delta y_j$
( $\Delta$ の階乗は非常に小さいので)
$= \frac{ y_{j+1} - y_{j} }{h}$

また、式(1-3)を2乗すると

$h^{2}D^{2} = (ln(1+\Delta))^2$

$\Delta^2 - \Delta^3 + \frac{11}{12}\Delta^{4} - \frac{5}{6}\Delta^5 + \frac{137}{180}\Delta^6 + \cdots$

となり、これを用いてyに関する2次の微分は次のように求められる。

$(\frac{d^{2}y}{dx^2})_j = D^{2}y = \frac{1}{h^2}\Bigl(\Delta^2 - \Delta^3 + \frac{11}{12}\Delta^{4} - \frac{5}{6}\Delta^5 + \frac{137}{180}\Delta^6 + \cdots)y_j$

$= \frac{\Delta^{2}y_j}{h^2}$

$= \frac{y_{j+2} - 2y_{j+1} + y_j }{h^2}$ (式1-4)

となる。

差分(微分)のまとめ

ある連続関数 y=f(x) , $a\lt x \lt b$ において、区間[a,b]をn等分し、その間隔を $h=\frac{b-a}{n}$ とすると
1階差分および2階差分は

$(\frac{dy}{dx})_j = \frac{ y_{j+1} - y_{j} }{h}$
$(\frac{d^{2}y}{dx^2})_j = \frac{y_{j+2} - 2y_{j+1} + y_j }{h^2}$

となる。
また、中心差分における1階差分は、

$(\frac{dy}{dx})_j = \frac{ y_{j+1} - y_{j-1} }{2h}$

となる。

[例題その1] 関数のにおける1次および2次微係数を求めよ。

手計算で
f'(x) = 3 x^2
f''(x) = 6x
なので、
x=3 のときは、
f'(3) = 3 x^2 = 27
f''(3) = 6x = 18

コード

■ コード

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

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

//y=x^3のx=2における微係数を求める
double func1(double x);
double bibun1a(double (*function)(double x),double x,double h);
double bibun1b(double (*function)(double x),double x,double h);
double bibun2a(double (*function)(double x),double x,double h);

int _tmain(int argc, _TCHAR* argv[])
{
  errno_t err;
  FILE *fp_output;
  double x  =3.;
  double diff1=0,ans1=0,diff2=0,ans2=0;

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

  for(float h=0.01;h<=1.;h=h+0.01){

     //1次微分その1
     double ans1=bibun1a(&func1,x,h);
     double ans1b=bibun1b(&func1,x,h);

     //2次微分
     double ans2=bibun2a(&func1,x,h);

     printf("x = %f\n",x);
     printf("微小区間(h)           = %f\n",h);
     printf("1階微分(前進差分)     = %f\n",ans1);
     printf("1階微分(中心差分)     = %f\n",ans1b);
     printf("2階微分(前進差分)     = %f\n",ans2);
     fprintf_s(fp_output,"%f %f %f %f\n",h,ans1,ans1b,ans2);
  }

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

double bibun1a(double (*function)(double x),double x,double h){

  double diff1a=(*function)(x+h)-(*function)(x);
  double ans=diff1a/h;
  return ans;
}

double bibun1b(double (*function)(double x),double x,double h){

  double diff1b=(*function)(x+h)-(*function)(x-h);
  double ans=diff1b/(2*h);
  return ans;
}


double bibun2a(double (*function)(double x),double x,double h){

  double diff2a=(*function)(x+2*h)-2.0*(*function)(x+h)+(*function)(x);
  double ans=diff2a/(h*h);
  return ans;
}

double func1(double x){
  double y;
  y=(x*x*x);
  return y;
}

■スクリプト

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 微分を実行する
% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
echo off
clear all
close all

load 'bibun_data.txt'

h=bibun_data(:,1);
b1a=bibun_data(:,2);
b1b=bibun_data(:,3);
b2a=bibun_data(:,4);

num=ones(size(h));
b1_real=num.*27.0;
b2_real=num.*18.0;

semilogx(h,b1a,'ro','linewidth',2);
hold on
semilogx(h,b1b,'mo','linewidth',2);
hold on
semilogx(h,b2a,'bo','linewidth',2);
hold on
semilogx(h,b1_real,'r-','linewidth',2);
hold on
semilogx(h,b2_real,'b-','linewidth',2);
hold on
grid on

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

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

xlabel('微小区間(h)','Fontsize',18)
ylabel('微係数','Fontsize',18)
title('x=3におけるx^3の1,2次微分','Fontsize',18) 

h=legend('1次微分(前進差分)',...,
'1次微分(中心差分)',...,
'2次微分(前進差分)',...,
'1次微分(理論)',...,
'2次微分(理論)',...,
3);
set(h,'FontSize',15)

axis([0.01 0.1 15 30]);
h=gca
get(h)
set(h,'LineWidth',2,...,
'FontSize',18)

print -djpeg suuti_bibun_byouga.jpg