07 漸近展開の有名な例
黒木玄2018-05-20～2019-04-03, 2023-06-22
Copyright 2018,2023 Gen Kuroki
License: MIT https://opensource.org/licenses/MIT (https://opensource.org/licenses/MIT)
Repository: https://github.com/genkuroki/Calculus (https://github.com/genkuroki/Calculus)
このファイルは次の場所できれいに閲覧できる:
http://nbviewer.jupyter.org/github/genkuroki/Calculus/blob/master/07%20example%20of%20asymptotic%20expansion.ipynb(http://nbviewer.jupyter.org/github/genkuroki/Calculus/blob/master/07%20example%20of%20asymptotic%20expansion.ipynb)
https://genkuroki.github.io/documents/Calculus/07%20example%20of%20asymptotic%20expansion.pdf(https://genkuroki.github.io/documents/Calculus/07%20example%20of%20asymptotic%20expansion.pdf)
このファイルは Julia Box (https://juliabox.com) で利用できる.
自分のパソコンにJulia言語 (https://julialang.org/)をインストールしたい場合にはWindowsへのJulia言語のインストール (http://nbviewer.jupyter.org/gist/genkuroki/81de23edcae631a995e19a2ecf946a4f)
Julia v1.1.0 の Windows 8.1 へのインストール(https://nbviewer.jupyter.org/github/genkuroki/msfd28/blob/master/install.ipynb)
を参照せよ. 前者は古く, 後者の方が新しい.
論理的に完璧な説明をするつもりはない. 細部のいい加減な部分は自分で訂正・修正せよ.

目次
 1 漸近展開の有名な例
1.1の定義1.2
1.3の漸近展開
1.4の漸近展開の別の導出の仕方
1.5 以上のまとめ
1.6の漸近展開を用いた数値計算
1.7 漸近展開の入門書
1.8 指数積分函数
1.9 連分数展開による数値計算

𝐹𝑛 (𝑥) 𝑛+1𝐹𝑛 (𝑥) = 𝑂(𝑥 )𝐹0 (𝑥)𝐹0 (𝑥)
𝐹0 (1/10)

In [1]:



12
3
4
5
6
7
8
910
1112131415
161718
1920
21
2223
24
25
26
27
28
2930
31
32
33

using Base.MathConstants
using Base64
using Printf
using Statistics
const e = ℯendof(a) = lastindex(a)linspace(start, stop, length) = range(start, stop, length=length)​using Plots#gr(); ENV["PLOTS_TEST"] = "true"#clibrary(:colorcet)#clibrary(:misc)
default(fmt=:png)​function pngplot(P...; kwargs...)sleep(0.1)
pngfile = tempname() * ".png"
savefig(plot(P...; kwargs...), pngfile)showimg("image/png", pngfile)end
pngplot(; kwargs...) = pngplot(plot!(; kwargs...))​
showimg(mime, fn) = open(fn) do fbase64 = base64encode(f)display("text/html", """<img src="data:$mime;base64,$base64">""")
end​using SymPy#sympy.init_printing(order="lex") # default#sympy.init_printing(order="rev-lex")​
using SpecialFunctions
using QuadGKIn [2]:

123456789101112131415161718# Override the Base.show definition of SymPy.jl:# https://github.com/JuliaPy/SymPy.jl/blob/29c5bfd1d10ac53014fa7fef468bc8deccadc2fc/src/​@eval SymPy function Base.show(io::IO, ::MIME"text/latex", x::SymbolicObject)
print(io, as_markdown("\\displaystyle " * sympy.latex(x, mode="plain", fold_short_fra
end@eval SymPy function Base.show(io::IO, ::MIME"text/latex", x::AbstractArray{Sym})
function toeqnarray(x::Vector{Sym})a = join(["\\displaystyle " * sympy.latex(x[i]) for i in 1:length(x)], "\\\\")"""\\left[ \\begin{array}{r}$a\\end{array} \\right]"""end
function toeqnarray(x::AbstractArray{Sym,2})sz = size(x)a = join([join("\\displaystyle " .* map(sympy.latex, x[i,:]), "&") for i in 1:sz"\\left[ \\begin{array}{" * repeat("r",sz[2]) * "}" * a * "\\end{array}\\right]"
end
print(io, as_markdown(toeqnarray(x)))
end

1 漸近展開の有名な例
1.1

𝐹𝑛 (𝑥) の定義

1.2

𝐹𝑛 (𝑥) = 𝑂(𝑥𝑛+1 )

𝑛 は非負の整数であるとする. 𝑥 > 0 の函数 𝐹𝑛 (𝑥) を∞𝐹𝑛 (𝑥) = 𝑛! 𝑒1/𝑥 ∫1/𝑥 𝑒−𝑡 𝑡−𝑛−1 𝑑𝑡と定める. 𝐹𝑛 (𝑥) の 𝑥 > 0 が 0 に近いときの様子を調べたい.
𝑥 > 0 のとき 𝑡 ≧ 1/𝑥 ならば 𝑡−𝑛−1 ≦ 𝑥𝑛+1 なので∞
0 < 𝐹𝑛 (𝑥) ≦ 𝑛! 𝑒1/𝑥 ∫1/𝑥 𝑒−𝑡 𝑥𝑛+1 𝑑𝑡 = 𝑛! 𝑥𝑛+1 .
ゆえに,𝐹𝑛 (𝑥) = 𝑂(𝑥𝑛+1 ) (𝑥 ↘ 0).
1.3

𝐹0 (𝑥) の漸近展開

1.4

𝐹0 (𝑥) の漸近展開の別の導出の仕方(−𝑒−𝑡 )′ = 𝑒−𝑡 を用いた部分積分によって∞𝐹𝑛 (𝑥) = 𝑛! 𝑒1/𝑥 ([−𝑒−𝑡 𝑡−𝑛−1 ]∞1/𝑥 + ∫1/𝑥 𝑒−𝑡 (−(𝑛 + 1)𝑡−𝑛−2 )𝑑𝑡)
= 𝑛! 𝑥𝑛+1 − 𝐹𝑛+1 (𝑥).
ゆえに 𝑥 > 0 において,
𝐹0 (𝑥) = 0!𝑥 − 𝐹1 (𝑥) = 0!𝑥 − 1! 𝑥2 + 𝐹2 (𝑥) = ⋯2 + ⋯ + (−1)𝑛 𝑛!𝑥𝑛+1 + (−1)𝑛+1 𝐹𝑛+1 (𝑥)= 0!𝑥−1!
𝑥𝑛= ∑(−1)𝑘𝑘! 𝑥𝑘+1 + (−1)𝑛+1 𝐹𝑛+1 (𝑥).
𝑘=0𝑥 > 0 であるとする. 𝐹0 (𝑥) は 𝑡 = 1/𝑥 + 𝑢 という置換によって次のように書き直される:∞ −𝑡∞ −𝑢𝐹0 (𝑥) = 𝑒1/𝑥 ∫1/𝑥 𝑒𝑡 𝑑𝑡 = 𝑥 ∫0 1 𝑒+ 𝑥𝑢 𝑑𝑢.
ゆえに

1 = 𝑛−1 (−1)𝑘𝑧𝑘 + (−1)𝑛 𝑧𝑛
1+𝑧 ∑1+𝑧
𝑘=0

を

𝑧 = 𝑥𝑢 に適用した結果を使うと,
𝑛−1∞∞ −𝑢 𝑛𝐹0 (𝑥) = ∑(−1)𝑘𝑥𝑘+1 ∫0 𝑒−𝑢 𝑢𝑘 𝑑𝑢 + (−1)𝑛 𝑥𝑛+1 ∫0 1𝑒+ 𝑢𝑥𝑢 𝑑𝑢
𝑘=0
𝑛−1∞ −𝑢 𝑛= ∑(−1)𝑘𝑘!𝑥𝑘+1 + (−1)𝑛 𝑥𝑛+1 ∫0 1𝑒+ 𝑢𝑥𝑢 𝑑𝑢.
𝑘=0

さらに∞ −𝑢 𝑛∞
0 < 𝑥𝑛+1 ∫0 1𝑒+ 𝑢𝑥𝑢 𝑑𝑢 ≦ 𝑥𝑛+1 ∫0 𝑒−𝑢 𝑢𝑛 𝑑𝑢 = 𝑛! 𝑥𝑛+1なので∞ −𝑢 𝑛𝑥𝑛+1 ∫0 1𝑒+ 𝑢𝑥𝑢 𝑑𝑢 = 𝑂(𝑥𝑛+1 )(𝑥 ↘ 0).

1.5 以上のまとめ
以上をまとめると

𝑛−1

𝑛−1

𝑘=0

𝑘=0

𝐹0 (𝑥) = ∑(−1)𝑘𝑘! 𝑥𝑘+1 + (−1)𝑛 𝐹𝑛 (𝑥) = ∑(−1)𝑘𝑘! 𝑥𝑘+1 + 𝑂(𝑥𝑛+1 ) (𝑥 ↘ 0).
ここで

しかし,

𝑥 > 0で∞∞ −𝑢 𝑛𝐹𝑛 (𝑥) = 𝑛!𝑒1/𝑥 ∫1/𝑥 𝑒−𝑡 𝑡−𝑛−1 𝑑𝑡 = 𝑥𝑛+1 ∫0 1𝑒+ 𝑢𝑥𝑢 𝑑𝑢 = 𝑂(𝑥𝑛+1 ) (𝑥 ↘ 0).
∞(−1)𝑘𝑘! 𝑥𝑘+1 = 0!𝑥 − 1! 𝑥2 + 2! 𝑥3 − 3! 𝑥4 + ⋯∑𝑘=0は決して収束しない.

このように全てを足し上げると発散する場合であっても, 有限項と剰余項の和の形式で表わせば, すべてがwell-definedな量だけ

1.6

𝐹0 (1/10) の漸近展開を用いた数値計算

𝑥 > 0 のとき 0 < 𝐹𝑛 (𝑥) ≦ 𝑛! 𝑥𝑛+1 となるのであった. 𝑥 = 1/10 のとき 𝑛!𝑥𝑛+1 が最小になる
𝑛 は 𝑛 = 9,10 であることを確
𝑛−1−5認できる. そのとき 𝑛! 𝑥𝑛+1 = 3.6288 × 10 である. したがって, 𝑛 = 9,10 のときの∑(−1)𝑘𝑘! 𝑥𝑘+1 による 𝐹0 (𝑥) の近似は最良になると予想される.
In [3]:

𝑘=0

1 [(n, factorial(n)*(1/10)^(n+1)) for n in 7:13]

Out[3]: 7-element Vector{Tuple{Int64, Float64}}:(7, 5.0400000000000026e-5)(8, 4.032000000000002e-5)(9, 3.628800000000002e-5)(10, 3.628800000000002e-5)(11, 3.991680000000003e-5)(12, 4.790016000000003e-5)(13, 6.227020800000005e-5)
In [4]:



12
345
6
7
8
9101112F₀(x) = exp(1/x)*quadgk(u->exp(-1/(x*u))/u, 0, 1)[1]F₀_ae(x, n=10) = sum(k->(-1)^k*factorial(k)*x^(k+1), 0:n-1)​@show Y = F₀(1/10)@show Y_ae9 = F₀_ae(1/10, 9)@show Y_ae10 = F₀_ae(1/10, 10)​n = 1:20
plot(size=(500,350))
plot!(ylims=(-0.0001,0.0001))
plot!(legend=:top)
plot!(n, F₀_ae.(1/10, n) .- Y, label="error")

Y = F₀(1 / 10) = 0.09156333393978829
Y_ae9 = F₀_ae(1 / 10, 9) = 0.09158192
Y_ae10 = F₀_ae(1 / 10, 10) = 0.091545632
Out[4]:

1.7 漸近展開の入門書

𝑥>0

このようにで決して収束しないべき級数であっても, 適切な解釈のもとで数学的に十分な意味を持つことがある. このような理由で数学では発散級数に関するたくさんの深い研究がある.
その方面の入門書としては次の文献がある. 非常に面白い本なのでおすすめできる.
大久保謙二郎・河野實彦共著『漸近展開』, 新しい応用の数学12, 教育出版, 1976, 1996(https://www.amazon.co.jp/dp/4316376306).

1.8 指数積分函数

𝐹0 (𝑥) は実用的にも重要な函数である. それと本質的に同じ

∞ −𝑡
1 −𝑧/𝑢𝐸1 (𝑧) = 𝑒−𝑧 𝐹0 (1/𝑧) = ∫𝑧 𝑒𝑡 𝑑𝑡 = ∫0 𝑒 𝑢 𝑑𝑢は指数積分函数と呼ばれる特殊函数の1つである(𝑡 = 𝑧/𝑢 ). この函数の数値計算については
http://nbviewer.jupyter.org/github/stevengj/18S096/blob/iap2017/pset3/pset3-solutions.ipynb(http://nbviewer.jupyter.org/github/stevengj/18S096/blob/iap2017/pset3/pset3-solutions.ipynb)が非常に面白い解説になっている 特殊函数の数値計算の最適化に興味がある人は是非とも閲覧して欲しい

1.9 連分数展開による数値計算
以下は

𝑥 = 1/𝑧 のときの𝐺(𝑧) = 𝐹0 (1/𝑧) = 𝑒𝑧 𝐸1 (𝑧) = 𝑒𝑧 ∫0 𝑒−𝑧/𝑢𝑢 𝑑𝑢1の連分数展開による数値計算である. この函数の連分数展開については
一松信著『特殊関数入門』数学選書, 森北出版, 1999 (https://www.amazon.co.jp/dp/4627038216)の第3章の最後のp.74にある例3.10および第6章のpp.132-133にある例6.4を参照せよ. この本も非常に面白い本なのでおすすめできる.
In [5]:

123
45678910function G_cf(z; n::Int=2)
cf = 1 + (n+1)/zfor i = n:-1:1cf = z + (1+i)/cf
cf = 1 + i/cfend
return 1 / (z + 1/cf)
end​G(z) = exp(z)*quadgk(u->exp(-z/u)/u, 0, 1)[1]Out[5]: G (generic function with 1 method)
In [6]:

1234
5
6sympy.init_printing(order="lex")z = symbols("z")
for n in 1:3
cf = G_cf(z, n=n)display(cf)end

1𝑧 + 1+ 1 1 2𝑧+ 21+ 𝑧

𝑧 + 1+𝑧 + 1+

1𝑧+

112
1+ 2 3𝑧+ 31+ 𝑧

1𝑧+1+

11𝑧+2

2

31+ 3 4𝑧+ 41+ 𝑧

In [7]:1234
5sympy.init_printing(order="rev-lex")
display(series(sum((-1)^n*factorial(n)*z^(n+1) for n in 0:9), n=11))
for n in 1:3
display(series(G_cf(1/z, n=n), n=11))
end𝑧 − 𝑧2 + 2𝑧3 − 6𝑧4 + 24𝑧5 − 120𝑧6 + 720𝑧7 − 5040𝑧8 + 40320𝑧9 − 362880𝑧10𝑧 − 𝑧2 + 2𝑧3 − 6𝑧4 + 24𝑧5 − 108𝑧6 + 504𝑧7 − 2376𝑧8 + 11232𝑧9 − 53136𝑧10 + 𝑂 (𝑧11 )𝑧 − 𝑧2 + 2𝑧3 − 6𝑧4 + 24𝑧5 − 120𝑧6 + 720𝑧7 − 4896𝑧8 + 35712𝑧9 − 269568𝑧10 + 𝑂 (𝑧11 )𝑧 − 𝑧2 + 2𝑧3 − 6𝑧4 + 24𝑧5 − 120𝑧6 + 720𝑧7 − 5040𝑧8 + 40320𝑧9 − 360000𝑧10 + 𝑂 (𝑧11 )𝑛 = 1,2,3 の連分数近似であっても, 大きめの 𝑧 における近似が相当にうまく行っていることを以下で確認する.
In [8]:



1 G_cf(4), G(4)Out[8]: (0.20636792452830188, 0.20634564990105575)
In [9]:



1
23
4
5
67
8
9
101112x = 0.05:0.05:4@time y = G.(x)@time y_cf1 = G_cf.(x, n=1)
@time y_cf2 = G_cf.(x, n=2)
@time y_cf3 = G_cf.(x, n=3)
plot(size=(500,350))
plot!(title="y = G(x) = e^x E_1(x)", titlefontsize=11)
plot!(xlims=(0,maximum(x)), ylims=(0,2.3))
plot!(x, y, label="numetical integration")plot!(x, y_cf1, label="continued fractional approximation n=1", ls=:dash)plot!(x, y_cf2, label="continued fractional approximation n=2", ls=:dash)plot!(x, y_cf3, label="continued fractional approximation n=3", ls=:dash)0.124564 seconds (148.57 k allocations: 10.049 MiB, 99.61% compilation time)
0.049663 seconds (98.27 k allocations: 6.637 MiB, 99.27% compilation time)
0.000048 seconds (6 allocations: 1008 bytes)
0.000023 seconds (6 allocations: 1008 bytes)
Out[9]: