数値微分と数値積分

関数f(x)のあるx=aでのf'(a),f"(a)を数値的に近似計算をする。

微分の定義は
　f'(a) = lim[h→0] (f(a+h) - f(a))/h
これを近似するのに
　lim[h→0] h = �凅
を代用すればよい、と考えるのが数値微分の基本的な発想である。

次のように前方差分、後方微分、中心微分を定義する。

前方差分(forward difference)
　�凉/�凅 = (f(x+�凅) - f(x))/�凅
後方差分または後退差分(backward difference)
　∇y/∇x = (f(x) - f(x-�凅))/�凅
中心差分(central difference)
　δy/δx = (f(x+�凅) - f(x-�凅))/2�凅
これは、�凉/�凅と∇y/∇xの平均となる。

テイラー展開をしてみる
　f(x+�凅) = f(x)+f'(x)�凅/1!+f"(x)(�凅)2/2!+ … +f(n)(�凅)n/n!+ …
これから
　{f(x+�凅)-f(x)}/�凅 = f'(x)/1!+f"(x)(�凅)1/2!+ … +f(n)(�凅)n-1/n!+ …
　∴�凅/�凉 = f'(x)+f"(x)(�凅)1/2!+ … +f(n)(�凅)n-1/n!+ …

同様に
　f(x-�凅) = f(x)-f'(x)�凅/1!+f"(x)(�凅)2/2!+ … -f(n)(�凅)n/n!+ …
これから
　{f(x)-f(x-�凅)}/�凅 = f'(x)/1!-f"(x)(�凅)1/2!+ … +f(n)(�凅)n-1/n!+ …
　∴∇x/∇y = f'(x)-f"(x)(�凅)1/2!+ … +f(n)(�凅)n-1/n!+ …


n階微係数の近似
上のテイラー展開した式を用いる。
読み飛ばした者は結果(太文字)だけでも見るように。
　�凉/�凅-∇y/∇x = f"(x)�凅+f""(x)(�凅)2/2･4!+ …
これを�凅で割れば
　(�凉/�凅-∇y/∇x)/�凅 = f"(x)+f""(x)�凅+ …
f""(x)�凅以上の項は小さいので無視すれば
(�凉/�凅-∇y/∇x)/�凅はf"(x)の近似とできる。

また、これを
　δ2y/δx2と表し、2階の中心差分という。

#include <stdio.h>
#include <math.h>
#define dx (0.1)	//�凅
/**
 * xの関数f(x)について
 * f"(a)の近似値を求める。
 * ----
 * ここでは
 * f(x)=x^3
 * f"(x)=6xのハズだが…
 */
float f(float x){
	return pow(x,3);	//f(x)=x^3
}
main(){
	float x,fpp;
	for(x=0;x<10;x++){
		fpp=(f(x+dx)-2*f(x)+f(x-dx))/dx/dx;
		printf("f\"(%1.0f) = %f\n",x,fpp);
	}
}
f(x)=x3という式について
f"(x)を求めている。
式は単純なのでプログラムに任せるまでもなく、
f"(x)=6xとわかる。
実際にプログラムで、
f"(0)、f"(1)、f"(2) … f"(9)
の計算結果を表示している。
0,6,12,18,,,54
と、出力されれば完全だが…

誤差は�凅(プログラムではdx)の二乗に比例する。
即ち、小さければ小さい程(理論上)精密になるが、それは桁数を
いくらでも取れる場合だけで、f(x+�凅)-f(x)の結果が0や無意味な
数値になってしまう。

一般の場合
δny/δxnについては(n∈Ｎ)
nが偶数の時
　δny/δxn = {(δn-2y/δ(x+�凅)n-2) - 2(δn-2y/δxn-2) + (δn-2y/δ(x-�凅)n-2)}/�凅2
nが奇数の時
　δny/δxn = {(δn-1y/δ(x+�凅)n-1) - (δn-1y/δ(x-�凅)n-1)}/2�凅
これをスターリング(Stirling)の方式という。


数値積分
定積分は次のように表せる
　∫[a〜b] f(x)･dx = lim[n→∞](b-a)/nΣ[k=1〜n]f(a+k(b-a)/n)
と書けてlim[n→∞]nの代わりに「十分大きなn」を使うことで近似する。

定積分は下の図の青色の斜線の面積を求めることになる。

xのa〜bをn等分し、端からx0(=a)、x1、x2、…xn-1、xn(=b)とする。
xkはa+k(b-a)/nである。
また、分割したxの幅を�凅とする。
　�凅 = (b-a)/n = xn-xn-1
よって、
　I = Σ�凅*{f(xn)+f(xn+1)}/2

●f(x)=4/(1+x2)のxの定義域[0〜1]でdxで積分する。

x=tanθとする。
　1/(1+x2)=cos2θ
　dx=1/cos2θ
であるから、
∫f(x)･dx = ∫1･dθ
　　　　　= π

#include <stdio.h>
/**
 * f(x) = 4/(1+x^2)の積分
 */
float f(float x){
	return 4/(1+x*x);
}
main(){
	float k,n,a,b;
	float sum=0.0,dx;
	
	printf("分割数 n=");scanf("%f",&n);
	a=0;b=1;
	dx=(b-a)/n;
	
	for(k=0;k<n;k++){
		sum += dx * (f(a+(b-a)*k/n)+f(a+(b-a)*(k+1)/n))/2;
	}
	printf("∫4/(1+x^2)･dx = %f",sum);
}
初めにnを指定する。
256とすると3.141590と表示され、大体πとなる。

ちなみに、k,n,a,bはint型でもいいような気がするが
dx=(b-a)/n
や(a+(b-a)*k/n)
が、勝手にintに縛られ、いちいち変数の前に(float)をつけるのは
面倒なので、初めにfloat型で宣言した。