sinのpade近似

準備

In [1]:
using Pkg
Pkg.add("Plots")
  Updating registry at `/opt/julia/registries/General`
  Updating git-repo `https://github.com/JuliaRegistries/General.git`
 Resolving package versions...
  Updating `/opt/julia/environments/v1.0/Project.toml`
 [no changes]
  Updating `/opt/julia/environments/v1.0/Manifest.toml`
 [no changes]
In [2]:
using Plots
ENV["PLOTS_TEST"] = false
┌ Info: Recompiling stale cache file /opt/julia/compiled/v1.0/Plots/ld3vC.ji for Plots [91a5bcdd-55d7-5caf-9e0b-520d859cae80]
└ @ Base loading.jl:1184
Out[2]:
false
In [3]:
Pkg.add("BenchmarkTools")
using BenchmarkTools
 Resolving package versions...
  Updating `/opt/julia/environments/v1.0/Project.toml`
 [no changes]
  Updating `/opt/julia/environments/v1.0/Manifest.toml`
 [no changes]
In [4]:
Pkg.add("SymPy")
using SymPy
 Resolving package versions...
  Updating `/opt/julia/environments/v1.0/Project.toml`
 [no changes]
  Updating `/opt/julia/environments/v1.0/Manifest.toml`
 [no changes]
In [5]:
x = Sym(:x)
Out[5]:
\begin{equation*}x\end{equation*}

普通にやってみる

sinをマクローリン展開

$$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)$$

近似式を定義

$$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}$$

等しいと仮定して

$$\begin{aligned} f(x) &= g(x) \\ 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} \\ \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 \end{aligned}$$

$P(x)$ の次数 $m$ , $Q(x)$ の次数 $n$ とすると, $x^{m+n}$ までの係数までの係数を比較。

$$\begin{aligned} a_0 &= 0 \\ a_1 &= 1 \\ a_2 &= b_1 \\ a_3 &= b_2 - \frac{1}{6} \\ b_3 - \frac{1}{6}b_1 &= 0 \\ -\frac{1}{6}b_2 + \frac{1}{120} &= 0 \\ -\frac{1}{6}b_3 + \frac{1}{120}b_1 &= 0 \end{aligned}$$

これを解いて

$$ g(x) = \frac{x - \frac{7}{60}x^3}{1+\frac{1}{20}x^2} $$

In [6]:
function padesin1(x)
    (x - 7/60*x^3)/(1 + 1/20*x^2)
end
Out[6]:
padesin1 (generic function with 1 method)
In [7]:
padesin1(x)
Out[7]:
\begin{equation*}\frac{- 0.116666666666667 x^{3} + x}{0.05 x^{2} + 1}\end{equation*}
In [8]:
simplify(padesin1(x))
Out[8]:
\begin{equation*}\frac{x \left(- 0.116666666666667 x^{2} + 1\right)}{0.05 x^{2} + 1}\end{equation*}
In [9]:
t = collect(-pi:0.01:pi)
y = [padesin1(i) for i in t]
y = hcat(y, [sin(i) for i in t])
plot(t, y)
Out[9]:
-3 -2 -1 0 1 2 3 -1.0 -0.5 0.0 0.5 1.0 y1 y2

ここまでで思ったこと

  • なんか項数へってね?
  • sin が奇関数であることを利用できないか.

有理関数が奇関数となるためには

有理関数

$$\begin{aligned} f(x) &= \frac{P(x)}{Q(x)} \\ P(x) &= p_0 + p_1x + p_2x^2 + p_3x^3 \cdots \\ Q(x) &= 1 + q_1x + q_2x^2 + q_3x^3 \cdots \end{aligned}$$

$f(x)$が奇関数であるとき$f(x) = -f(-x)$であるから

$$\begin{aligned} \frac{P(x)}{Q(x)} &= - \frac{P(-x)}{Q(-x)} \\ \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} \\ \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} \end{aligned}$$

この等式を見ると$p_{2n}$と$q_{2n+1}$を0と置くと成り立つ. 項数が減っていたのはこのせい. 

もう一度チャレンジ

奇関数となるよう近似式を定義し、項数を増やしてみる

$$ g(x) = \frac{a_0x+a_1x^3+a_2x^5}{1+q_0x^2+q_1x^4} $$

さっきと同じように係数比較して頑張ると

$$ 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} $$

In [10]:
function padesin2(x)
    (x - 426/3024*x^3 + 25/6048*x^5)/(1 + 13/504*x^2 + 1/10080*x^4)
end
Out[10]:
padesin2 (generic function with 1 method)
In [11]:
padesin2(x)
Out[11]:
\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*}
In [12]:
simplify(padesin2(x))
Out[12]:
\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*}
In [13]:
t = collect(-pi:0.01:pi)
y = [padesin2(i) for i in t]
y = hcat(y, [sin(i) for i in t])
plot(t, y)
Out[13]:
-3 -2 -1 0 1 2 3 -1.0 -0.5 0.0 0.5 1.0 y1 y2

まだ気になる

