微分を使わずテイラー展開

#### テイラー展開

テイラー展開の公式は次のとおりです。

```math
f(x)=\sum_{n=0}^{\infty}\frac{f^{(n)}(x_0)}{n!}(x-x_0)^n=\sum_{n=0}^{\infty}a_n(x-x_0)^n
```

　したがって、微分（$f^{(n)}$）を使わなければ各次の係数 $a_n$ は求まりそうもありません。微分といっても、展開すべき関数が $\sin x$、$\cos x$、$e^x$ 程度なら高次導関数もシンプルであり、あまり抵抗は感じません。しかし、例えば $\tan x$ の場合はどうでしょうか？

　　$(\tan x)'=\tan^2 x + 1 $
　　$(\tan x)''=2\tan x\,(\tan^2x + 1)$
　　$(\tan x)'''=2(\tan^2x + 1)^2 + 4\tan^2x\,(\tan^2x + 1)$
　　$(\tan x)''''=16\tan x\,(\tan^2x + 1)^2 + 8\tan^3x(\tan^2x + 1)$
　　　　　　　　　$\vdots$

　高次になるに従ってどんどん複雑になっていきます。とても手計算ではやっていられません。上記も自分で導出したものではなく、MATLAB の Symbolic Math Toolbox が教えてくれた答です。これが無ければ全くのお手上げです。係数 $a_n$ を求めるためのもっと簡単な方法はないものでしょうか？

#### まずは正攻法（微分を使う）による展開
```math
f_{\small{N}}(x)=\sum_{n=0}^{N}\frac{f^{(n)}(x_0)}{n!}(x-x_0)^n=\sum_{n=0}^{N}a_n(x-x_0)^n
```

　少々面倒ではありますが、まずは比較の基準を作るために、展開の中心 $x_0$ を $0.3$ 、最高次数 $N$ を $18$ として、微分を使う方法で $\tan x$ のテイラー展開 $f_{\small{N}}(x)$ を求めた結果が**図１**のピンク色の曲線です。青色の曲線は $\tan x$ の真値です。
 
![taylor_1_2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/932484/de98bd03-5999-7ff9-6f45-8a845a48e73a.png)

#### 連立方程式による係数の算出

　上の例のように次数 $n$ が $N$ までの有限の値であれば、各次の係数 $a_0$～$a_{\small{N}}$ は微分を使わなくても求まりそうな気もします。例えば、元の関数 $f(x)$ の$x_1$～$x_{\small{N}+1}$における値 $f(x_1)$～$f(x_{\small{N}+1})$ を使用して、$a_0$～$a_{\small{N}}$を未知数とした次の連立方程式を解けばどうでしょうか？

```math
\left[
\begin{array}{c}
\,1 & x_1{-}x_0 & (x_1{-}x_0)^2 & \cdots & (x_1{-}x_0)^{\small{N}}\\
\,1 & x_2{-}x_0 & (x_2{-}x_0)^2 & \cdots & (x_2{-}x_0)^{\small{N}}\\
\,1 & x_3{-}x_0 & (x_3{-}x_0)^2 & \cdots & (x_3{-}x_0)^{\small{N}}\\
\,\vdots & \vdots & \vdots &  & \vdots\\
\,1 & x_{\small{N}+1}{-}x_0 & (x_{\small{N}+1}{-}x_0)^2 & \cdots & (x_{\small{N}+1}{-}x_0)^{\small{N}}\\
\end{array}
\right]
\cdot
\left[\begin{array}{c}
a_0\strut\\
a_1\\
a_2\\
\vdots\\
a_{\small{N}}\strut
\end{array}
\right]=
\left[\begin{array}{c}
f(x_1)\strut\\
f(x_2)\\
f(x_3)\\
\vdots\\
f(x_{\small{N}+1})\strut
\end{array}
\right]
```

　しかし、$f(x){=}\tan x$、$x_1\cdots x_{\small{N}+1}{=}x_0{\mp}2.5/2$ として**やってみましたがダメでした**。求められた $a_0$～$a_{\small{N}}$ を使って描いた曲線が**図２**の緑色の線です。$f(x)$ の近似曲線にはなっていますが、テイラー展開によるカーブとは若干異なります。

#### 再挑戦：サンプリング点を局所に集中

