実数だけでローラン展開

はじめに

　テイラー展開は高校の数学でも習いましたが、ローラン展開はもっと上級の複素関数論の教科書の中で初めて出てきます。テイラー展開の方は、実数の世界で体験した公式が、複素数の世界でも成り立つという自然な流れで解釈できるので、それほど悩ましいことはありません。しかし、ローラン展開は複素数の世界で突然現れ、難解で理解し難く、なかなか馴染みにくい存在です。

　そこで、ここではローラン展開を一足飛びに複素数の領域で考えるのではなく、まずは実数の領域だけで疑似体験することで、少しでも拒否反応を和らげてみたいと思います。

ローラン展開

　ローラン展開の公式を示せば次のとおりです。

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}

　(1)式だけを見れば、下記のテイラー展開の $n{=}0$ を $n{=}{-}\infty$ に変えたに過ぎません。

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

代用ローラン展開

　係数 $a_n$ の式の中には微分どころか複素積分が入っており、テイラー展開に比べると、その係数を求めるのは遥かに難しそうです。実数だけでローラン展開を考えるにはどうすれば良いでしょうか？

　そこで、「微分を使わずテイラー展開」で示したのと同様に、ここでも、狭域多点サンプリングから得られる連立方程式を高精度で解くことで求めます。この方法が正しいという確証はありませんが、理論は度外視した実験数学です（しかし、結果を見ると、まんざら間違いでもなさそうです）。

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

　ここでは、展開によって近似可能な $x$ の範囲の限界も知りたいので、$N$ の値は「微分を使わずテイラー展開」のときよりもさらに $\infty$ に近い値の $50$ とします。すると、計算精度もさらに必要になるので、より高く$10$進$200$桁とします。したがって、「微分を使わずテイラー展開」と同様、やはりMATLAB 本体だけではなく、Symbolic Math Toolbox が必要です。

ローラン展開は役に立たない？

　(4)式を見ると、次数が負になる項もあるので、例えば $f(x)=1/x^2$ のような関数であれば、展開の結果は、当然、係数 $a_{\text{-}2}$ だけが $1.0$ で、残りは全て $0$ となるはずです。

　ローラン展開の優位性を示すための基準として、この関数を $x_0{=}0.5$ としてテイラー展開した結果を図１に示しておきます。

　同じ条件で代用ローラン展開した結果を図２に示します。当然の期待として、その結果は真値のカーブと一致するはずでした。しかし、何と、テイラー展開の結果と寸分変わりません。折角準備していた指定席$a_{\text{-}2}$は使われないままです。ローラン展開は何の役に立つのでしょうか？

ローラン展開はやはり役立つ！

　今度は、展開の中心だけを $x_0{=}0$ に変え、他は先と同一条件で再トライした結果を図３に示します。これは見事に真値と一致しています。指定席は無駄ではありませんでした。どうも、$x_0$ を、$f(x){=}\pm\infty$ となるような特異点に置くと真価を発揮するようです。

　期待に沿うような簡単な関数だけで誤魔化そうとしているのだろうと疑われては困ります。最後に、もっと複雑で、さらに特異点を２つ持つ関数 $f(x)=1/x^2+0.2/{(x+1)^2}+10\sin(10x)+10$ を、$x_0{=}0$ を中心に展開した結果を図４に示しておきます。

まとめ

　以上のカーブを観察して、実数の世界だけでローラン展開の次の性質を体験することができました。

①　展開の中心が特異点以外の場合には、テイラー展開の結果と変わらない。
②　特異点を中心にローラン展開すると、特異点の両側の、テイラー展開よりも広い領域で良好な近似が得られる。
③　$x_0$以外にも特異点がある場合、近似が成立するのは $x_0\pm|,x_0{-}(x_0$に至近の他の特異点$),|$ の範囲である。

プログラム

Laurent_exp_subst.m

% 【動作させるには 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)