実数だけでローラン展開

f(z)=\sum_{n=-\infty}^{\infty}a_n(z-z_0)^n\tag{1}

a_n=\frac{1}{2\pi i}\oint_C \frac{f(z)}{(z-z_0)^{n+1}}dz\tag{2}

f(z)=\sum_{n=0}^{\infty}a_n(z-z_0)^n\tag{3}

f_{\small{N}}(x)=\sum_{n=-N}^{N}a_n(x-x_0)^n\tag{4}

% 【動作させるには Symbolic Math Toolbox が必要】
% 
% ■■■　実数だけでローラン展開してみる実験　■■■

% 「微分を使わずテイラー展開」（Taylor_exp_subst.m）の手法を適用。
% 　　　（サンプリングデータを局所に集中させて、　 ）
% 　　　（連立方程式を高精度で解いて係数a_nを求める）

% 正規のローラン展開の公式は下記。
% 
%        +∞
% f(z) =  Σ  a_n (z-z0)^n       ---- (1)
%       n=-∞
% 
%        1         f(z)
% a_n = ---- ∫------------ dz   ---- (2)
%       2πi   (z-z0)^(n+1)
% 
% しかし、先の代用テイラー展開において正規の公式による a_n を使わなかったのと同様に、
% ここでも(2)式は使わずに、連立方程式を高精度で解いて a_n を求めてみようと思う。
% また、ここでは変数は実数のみに限定するので、(1)式の z, z0 は x, x0 と置き換える。
% すると、使用する式は下式となる。
% 
%         +∞                 +N              
% f(x) =  Σ  a_n (x-x0)^n ≒ Σ  a_n (x-x0)^n       ---- (3)
%        n=-∞               n=-N             
% 
% 結局、下記のテイラー展開の公式の n=0 が n=-∞ に変わるだけの違いとなる。
% 
%         +∞                 +N              
% f(x) =  Σ  a_n (x-x0)^n ≒ Σ  a_n (x-x0)^n       ---- (4)
%         n=0                 n=0             
% 
% 代用テイラー展開の適用例においては N を 18 としていた。しかしここでは、展開に
% よって近似可能なxの範囲の限界も知りたいので、N をより∞に近い 50 とする。
% これにより、(3)式の未知数 a_n の個数は 2*50+1=101 個となり、101個のサンブリング
% 点と、101個の方程式が必要になる。
% 
% これらのサンプリング点は、x0 を中心とした狭い範囲（Δx=0.05）の中に、均等に分散
% させることにする。
% 
% 計算に使う有効数字の桁数は必要最小限に抑えたいところだが、N を 50 にして試して
% みると、10進100桁でも不足ぎみとなる。そこで安全サイドに大奮発して、10進200桁で
% 計算することにした。（そのため、本プログラムの結果が出るまで30秒近くかかる）
% 
% ①　まず、比較の基準を作るため、微分を使う正規のテイラー展開によって
% 　　f(x)=1/(x^2) の x0=0.5、N=50 での展開結果を求めて、そのグラフを描く。
% 
% ②　次に、同一の関数 f(x)=1/(x^2) を x0=0.5 として代用ローラン展開してみる。
% 　　当然、a_(-50)～a_50 のうちa_(-2) だけが 1.0、他は 0 になることが期待される。
% 　　
% 　　結果：
% 　　しかし、期待は外れて a_(-50)～a_(-1) は全て 0。 a_n としては n>=0 の項しか
% 　　残らない。これによるカーブは、正規のテイラー展開と完全に重なってしまう。
% 　　ローラン展開でせっかく導入した n<0 の項は何の役にも立っていない。
% 　　
% ③　今度は、Δx を、その中心に f(x)=1/(x^2) の特異点 x=0 がくるように再配置（す
% 　　なわち、x0=0）して、②と同様の計算をしてみる。
% 　　ただし、サンプリング点 x がぴったりと x0 と重なると x-x0=0 となってゼロ割エ
% 　　ラーが出るので、Δxの配置は、サンプリング間隔の半分（ 0.5*Δx/(2*N) ）だけ
% 　　ずらしておく。
% 　　
% 　　結果：
% 　　代用ローラン展開の結果は、全域が見事に真値と一致する。 すなわち、a_(-2) だけ
% 　　が 1.0、他は 0 になる。ローラン展開の面目躍如である。
% 　　
% ④　次は、もっと複雑な関数で、さらに特異点を２つ（x=0,-1）持つ
% 　　fx=1/(x^2)+0.2/((x+1)^2)+10*sin(10*x)+10 について、③と同様の計算をする。
% 　　
% 　　結果：
% 　　ここでも、代用ローラン展開は真値を良く近似できている。
% 　　さらに、近似可能なxの範囲の限界についても、次の性質があることが把握できる。
% 　　
% 　　特異点が無い関数（sinなど）では、N を大きくして行けばテイラー展開で近似可能
% 　　な x の範囲はどこまでも広がっていくように見えた。
% 　　しかし、特異点がある場合には、テイラー展開においても近似範囲には限界があり
% 　　そうなことが、①や②のグラフからも予想できる。
% 　　さらに、④の結果から、ローラン展開でも同様の限界があることが分かる。
% 　　①②④に共通の傾向として、近似できる x の幅は x0 ± |x0-(x0に至近の特異点)| 
% 　　となりそうである。|x0-(x0に至近の特異点)| は、複素数でのローラン展開の収束
% 　　半径に相当する。