　それもそのはず、テイラー展開の係数を求めるのであれば、前式の右辺は $f(x)$ ではなく、$f_{\small{N}}(x)$ としなければなりません。しかし、$f_{\small{N}}(x)$ は未知なので 仕方なく $f(x)$ で代用していたのです。

　$f(x){\neq}f_{\small{N}}(x)$ の程度は、$x$ が展開の中心 $x_0$ から離れるほど強くなるのですから、$x_1$～$x_{\small{N}+1}$ として $x_0$ の至近の点だけ選べばどうでしょうか？

しかし、$x_1\cdots x_{\small{N}+1}{=}x_0{\mp}0.1/2$ として**やってみましたがダメでした**。結果は**図３**の緑色の曲線です。図２よりさらに悪い方向に行っているように見えます。

![taylor_3_4.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/932484/51917b3a-9355-2683-45a1-7ea452155f21.png)

#### 再々挑戦：計算精度を上げる

　しかし、図３のズレの原因は、図２のそれとは異なるように思えます。上の方法では、$f(x_1)$～$f(x_{\small{N}+1})$ が殆ど同一の値になるため、丸め誤差の影響がモロに現れているのではないでしょうか。では、有効数字の桁数を増やしてみたらどうでしょうか？

　**やってみたところ**、遂に、正規のテイラー展開によるカーブと 代用処方によるカーブを**ピッタリと一致**させることができました。結果は**図４**の緑色の曲線です。MATLABの通常の有効桁数は10進16桁ですが、**10進28桁**まで増やすことでやっと解決しました。

#### まとめ

　計算精度さえ上げれば、微分を使わずにテイラー展開の係数が求められることが分かりました。しかし、MATLABの本体には有効桁数を調整する機能はなさそうです。これを実現するにはSymbolic Math Toolbox が必要です。結局、難しい関数を解析的に微分するにしても、微分を避けて連立方程式で解くにしても、Symbolic Math Toolboxを使うことになってしまいました。

　Symbolic Math Toolboxの機能として、数式のままでの解析は良く知られていますが、数値計算の有効数字の桁数の調整も見逃せないお役立ち機能です。

#### プログラム

　以上、$\tan x$ の例だけについて述べましたが、添付のリストでtype=1,2,3と変更して実行すれば、$\sin x$、$\cos x$、$e^x$ での結果も確認することができます。

