2021/11/5投稿
0.メニュー
1.対数関数logと積分の関係
2.数値積分、近似公式
3.Cでの計算コード
4.計算結果
logについて
自然対数を $\log$ としその底は $e=2.7182818 \cdot\cdot\cdot $ とする
1.対数関数logと積分の関係
$\log(x)$の微分は次のように与えられる
\frac{d}{dx}\log(x) = \frac{1}{x}
従って次のような関係式が得られる
\log(a) = \int_{1}^{a} \frac{dx}{x}
上式の右辺を数値積分で計算することを考える
2.数値積分、近似公式
積分区間を離散化
\space\space\space x_k := x_0 + k\Delta x\\
図
シンプソンの公式
\int_{x_{k-1}}^{x_{k}} f(x) dx \simeq \frac{\Delta x}{6}
\Bigl(
f(x_{k-1})
+4f\Bigl( x_{k-1}+\frac{\Delta x}{2} \Bigl)
+f(x_k)
\Bigl) \\
上の公式を区間$x=[x_0,x_N]$に適用する
\int_{x_0}^{x_N} f(x) dx = \sum_{k=1}^N \int_{x_{k-1}}^{x_{k}} f(x) dx \\
\simeq
\sum_{k=1}^N \frac{\Delta x}{6}
\Bigl( f(x_{k-1})+4f\Bigl(x_{k-1}+\frac{\Delta x}{2}\Bigl)+f(x_k) \Bigl)
$f(x) = 1/x$を適用すると次のような近似公式を得る
自然対数の近似公式
x_k := 1 + k\Delta x \space\space\space N = \frac{a-1}{\Delta x}\\
\log(a) \simeq \sum_{k=1}^N \frac{\Delta x}{6}
\Bigl( \frac{1}{x_{k-1}} + \frac{4}{x_{k-1}+ \frac{\Delta x}{2}} + \frac{1}{x_k} \Bigl)
3.Cでの計算コード
mylog.c
#include <stdio.h>
#include <math.h>
// logの近似値を計算する
// xの整数n乗を返す関数
double zpow(double x, int n){
double re = 1.0;
int i;
for (i = 1; i <= n; i++){
re *= x;
}
return re;
}
// nの階乗
double fac(int n){
double re = 1.0;
int i;
for (i = 1; i <= n; i++){
re *= (double)i;
}
return re;
}
// 自作対数関数 シンプソン積分
double mylog(double a){
int N = 1000, k;
double re = 0.0, h = (a - 1.0)/((double)N), x_k, x_k_1;
for(k = 1; k <= N; k++){
x_k_1 = 1.0 + ((double)k - 1.0)*h;
x_k = 1.0 + ((double)k)*h;
re += (h/6.0)*(1.0/x_k_1 + 4.0/(x_k_1 + 0.5*h) + 1/x_k);
}
return re;
}
// main関数
int main(void){
// 変数宣言
double val, a = 30.0, mylog_val;
//
// 数値計算
mylog_val = mylog(a);
printf(" a := %lf \n", a);
printf(" log(a) = %.16f\n", mylog_val);
// 絶対誤差
printf(" Error: %.16f", fabs(log(a) - mylog_val));
//
return 0;
}
4.計算結果
今回は$\log(30)$の近似値を計算してみた
cのmathライブラリーの組み込み関数logで計算した値
との絶対誤差も計算してみたよ
精度上の注意
因みにこの記事には載せていないがほかの値で計算してみた結果 $\log(a)$の$a$の値を大きくするほど精度が悪くなる $a=30$で計算したが、この精度はあくまでも一例に過ぎない。
積分区間を100分割したとき
a := 30.000000
log(a) = 3.4012114387474641
Error: 0.0000140570853087
およそ4ケタぐらい一致
もっと分割数を増やしてみよう
積分区間を1000分割したとき
a := 30.000000
log(a) = 3.4011973831349214
Error: 0.0000000014727659
今度は8ケタぐらい一致!!
さらに分割数を増やしてみよう
積分区間を10000分割したとき
a := 30.000000
log(a) = 3.4011973816622989
Error: 0.0000000000001434
oh...12ケタも一致している!!!
ただし、計算に数秒くらいかかった。
それにしても組み込み関数のlogではどんなアルゴリズムで計算しているのか気になる......