よく見ると$\frac{1}{10080}x^4$は影響少なそう. なので$x^4$の項を省いて

$$ g(x) = \frac{a_0x+a_1x^3+a_2x^5}{1+q_0x^2} $$

で近似すると

$$ g(x) = \frac{x-\frac{1}{7}x^3+\frac{11}{2520}x^5}{1+\frac{1}{42}x^2} $$

In [14]:
function padesin3(x)
    # (x - 1/7*x^3 + 11/2520*x^5)/(1 + 1/42*x^2)
    # 計算が少なくなりそうに式変形
    x2 = x^2
    x*(11.0/60.0*x2 - 137.0/10.0 + 37044.0/60.0/(x2+42.0));
end
Out[14]:
padesin3 (generic function with 1 method)
In [15]:
padesin3(x)
Out[15]:
\begin{equation*}x \left(0.183333333333333 x^{2} - 13.7 + \frac{617.4}{x^{2} + 42.0}\right)\end{equation*}
In [16]:
simplify(padesin3(x))
Out[16]:
\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*}
In [17]:
t = collect(-pi:0.01:pi)
y = [padesin3(i) for i in t]
y = hcat(y, [sin(i) for i in t])
plot(t, y)
Out[17]:
-3 -2 -1 0 1 2 3 -1.0 -0.5 0.0 0.5 1.0 y1 y2

さらに式変形

近似的に因数分解してみた

In [18]:
function padesin4(angle)
    angle2 = angle^2
    return (11*angle2-252)*(angle2-10)*angle / (2520 + 60 * angle2)
end
Out[18]:
padesin4 (generic function with 1 method)
In [19]:
padesin4(x)
Out[19]:
\begin{equation*}\frac{x \left(x^{2} - 10\right) \left(11 x^{2} - 252\right)}{60 x^{2} + 2520}\end{equation*}
In [20]:
t = collect(-pi:0.01:pi)
y = [padesin4(i) for i in t]
y = hcat(y, [sin(i) for i in t])
plot(t, y)
Out[20]:
-3 -2 -1 0 1 2 3 -1.0 -0.5 0.0 0.5 1.0 y1 y2
In [21]:
function padesin5(angle)
    angle2 = angle^2
    return (10.8*angle2-252)*(angle2-10)*angle / (2520 + 60 * angle2)
end
Out[21]:
padesin5 (generic function with 1 method)
In [22]:
padesin5(x)
Out[22]:
\begin{equation*}\frac{x \left(x^{2} - 10\right) \left(10.8 x^{2} - 252\right)}{60 x^{2} + 2520}\end{equation*}
In [23]:
t = collect(-pi:0.01:pi)
y = [padesin5(i) for i in t]
y = hcat(y, [sin(i) for i in t])
plot(t, y)
Out[23]:
-3 -2 -1 0 1 2 3 -1.0 -0.5 0.0 0.5 1.0 y1 y2

$\pm\pi$ で連続的につながるような関数

$\pm\pi$で0かつ傾きが同じだったらきれいになりそうな感じがしたので、条件に加える。

$$ \begin{aligned} g(x) &= \frac{p_0x+p_1x^3+p_2x^5}{1+q_0x^2+q_1x^4} \\ &= x-\frac{x^3}{6}+\frac{x^5}{120} \\ g(\pi) &= g(-\pi) = 0 \\ g'(\pi) &= g'(-\pi) = \sin'(\pi) \end{aligned} $$

めんどくさくなりそうな気がしたので、juliaに連立方程式を解かせる。

In [24]:
a = [1 0      0      0      0     ;
     0 1      0      -1     0     ;
     0 0      1      1/6    -1    ;
     1 2*pi^2 5*pi^4 0      0     ;
     1 pi^2   pi^4   0      0      ]
Out[24]:
5×5 Array{Float64,2}:
 1.0   0.0       0.0      0.0        0.0
 0.0   1.0       0.0     -1.0        0.0
 0.0   0.0       1.0      0.166667  -1.0
 1.0  19.7392  487.045    0.0        0.0
 1.0   9.8696   97.4091   0.0        0.0
In [25]:
b = [1, -1/6, 1/120, -1, 0]
Out[25]:
5-element Array{Float64,1}:
  1.0                 
 -0.16666666666666666 
  0.008333333333333333
 -1.0                 
  0.0                 
In [26]:
xs = a\b
Out[26]:
5-element Array{Float64,1}:
  1.0                  
 -0.10132118364233778  
 -0.0                  
  0.06534548302432888  
  0.0025575805040548138
In [27]:
function padesin6(x)
    global xs
    x2 = x^2
    x * (xs[1] + x2*xs[2])/(1 + x2*(xs[4]+ xs[5]*x2))
end
Out[27]:
padesin6 (generic function with 1 method)
In [28]:
t = collect(-pi:0.01:pi)
y = [padesin6(i) for i in t]
y = hcat(y, [sin(i) for i in t])
plot(t, y)
Out[28]:
-3 -2 -1 0 1 2 3 -1.0 -0.5 0.0 0.5 1.0 y1 y2