```Taylor_exp_subst.m
% 【動作させるには Symbolic Math Toolbox が必要】
% 
% ■■■　微分を使わずにテイラー展開の係数を求めてみる実験　■■■
% 
% ①　まずは正攻法のテイラー展開で係数を求めてグラフを描く。
% ②　次に真値のサンプリングから連立方程式を作り
% 　　　それを解くことで係数を求めてグラフを描く。（しかし、期待外れ）
% ③　同上の方法を、サンプリングデータを局所に集中させたうえで実行。
% 　　　（しかし、またまた期待外れ。かえって悪化。）
% ④　計算精度を上げるために有効数字の桁数を上げたうえ、同様に実行。
% 　　　（めでたく期待どおりの結果）

clear
close all

syms x          % シンボリック変数の設定
                % 　このシンボルを含む記述は、数値ではなくシンボル式になる。

digits(16);     % 無指定の部分での計算精度は、MATLAB本体標準の10進16桁とする。
                % 　これを指定しておかないと、Symbolic Math Toolbox の規定値の
                % 　10進32桁となり、本体と同一レベルでの比較ができなくなる。

% テイラー展開する関数の選択
type=4;         % 1: sin,  2: cos,  3: exp,  4: tan

switch type     % 展開される関数の種類に応じて、グラフ映えする条件を設定。
  case 1                        % ■ sin ■
    fx=sin(x);                  % 展開される関数の式
    xx=linspace(-10,10,200);    % グラフ化する横軸の範囲
    x0=0.3;                     % テイラー展開の中心座標
    N=18;                       % テイラー展開の最高次数
    y_zone=[-8,8];              % グラフの縦軸範囲の制限
    ww=8;                       % 連立方程式のサンプリング点x1～x(N+1)を
                                % 　x0を中心とした 広めの幅ww内に均等配置する。
    wn=0.2;                     % 同じく、狭めの幅wn
    Q=36;                       % 高精度計算するときの10進数の有効桁数
  case 2                        % ■ cos ■
    fx=cos(x);
    xx=linspace(-10,10,200);
    x0=0.3;
    N=18;
    y_zone=[-8,8];
    ww=8;
    wn=0.2;
    Q=37;
  case 3                        % ■ exp ■
    fx=exp(x);
    xx=linspace(-15,10,200);
    x0=0.3;
    N=18;
    y_zone=[-4000,22000];
    ww=12;
    wn=0.2;
    Q=37;
  otherwise                     % ■ tan ■
    fx=tan(x);
    xx=linspace(-2,2,200);
    x0=0.3;
    N=18;
    y_zone=[-4,8];
    ww=2.5;
    wn=0.1;
    Q=28;
end

% 真値 f(x) のグラフ用のデータ

Real=double(subs(fx,xx));    % 関数の真値のグラフ用のデータ。
                             % 　シンボリック関数fxに横軸のデータxxを代入して数値化。
                             % 　Symbolic Math Toolbox の既定10進32桁は過剰品質なので
                             % 　doubleで、MATLAB標準の10進16桁に落とす。

% 正規テイラー展開結果 fN(x) のグラフ用のデータ
%        N
% fN(x)=Σ f^(n)(x0)/n! (x-x0)^n
%       n=0

F=fx;
a0=subs(F,x0);                    % テイラー展開の定数項（関数f(x)にx=x0を代入して求める）。
Tay=a0+0*xx;                      % そのグラフ用のデータ（これに高次の項を加えていく）。
disp(strcat('f(x)=', string(F)))  % 展開される関数の式をコマンドラインに表示。
for n=1:N
  F=diff(F);                                            % 各次の導関数の式を求める。
  disp(strcat('f(', num2str(n), ')(x)=', string(F)))    % 導関数の式をコマンドラインに表示。
  Tay=Tay+(1/prod(1:n))*double(subs(F,x0))*(xx-x0).^n;  % テイラー展開の各項を順次合成。
                                                        % 　prod(1:n)=n!
end

% ■■■■■　f(x)の真値のグラフと、正規テイラー展開結果fN(x)のグラフを描画

figure(1)
hold on
if string(fx)=='tan(x)'      % tan(x)の場合、Realに邪魔な線が出るのでそれを消す。
  V=[Real Real(length(Real))];      % Realの左端の要素を削除し、他の要素を１つずつ左に
  V(1)=[];                          % 　ずらし、右端にはそれと同一要素を追加してVとおく。
  nn=find((Real-V)>10);             % +∞ → -∞ の切り替わり点を探す。
  nnn=length(nn);
  if nnn~=0
    for k=1:nnn
      Real(nn(k))=NaN;
    end
  end
  % 消した代わりに、π/2±nπ の部分に細い縦線を描く。
  V=-10*pi-pi/2:pi:10*pi+pi/2;
  nn=find(V>=xx(1) & V<=xx(length(xx)));
  nnn=length(nn);
  if nnn~=0
    h0=plot([V(nn);V(nn)],y_zone','-b');
  end
end
h1=plot(xx,Real,'-b','LineWidth',6);               % f(x)の真値のグラフ
h2=plot(x0,a0,'ob','MarkerSize',10,'LineWidth',2); % 展開の中心点(x0,f(x0))を〇で表示。
h3=plot(xx,Tay,'-m','LineWidth',4);                % テイラー展開結果のfN(x)のグラフ
legend([h1 h2 h3],{'真値','展開の中心','正規テイラー展開'},'Location','north')
title(strcat(string(fx), ' の ', num2str(N), '次までのテイラー展開'))
ylim(y_zone);
grid on
AX=axis;
xlabel('x')
ylabel('f(x), f_N(x)')
title(strcat('f(x)=',string(fx),'の真値と',num2str(N),'次までの正規テイラー展開'))
ax=axis;
tx=ax(1)-(ax(2)-ax(1))*0.13;
ty=ax(4)+(ax(4)-ax(3))*0;
text(tx,ty,'図１','FontSize',20)

% ■■■■■　代用テイラー展開（広域サンプリング連立方程式）結果のグラフを追加

% 　下記を解いて係数a0～aNを求める。ただし、M=N+1。
% 
% 　f(x1) = a0 + a1*(x1-x0)^1 + a2*(x1-x0)^2 + a3*(x1-x0)^3 + ... + aN*(x1-x0)^N
% 　f(x2) = a0 + a1*(x2-x0)^1 + a2*(x2-x0)^2 + a3*(x2-x0)^3 + ... + aN*(x2-x0)^N
% 　f(x3) = a0 + a1*(x3-x0)^1 + a2*(x3-x0)^2 + a3*(x3-x0)^3 + ... + aN*(x3-x0)^N
% 　　｜　　｜　　　　｜　　　　　　｜　　　　　　　｜　　　　　　　　　｜
% 　f(xM) = a0 + a1*(xM-x0)^1 + a2*(xM-x0)^2 + a3*(xM-x0)^3 + ... + aN*(xM-x0)^N
% 
% 以下では、xs = [x1-x0; x2-x0; x3-x0; ... ; xM-x0] に相当する。

xs=linspace(-0.5,0.5,N+1)'*vpa(ww);  % vpaの機能：( )内の変数をdigitsで設定されて
                                     % 　　　　　　いる精度（10進の桁数）に置き換える。
                                     % 以降、この値を含む式の結果は全て同一の精度となる。
                                     % 　xs はもちろん、ff、A、aa なども同一精度となる。
                                     % この段階では、先に設定されたdigits(16)が有効。

ff=subs(fx,xs+x0);    % サンプリング点x1,x2,x3,...における fx の値

A=[];                 % 行列Aの作成
                      %     [ 1   (x1-x0)^1   (x1-x0)^2   (x1-x0)^3  --- (x1-x0)^N ]
                      %     [ 1   (x2-x0)^1   (x2-x0)^2   (x2-x0)^3  --- (x2-x0)^N ]
for n=0:N             %   A=[ 1   (x3-x0)^1   (x3-x0)^2   (x3-x0)^3  --- (x3-x0)^N ]
  A=[A xs.^n];        %     [ |       |           |           |              |     ]
end                   %     [ 1   (xM-x0)^1   (xM-x0)^2   (xM-x0)^3  --- (xM-x0)^N ]

aa=inv(A)*ff;         % 係数a0,a1,a2,...の算出

Subst=xx*0;           % 代用テイラー展開のカーブ（初期化）
for n=0:N             % 初期化カーブに高次の曲線を順次重ねていく。
  Subst = Subst + aa(n+1)*((xx-x0).^n);   % 注：　aa(1)=a0, aa(N+1)=aN
end

figure(2)
hold on
if string(fx)=='tan(x)'      % tan(x)の場合、π/2±nπ の部分に細い縦線を貼り付ける。
  h0o=copyobj(h0,gca);
end
h1o=copyobj(h1,gca);         % figure(1)の真値カーブを貼り付ける。
h2o=copyobj(h2,gca);         % 　〃　　 の展開の中心点を貼り付ける。
h3o=copyobj(h3,gca);         % 　〃　　 のfN(x)のカーブを貼り付ける。
axis(AX);
h4o=plot(xx,Subst,'-g','LineWidth',2);
legend([h1o h2o h3o h4o],{'真値','展開の中心','正規テイラー展開','代用展開'}, ...
                          'Location','north')
grid on
xlabel('x')
ylabel('f(x), f_N(x)')
title({strcat('f(x)=',string(fx),'の真値と',num2str(N),'次までの正規テイラー展開'); ...
       strcat('および　広域(Δx=',num2str(ww),')サンプル連立方程式による代用展開')})
text(tx,ty,'図２','FontSize',20)

% ■■■■■　代用テイラー展開（狭域サンプリング連立方程式・標準精度）結果のグラフを追加
% 上記の「広域サンプリング連立方程式」と殆ど同一の処理（wwがwnに変わるだけ）

xs=linspace(-0.5,0.5,N+1)'*vpa(wn);   % この段階では、先に設定されたdigits(16)が有効。

ff=subs(fx,xs+x0);    % サンプリング点x1,x2,x3,...における fx の値

A=[];                 % 行列Aの作成
                      %     [ 1   (x1-x0)^1   (x1-x0)^2   (x1-x0)^3  --- (x1-x0)^N ]
                      %     [ 1   (x2-x0)^1   (x2-x0)^2   (x2-x0)^3  --- (x2-x0)^N ]
for n=0:N             %   A=[ 1   (x3-x0)^1   (x3-x0)^2   (x3-x0)^3  --- (x3-x0)^N ]
  A=[A xs.^n];        %     [ |       |           |           |              |     ]
end                   %     [ 1   (xM-x0)^1   (xM-x0)^2   (xM-x0)^3  --- (xM-x0)^N ]

aa=inv(A)*ff;         % 係数a0,a1,a2,...の算出

Subst=xx*0;           % 代用テイラー展開のカーブ（初期化）
for n=0:N             % 初期化カーブに高次の曲線を順次重ねていく。
  Subst = Subst + aa(n+1)*((xx-x0).^n);   % 注：　aa(1)=a0, aa(N+1)=aN
end

figure(3)
hold on
if string(fx)=='tan(x)'      % tan(x)の場合、π/2±nπ の部分に細い縦線を貼り付ける。
  h0o=copyobj(h0,gca);
end
h1o=copyobj(h1,gca);         % figure(1)の真値カーブを貼り付ける。
h2o=copyobj(h2,gca);         % 　〃　　 の展開の中心点を貼り付ける。
h3o=copyobj(h3,gca);         % 　〃　　 のfN(x)のカーブを貼り付ける。
axis(AX);
h4o=plot(xx,Subst,'-g','LineWidth',2);
legend([h1o h2o h3o h4o],{'真値','展開の中心','正規テイラー展開','代用展開'}, ...
                          'Location','north')
grid on
xlabel('x')
ylabel('f(x), f_N(x)')
title({strcat('f(x)=',string(fx),'の真値と',num2str(N),'次までの正規テイラー展開'); ...
       strcat('および　狭域(Δx=',num2str(wn), ...
              ')サンプル連立方程式（標準精度）による代用展開')})
text(tx,ty,'図３','FontSize',20)

% ■■■■■　代用テイラー展開（狭域サンプリング連立方程式・高精度）結果のグラフを追加
% 上記の「広域サンプリング連立方程式・標準精度」と殆ど同一の処理（digits(Q)の追加のみ）

digits(Q)

xs=linspace(-0.5,0.5,N+1)'*vpa(wn);   % 新しく設定した計算精度を反映させる

ff=subs(fx,xs+x0);    % サンプリング点x1,x2,x3,...における fx の値

A=[];                 % 行列Aの作成
                      %     [ 1   (x1-x0)^1   (x1-x0)^2   (x1-x0)^3  --- (x1-x0)^N ]
                      %     [ 1   (x2-x0)^1   (x2-x0)^2   (x2-x0)^3  --- (x2-x0)^N ]
for n=0:N             %   A=[ 1   (x3-x0)^1   (x3-x0)^2   (x3-x0)^3  --- (x3-x0)^N ]
  A=[A xs.^n];        %     [ |       |           |           |              |     ]
end                   %     [ 1   (xM-x0)^1   (xM-x0)^2   (xM-x0)^3  --- (xM-x0)^N ]

aa=inv(A)*ff;         % 係数a0,a1,a2,...の算出

Subst=xx*0;           % 代用テイラー展開のカーブ（初期化）
for n=0:N             % 初期化カーブに高次の曲線を順次重ねていく。
  Subst = Subst + aa(n+1)*((xx-x0).^n);   % 注：　aa(1)=a0, aa(N+1)=aN
end

figure(4)
hold on
if string(fx)=='tan(x)'      % tan(x)の場合、π/2±nπ の部分に細い縦線を貼り付ける。
  h0o=copyobj(h0,gca);
end
h1o=copyobj(h1,gca);         % figure(1)の真値カーブを貼り付ける。
h2o=copyobj(h2,gca);         % 　〃　　 の展開の中心点を貼り付ける。
h3o=copyobj(h3,gca);         % 　〃　　 のfN(x)のカーブを貼り付ける。
axis(AX);
h4o=plot(xx,Subst,'-g','LineWidth',2);
legend([h1o h2o h3o h4o],{'真値','展開の中心','正規テイラー展開','代用展開'}, ...
                          'Location','north')
grid on
xlabel('x')
ylabel('f(x), f_N(x)')
title({strcat('f(x)=',string(fx),'の真値と',num2str(N),'次までの正規テイラー展開'); ...
       strcat('および　狭域(Δx=',num2str(wn), ...
              ')サンプル連立方程式（10進',num2str(Q),'桁精度）による代用展開')})
text(tx,ty,'図４','FontSize',20)
```