clear
close all

syms x          % シンボリック変数の設定
                % 　このシンボルを含む記述は、数値ではなくシンボル式になる。
digits(200);    % 計算精度は大盤振る舞いの10進200桁とする。
N=50;           % 展開の次数は±50次とする。
wn=vpa(0.05);   % Δxの値（10進200桁精度）
xx=vpa(linspace(-2,2,200));  % グラフ化する横軸の範囲（10進200桁精度）

% ■■■■■■　ステップ１　■■■■■■
% 関数 fx=1/(x^2) の真値と x0=0.5 を中心とした第50次までの正規テイラー展開
% ■■■■■■■■■■■■■■■■■■■

% 真値のグラフ用のデータ
fx=1/(x^2);
Real=subs(fx,xx);

% 正規のテイラー展開のデータ
x0=vpa(0.5);                      % 展開の中心の x 座標（10進200桁精度）
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))*(subs(F,x0))*(xx-x0).^n;        % テイラー展開の各項を順次合成。
                                                        % 　prod(1:n)=n!
end

figure(1)
hold on
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)のグラフ
h5=plot([0 0],[0 50],'-r','LineWidth',1);          % 特異点の表示線
legend([h1 h2 h3 h5],{'真値','展開の中心','正規テイラー展開','特異点'}, ...
                      'Location','northeast')
title(strcat(string(fx), ' の ', num2str(N), '次までのテイラー展開'))
ylim([0 50]);
grid on
AX=axis;
xlabel('x')
ylabel('f(x)')
title(strcat('f(x)=',string(fx),'の真値と',num2str(N),'次までの正規テイラー展開'))
ax=axis;
tx=ax(1)-(ax(2)-ax(1))*0.15;
ty=ax(4)+(ax(4)-ax(3))*0;
text(tx,ty,'図１','FontSize',20)

% ■■■■■■　ステップ２　■■■■■■
% 関数 fx=1/(x^2) の真値と x0=0.5 を中心とした第50次までの正規テイラー展開
% および　x0=0.5 を中心とした ±50 次までの代用ローラン展開
% ■■■■■■■■■■■■■■■■■■■

