sinのpade近似

準備

In [1]:
Pkg.add("Plots") 
Pkg.add("GR")
using Plots

ENV["PLOTS_TEST"] = false

Plots.reset_defaults()
Plots.gr(
    legend=false,
    titlefont=Plots.font("sans-serif", 12),
    legendfont=Plots.font("sans-serif", 8),
    guidefont=Plots.font("sans-serif", 10),
    tickfont=Plots.font("sans-serif", 8),
    markerstrokewidth=0
)
Out[1]:
Plots.GRBackend()
In [2]:
Pkg.add("BenchmarkTools")
using BenchmarkTools
INFO: Package BenchmarkTools is already installed
INFO: METADATA is out-of-date — you may not have the latest version of BenchmarkTools
INFO: Use `Pkg.update()` to get the latest versions of your packages
In [3]:
Pkg.add("SymPy")
using SymPy
INFO: Package SymPy is already installed
INFO: METADATA is out-of-date — you may not have the latest version of SymPy
INFO: Use `Pkg.update()` to get the latest versions of your packages
In [4]:
x = Sym(:x)
Out[4]:
$$x$$

普通にやってみる

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 [5]:
function padesin1(x)
    (x - 7/60*x^3)/(1 + 1/20*x^2)
end
Out[5]:
padesin1 (generic function with 1 method)
In [6]:
padesin1(x)
Out[6]:
$$\frac{- 0.116666666666667 x^{3} + x}{0.05 x^{2} + 1}$$
In [7]:
simplify(padesin1(x))
Out[7]:
$$\frac{x \left(- 0.116666666666667 x^{2} + 1\right)}{0.05 x^{2} + 1}$$
In [8]:
t = collect(-pi:0.01:pi)
y = [padesin1(i) for i in t]
y = hcat(y, [sin(i) for i in t])
Plots.plot(t, y, fmt = :png)
Out[8]:

ここまでで思ったこと

  • なんか項数へってね?
  • 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 [9]:
function padesin2(x)
    (x - 426/3024*x^3 + 25/6048*x^5)/(1 + 13/504*x^2 + 1/10080*x^4)
end
Out[9]:
padesin2 (generic function with 1 method)
In [10]:
padesin2(x)
Out[10]:
$$\frac{0.00413359788359788 x^{5} - 0.140873015873016 x^{3} + x}{9.92063492063492 \cdot 10^{-5} x^{4} + 0.0257936507936508 x^{2} + 1}$$
In [11]:
simplify(padesin2(x))
Out[11]:
$$\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}$$
In [12]:
t = collect(-pi:0.01:pi)
y = [padesin2(i) for i in t]
y = hcat(y, [sin(i) for i in t])
Plots.plot(t, y, fmt = :png)
Out[12]:

まだ気になる

よく見ると$\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 [13]:
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[13]:
padesin3 (generic function with 1 method)
In [14]:
padesin3(x)
Out[14]:
$$x \left(0.183333333333333 x^{2} - 13.7 + \frac{617.4}{x^{2} + 42.0}\right)$$
In [15]:
simplify(padesin3(x))
Out[15]:
$$\frac{x}{x^{2} + 42.0} \left(\left(0.183333333333333 x^{2} - 13.7\right) \left(x^{2} + 42.0\right) + 617.4\right)$$
In [16]:
t = collect(-pi:0.01:pi)
y = [padesin3(i) for i in t]
y = hcat(y, [sin(i) for i in t])
Plots.plot(t, y, fmt = :png)
Out[16]:

さらに式変形

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

In [17]:
function padesin4(angle)
    angle2 = angle^2
    return (11*angle2-252)*(angle2-10)*angle / (2520 + 60 * angle2)
end
Out[17]:
padesin4 (generic function with 1 method)
In [18]:
t = collect(-pi:0.01:pi)
y = [padesin4(i) for i in t]
y = hcat(y, [sin(i) for i in t])
Plots.plot(t, y, fmt = :png)
Out[18]:
In [19]:
function padesin5(angle)
    angle2 = angle^2
    return (10.8*angle2-252)*(angle2-10)*angle / (2520 + 60 * angle2)
