{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# sinのpade近似"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 準備"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\u001b[32m\u001b[1m Updating\u001b[22m\u001b[39m registry at `/opt/julia/registries/General`\n",
"\u001b[32m\u001b[1m Updating\u001b[22m\u001b[39m git-repo `https://github.com/JuliaRegistries/General.git`\n",
"\u001b[?25l\u001b[2K\u001b[?25h\u001b[32m\u001b[1m Resolving\u001b[22m\u001b[39m package versions...\n",
"\u001b[32m\u001b[1m Updating\u001b[22m\u001b[39m `/opt/julia/environments/v1.0/Project.toml`\n",
"\u001b[90m [no changes]\u001b[39m\n",
"\u001b[32m\u001b[1m Updating\u001b[22m\u001b[39m `/opt/julia/environments/v1.0/Manifest.toml`\n",
"\u001b[90m [no changes]\u001b[39m\n"
]
}
],
"source": [
"using Pkg\n",
"Pkg.add(\"Plots\") "
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"scrolled": true
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"┌ Info: Recompiling stale cache file /opt/julia/compiled/v1.0/Plots/ld3vC.ji for Plots [91a5bcdd-55d7-5caf-9e0b-520d859cae80]\n",
"└ @ Base loading.jl:1184\n"
]
},
{
"data": {
"text/plain": [
"false"
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"using Plots\n",
"ENV[\"PLOTS_TEST\"] = false"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\u001b[32m\u001b[1m Resolving\u001b[22m\u001b[39m package versions...\n",
"\u001b[32m\u001b[1m Updating\u001b[22m\u001b[39m `/opt/julia/environments/v1.0/Project.toml`\n",
"\u001b[90m [no changes]\u001b[39m\n",
"\u001b[32m\u001b[1m Updating\u001b[22m\u001b[39m `/opt/julia/environments/v1.0/Manifest.toml`\n",
"\u001b[90m [no changes]\u001b[39m\n"
]
}
],
"source": [
"Pkg.add(\"BenchmarkTools\")\n",
"using BenchmarkTools"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\u001b[32m\u001b[1m Resolving\u001b[22m\u001b[39m package versions...\n",
"\u001b[32m\u001b[1m Updating\u001b[22m\u001b[39m `/opt/julia/environments/v1.0/Project.toml`\n",
"\u001b[90m [no changes]\u001b[39m\n",
"\u001b[32m\u001b[1m Updating\u001b[22m\u001b[39m `/opt/julia/environments/v1.0/Manifest.toml`\n",
"\u001b[90m [no changes]\u001b[39m\n"
]
}
],
"source": [
"Pkg.add(\"SymPy\")\n",
"using SymPy"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"\\begin{equation*}x\\end{equation*}"
],
"text/plain": [
"x"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"x = Sym(:x)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 普通にやってみる\n",
"\n",
"sinをマクローリン展開\n",
"\n",
"$$f(x) = sin(x) = \\sum^{\\infty}_{n=0}{\\frac{(-1)^{n-1}}{(2n+1)!}x^{2n+1}} = x - \\frac{x^3}{6} + \\frac{x^5}{120} + O(x^7)$$\n",
"\n",
"近似式を定義\n",
"\n",
"$$g(x) = \\frac{P(x)}{Q(x)} = \\frac{a_0 + a_1x + a_2x^2 + a_3x^3}{1 + b_1x + b_2x^2 + b_3x^3}$$\n",
"\n",
"等しいと仮定して\n",
"\n",
"$$\\begin{aligned}\n",
" f(x) &= g(x) \\\\\n",
" x - \\frac{x^3}{6} + \\frac{x^5}{120} &= \\frac{a_0 + a_1x + a_2x^2 + a_3x^3}{1 + b_1x + b_2x^2 + b_3x^3} \\\\\n",
" \\left(x - \\frac{x^3}{6} + \\frac{x^5}{120}\\right)(1 + b_1x + b_2x^2 + b_3x^3) &= a_0 + a_1x + a_2x^2 + a_3x^3\n",
"\\end{aligned}$$\n",
"\n",
"$P(x)$ の次数 $m$ , $Q(x)$ の次数 $n$ とすると, $x^{m+n}$ までの係数までの係数を比較。\n",
"\n",
"$$\\begin{aligned}\n",
"a_0 &= 0 \\\\\n",
"a_1 &= 1 \\\\\n",
"a_2 &= b_1 \\\\\n",
"a_3 &= b_2 - \\frac{1}{6} \\\\\n",
"b_3 - \\frac{1}{6}b_1 &= 0 \\\\\n",
"-\\frac{1}{6}b_2 + \\frac{1}{120} &= 0 \\\\\n",
"-\\frac{1}{6}b_3 + \\frac{1}{120}b_1 &= 0\n",
"\\end{aligned}$$\n",
"\n",
"これを解いて\n",
"\n",
"$$\n",
"g(x) = \\frac{x - \\frac{7}{60}x^3}{1+\\frac{1}{20}x^2}\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"padesin1 (generic function with 1 method)"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function padesin1(x)\n",
" (x - 7/60*x^3)/(1 + 1/20*x^2)\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"\\begin{equation*}\\frac{- 0.116666666666667 x^{3} + x}{0.05 x^{2} + 1}\\end{equation*}"
],
"text/plain": [
" 3 \n",
"- 0.116666666666667⋅x + x\n",
"──────────────────────────\n",
" 2 \n",
" 0.05⋅x + 1 "
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"padesin1(x)"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"\\begin{equation*}\\frac{x \\left(- 0.116666666666667 x^{2} + 1\\right)}{0.05 x^{2} + 1}\\end{equation*}"
],
"text/plain": [
" ⎛ 2 ⎞\n",
"x⋅⎝- 0.116666666666667⋅x + 1⎠\n",
"──────────────────────────────\n",
" 2 \n",
" 0.05⋅x + 1 "
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"simplify(padesin1(x))"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"t = collect(-pi:0.01:pi)\n",
"y = [padesin1(i) for i in t]\n",
"y = hcat(y, [sin(i) for i in t])\n",
"plot(t, y)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## ここまでで思ったこと\n",
"\n",
"- なんか項数へってね？\n",
"- sin が奇関数であることを利用できないか.\n",
"\n",
"## 有理関数が奇関数となるためには\n",
"\n",
"有理関数\n",
"\n",
"$$\\begin{aligned}\n",
"f(x) &= \\frac{P(x)}{Q(x)} \\\\\n",
"P(x) &= p_0 + p_1x + p_2x^2 + p_3x^3 \\cdots \\\\\n",
"Q(x) &= 1 + q_1x + q_2x^2 + q_3x^3 \\cdots\n",
"\\end{aligned}$$\n",
"\n",
"$f(x)$が奇関数であるとき$f(x) = -f(-x)$であるから\n",
"\n",
"$$\\begin{aligned}\n",
"\\frac{P(x)}{Q(x)} &= - \\frac{P(-x)}{Q(-x)} \\\\\n",
"\\frac{p_0 + p_1x + p_2x^2 + p_3x^3 \\cdots}{1 + q_1x + q_2x^2 + q_3x^3 \\cdots} &= - \\frac{p_0 - p_1x + p_2x^2 - p_3x^3 \\cdots}{1 - q_1x + q_2x^2 - q_3x^3 \\cdots} \\\\\n",
"\\frac{p_0 + p_1x + p_2x^2 + p_3x^3 \\cdots}{1 + q_1x + q_2x^2 + q_3x^3 \\cdots} &= \\frac{- p_0 + p_1x - p_2x^2 + p_3x^3 \\cdots}{1 - q_1x + q_2x^2 - q_3x^3 \\cdots}\n",
"\\end{aligned}$$\n",
"\n",
"この等式を見ると$p_{2n}$と$q_{2n+1}$を0と置くと成り立つ. 項数が減っていたのはこのせい. \n",
"\n",
"## もう一度チャレンジ\n",
"\n",
"奇関数となるよう近似式を定義し、項数を増やしてみる\n",
"\n",
"$$\n",
"g(x) = \\frac{a_0x+a_1x^3+a_2x^5}{1+q_0x^2+q_1x^4}\n",
"$$\n",
"\n",
"さっきと同じように係数比較して頑張ると\n",
"\n",
"$$\n",
"g(x) = \\frac{x-\\frac{426}{3024}x^3+\\frac{25}{6048}x^5}{1+\\frac{13}{504}x^2+\\frac{1}{10080}x^4}\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"padesin2 (generic function with 1 method)"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function padesin2(x)\n",
" (x - 426/3024*x^3 + 25/6048*x^5)/(1 + 13/504*x^2 + 1/10080*x^4)\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"\\begin{equation*}\\frac{0.00413359788359788 x^{5} - 0.140873015873016 x^{3} + x}{9.92063492063492 \\cdot 10^{-5} x^{4} + 0.0257936507936508 x^{2} + 1}\\end{equation*}"
],
"text/plain": [
" 5 3 \n",
"0.00413359788359788⋅x - 0.140873015873016⋅x + x \n",
"──────────────────────────────────────────────────\n",
" 4 2 \n",
"9.92063492063492e-5⋅x + 0.0257936507936508⋅x + 1"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"padesin2(x)"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"\\begin{equation*}\\frac{x \\left(0.00413359788359788 x^{4} - 0.140873015873016 x^{2} + 1\\right)}{9.92063492063492 \\cdot 10^{-5} x^{4} + 0.0257936507936508 x^{2} + 1}\\end{equation*}"
],
"text/plain": [
" ⎛ 4 2 ⎞\n",
"x⋅⎝0.00413359788359788⋅x - 0.140873015873016⋅x + 1⎠\n",
"─────────────────────────────────────────────────────\n",
" 4 2 \n",
" 9.92063492063492e-5⋅x + 0.0257936507936508⋅x + 1 "
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"simplify(padesin2(x))"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"t = collect(-pi:0.01:pi)\n",
"y = [padesin2(i) for i in t]\n",
"y = hcat(y, [sin(i) for i in t])\n",
"plot(t, y)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## まだ気になる\n",
"\n",
"よく見ると$\\frac{1}{10080}x^4$は影響少なそう. なので$x^4$の項を省いて\n",
"\n",
"$$\n",
"g(x) = \\frac{a_0x+a_1x^3+a_2x^5}{1+q_0x^2}\n",
"$$\n",
"\n",
"で近似すると\n",
"\n",
"$$\n",
"g(x) = \\frac{x-\\frac{1}{7}x^3+\\frac{11}{2520}x^5}{1+\\frac{1}{42}x^2}\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"padesin3 (generic function with 1 method)"
]
},
"execution_count": 14,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function padesin3(x)\n",
" # (x - 1/7*x^3 + 11/2520*x^5)/(1 + 1/42*x^2)\n",
" # 計算が少なくなりそうに式変形\n",
" x2 = x^2\n",
" x*(11.0/60.0*x2 - 137.0/10.0 + 37044.0/60.0/(x2+42.0));\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"\\begin{equation*}x \\left(0.183333333333333 x^{2} - 13.7 + \\frac{617.4}{x^{2} + 42.0}\\right)\\end{equation*}"
],
"text/plain": [
" ⎛ 2 617.4 ⎞\n",
"x⋅⎜0.183333333333333⋅x - 13.7 + ─────────⎟\n",
" ⎜ 2 ⎟\n",
" ⎝ x + 42.0⎠"
]
},
"execution_count": 15,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"padesin3(x)"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"\\begin{equation*}\\frac{x}{x^{2} + 42.0} \\left(\\left(0.183333333333333 x^{2} - 13.7\\right) \\left(x^{2} + 42.0\\right) + 617.4\\right)\\end{equation*}"
],
"text/plain": [
" ⎛⎛ 2 ⎞ ⎛ 2 ⎞ ⎞\n",
"x⋅⎝⎝0.183333333333333⋅x - 13.7⎠⋅⎝x + 42.0⎠ + 617.4⎠\n",
"─────────────────────────────────────────────────────\n",
" 2 \n",
" x + 42.0 "
]
},
"execution_count": 16,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"simplify(padesin3(x))"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 17,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"t = collect(-pi:0.01:pi)\n",
"y = [padesin3(i) for i in t]\n",
"y = hcat(y, [sin(i) for i in t])\n",
"plot(t, y)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## さらに式変形\n",
"\n",
"近似的に因数分解してみた"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"padesin4 (generic function with 1 method)"
]
},
"execution_count": 18,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function padesin4(angle)\n",
" angle2 = angle^2\n",
" return (11*angle2-252)*(angle2-10)*angle / (2520 + 60 * angle2)\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"\\begin{equation*}\\frac{x \\left(x^{2} - 10\\right) \\left(11 x^{2} - 252\\right)}{60 x^{2} + 2520}\\end{equation*}"
],
"text/plain": [
" ⎛ 2 ⎞ ⎛ 2 ⎞\n",
"x⋅⎝x - 10⎠⋅⎝11⋅x - 252⎠\n",
"─────────────────────────\n",
" 2 \n",
" 60⋅x + 2520 "
]
},
"execution_count": 19,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"padesin4(x)"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 20,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"t = collect(-pi:0.01:pi)\n",
"y = [padesin4(i) for i in t]\n",
"y = hcat(y, [sin(i) for i in t])\n",
"plot(t, y)"
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"padesin5 (generic function with 1 method)"
]
},
"execution_count": 21,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function padesin5(angle)\n",
" angle2 = angle^2\n",
" return (10.8*angle2-252)*(angle2-10)*angle / (2520 + 60 * angle2)\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"\\begin{equation*}\\frac{x \\left(x^{2} - 10\\right) \\left(10.8 x^{2} - 252\\right)}{60 x^{2} + 2520}\\end{equation*}"
],
"text/plain": [
" ⎛ 2 ⎞ ⎛ 2 ⎞\n",
"x⋅⎝x - 10⎠⋅⎝10.8⋅x - 252⎠\n",
"───────────────────────────\n",
" 2 \n",
" 60⋅x + 2520 "
]
},
"execution_count": 22,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"padesin5(x)"
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 23,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"t = collect(-pi:0.01:pi)\n",
"y = [padesin5(i) for i in t]\n",
"y = hcat(y, [sin(i) for i in t])\n",
"plot(t, y)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## $\\pm\\pi$ で連続的につながるような関数\n",
"\n",
"$\\pm\\pi$で0かつ傾きが同じだったらきれいになりそうな感じがしたので、条件に加える。\n",
"\n",
"$$\n",
"\\begin{aligned}\n",
"g(x) &= \\frac{p_0x+p_1x^3+p_2x^5}{1+q_0x^2+q_1x^4} \\\\\n",
" &= x-\\frac{x^3}{6}+\\frac{x^5}{120} \\\\\n",
"g(\\pi) &= g(-\\pi) = 0 \\\\\n",
"g'(\\pi) &= g'(-\\pi) = \\sin'(\\pi)\n",
"\\end{aligned}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"めんどくさくなりそうな気がしたので、juliaに連立方程式を解かせる。"
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"5×5 Array{Float64,2}:\n",
" 1.0 0.0 0.0 0.0 0.0\n",
" 0.0 1.0 0.0 -1.0 0.0\n",
" 0.0 0.0 1.0 0.166667 -1.0\n",
" 1.0 19.7392 487.045 0.0 0.0\n",
" 1.0 9.8696 97.4091 0.0 0.0"
]
},
"execution_count": 24,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"a = [1 0 0 0 0 ;\n",
" 0 1 0 -1 0 ;\n",
" 0 0 1 1/6 -1 ;\n",
" 1 2*pi^2 5*pi^4 0 0 ;\n",
" 1 pi^2 pi^4 0 0 ]"
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"5-element Array{Float64,1}:\n",
" 1.0 \n",
" -0.16666666666666666 \n",
" 0.008333333333333333\n",
" -1.0 \n",
" 0.0 "
]
},
"execution_count": 25,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"b = [1, -1/6, 1/120, -1, 0]"
]
},
{
"cell_type": "code",
"execution_count": 26,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"5-element Array{Float64,1}:\n",
" 1.0 \n",
" -0.10132118364233778 \n",
" -0.0 \n",
" 0.06534548302432888 \n",
" 0.0025575805040548138"
]
},
"execution_count": 26,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"xs = a\\b"
]
},
{
"cell_type": "code",
"execution_count": 27,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"padesin6 (generic function with 1 method)"
]
},
"execution_count": 27,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function padesin6(x)\n",
" global xs\n",
" x2 = x^2\n",
" x * (xs[1] + x2*xs[2])/(1 + x2*(xs[4]+ xs[5]*x2))\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 28,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 28,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"t = collect(-pi:0.01:pi)\n",
"y = [padesin6(i) for i in t]\n",
"y = hcat(y, [sin(i) for i in t])\n",
"plot(t, y)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 参考\n",
"- [Padé Approximant - WolframMathWorld](http://mathworld.wolfram.com/PadeApproximant.html)\n",
"- [Cosecant - WolframMathWorld](http://mathworld.wolfram.com/Cosecant.html)\n",
"- [pade.mws](http://homepages.math.uic.edu/~jan/MCS471/Lec25/pade.html)"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Julia 1.0.0",
"language": "julia",
"name": "julia-1.0"
},
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "1.0.0"
}
},
"nbformat": 4,
"nbformat_minor": 2
}