% 代用ローラン展開のデータ
xs=linspace(-0.5,0.5,2*N+1)'*wn;
xs=xs+0.5*(xs(2)-xs(1));           % 半ピッチだけずらす
ff=subs(fx,xs+x0); % サンプリング点x1,x2,...,x(2*N+1) における fx の値
A=[];              % 行列Aの作成（ M=2*N+1 ）
                   %     [ 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=-N: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;      % 係数 a_(-50),...,a50 の算出
%figure(5)
%bar([-N:N]',aa(1:2*N+1))
Subst=xx*0;        % 代用ローラン展開のカーブ（初期化）
for n=-N:N         % 初期化カーブに各次の曲線を順次重ねていく。
  Subst = Subst + aa(n+1+N)*((xx-x0).^n);   % 注：　aa(1)=a_(-N), aa(M)=aN
end

figure(2)
hold on
h1o=copyobj(h1,gca);         % figure(1)の真値カーブを貼り付ける。
h2o=copyobj(h2,gca);         % 　〃　　 の展開の中心点を貼り付ける。
h3o=copyobj(h3,gca);         % 　〃　　 のfN(x)のカーブを貼り付ける。
h5o=copyobj(h5,gca);         % 　〃　　 の特異点の位置を貼り付ける。
axis(AX);
h4o=plot(xx,Subst,'-g','LineWidth',2);   % 代用ローラン展開のカーブ
legend([h1o h2o h3o h4o h5o],{'真値','展開の中心','正規テイラー展開', ...
                              '代用ローラン展開','特異点'},'Location','northeast')
grid on
xlabel('x')
ylabel('f(x)')
title({strcat('f(x)=',string(fx),'の真値と',num2str(N),'次までの正規テイラー展開'); ...
       strcat('および　狭域(Δx=',num2str(double(wn)), ...
       ')サンプリング高精度連立方程式による代用ローラン展開')})
text(tx,ty,'図２','FontSize',20)

% ■■■■■■　ステップ３　■■■■■■
% 関数 fx=1/(x^2) の真値
% および　x0=0（特異点）を中心とした ±50 次までの代用ローラン展開
% ■■■■■■■■■■■■■■■■■■■

% 代用ローラン展開のデータ
x0=vpa(0);
xs=linspace(-0.5,0.5,2*N+1)'*wn;
xs=xs+0.5*(xs(2)-xs(1));
ff=subs(fx,xs+x0);
A=[];
for n=-N:N
  A=[A xs.^n];
end
aa=inv(A)*ff;
%figure(6)
%bar([-N:N]',aa(1:2*N+1))
Subst=xx*0;
for n=-N:N
  Subst = Subst + aa(n+1+N)*((xx-x0).^n);
end

figure(3)
hold on
h1o=copyobj(h1,gca);
h2o=plot(x0,0,'ob','MarkerSize',10,'LineWidth',2);
h5o=copyobj(h5,gca);
axis(AX);
h4o=plot(xx,Subst,'-g','LineWidth',2);
legend([h1o h2o h4o h5o],{'真値','展開の中心','代用ローラン展開','特異点'}, ...
                          'Location','northeast')
grid on
xlabel('x')
ylabel('f(x)')
title({strcat('f(x)=',string(fx),'の真値'); ...
       strcat('および　狭域(Δx=',num2str(double(wn)), ...
       ')サンプリング高精度連立方程式による代用ローラン展開')})
text(tx,ty,'図３','FontSize',20)

% ■■■■■■　ステップ４　■■■■■■
% 関数 fx=1/(x^2)+0.2/((x+1)^2)+10*sin(10*x)+10 の真値
% および　x0=0（特異点）を中心とした ±50 次までの代用ローラン展開
% ■■■■■■■■■■■■■■■■■■■

% 真値のグラフ用のデータ
fx=1/(x^2)+0.2/((x+1)^2)+10*sin(10*x)+10;
Real=subs(fx,xx);

% 代用ローラン展開のデータ
x0=vpa(0);
xs=linspace(-0.5,0.5,2*N+1)'*wn;
xs=xs+0.5*(xs(2)-xs(1));
ff=subs(fx,xs+x0);
A=[];
for n=-N:N
  A=[A xs.^n];
end
aa=inv(A)*ff;
%figure(7)
%bar([-N:N]',aa(1:2*N+1))
Subst=xx*0;
for n=-N:N
  Subst = Subst + aa(n+1+N)*((xx-x0).^n);
end

figure(4)
hold on
h1o=plot(xx,Real,'-b','LineWidth',6);
h2o=plot(x0,0,'ob','MarkerSize',10,'LineWidth',2);
h5o=plot([-1 -1 0 0 0],[0 50 NaN 0 50],'-r','LineWidth',1);
axis(AX);
h4o=plot(xx,Subst,'-g','LineWidth',2);
legend([h1o h2o h4o h5o],{'真値','展開の中心','代用ローラン展開','特異点'}, ...
                          'Location','northeast')
grid on
xlabel('x')
ylabel('f(x)')
title({strcat('f(x)=',string(fx),'の真値　　および'); ...
       strcat('狭域(Δx=',num2str(double(wn)), ...
       ')サンプリング高精度連立方程式による±',num2str(N),'次代用ローラン展開')})
text(tx,ty,'図４','FontSize',20)