end
Out[19]:
padesin5 (generic function with 1 method)
In [20]:
t = collect(-pi:0.01:pi)
y = [padesin5(i) for i in t]
y = hcat(y, [sin(i) for i in t])
Plots.plot(t, y, fmt = :png)
Out[20]:

$\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 [21]:
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[21]:
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 [22]:
b = [1, -1/6, 1/120, -1, 0]
Out[22]:
5-element Array{Float64,1}:
  1.0       
 -0.166667  
  0.00833333
 -1.0       
  0.0       
In [23]:
xs = a\b
Out[23]:
5-element Array{Float64,1}:
  1.0       
 -0.101321  
 -0.0       
  0.0653455 
  0.00255758
In [24]:
function padesin6(x)
    x2 = x^2
    x * (xs[1] + x2*xs[2])/(1 + x2*(xs[4]+ xs[5]*x2))
end
Out[24]:
padesin6 (generic function with 1 method)
In [25]:
t = collect(-pi:0.01:pi)
y = [padesin6(i) for i in t]
y = hcat(y, [sin(i) for i in t])
Plots.plot(t, y, fmt = :png)
Out[25]:

ベンチマーク

In [26]:
function bench(f)
    @time for i = 1:10^7
        f(pi/i)
    end
end

bench(sin)
bench(padesin1)
bench(padesin2)
bench(padesin3)
bench(padesin4)
bench(padesin5)
bench(padesin6)
  0.078229 seconds
  0.000000 seconds
  1.411019 seconds
  0.000000 seconds
  0.000000 seconds
  0.000000 seconds
  3.708497 seconds (160.00 M allocations: 2.384 GiB, 10.52% gc time)

juliaの最適化とかgcの気持ちがわからない…

SincMの近似

そのままやってみる

sinのパデ近似の結果を使って、SincM関数を何とかしてみる.

$$\begin{aligned} SincM(x) &= \frac{sin(Mx)}{sin(x)} \\ &= \frac{60Mx-7M^3x^3}{60 + 3M^2x^2}\times\frac{60 + 3x^2}{60x-7x^3} \\ &= \frac{Mx(60-7M^2x^2)(60+3x^2)}{x(60+3M^2x^2)(60-7x^2)}\\ &= \frac{M(60-7M^2x^2)(60+3x^2)}{(60+3M^2x^2)(60-7x^2)} \end{aligned}$$

In [27]:
function sincm(m, x)
    (m*(60-7*m^2*x^2)*(60+3*x^2))/((60+3*m^2*x^2)*(60-7*x^2))
end
Out[27]:
sincm (generic function with 1 method)
In [28]:
m = Sym(:m)
sincm(m, x)
Out[28]:
$$\frac{m \left(3 x^{2} + 60\right) \left(- 7 m^{2} x^{2} + 60\right)}{\left(- 7 x^{2} + 60\right) \left(3 m^{2} x^{2} + 60\right)}$$

↓ 現実

In [29]:
Plots.plot(sincm(11,x), fmt = :png,)
Out[29]:
In [30]:
sincm(5,0)
Out[30]:
5.0

↓ 理想

In [31]:
Plots.plot(sin(11*x)/sin(x), fmt = :png)
Out[31]:

なんかグラフおかしくね?

うーん. たぶんsin(Mx)がだめ. sin(x)が$[-\pi, \pi]$での近似なのでsin(Mx)は$[-\pi/M, \pi/M]$でしか有効じゃない.

ローラン展開を使う(意味ないが)

なんか特異点周りで近似するローラン展開とかいう式を発見したので使ってみるが、問題が解消されるわけではない.

$$\begin{aligned} SincM(x) &= \frac{sin(Mx)}{sin(x)} \\ &= \frac{60Mx-7M^3x^3}{60 + 3M^2x^2}\times(\frac{1}{x}+\frac{1}{6}x+\frac{7}{360}x^3) \\ &= \frac{Mx(60-7M^2x^2)}{60 + 3M^2x^2}\times(\frac{1}{x}+\frac{1}{6}x+\frac{7}{360}x^3) \\ & = \frac{M(60-7M^2x^2)}{60 + 3M^2x^2}\times(1+\frac{1}{6}x^2+\frac{7}{360}x^4) \end{aligned}$$