自然対数と指数関数の値をテイラー展開で求める関数を自作する！

#概要
自然対数関数$\log_ex$と指数関数$e^x$の値をテイラー展開を用いてコンピュータ上で近似していきたい.
具体的にはmath.hで宣言されている関数とほぼ遜色ないレベルで近似できることが目標である.
以下の順序で話していく.

1. テイラー(マクローリン)展開の定義
2. 自然対数関数と指数関数のテイラー(マクローリン)展開式
3. とりあえずコードに落とし込んでみる（C/C++）
4. 問題点(ベキ級数の収束半径,収束速度)
5. 問題点の解決策
6. 解決策を踏まえたコード

３まで読み飛ばしてもよい.

• 想定読者：高校数学がある程度理解できていれば何とかなると思う.(論理ギャップは埋めた)
• 参考図書：『共立講座 21世紀の数学 微分積分』(黒田成俊)

#0.級数とは
まず,級数という言葉だけ確認しておく.
有限個の数,$a_0,a_1,a_2,\cdots,a_n$の和を

a_0+a_1+a_2+\cdots+a_n=\sum_{k=0}^n a_k

とかく.数列{$a_n$}に対して,その形式的な和$a_0+a_1+a_2+\cdots+a_n+\cdots$を級数といい,

a_0+a_1+a_2+\cdots+a_n+\cdots=\sum_{k=0}^\infty a_k

とかく.$\sum_{k=0}^\infty a_k$が意味をもつのは,収束するときだけであることを忘れてはならない.

\sum_{k=0}^\infty a_k=\lim_{n\to\infty}\sum_{k=0}^n a_k

である.つまり$n$項までの和を計算してから$n\to\infty$とすればよい.例えば,

\begin{eqnarray}
\sum_{k=0}^\infty k &=& \lim_{n\to\infty}\sum_{k=0}^n k\\
&=&\lim_{n\to\infty}\frac{1}{2}n(n+1)=\infty\\
\end{eqnarray}

$\sum_{k=0}^\infty k$は意味をもたないことが分かる.
#1.テイラー展開とは
###定義

関数$f(x)$は区間$J$で定義されているとし、$a\in J$ とする.いま, $a\in I$ $\subset J$を満たすある開区間$I$において,$f(x)$が$x-a$のベキ級数で表されているとする.つまり

f(x)=\sum_{n=0}^\infty c_n(x-a)^n,x\in I. \\

そのとき$f(x)$は$I$上で無限回微分可能で,

c_n=\frac{f^{(n)}(a)}{n!}\\

となる(証明略).したがって

f(x)=\sum_{n=0}^\infty \frac{f^{(n)}(a)}{n!}(x-a)^n,x\in I \tag{1}\\

である.$(1)$式を$f(x)$の$x=a$におけるテイラー展開という.とくに$a=0$の場合,すなわち$x=0$におけるテイラー展開をマクローリン展開という.ここで,$f^{(n)}(x)$は$f(x)$を$n$回微分したものとする.　$f^{(1)}(x)=f'(x),f^{(2)}(x)=f''(x)$である.
関数$f(x)$が$a$を含むある開区間で無限回微分可能ならば,形式的には$(1)$の右辺の無限級数が考えられる.この無限級数をテイラー級数という.

$f(x)$が$x=a$でテイラー展開できるとは,

1. $x=a$におけるテイラー級数の収束半径が0でない.
2. $x=a$の近くで$f(x)$とテイラー級数が一致する.
   の２条件が成り立つことである.

この項で扱う自然対数関数と指数関数は無限回微分可能で,自然対数関数については$x=1$でテイラー展開可能.指数関数についてはマクローリン展開可能なのでこの条件は今は特に気にする必要はない.
###テイラー展開のポイント
テイラー級数が意味をもてば,よく分からない関数値が実数係数多項式という非常に簡単なもので近似出来るのが,テイラー展開の醍醐味である.
かなり野暮ったい級数が出てきたが,まずやることはひたすら微分することである.(暴論)

