Old_Ver_for_Julia_v0.6/07 example of asymptotic expansion.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 07 漸近展開の有名な例\n",
    "\n",
    "黒木玄\n",
    "\n",
    "2018-05-20\n",
    "\n",
    "* Copyright 2018 Gen Kuroki\n",
    "* License: MIT https://opensource.org/licenses/MIT\n",
    "* Repository: https://github.com/genkuroki/Calculus\n",
    "\n",
    "このファイルは次の場所できれいに閲覧できる:\n",
    "\n",
    "* http://nbviewer.jupyter.org/github/genkuroki/Calculus/blob/master/07%20example%20of%20asymptotic%20expansion.ipynb\n",
    "\n",
    "* https://genkuroki.github.io/documents/Calculus/07%20example%20of%20asymptotic%20expansion.pdf\n",
    "\n",
    "このファイルは <a href=\"https://juliabox.com\">Julia Box</a> で利用できる.\n",
    "\n",
    "自分のパソコンに<a href=\"https://julialang.org/\">Julia言語</a>をインストールしたい場合には\n",
    "\n",
    "* <a href=\"http://nbviewer.jupyter.org/gist/genkuroki/81de23edcae631a995e19a2ecf946a4f\">WindowsへのJulia言語のインストール</a>\n",
    "\n",
    "を参照せよ.\n",
    "\n",
    "論理的に完璧な説明をするつもりはない. 細部のいい加減な部分は自分で訂正・修正せよ.\n",
    "\n",
    "$\n",
    "\\newcommand\\eps{\\varepsilon}\n",
    "\\newcommand\\ds{\\displaystyle}\n",
    "\\newcommand\\Z{{\\mathbb Z}}\n",
    "\\newcommand\\R{{\\mathbb R}}\n",
    "\\newcommand\\C{{\\mathbb C}}\n",
    "\\newcommand\\QED{\\text{□}}\n",
    "\\newcommand\\root{\\sqrt}\n",
    "$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "toc": true
   },
   "source": [
    "<h1>目次<span class=\"tocSkip\"></span></h1>\n",
    "<div class=\"toc\"><ul class=\"toc-item\"><li><span><a href=\"#漸近展開の有名な例\" data-toc-modified-id=\"漸近展開の有名な例-1\"><span class=\"toc-item-num\">1&nbsp;&nbsp;</span>漸近展開の有名な例</a></span><ul class=\"toc-item\"><li><span><a href=\"#$F_n(x)$-の定義\" data-toc-modified-id=\"$F_n(x)$-の定義-1.1\"><span class=\"toc-item-num\">1.1&nbsp;&nbsp;</span>$F_n(x)$ の定義</a></span></li><li><span><a href=\"#$F_n(x)-=-O(x^{n+1})$\" data-toc-modified-id=\"$F_n(x)-=-O(x^{n+1})$-1.2\"><span class=\"toc-item-num\">1.2&nbsp;&nbsp;</span>$F_n(x) = O(x^{n+1})$</a></span></li><li><span><a href=\"#$F_0(x)$-の漸近展開\" data-toc-modified-id=\"$F_0(x)$-の漸近展開-1.3\"><span class=\"toc-item-num\">1.3&nbsp;&nbsp;</span>$F_0(x)$ の漸近展開</a></span></li><li><span><a href=\"#$F_0(x)$-の漸近展開の別の導出の仕方\" data-toc-modified-id=\"$F_0(x)$-の漸近展開の別の導出の仕方-1.4\"><span class=\"toc-item-num\">1.4&nbsp;&nbsp;</span>$F_0(x)$ の漸近展開の別の導出の仕方</a></span></li><li><span><a href=\"#以上のまとめ\" data-toc-modified-id=\"以上のまとめ-1.5\"><span class=\"toc-item-num\">1.5&nbsp;&nbsp;</span>以上のまとめ</a></span></li><li><span><a href=\"#$F_0(1/10)$-の漸近展開を用いた数値計算\" data-toc-modified-id=\"$F_0(1/10)$-の漸近展開を用いた数値計算-1.6\"><span class=\"toc-item-num\">1.6&nbsp;&nbsp;</span>$F_0(1/10)$ の漸近展開を用いた数値計算</a></span></li><li><span><a href=\"#漸近展開の入門書\" data-toc-modified-id=\"漸近展開の入門書-1.7\"><span class=\"toc-item-num\">1.7&nbsp;&nbsp;</span>漸近展開の入門書</a></span></li><li><span><a href=\"#指数積分函数\" data-toc-modified-id=\"指数積分函数-1.8\"><span class=\"toc-item-num\">1.8&nbsp;&nbsp;</span>指数積分函数</a></span></li><li><span><a href=\"#連分数展開による数値計算\" data-toc-modified-id=\"連分数展開による数値計算-1.9\"><span class=\"toc-item-num\">1.9&nbsp;&nbsp;</span>連分数展開による数値計算</a></span></li></ul></li></ul></div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "using Plots\n",
    "gr(); ENV[\"PLOTS_TEST\"] = \"true\"\n",
    "#clibrary(:colorcet)\n",
    "clibrary(:misc)\n",
    "\n",
    "function pngplot(P...; kwargs...)\n",
    "    sleep(0.1)\n",
    "    pngfile = tempname() * \".png\"\n",
    "    savefig(plot(P...; kwargs...), pngfile)\n",
    "    showimg(\"image/png\", pngfile)\n",
    "end\n",
    "pngplot(; kwargs...) = pngplot(plot!(; kwargs...))\n",
    "\n",
    "showimg(mime, fn) = open(fn) do f\n",
    "    base64 = base64encode(f)\n",
    "    display(\"text/html\", \"\"\"<img src=\"data:$mime;base64,$base64\">\"\"\")\n",
    "end\n",
    "\n",
    "using SymPy\n",
    "#sympy[:init_printing](order=\"lex\") # default\n",
    "#sympy[:init_printing](order=\"rev-lex\")\n",
    "\n",
    "using SpecialFunctions\n",
    "using QuadGK"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 漸近展開の有名な例"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### $F_n(x)$ の定義\n",
    "\n",
    "$n$ は非負の整数であるとする. $x>0$ の函数 $F_n(x)$ を\n",
    "\n",
    "$$\n",
    "F_n(x) = n!\\,e^{1/x}\\int_{1/x}^\\infty e^{-t} t^{-n-1}\\,dt\n",
    "$$\n",
    "\n",
    "と定める.  $F_n(x)$ の $x>0$ が $0$ に近いときの様子を調べたい."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### $F_n(x) = O(x^{n+1})$\n",
    "\n",
    "$x>0$ のとき $t\\geqq 1/x$ ならば $t^{-n-1}\\leqq x^{n+1}$ なので\n",
    "\n",
    "$$\n",
    "0 < F_n(x) \\leqq n!\\,e^{1/x}\\int_{1/x}^\\infty e^{-t}x^{n+1}\\,dt = n!\\, x^{n+1}.\n",
    "$$\n",
    "\n",
    "ゆえに,\n",
    "\n",
    "$$\n",
    "F_n(x) = O(x^{n+1}) \\quad (x\\searrow 0).\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### $F_0(x)$ の漸近展開\n",
    "\n",
    "$(-e^{-t})'=e^{-t}$ を用いた部分積分によって\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "F_n(x) &= n!\\,e^{1/x}\\left(\n",
    "\\left[-e^{-t}t^{-n-1}\\right]_{1/x}^\\infty + \\int_{1/x}^\\infty e^{-t}(-(n+1)t^{-n-2})dt\n",
    "\\right)\n",
    "\\\\ &=\n",
    "n!\\,x^{n+1} - F_{n+1}(x).\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "ゆえに $x>0$ において, \n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "F_0(x) &= 0!\\,x - F_1(x) =\n",
    "0!\\,x - 1!\\,x^2 + F_2(x) = \\cdots \n",
    "\\\\ &=\n",
    "0!\\,x - 1!\\,x^2 + \\cdots + (-1)^n n! x^{n+1} + (-1)^{n+1} F_{n+1}(x) \n",
    "\\\\ &=\n",
    "\\sum_{k=0}^n (-1)^k k!\\, x^{k+1} + (-1)^{n+1} F_{n+1}(x).\n",
    "\\end{aligned}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### $F_0(x)$ の漸近展開の別の導出の仕方\n",
    "\n",
    "$x>0$ であるとする. $F_0(x)$ は $t=1/x+u$ という置換によって次のように書き直される:\n",
    "\n",
    "$$\n",
    "F_0(x) = e^{1/x}\\int_{1/x}^\\infty \\frac{e^{-t}}{t}\\,dt =\n",
    "x\\int_0^\\infty \\frac{e^{-u}}{1+xu}\\,du.\n",
    "$$\n",
    "\n",
    "ゆえに\n",
    "\n",
    "$$\n",
    "\\frac{1}{1+z} = \\sum_{k=0}^{n-1}(-1)^k z^k + (-1)^n\\frac{z^n}{1+z}\n",
    "$$\n",
    "\n",
    "を $z=xu$ に適用した結果を使うと, \n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "F_0(x) &= \\sum_{k=0}^{n-1} (-1)^k x^{k+1}\\int_0^\\infty e^{-u}u^k\\,du + \n",
    "(-1)^n x^{n+1}\\int_0^\\infty \\frac{e^{-u}u^n}{1+xu}\\,du\n",
    "\\\\ &=\n",
    "\\sum_{k=0}^{n-1} (-1)^k k! x^{k+1} + \n",
    "(-1)^n x^{n+1}\\int_0^\\infty \\frac{e^{-u}u^n}{1+xu}\\,du.\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "さらに\n",
    "\n",
    "$$\n",
    "0 < x^{n+1}\\int_0^\\infty \\frac{e^{-u}u^n}{1+xu}\\,du \\leqq \n",
    "x^{n+1}\\int_0^\\infty e^{-u}u^n\\,du = n!\\,x^{n+1}\n",
    "$$\n",
    "\n",
    "なので\n",
    "\n",
    "$$\n",
    "x^{n+1}\\int_0^\\infty \\frac{e^{-u}u^n}{1+xu}\\,du = O(x^{n+1}) \\qquad(x\\searrow 0).\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 以上のまとめ\n",
    "\n",
    "以上をまとめると\n",
    "\n",
    "$$\n",
    "F_0(x) = \\sum_{k=0}^{n-1} (-1)^k k!\\, x^{k+1} + (-1)^n F_n(x) =\n",
    "\\sum_{k=0}^{n-1}(-1)^k k!\\,x^{k+1} + O(x^{n+1}) \\quad (x\\searrow 0).\n",
    "$$\n",
    "\n",
    "ここで\n",
    "\n",
    "$$\n",
    "F_n(x) = n!e^{1/x}\\int_{1/x}^\\infty e^{-t}t^{-n-1}\\,dt =\n",
    "x^{n+1}\\int_0^\\infty \\frac{e^{-u}u^n}{1+xu}\\,du = O(x^{n+1})\n",
    "\\quad(x\\searrow 0).\n",
    "$$\n",
    "\n",
    "しかし, $x>0$ で\n",
    "\n",
    "$$\n",
    "\\sum_{k=0}^\\infty (-1)^k k!\\, x^{k+1} =\n",
    "0!\\,x - 1!\\,x^2 + 2!\\,x^3 - 3!\\,x^4 + \\cdots\n",
    "$$\n",
    "\n",
    "は決して収束しない. \n",
    "\n",
    "このように全てを足し上げると発散する場合であっても, 有限項と剰余項の和の形式で表わせば, すべてがwell-definedな量だけを使って議論を進めることができる場合がある. 上の場合には剰余項は $(-1)^n F_n(x)$ である."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### $F_0(1/10)$ の漸近展開を用いた数値計算\n",
    "\n",
    "$x>0$ のとき $0<F_n(x)\\leqq n!\\,x^{n+1}$ となるのであった. $x=1/10$ のとき $n!\\,x^{n+1}$ が最小になる $n$ は $n=9,10$ であることを確認できる. そのとき $n!\\,x^{n+1} = 3.6288\\times10^{-5}$ である.  したがって, $n=9,10$ のときの $\\ds\\sum_{k=0}^{n-1} (-1)^k k!\\, x^{k+1}$ による $F_0(x)$ の近似は最良になると予想される."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "7-element Array{Tuple{Int64,Float64},1}:\n",
       " (7, 5.04e-5)    \n",
       " (8, 4.032e-5)   \n",
       " (9, 3.6288e-5)  \n",
       " (10, 3.6288e-5) \n",
       " (11, 3.99168e-5)\n",
       " (12, 4.79002e-5)\n",
       " (13, 6.22702e-5)"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "[(n, factorial(n)*(1/10)^(n+1)) for n in 7:13]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Y = F₀(1 / 10) = 0.09156333393978826\n",
      "Y_ae9 = F₀_ae(1 / 10, 9) = 0.09158192\n",
      "Y_ae10 = F₀_ae(1 / 10, 10) = 0.091545632\n"
     ]
    },
    {
     "data": {
      "image/svg+xml": [
       "<?xml version=\"1.0\" encoding=\"utf-8\"?>\n",
       "<svg xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\" width=\"500\" height=\"350\" viewBox=\"0 0 500 350\">\n",
       "<defs>\n",
       "  <clipPath id=\"clip8300\">\n",
       "    <rect x=\"0\" y=\"0\" width=\"500\" height=\"350\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<polygon clip-path=\"url(#clip8300)\" points=\"\n",
       "0,350 500,350 500,0 0,0 \n",
       "  \" fill=\"#ffffff\" fill-opacity=\"1\"/>\n",
       "<defs>\n",
       "  <clipPath id=\"clip8301\">\n",
       "    <rect x=\"100\" y=\"0\" width=\"351\" height=\"350\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<polygon clip-path=\"url(#clip8300)\" points=\"\n",
       "66.2644,325.869 480.315,325.869 480.315,11.811 66.2644,11.811 \n",
       "  \" fill=\"#ffffff\" fill-opacity=\"1\"/>\n",
       "<defs>\n",
       "  <clipPath id=\"clip8302\">\n",
       "    <rect x=\"66\" y=\"11\" width=\"415\" height=\"315\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<polyline clip-path=\"url(#clip8302)\" style=\"stroke:#000000; stroke-width:0.5; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  153.433,325.869 153.433,11.811 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip8302)\" style=\"stroke:#000000; stroke-width:0.5; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  262.394,325.869 262.394,11.811 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip8302)\" style=\"stroke:#000000; stroke-width:0.5; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  371.354,325.869 371.354,11.811 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip8302)\" style=\"stroke:#000000; stroke-width:0.5; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  480.315,325.869 480.315,11.811 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip8302)\" style=\"stroke:#000000; stroke-width:0.5; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  66.2644,325.869 480.315,325.869 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip8302)\" style=\"stroke:#000000; stroke-width:0.5; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  66.2644,247.354 480.315,247.354 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip8302)\" style=\"stroke:#000000; stroke-width:0.5; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  66.2644,168.84 480.315,168.84 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip8302)\" style=\"stroke:#000000; stroke-width:0.5; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  66.2644,90.3254 480.315,90.3254 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip8302)\" style=\"stroke:#000000; stroke-width:0.5; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  66.2644,11.811 480.315,11.811 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip8300)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  66.2644,325.869 480.315,325.869 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip8300)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  66.2644,325.869 66.2644,11.811 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip8300)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  153.433,325.869 153.433,321.158 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip8300)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  262.394,325.869 262.394,321.158 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip8300)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  371.354,325.869 371.354,321.158 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip8300)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  480.315,325.869 480.315,321.158 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip8300)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  66.2644,325.869 72.4751,325.869 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip8300)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  66.2644,247.354 72.4751,247.354 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip8300)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  66.2644,168.84 72.4751,168.84 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip8300)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  66.2644,90.3254 72.4751,90.3254 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip8300)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  66.2644,11.811 72.4751,11.811 \n",
       "  \"/>\n",
       "<g clip-path=\"url(#clip8300)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:middle;\" transform=\"rotate(0, 153.433, 339.169)\" x=\"153.433\" y=\"339.169\">5</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip8300)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:middle;\" transform=\"rotate(0, 262.394, 339.169)\" x=\"262.394\" y=\"339.169\">10</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip8300)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:middle;\" transform=\"rotate(0, 371.354, 339.169)\" x=\"371.354\" y=\"339.169\">15</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip8300)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:middle;\" transform=\"rotate(0, 480.315, 339.169)\" x=\"480.315\" y=\"339.169\">20</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip8300)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:end;\" transform=\"rotate(0, 61.2644, 330.369)\" x=\"61.2644\" y=\"330.369\">-0.00010</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip8300)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:end;\" transform=\"rotate(0, 61.2644, 251.854)\" x=\"61.2644\" y=\"251.854\">-0.00005</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip8300)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:end;\" transform=\"rotate(0, 61.2644, 173.34)\" x=\"61.2644\" y=\"173.34\">0.00000</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip8300)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:end;\" transform=\"rotate(0, 61.2644, 94.8254)\" x=\"61.2644\" y=\"94.8254\">0.00005</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip8300)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:end;\" transform=\"rotate(0, 61.2644, 16.311)\" x=\"61.2644\" y=\"16.311\">0.00010</text>\n",
       "</g>\n",
       "<polyline clip-path=\"url(#clip8302)\" style=\"stroke:#009af9; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  66.2644,-13079.2 88.0565,2623.72 109.849,-516.852 131.641,425.321 153.433,48.452 175.225,236.887 197.017,123.826 218.809,202.968 240.601,139.654 262.394,196.637 \n",
       "  284.186,139.654 305.978,202.335 327.77,127.118 349.562,224.9 371.354,88.0053 393.146,293.348 414.939,-35.2002 436.731,523.332 458.523,-482.026 480.315,1428.15 \n",
       "  \n",
       "  \"/>\n",
       "<polygon clip-path=\"url(#clip8300)\" points=\"\n",
       "174.324,56.931 249.727,56.931 249.727,26.691 174.324,26.691 \n",
       "  \" fill=\"#ffffff\" fill-opacity=\"1\"/>\n",
       "<polyline clip-path=\"url(#clip8300)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  174.324,56.931 249.727,56.931 249.727,26.691 174.324,26.691 174.324,56.931 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip8300)\" style=\"stroke:#009af9; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  179.324,41.811 209.324,41.811 \n",
       "  \"/>\n",
       "<g clip-path=\"url(#clip8300)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:start;\" transform=\"rotate(0, 214.324, 46.311)\" x=\"214.324\" y=\"46.311\">error</text>\n",
       "</g>\n",
       "</svg>\n"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "F₀(x) = exp(1/x)*quadgk(u->exp(-1/(x*u))/u, 0, 1)[1]\n",
    "F₀_ae(x, n=10) = sum(k->(-1)^k*factorial(k)*x^(k+1), 0:n-1)\n",
    "\n",
    "@show Y = F₀(1/10)\n",
    "@show Y_ae9 = F₀_ae(1/10, 9)\n",
    "@show Y_ae10 = F₀_ae(1/10, 10)\n",
    "\n",
    "n = 1:20\n",
    "plot(size=(500,350))\n",
    "plot!(ylims=(-0.0001,0.0001))\n",
    "plot!(legend=:top)\n",
    "plot!(n, F₀_ae.(1/10, n)-Y, label=\"error\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 漸近展開の入門書\n",
    "\n",
    "このように $x > 0$ で決して収束しないべき級数であっても, 適切な解釈のもとで数学的に十分な意味を持つことがある.  このような理由で数学では発散級数に関するたくさんの深い研究がある.\n",
    "\n",
    "その方面の入門書としては次の文献がある. 非常に面白い本なのでおすすめできる.\n",
    "\n",
    "* <a href=\"https://www.amazon.co.jp/dp/4316376306\">大久保謙二郎・河野實彦共著『漸近展開』, 新しい応用の数学12, 教育出版, 1976, 1996</a>."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 指数積分函数\n",
    "\n",
    "$F_0(x)$ は実用的にも重要な函数である. それと本質的に同じ\n",
    "\n",
    "$$\n",
    "E_1(z) = e^{-z}F_0(1/z) = \\int_z^\\infty \\frac{e^{-t}}{t}\\,dt =\n",
    "\\int_0^1 \\frac{e^{-z/u}}{u} du\n",
    "$$\n",
    "\n",
    "は指数積分函数と呼ばれる特殊函数の1つである($t=z/u$). この函数の数値計算については\n",
    "\n",
    "* http://nbviewer.jupyter.org/github/stevengj/18S096/blob/iap2017/pset3/pset3-solutions.ipynb\n",
    "\n",
    "が非常に面白い解説になっている. 特殊函数の数値計算の最適化に興味がある人は是非とも閲覧して欲しい."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 連分数展開による数値計算\n",
    "\n",
    "以下は $x=1/z$ のときの\n",
    "\n",
    "$$\n",
    "G(z) = F_0(1/z) = e^z E_1(z) = e^z \\int_0^1 \\frac{e^{-z/u}}{u}du\n",
    "$$\n",
    "\n",
    "の連分数展開による数値計算である. この函数の連分数展開については\n",
    "\n",
    "* <a href=\"https://www.amazon.co.jp/dp/4627038216\">一松信著『特殊関数入門』数学選書, 森北出版, 1999</a>\n",
    "\n",
    "の第3章の最後のp.74にある例3.10および第6章のpp.132-133にある例6.4を参照せよ. この本も非常に面白い本なのでおすすめできる."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "G (generic function with 1 method)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function G_cf(z; n::Int=2)\n",
    "    cf = 1 + (n+1)/z\n",
    "    for i = n:-1:1\n",
    "        cf = z + (1+i)/cf\n",
    "        cf = 1 + i/cf\n",
    "    end\n",
    "    return 1 / (z + 1/cf)\n",
    "end\n",
    "\n",
    "G(z) = exp(z)*quadgk(u->exp(-z/u)/u, 0, 1)[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$\\frac{1}{z + \\frac{1}{1 + \\frac{1}{z + \\frac{2}{1 + \\frac{2}{z}}}}}$$"
      ],
      "text/plain": [
       "        1        \n",
       "-----------------\n",
       "          1      \n",
       "z + -------------\n",
       "            1    \n",
       "    1 + ---------\n",
       "              2  \n",
       "        z + -----\n",
       "                2\n",
       "            1 + -\n",
       "                z"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$$\\frac{1}{z + \\frac{1}{1 + \\frac{1}{z + \\frac{2}{1 + \\frac{2}{z + \\frac{3}{1 + \\frac{3}{z}}}}}}}$$"
      ],
      "text/plain": [
       "            1            \n",
       "-------------------------\n",
       "              1          \n",
       "z + ---------------------\n",
       "                1        \n",
       "    1 + -----------------\n",
       "                  2      \n",
       "        z + -------------\n",
       "                    2    \n",
       "            1 + ---------\n",
       "                      3  \n",
       "                z + -----\n",
       "                        3\n",
       "                    1 + -\n",
       "                        z"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$$\\frac{1}{z + \\frac{1}{1 + \\frac{1}{z + \\frac{2}{1 + \\frac{2}{z + \\frac{3}{1 + \\frac{3}{z + \\frac{4}{1 + \\frac{4}{z}}}}}}}}}$$"
      ],
      "text/plain": [
       "                1                \n",
       "---------------------------------\n",
       "                  1              \n",
       "z + -----------------------------\n",
       "                    1            \n",
       "    1 + -------------------------\n",
       "                      2          \n",
       "        z + ---------------------\n",
       "                        2        \n",
       "            1 + -----------------\n",
       "                          3      \n",
       "                z + -------------\n",
       "                            3    \n",
       "                    1 + ---------\n",
       "                              4  \n",
       "                        z + -----\n",
       "                                4\n",
       "                            1 + -\n",
       "                                z"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "sympy[:init_printing](order=\"lex\")\n",
    "z = symbols(\"z\")\n",
    "for n in 1:3\n",
    "    cf = G_cf(z, n=n)\n",
    "    display(cf)\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$z - z^{2} + 2 z^{3} - 6 z^{4} + 24 z^{5} - 120 z^{6} + 720 z^{7} - 5040 z^{8} + 40320 z^{9} - 362880 z^{10}$$"
      ],
      "text/plain": [
       "     2      3      4       5        6        7         8          9           \n",
       "z - z  + 2*z  - 6*z  + 24*z  - 120*z  + 720*z  - 5040*z  + 40320*z  - 362880*z\n",
       "\n",
       "10\n",
       "  "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$$z - z^{2} + 2 z^{3} - 6 z^{4} + 24 z^{5} - 108 z^{6} + 504 z^{7} - 2376 z^{8} + 11232 z^{9} - 53136 z^{10} + \\mathcal{O}\\left(z^{11}\\right)$$"
      ],
      "text/plain": [
       "     2      3      4       5        6        7         8          9          1\n",
       "z - z  + 2*z  - 6*z  + 24*z  - 108*z  + 504*z  - 2376*z  + 11232*z  - 53136*z \n",
       "\n",
       "0    / 11\\\n",
       "  + O\\z  /"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$$z - z^{2} + 2 z^{3} - 6 z^{4} + 24 z^{5} - 120 z^{6} + 720 z^{7} - 4896 z^{8} + 35712 z^{9} - 269568 z^{10} + \\mathcal{O}\\left(z^{11}\\right)$$"
      ],
      "text/plain": [
       "     2      3      4       5        6        7         8          9           \n",
       "z - z  + 2*z  - 6*z  + 24*z  - 120*z  + 720*z  - 4896*z  + 35712*z  - 269568*z\n",
       "\n",
       "10    / 11\\\n",
       "   + O\\z  /"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$$z - z^{2} + 2 z^{3} - 6 z^{4} + 24 z^{5} - 120 z^{6} + 720 z^{7} - 5040 z^{8} + 40320 z^{9} - 360000 z^{10} + \\mathcal{O}\\left(z^{11}\\right)$$"
      ],
      "text/plain": [
       "     2      3      4       5        6        7         8          9           \n",
       "z - z  + 2*z  - 6*z  + 24*z  - 120*z  + 720*z  - 5040*z  + 40320*z  - 360000*z\n",
       "\n",
       "10    / 11\\\n",
       "   + O\\z  /"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "sympy[:init_printing](order=\"rev-lex\")\n",
    "display(series(sum((-1)^n*factorial(n)*z^(n+1) for n in 0:9), n=11))\n",
    "for n in 1:3\n",
    "    display(series(G_cf(1/z, n=n), n=11))\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$n=1,2,3$ の連分数近似であっても, 大きめの $z$ における近似が相当にうまく行っていることを以下で確認する."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.20636792452830188, 0.20634564990105575)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "G_cf(4), G(4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  2.022534 seconds (266.78 k allocations: 13.742 MiB)\n",
      "  0.047182 seconds (16.07 k allocations: 832.076 KiB)\n",
      "  0.016745 seconds (5.59 k allocations: 284.821 KiB)\n",
      "  0.016846 seconds (5.59 k allocations: 285.634 KiB)\n"
     ]
    },
    {
     "data": {
      "image/svg+xml": [
       "<?xml version=\"1.0\" encoding=\"utf-8\"?>\n",
       "<svg xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\" width=\"500\" height=\"350\" viewBox=\"0 0 500 350\">\n",
       "<defs>\n",
       "  <clipPath id=\"clip8500\">\n",
       "    <rect x=\"0\" y=\"0\" width=\"500\" height=\"350\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<polygon clip-path=\"url(#clip8500)\" points=\"\n",
       "0,350 500,350 500,0 0,0 \n",
       "  \" fill=\"#ffffff\" fill-opacity=\"1\"/>\n",
       "<defs>\n",
       "  <clipPath id=\"clip8501\">\n",
       "    <rect x=\"100\" y=\"0\" width=\"351\" height=\"350\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<polygon clip-path=\"url(#clip8500)\" points=\"\n",
       "32.4745,325.45 480.315,325.45 480.315,31.4961 32.4745,31.4961 \n",
       "  \" fill=\"#ffffff\" fill-opacity=\"1\"/>\n",
       "<defs>\n",
       "  <clipPath id=\"clip8502\">\n",
       "    <rect x=\"32\" y=\"31\" width=\"449\" height=\"295\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<polyline clip-path=\"url(#clip8502)\" style=\"stroke:#000000; stroke-width:0.5; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  32.4745,325.45 32.4745,31.4961 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip8502)\" style=\"stroke:#000000; stroke-width:0.5; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  144.435,325.45 144.435,31.4961 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip8502)\" style=\"stroke:#000000; stroke-width:0.5; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  256.395,325.45 256.395,31.4961 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip8502)\" style=\"stroke:#000000; stroke-width:0.5; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  368.355,325.45 368.355,31.4961 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip8502)\" style=\"stroke:#000000; stroke-width:0.5; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  480.315,325.45 480.315,31.4961 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip8502)\" style=\"stroke:#000000; stroke-width:0.5; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  32.4745,325.45 480.315,325.45 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip8502)\" style=\"stroke:#000000; stroke-width:0.5; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  32.4745,261.547 480.315,261.547 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip8502)\" style=\"stroke:#000000; stroke-width:0.5; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  32.4745,197.644 480.315,197.644 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip8502)\" style=\"stroke:#000000; stroke-width:0.5; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  32.4745,133.741 480.315,133.741 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip8502)\" style=\"stroke:#000000; stroke-width:0.5; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  32.4745,69.8378 480.315,69.8378 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip8500)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  32.4745,325.45 480.315,325.45 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip8500)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  32.4745,325.45 32.4745,31.4961 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip8500)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  32.4745,325.45 32.4745,321.04 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip8500)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  144.435,325.45 144.435,321.04 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip8500)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  256.395,325.45 256.395,321.04 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip8500)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  368.355,325.45 368.355,321.04 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip8500)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  480.315,325.45 480.315,321.04 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip8500)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  32.4745,325.45 39.1921,325.45 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip8500)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  32.4745,261.547 39.1921,261.547 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip8500)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  32.4745,197.644 39.1921,197.644 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip8500)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  32.4745,133.741 39.1921,133.741 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip8500)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  32.4745,69.8378 39.1921,69.8378 \n",
       "  \"/>\n",
       "<g clip-path=\"url(#clip8500)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:middle;\" transform=\"rotate(0, 32.4745, 338.75)\" x=\"32.4745\" y=\"338.75\">0</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip8500)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:middle;\" transform=\"rotate(0, 144.435, 338.75)\" x=\"144.435\" y=\"338.75\">1</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip8500)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:middle;\" transform=\"rotate(0, 256.395, 338.75)\" x=\"256.395\" y=\"338.75\">2</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip8500)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:middle;\" transform=\"rotate(0, 368.355, 338.75)\" x=\"368.355\" y=\"338.75\">3</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip8500)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:middle;\" transform=\"rotate(0, 480.315, 338.75)\" x=\"480.315\" y=\"338.75\">4</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip8500)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:end;\" transform=\"rotate(0, 27.4745, 329.95)\" x=\"27.4745\" y=\"329.95\">0.0</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip8500)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:end;\" transform=\"rotate(0, 27.4745, 266.047)\" x=\"27.4745\" y=\"266.047\">0.5</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip8500)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:end;\" transform=\"rotate(0, 27.4745, 202.144)\" x=\"27.4745\" y=\"202.144\">1.0</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip8500)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:end;\" transform=\"rotate(0, 27.4745, 138.241)\" x=\"27.4745\" y=\"138.241\">1.5</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip8500)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:end;\" transform=\"rotate(0, 27.4745, 74.3378)\" x=\"27.4745\" y=\"74.3378\">2.0</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip8500)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:16; text-anchor:middle;\" transform=\"rotate(0, 256.395, 14.4)\" x=\"256.395\" y=\"14.4\">y = G(x) = e^x E_1(x)</text>\n",
       "</g>\n",
       "<polyline clip-path=\"url(#clip8502)\" style=\"stroke:#009af9; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  38.0725,-6.13386 43.6705,67.9664 49.2685,107.993 54.8665,134.591 60.4645,154.077 66.0625,169.202 71.6605,181.407 77.2585,191.531 82.8565,200.108 88.4546,207.496 \n",
       "  94.0526,213.944 99.6506,219.635 105.249,224.703 110.847,229.253 116.445,233.365 122.043,237.104 127.641,240.522 133.239,243.66 138.837,246.554 144.435,249.233 \n",
       "  150.033,251.721 155.631,254.038 161.229,256.204 166.827,258.232 172.425,260.137 178.023,261.929 183.621,263.619 189.219,265.216 194.817,266.727 200.415,268.16 \n",
       "  206.013,269.52 211.611,270.814 217.209,272.046 222.807,273.221 228.405,274.343 234.003,275.415 239.601,276.44 245.199,277.423 250.797,278.365 256.395,279.27 \n",
       "  261.993,280.138 267.591,280.974 273.189,281.777 278.787,282.551 284.385,283.297 289.983,284.017 295.581,284.711 301.179,285.382 306.777,286.03 312.375,286.657 \n",
       "  317.973,287.264 323.571,287.851 329.169,288.42 334.767,288.971 340.365,289.506 345.963,290.024 351.561,290.528 357.159,291.017 362.757,291.492 368.355,291.954 \n",
       "  373.953,292.403 379.551,292.84 385.149,293.265 390.747,293.679 396.345,294.082 401.943,294.474 407.541,294.857 413.139,295.23 418.737,295.594 424.335,295.949 \n",
       "  429.933,296.296 435.531,296.634 441.129,296.965 446.727,297.288 452.325,297.603 457.923,297.911 463.521,298.212 469.119,298.507 474.717,298.795 480.315,299.077 \n",
       "  \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip8502)\" style=\"stroke:#e26f46; stroke-width:1; stroke-opacity:1; fill:none\" stroke-dasharray=\"8, 5\" points=\"\n",
       "  38.0725,-588.101 43.6705,-159.865 49.2685,-15.7969 54.8665,57.1278 60.4645,101.525 66.0625,131.608 71.6605,153.473 77.2585,170.171 82.8565,183.403 88.4546,194.189 \n",
       "  94.0526,203.182 99.6506,210.818 105.249,217.399 110.847,223.143 116.445,228.211 122.043,232.723 127.641,236.773 133.239,240.432 138.837,243.759 144.435,246.8 \n",
       "  150.033,249.593 155.631,252.17 161.229,254.556 166.827,256.774 172.425,258.842 178.023,260.775 183.621,262.588 189.219,264.292 194.817,265.897 200.415,267.412 \n",
       "  206.013,268.845 211.611,270.203 217.209,271.492 222.807,272.717 228.405,273.884 234.003,274.996 239.601,276.058 245.199,277.073 250.797,278.044 256.395,278.975 \n",
       "  261.993,279.867 267.591,280.724 273.189,281.547 278.787,282.338 284.385,283.101 289.983,283.835 295.581,284.542 301.179,285.225 306.777,285.884 312.375,286.521 \n",
       "  317.973,287.137 323.571,287.733 329.169,288.31 334.767,288.868 340.365,289.41 345.963,289.935 351.561,290.444 357.159,290.938 362.757,291.418 368.355,291.884 \n",
       "  373.953,292.338 379.551,292.778 385.149,293.207 390.747,293.624 396.345,294.031 401.943,294.426 407.541,294.812 413.139,295.187 418.737,295.554 424.335,295.911 \n",
       "  429.933,296.26 435.531,296.6 441.129,296.932 446.727,297.257 452.325,297.574 457.923,297.883 463.521,298.186 469.119,298.482 474.717,298.772 480.315,299.055 \n",
       "  \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip8502)\" style=\"stroke:#3da44d; stroke-width:1; stroke-opacity:1; fill:none\" stroke-dasharray=\"8, 5\" points=\"\n",
       "  38.0725,-399.683 43.6705,-76.1654 49.2685,33.9781 54.8665,90.5526 60.4645,125.55 66.0625,149.653 71.6605,167.455 77.2585,181.263 82.8565,192.364 88.4546,201.538 \n",
       "  94.0526,209.283 99.6506,215.937 105.249,221.733 110.847,226.841 116.445,231.389 122.043,235.471 127.641,239.161 133.239,242.519 138.837,245.592 144.435,248.416 \n",
       "  150.033,251.024 155.631,253.441 161.229,255.689 166.827,257.787 172.425,259.75 178.023,261.592 183.621,263.324 189.219,264.957 194.817,266.499 200.415,267.959 \n",
       "  206.013,269.343 211.611,270.656 217.209,271.906 222.807,273.096 228.405,274.231 234.003,275.315 239.601,276.351 245.199,277.342 250.797,278.293 256.395,279.204 \n",
       "  261.993,280.079 267.591,280.92 273.189,281.729 278.787,282.507 284.385,283.257 289.983,283.98 295.581,284.678 301.179,285.352 306.777,286.002 312.375,286.632 \n",
       "  317.973,287.24 323.571,287.829 329.169,288.4 334.767,288.953 340.365,289.489 345.963,290.009 351.561,290.514 357.159,291.004 362.757,291.48 368.355,291.942 \n",
       "  373.953,292.392 379.551,292.83 385.149,293.256 390.747,293.67 396.345,294.074 401.943,294.467 407.541,294.85 413.139,295.224 418.737,295.588 424.335,295.944 \n",
       "  429.933,296.291 435.531,296.63 441.129,296.96 446.727,297.283 452.325,297.599 457.923,297.907 463.521,298.209 469.119,298.504 474.717,298.792 480.315,299.074 \n",
       "  \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip8502)\" style=\"stroke:#c271d2; stroke-width:1; stroke-opacity:1; fill:none\" stroke-dasharray=\"8, 5\" points=\"\n",
       "  38.0725,-291.975 43.6705,-30.4484 49.2685,60.0211 54.8665,107.342 60.4645,137.158 66.0625,158.055 71.6605,173.739 77.2585,186.08 82.8565,196.131 88.4546,204.531 \n",
       "  94.0526,211.693 99.6506,217.9 105.249,223.348 110.847,228.182 116.445,232.51 122.043,236.415 127.641,239.962 133.239,243.202 138.837,246.176 144.435,248.919 \n",
       "  150.033,251.459 155.631,253.819 161.229,256.018 166.827,258.075 172.425,260.003 178.023,261.815 183.621,263.521 189.219,265.131 194.817,266.654 200.415,268.097 \n",
       "  206.013,269.465 211.611,270.766 217.209,272.004 222.807,273.184 228.405,274.31 234.003,275.386 239.601,276.415 245.199,277.401 250.797,278.346 256.395,279.252 \n",
       "  261.993,280.123 267.591,280.96 273.189,281.765 278.787,282.54 284.385,283.288 289.983,284.008 295.581,284.703 301.179,285.375 306.777,286.024 312.375,286.651 \n",
       "  317.973,287.258 323.571,287.846 329.169,288.415 334.767,288.967 340.365,289.502 345.963,290.021 351.561,290.525 357.159,291.014 362.757,291.489 368.355,291.951 \n",
       "  373.953,292.401 379.551,292.838 385.149,293.263 390.747,293.677 396.345,294.08 401.943,294.473 407.541,294.856 413.139,295.229 418.737,295.593 424.335,295.948 \n",
       "  429.933,296.295 435.531,296.634 441.129,296.964 446.727,297.287 452.325,297.602 457.923,297.91 463.521,298.212 469.119,298.507 474.717,298.795 480.315,299.077 \n",
       "  \n",
       "  \"/>\n",
       "<polygon clip-path=\"url(#clip8500)\" points=\"\n",
       "208.304,121.976 465.315,121.976 465.315,46.3761 208.304,46.3761 \n",
       "  \" fill=\"#ffffff\" fill-opacity=\"1\"/>\n",
       "<polyline clip-path=\"url(#clip8500)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  208.304,121.976 465.315,121.976 465.315,46.3761 208.304,46.3761 208.304,121.976 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip8500)\" style=\"stroke:#009af9; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  213.304,61.4961 243.304,61.4961 \n",
       "  \"/>\n",
       "<g clip-path=\"url(#clip8500)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:start;\" transform=\"rotate(0, 248.304, 65.9961)\" x=\"248.304\" y=\"65.9961\">numetical integration</text>\n",
       "</g>\n",
       "<polyline clip-path=\"url(#clip8500)\" style=\"stroke:#e26f46; stroke-width:1; stroke-opacity:1; fill:none\" stroke-dasharray=\"8, 5\" points=\"\n",
       "  213.304,76.6161 243.304,76.6161 \n",
       "  \"/>\n",
       "<g clip-path=\"url(#clip8500)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:start;\" transform=\"rotate(0, 248.304, 81.1161)\" x=\"248.304\" y=\"81.1161\">continued fractional approximation n=1</text>\n",
       "</g>\n",
       "<polyline clip-path=\"url(#clip8500)\" style=\"stroke:#3da44d; stroke-width:1; stroke-opacity:1; fill:none\" stroke-dasharray=\"8, 5\" points=\"\n",
       "  213.304,91.7361 243.304,91.7361 \n",
       "  \"/>\n",
       "<g clip-path=\"url(#clip8500)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:start;\" transform=\"rotate(0, 248.304, 96.2361)\" x=\"248.304\" y=\"96.2361\">continued fractional approximation n=2</text>\n",
       "</g>\n",
       "<polyline clip-path=\"url(#clip8500)\" style=\"stroke:#c271d2; stroke-width:1; stroke-opacity:1; fill:none\" stroke-dasharray=\"8, 5\" points=\"\n",
       "  213.304,106.856 243.304,106.856 \n",
       "  \"/>\n",
       "<g clip-path=\"url(#clip8500)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:start;\" transform=\"rotate(0, 248.304, 111.356)\" x=\"248.304\" y=\"111.356\">continued fractional approximation n=3</text>\n",
       "</g>\n",
       "</svg>\n"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x = 0.05:0.05:4\n",
    "@time y = G.(x)\n",
    "@time y_cf1 = G_cf.(x, n=1)\n",
    "@time y_cf2 = G_cf.(x, n=2)\n",
    "@time y_cf3 = G_cf.(x, n=3)\n",
    "plot(size=(500,350))\n",
    "plot!(title=\"y = G(x) = e^x E_1(x)\", titlefontsize=11)\n",
    "plot!(xlims=(0,maximum(x)), ylims=(0,2.3))\n",
    "plot!(x, y, label=\"numetical integration\")\n",
    "plot!(x, y_cf1, label=\"continued fractional approximation n=1\", ls=:dash)\n",
    "plot!(x, y_cf2, label=\"continued fractional approximation n=2\", ls=:dash)\n",
    "plot!(x, y_cf3, label=\"continued fractional approximation n=3\", ls=:dash)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "_draft": {
   "nbviewer_url": "https://gist.github.com/4177ac45105dbd51f53a91555f8b48e8"
  },
  "gist": {
   "data": {
    "description": "07 example of asymptotic expansion",
    "public": true
   },
   "id": "4177ac45105dbd51f53a91555f8b48e8"
  },
  "kernelspec": {
   "display_name": "Julia 0.6.4",
   "language": "julia",
   "name": "julia-0.6"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.6.4"
  },
  "toc": {
   "nav_menu": {
    "height": "207px",
    "width": "242px"
   },
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": true,
   "title_cell": "目次",
   "title_sidebar": "目次",
   "toc_cell": true,
   "toc_position": {
    "height": "calc(100% - 180px)",
    "left": "10px",
    "top": "150px",
    "width": "214px"
   },
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}