私の理解では、精度を上げるには、テイラー級数のより多くの項を考慮する必要があります。たとえば、テイラー級数で e(1) を計算しようとするとどうなるかを考えてみましょう。
$e(x) = \sum\limits_{n=0}^{\infty} frac{x^n}{n!}$
e(1) の展開の最初のいくつかの項を考慮することができます。
n value of nth term sum
0 x^0/0! = 1 1
1 x^1/1! = 1 2
2 x^2/2! = 0.5 2.5
3 x^3/3! = 0.16667 2.66667
4 x^4/4! = 0.04167 2.70834
最初に、より多くの項を追加するにつれて、e(1) の正確な値に近づいていることと、連続する合計の差が小さくなっていることに注意してください。
したがって、e(x) の実装は次のように記述できます。
#include <stdbool.h>
#include <stdio.h>
#include <math.h>
typedef float (*term)(int, int);
float evalSum(int, int, int, term);
float expTerm(int, int);
int fact(int);
int mypow(int, int);
bool sgn(float);
const int maxTerm = 10; // number of terms to evaluate in series
const float epsilon = 0.001; // the accepted error
int main(void)
{
// change these values to modify the range and increment
float start = -2;
float end = 2;
float inc = 1;
for(int x = start; x <= end; x += inc)
{
float value = 0;
float prev = 0;
for(int ndx = 0; ndx < maxTerm; ndx++)
{
value = evalSum(0, ndx, x, expTerm);
float diff = fabs(value-prev);
if((sgn(value) && sgn(prev)) && (diff < epsilon))
break;
else
prev = value;
}
printf("the approximate value of exp(%d) is %f\n", x, value);
}
return 0;
}
私は、目的の精度を得るために展開で 10 を超える項を使用する必要がないことを推測として使用しました。したがって、内側の for ループはn、範囲 [0,10] の値をループする場所です。
また、必要な精度に達しているかどうかを確認する専用の行がいくつかあります。まず、今回の評価と前回の評価の差の絶対値を計算し、差の絶対値をとります。差がイプシロン値 (1E-3) より小さいかどうかを確認することは、ループを早期に終了する基準の 1 つです。また、e(-1) の値を計算する際の変動により、現在の値と前の値の符号が同じであることを確認する必要がありました。これは、条件の最初の節が行っていることです。
float evalSum(int start, int end, int val, term fnct)
{
float sum = 0;
for(int n = start; n <= end; n++)
{
sum += fnct(n, val);
}
return sum;
}
これは、級数の最初の n 項を評価するために私が書いた効用関数です。 startは開始値 (このコードは常に 0) であり、endは終了値です。最後のパラメーターは、特定の項の計算方法を表す関数へのポインターです。このコードでは、fnct は、整数パラメーターを受け取り、float を返す任意の関数へのポインターにすることができます。
float expTerm(int n, int x)
{
return (float)mypow(x,n)/(float)fact(n);
}
この 1 行の関数に埋もれている場所で、ほとんどの作業が行われます。この関数は、e(n) のテイラー展開の閉形式を表します。上記を注意深く見ると、与えられた x と n の値に対して $\fract{x^n}{n!}$ を計算していることがわかります。ヒントとして、コサイン部分を実行するには、cos のテイラー展開の項の閉を評価する関数を作成する必要があります。これは $(-1)^n\fact{x^{2n}}{(2n)!}$ で与えられます。
int fact(int n)
{
if(0 == n)
return 1; // by defination
else if(1 == n)
return 1;
else
return n*fact(n-1);
}
これは階乗関数の標準的な実装です。ここで見る特別なものは何もありません。
int mypow(int base, int exp)
{
int result = 1;
while(exp)
{
if(exp&1) // b&1 quick check for odd power
{
result *= base;
}
exp >>=1; // exp >>= 1 quick division by 2
base *= base;
}
return result;
}
累乗を行うためのカスタム関数。確かに のバージョンを使用することも<math.h>できましたが、整数のべき乗のみを行うことがわかっていたので、最適化されたバージョンを作成できました。<math.h>ヒント: cosine を実行する場合、浮動小数点ベースを扱うには、おそらく からのバージョンを使用する必要があります。
bool sgn(float x)
{
if(x < 0) return false;
else return true;
}
浮動小数点値の符号を決定するための信じられないほど単純な関数で、true を返すと正で、それ以外の場合は false を返します。
このコードは、gcc バージョン 4.8.4 を使用して Ubuntu-14.04 でコンパイルされました。
******@crossbow:~/personal/projects$ gcc -std=c99 -pedantic -Wall series.c -o series
******@crossbow:~/personal/projects$ ./series
the approximate value of exp(-2) is 0.135097
the approximate value of exp(-1) is 0.367857
the approximate value of exp(0) is 1.000000
the approximate value of exp(1) is 2.718254
the approximate value of exp(2) is 7.388713
bcを使用して指定される期待値は次のとおりです。
******@crossbow:~$ bc -l
bc 1.06.95
Copyright 1991-1994, 1997, 1998, 2000, 2004, 2006 Free Software Foundation, Inc.
This is free software with ABSOLUTELY NO WARRANTY.
For details type `warranty'.
e(-2)
.13533528323661269189
e(-1)
.36787944117144232159
e(0)
1.00000000000000000000
e(1)
2.71828182845904523536
e(2)
7.38905609893065022723
ご覧のとおり、値は要求した許容範囲内に収まっています。コサイン部分を行うための演習として残します。
これが役に立てば幸いです、
-T