#2.自然対数関数と指数関数の級数展開式
###Ⅰ.指数関数
まずは指数関数の微分について考えよう.($(e^x)'=e^x$が分かっているなら,中段の本題まで読み飛ばしてよい)
指数関数$e^x$は全ての$x\in\mathbb{R}$で微分可能である. 実際,$f(x)=e^x$とおくと,$f(x)$は全ての実数で定義されていて,微分の定義より

\begin{eqnarray}
f'(x)&=&\lim_{h\to0}\frac{f(x+h)-f(x)}{h}\\
&=&\lim_{h\to0}\frac{e^{x+h}-e^x}{h}\\
&=&e^x\lim_{h\to0}\frac{e^h-1}{h}\\
\end{eqnarray}

$p(h)=e^h-1,q(h)=h$ とおくと $h=0$で$p(0),q(0)$は明らかに連続で, $p(0)=q(0)=0$ , $q'(0)=1$より ロピタルの法則 が使えて

\lim_{h\to0}\frac{e^h-1}{h}=\lim_{h\to0}\frac{p(h)}{q(h)}=\frac{p'(0)}{q'(0)}=\frac{e^0}{1}=1.

が成り立つ.したがって$f'(x)=e^x$を得た.
$f(x)=f'(x)=e^x$だからすべての$k\in\mathbb{N}$に対して $f^{(k)}=e^x$が成り立つ.すなわち何回微分しても$e^x$である.

ここからが本題である.$e^x$の$x=0$でのテイラー(マクローリン)展開を考えよう.$0!=1 , f^{(0)}(x)=f(x)$に注意して$(1)$式に$f(x)=e^x,a=0$を代入すると

\begin{eqnarray}
e^x&=&\sum_{n=0}^\infty\frac{(e^0)^{(n)}}{n!}x^n\\
&=&\frac{e^0}{0!}+\frac{(e^0)^{(1)}}{1!}x+\frac{(e^0)^{(2)}}{2!}x^2+\frac{(e^0)^{(3)}}{3!}x^3+\cdots+\frac{(e^0)^{(k)}}{k!}x^k+\cdots\\
&=&1+x+\frac{x^2}{2!}+\frac{x^3}{3!}+\cdots+\frac{x^k}{k!}+\cdots\\
&=&\sum_{n=0}^\infty\frac{x^n}{n!}\tag{2}
\end{eqnarray}

を得る.
###Ⅱ.自然対数関数
次に自然対数関数の微分について考えよう.自然対数関数は指数関数の逆関数として定義される.つまり$e^{\log_ex}=\log_ee^x=x$である.
まずは$\log_ex$の基本性質から見ていく. $a,b>0$で

\begin{eqnarray}
\log_ea^p&=&p\log_ea\tag{a}\\
\log_e\frac{1}{a}&=&-\log_ea\tag{b}\\
\log_eab&=&\log_ea+\log_eb\tag{c}\\
\log_e\frac{a}{b}&=&\log_ea-\log_eb\tag{d}\\
\log_e1&=&0\tag{e}\\
\end{eqnarray}

(証明)
自然対数関数が指数関数の逆関数であることを用いて証明する.

$(a):$ $\log_ea=t$とおく.両辺に指数関数をとって,$a=e^t$.両辺を$p$乗して$a^p=e^{tp}$.両辺に対数を取れば,$\log_ea^p=tp=p\log_ea$を得る.

$(b):$ $(a)$で$p=-1$ととればよい.

$(c):$ $\log_ea=s,\log_eb=t$とおく.この２式それぞれの両辺に指数関数をとって,$a=e^s,b=e^t$となる.$ab=e^se^t=e^{s+t}$この式の一番左と一番右の辺に自然対数をとって,$\log_eab=s+t=\log_ea+\log_eb$を得る.

$(d):$ まず$\log_e\frac{a}{b}=\log_eab^{-1}$である. $(c)$から,$\log_eab^{-1}=\log_ea+\log_eb^{-1}.$また,$(a)$より$\log_ea+\log_eb^{-1}=\log_ea-\log_eb$を得る.

$(e):$ $(a)$において$p=0$とすると$\log_ea^0=\log_e1=0$を得る.

(証明終)

対数関数$\log_ex$は$x>0$で微分可能である. 実際,$f(x)=\log_ex$とおくと,$f(x)$は$x>0$で定義されていて,微分の定義より

\begin{eqnarray}
f'(x)&=&\lim_{h\to0}\frac{f(x+h)-f(x)}{h}\\
&=&\lim_{h\to0}\frac{\log_e(x+h)-\log_ex}{h}\\
&=&\lim_{h\to0}\frac{1}{h}\log_e(1+\frac{h}{x})\\
&=&\lim_{h\to0}\frac{1}{h}\log_e(1+\frac{h}{x})^{\frac{h}{x}\frac{x}{h}}\\
&=&\lim_{h\to0}\frac{1}{h}\frac{h}{x}\log_e(1+\frac{h}{x})^{\frac{x}{h}}\\
&=&\lim_{h\to0}\frac{1}{x}\log_e(1+\frac{h}{x})^{\frac{x}{h}}\\
&=&\frac{1}{x}
\end{eqnarray}

ここで最後の等式は,$e=\lim_{n\to\infty}(1+\frac{1}{n})^n$を用いた.

$f''(x)=-\frac{1}{x^2},f^{(3)}(x)=\frac{2}{x^3},f^{(4)}(x)=-\frac{6}{x^4},f^{(5)}(x)=\frac{24}{x^5},\cdots,f^{(n)}(x)=(-1)^{n-1}\frac{(n-1)!}{x^n},\cdots$
となる.したがって,$\log_ex$の$x=1$におけるテイラー展開は,

\begin{eqnarray}
\log_ex&=&\sum_{n=0}^\infty\frac{f^{(n)}(1)}{n!}(x-1)^n\\
&=&\log_e1+(x-1)-\frac{(x-1)^2}{2}+\frac{(x-1)^3}{3}-\frac{(x-1)^4}{4}+\frac{(x-1)^5}{5}+\cdots\\
&=&\sum_{n=1}^\infty(-1)^{n-1}\frac{(x-1)^n}{n}\tag{3}
\end{eqnarray}

となる.
#3.とりあえずコードにする
級数展開出来ればあとは任意項まで計算することで関数値を近似することができる(かもしれない).
まずは指数関数の方から.$(2)$式を使って,

exp.cpp

double myexp(double x)
{
    double total = 1.0;    //級数和
    double a = 1.0;        //各項
    for (int i = 1; i < 20; i++)
    {
        a = a / i * x; //前の項に(x/i)をかければ次の項になる.
        total += a;    //現在の項を級数和として加える
    }
    return total;
}

次に対数関数.$(3)$式を使って,

log.cpp

double mylog(double x)
{
    if (x <= 0.0) { return 0; }    //入力が0以下なら関数未定義なので0を返すことにする.
    double total = 0.0;    //級数和
    double a = 1.0;        //分子(x-1)^nを保存
    for (int i = 1; i < 20; i++)
    {
        a = a * (x - 1.0);    //各項の分子部分を計算
        /* 奇数項のとき(+) */
        if (i % 2 == 1)
            total += a / i;   
        /* 偶数項のとき(-) */
        else
            total -= a / i;
    }
    return total;
}

########
#4.問題点(収束半径,収束速度)
上記のコードでは指数関数,対数関数ともにxが少しでも大きくなると収束値が全く異なる値になってしまう問題が起きる.
実際に実行してみると,

main.cpp

#include<iostream>
int main()
{
    std::cout << "exp3.0 = "  << myexp(3.0)  << std::endl;
    std::cout << "log1.5 = "  << mylog(1.5)  << std::endl;
    std::cout << "exp12.0 = " << myexp(12.0) << std::endl;
    std::cout << "log4.0 = "  << mylog(4.0)  << std::endl;
    return 0;
}

実行結果 （括弧内はmath.hのexp,log関数による計算結果）

exp3.0  = 20.0855     //(20.0855)
log1.5  = 0.405465    //(0.405465)
exp12.0 = 159291      //(162755)
log4.0  = 4.52595e+07 //(1.38629)

特に対数関数の方は値の差が顕著であることが分かる.
これは級数の収束と関わりがある.

###ベキ級数
定義

\sum_{n=0}^\infty c_nx^n=c_0+c_1x+c_2x^2+c_3x^3+\cdots+c_nx^n+\cdots

という形の級数を$x$のベキ級数という.ベキ級数がある$x_0\neq0$で発散すれば,$|x|>|x_0|$を満たす全ての$x$で発散する.(証明略)
ベキ級数の収束に関して次の3つの場合しか起こらない.

$(1):\sum_{n=0}^\infty c_nx^nは,x=0でしか収束しない.$
$(2):あるr,(0 < r <\infty)が存在して,\sum_{n=0}^\infty c_nx^nは,|x| < rならば収束し,|x| > rならば発散する.$
$(3):\sum_{n=0}^\infty c_nx^nはすべてのxで収束する.$
(2)の$r$のことを,ベキ級数$\sum_{n=0}^\infty c_nx^n$の収束半径という.(1)のとき収束半径は0, (3)のとき収束半径は$\infty$とする.

テイラー級数はベキ級数であり,$e^x$と$\log_ex$の収束半径はそれぞれ$\infty$,$|x-1|<1$である.$e^x$はすべての実数でテイラー級数が収束するが,$\log_ex$は$0 < x < 2$でしか収束しないのでより考察する必要がある.
また$|x|$が大きくなれば級数の各項の値も大きくなり,収束速度が落ちる.$\log_ex$に関しても$x=0$の十分近くでは同様に級数の収束速度が落ちる.
#5.問題点の解決策
まず問題点を整理しよう.
###問題点

1. $\log_ex$は$x > 2$のとき級数展開が使えない.
2. 級数の収束速度.極端な値をとると多くの項を計算しなければ正確な近似値が得られない.

###解決策
１：いま,$x>2$とする.対数関数の基本性質$(b)$を使うことによって次のようにできる.

\begin{eqnarray}
\log_ex &=& \log_e\frac{1}{\frac{1}{x}}\\
&=&\log_e1-\log_e\frac{1}{x}\\
&=&-\log_e\frac{1}{x}
\end{eqnarray}

$x > 2$であったから$\frac{1}{x} < \frac{1}{2}$.これにより$x>2$のときは$-\log_e\frac{1}{x}$のテイラー展開を考えればよいことが分かった.

２：級数が安定して収束する$x$を探してみる.(ここでの安定とは,20項目までの計算で級数の値がほぼ変化しなくなることを指す.)
結果,$e^2,e^{-2},\log_e0.5$は安定していることが分かった.分割統治法のような実装にする.

まず指数関数から考えよう.$x>2$とする.

\begin{eqnarray}
e^x&=&e^{(x-2)}e^2
\end{eqnarray}

であることを使う.$x<-2$のときも$e^{-2}$を考えて符号を変えて同様にできる.

対数関数は$x<0.5$とする.対数関数の基本性質$(C)$より,

\begin{eqnarray}
\log_ex &=& \log_e(0.5\times2x)\\
&=&\log_e0.5+\log_e2x
\end{eqnarray}

であることを使う.
#6.解決策を踏まえたコード

exp.cpp

double myexpv2(double x)
{
	double total = 1.0;
	double a = 1.0;
	/*  x >= 5以上のとき級数の収束速度が落ちるので、
	*	e^x = e^2 * e^(x-2) を考える.
	*/
	if (x >= 5.0)
	{
		total = myexpv2(2.0) * myexpv2(x - 2.0);	//再起的に計算(|x|を小さくして安定範囲までもっていく)
	}
	else if(x <= -5.0)
	{
		total = myexpv2(-2.0) * myexpv2(x + 2.0);	//負でも同様
	}
	else
	{
		/* -5 < x < 5のときは級数展開で値を求める. */
		for (int i = 1; i < 20; i++)
		{
			a = a / i * x;
			total += a;
		}
	}
	return total;
}

log.cpp

double mylogv2(double x)
{
	if (x <= 0.0) { return 0; }
	double a = 1.0;
	double total = 0.0;
	/*	交代級数の収束範囲が(0 < x < 2)なのでx>=2のとき
	*	logx = log1/(1/x) = log1-log1/xを使って0 < x < 2の範囲に落とし込む.(0 < 1/x < 2)
	*/
	if(x >= 2.0)
	{
		total = -mylogv2(1.0 / x);
	}
	else
	{
		/*	x < 0.5のとき収束速度が落ちるので、
		*	logx = log(0.5 * x / 0.5) = log0.5 + log(x / 0.5)として安定しているlog0.5を利用しつつ第２項を大きくしていく.
		*/
		if (x < 0.5)
		{
			total = mylogv2(0.5) + mylogv2(x / 0.5);
		}
		else
		{
			//x = 1の周りのテイラー展開を考える
			for (int i = 1; i < 20; i++)
			{
				a = a * (x - 1.0);
				/* 奇数項のとき+ */
				if (i % 2 == 1)
					total += a / i;
				/* 偶数項のとき- */
				else
					total -= a / i;
			}
		}
	}
	return total;
}

main.cpp

#include<iostream>
int main()
{
    std::cout << "exp3.0 = "  << myexpv2(3.0)  << std::endl;
    std::cout << "log1.5 = "  << mylogv2(1.5)  << std::endl;
    std::cout << "exp12.0 = " << myexpv2(12.0) << std::endl;
    std::cout << "log4.0 = "  << mylogv2(4.0)  << std::endl;
    //追加実験
    std::cout << "exp200.0 = " << myexpv2(200.0) << std::endl;
    std::cout << "exp-120.0 = " << myexpv2(-120.0) << std::endl;
    std::cout << "log400.0 = "  << mylogv2(400.0)  << std::endl;
    std::cout << "log0.0000001 = "  << mylogv2(0.0000001)  << std::endl;
    return 0;
}

実行結果 (括弧内はmath.hのlog,exp関数の結果)

exp3.0   = 20.0855      //(20.0855)
log1.5   = 0.405465     //(0.405465)
exp12.0  = 162755       //(162755)
log4.0   = 1.38629      //(1.38629)

exp200.0 = 7.22597e+86  //(7.22597e+86)
exp-120.0= 7.66749e-53  //(7.66765e-53)
log400.0 = 5.99146      //(5.99146)
log0.0000001 = -16.1181 //(-16.1181)

#終わりに
誤字脱字,間違っているところがあるかもしれないので気を付けてくださいね！