小清水 (@curekoshimizu) です。

以前の記事に Excelで計算が破綻する例 を紹介しました。

紹介した数式は次でしたが、

$a=77617, b=33096$

$\displaystyle 333.75 b^6 + a^2 (11 a^2 b^2 - b^6 - 121 b^4 - 2) +5.5 b^8 + \frac{a}{2b}$

実に 人工的 な感じがします。

Abstract

今回紹介する話は 数学上自然な定義 なのに

それが うまくいかない という話になります。

お題は皆さん馴染みのある 三角函数 (sin) です。

三角函数の定義 といえばこの問題ですよね！

f:id:curekoshimizu:20161206083612p:plain:w450

‘99年東京大学入試問題。意表を突かれる入試問題です。

この問題への回答ではありませんが、

sin 函数の定義 ってなんでしょう。

大学以降、sin 函数は次で定義されることが多いです。

sin 函数の定義

$\displaystyle \mathrm{sin}(x) = x - \frac{x^3}{3!} + \frac{x^5}{5!} - \frac{x^7}{7!} +\cdots + (-1)^n \frac{x^{2n+1}}{(2n+1)!} + \cdots \quad (x \in \mathrm{R})$

この式は 0周りの Talylor 展開ですが、

x : 実数全体 で (広義一様) 収束する式になります。 (実は複素数全体まで拡張できる)

非常に非常に自然な定義 です！

sin 函数を定義通りにコンピュータで計算してみよう

この式を使って計算をしてみると、

例えば次のようなプログラムになるかと思います。

(剰余項に関する評価は次の節で)

#include <stdio.h>
#include <stdint.h>
#include <math.h>

double sin_taylor(double x)
{
    const double minus_xx = -x*x;

    const double term_tho = 1.0e-16;

    /* Taylor展開第1項までは直接計算 (x) */
    double term = x;
    double ret = x;

    /* Taylor展開第3項以降をループで計算 */
    for (uint32_t k = 1; ; k++) {

        uint32_t n = 2*k + 1;

        /* Taylor 展開第x^n=x^{2k+1}項の計算式: (-1)^k x^n/n! */
        term *= minus_xx / ( (n - 1.0) *n );

        if (fabs(term) <= term_tho) {
            break; /* 充分収束したことを示す条件 (次の節を参照) */
        }

        ret += term;
    } 

    return ret;
}

int main(int argc, char** argv)
{

    for (uint32_t ix = 0; ix <= 50; ix += 5) {

        double ret = sin_taylor((double)ix);

        printf("%d, %lf\n", ix, ret);
    }

    return 0;
}

結果をまとめると

x	sin 計算結果
0	0.000000
5	-0.958924
10	-0.544021
15	0.650288
20	-0.988004
25	-0.132351
30	-0.988004
35	-0.428035
40	1.069746
45	111.576759
50	-25547.371714

0 から離れると非常におかしい結果に変わっていきます。

実数上では $| \mathrm{sin}(x) | \leq 1$ 満たさないといけません。

しかし $\mathrm{sin}(40.0), \mathrm{sin}(45.0), \mathrm{sin}(50.0)$ は明らかにおかしいです。

評価を間違えたのでしょうか？

評価方法を数学的に考える

$x > 0$ : fixed としたとき、

Taylor の定理より

$\displaystyle \mathrm{sin}(x) = x - \frac{x^3}{3!} + \frac{x^5}{5!} - \frac{x^7}{7!} +\cdots + (-1)^{n-1} \frac{x^{2n-1}}{(2n-1)!} + (-1)^{n} \frac{\mathrm{sin}^{(2n+1)}(\theta)}{(2n+1)!} x^{2n+1}$

となる $\theta \in (0, x)$ が存在することになります。

ここでは記法として $\mathrm{sin}^{(2n+1)}$ を $\mathrm{sin}$ の $(2n+1)$ 階導函数としました。

$\displaystyle \left| \mathrm{sin}(x) - \left( x - \frac{x^3}{3!} + \frac{x^5}{5!} - \frac{x^7}{7!} +\cdots + (-1)^{n-1} \right) \right| = \left| (-1)^{n} \frac{\mathrm{sin}^{(2n+1)}(\theta)}{(2n+1)!} x^{2n+1} \right| \leq \frac{x^{2n+1}}{(2n+1)!}$