×

[PR]この広告は3ヶ月以上更新がないため表示されています。
ホームページを更新後24時間以内に表示されなくなります。

数値微分と数値積分

関数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)
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型で宣言した。


参考:岩波出版 「情報処理入門コース7 数値計算」戸川隼人著
goto index