テイラー展開
テイラー展開の公式は次のとおりです。
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$ を求めるためのもっと簡単な方法はないものでしょうか？
まずは正攻法（微分を使う）による展開
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$ の真値です。
連立方程式による係数の算出
上の例のように次数 $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}}$を未知数とした次の連立方程式を解けばどうでしょうか？
\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$ としてやってみましたがダメでした。結果は図３の緑色の曲線です。図２よりさらに悪い方向に行っているように見えます。
再々挑戦：計算精度を上げる
しかし、図３のズレの原因は、図２のそれとは異なるように思えます。上の方法では、$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$ での結果も確認することができます。
% 【動作させるには 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)