Chebyshev多項式の有限和で表される関数を本稿ではChebyshev補間という。Chebyshev補間
$$p(x)=\sum^{M-1}_{n=0}a_nT_n(x)$$に対して、その微分
$$p'(x)=\sum_{n=0}^{M-2}b_nT_n(x),\quad {}^\prime = \frac d{dx}$$を求める。係数 $a_n$ と $b_n$ の関係を調べる。
$$p(x)=\int p'(x)dx=\int \sum_{n=0}^{M-2}b_nT_n(x)dx $$であるから、$T_n(x)=\cos n\theta$ $(x=\cos\theta)$ という事実から
$$ \begin{align*} \int T_n(x)dx&=\int -\cos n\theta\sin\theta d\theta \\ &=-\frac{1}{2}\int(\sin (n+1)\theta-\sin (n-1)\theta)d\theta \\ &=\frac{1}{2}\left[\frac{\cos (n+1)\theta}{n+1}-\frac{\cos (n-1)\theta}{n-1}\right]\\ &=\frac{1}{2}\left[\frac{T_{n+1}(x)}{n+1}-\frac{T_{n-1}(x)}{n-1}\right]. \end{align*} $$よって積分による定数倍を除けば、
$$ \begin{align*} \int T_n(x)dx=\begin{cases} \frac{1}{2}\left[\frac{T_{n+1}(x)}{n+1}-\frac{T_{n-1}(x)}{n-1}\right] & (n\ge 2)\\ \frac14 T_2(x) & (n=1)\\ T_1(x) & (n=0). \end{cases} \end{align*} $$従って、 $$ \begin{align*} p(x)=\int p(x)'dx&=C +b_0T_1(x)+\frac{1}{4}b_1T_2(x)+\sum^{M-2}_{n=2}\frac{b_n}{2}\left[\frac{T_{n+1}(x)}{n+1}-\frac{T_{n-1}(x)}{n-1}\right] \\ &=\sum^{M-1}_{n=0}a_nT_n(x). \end{align*} $$
つまり $a_n$ と $b_n$ の関係は以下の漸化式で表すことができる。
$$ \begin{align*} b_{n-1}&=b_{n+1}+2na_n,\quad(n=M-1,M-2,\dots,2),\quad b_{M+1}=b_{M}=0,\\ 2b_0&=b_2+2a_1. \end{align*} $$注 最後の $b_0$ の項は $b_0'=2b_0$ とおいて漸化式を $n=M-1,M-2,\dots,1$ として計算し、最後に $b_0=b_0'/2$ と求めても同じである（以下の実装ではそのようにしてる。
以上の操作を関数 $p$ のChebyshev係数 $a_n$ から $p'$ のChebyshev係数 $b_n$ を求めるchebdiff関数としてまとめる。
versioninfo()
Julia Version 1.8.3 Commit 0434deb161e (2022-11-14 20:14 UTC) Platform Info: OS: macOS (arm64-apple-darwin21.3.0) CPU: 10 × Apple M1 Max WORD_SIZE: 64 LIBM: libopenlibm LLVM: libLLVM-13.0.1 (ORCJIT, apple-m1) Threads: 1 on 8 virtual cores
function chebdiff(a; I=[-1,1])# Input is Two-sided
M = length(a)
b = zeros(M+1)
for r = M-1:-1:1
b[r] = b[r+2] + 2*r*a[r+1]
end
b[1] /= 2.0
return b[setdiff(1:end,end)]*(2/(I[2]-I[1])) # Output is Two-sided
end
chebdiff (generic function with 1 method)
この関数を用いて、$e^{x}$ のChebyshev補間を何階か微分し、各係数を比較すると
include("FourierChebyshev.jl")
# f(x) = sinpi(x)
# reshape([-pi^2*a;chebdiff(chebdiff(a))],20,2)
f(x) = exp(x)
a = cheb(f)
reshape([a;chebdiff(a);chebdiff(chebdiff(a))],length(a),3)
14×3 Matrix{Float64}:
1.26607 1.26607 1.26607
1.13032 1.13032 1.13032
0.271495 0.271495 0.271495
0.0443368 0.0443368 0.0443368
0.00547424 0.00547424 0.00547424
0.000542926 0.000542926 0.000542926
4.49773e-5 4.49773e-5 4.49773e-5
3.19844e-6 3.19844e-6 3.19844e-6
1.99212e-7 1.99212e-7 1.99211e-7
1.10368e-8 1.10368e-8 1.10367e-8
5.5059e-10 5.5059e-10 5.4933e-10
2.49796e-11 2.49696e-11 2.48974e-11
1.0404e-12 1.03739e-12 0.0
3.98997e-14 0.0 0.0
maximum(abs,a-chebdiff(chebdiff(a)))
1.7396084572851578e-12
上の係数比較から分かるように、1回微分操作を行うと最後の係数が0になる。そのため微分操作は補間の精度を低下させる操作であることに注意する。
さらに定義域が $x\in [a,b]$ であるときは、変数変換 $[a,b]\to[-1,1]$ のスケール定数がかかる。すなわち
$$ x\mapsto \xi = 2\frac{x-a}{h}-1\quad\Rightarrow\quad d\xi=\frac2h dx,\quad h = b-a $$から
$$ p'(x) = \frac{d}{dx}\sum_{n=0}^{M-1}a_n T_n(x) = \frac{d}{d\xi}\sum_{n=0}^{M-1}a_n T_n(\xi) \frac{d\xi}{dx} = \sum_{n=0}^{M-2}b_n T_n(\xi) \frac2h = \sum_{n=0}^{M-2}\left(b_n \frac2h\right) T_n(x) $$として、係数を計算する。
using SpecialFunctions
f(x) = 2*exp(-x^2)/sqrt(π); dom = [-5,5]
a = cheb(x->erf(x),dom)
plot_cheb(cheb(f,dom),I=dom,label="\$2e^{-x^2}/\\sqrt{\\pi}\$")
plot_cheb!(chebdiff(a,I=dom),I=dom,label="\$\\frac{d}{dx}\\mathrm{erf}(x)\$")
もしもChebyshev補間 $p(x)$ がOne-sided Chebyshev補間
$$ p(x) = c_0 + 2\sum_{n=1}^{M-1}c_n T_n(x),\quad\begin{cases}c_0=a_0 & (n=0)\\ c_n=\frac{a_n}2 & (n\ge 1)\end{cases} $$で表されている場合は、上記の漸化式を $n=1$ まで計算すればよい。
$$ b_{n-1}=b_{n+1}+2nc_n,\quad(n=M-1,M-2,\dots,1),\quad b_{M+1}=b_{M}=0. $$この操作を実装する。
function chebdiff_oneside(a; I=[-1,1])# Input is One-sided
M = length(a)
b = zeros(M+1)
for r = M-1:-1:1
b[r] = b[r+2] + 2*r*a[r+1]
end
return b[setdiff(1:end,end)]*(2/(I[2]-I[1])) # Output is One-sided
end
chebdiff_oneside (generic function with 1 method)
f(x) = exp(x)
a = cheb(f) # "cheb" function returns Two-sided coefficients
c = a[:]
c[2:end] /= 2.
reshape([c;chebdiff_oneside(c);chebdiff_oneside(chebdiff_oneside(c))],length(a),3)
14×3 Matrix{Float64}:
1.26607 1.26607 1.26607
0.565159 0.565159 0.565159
0.135748 0.135748 0.135748
0.0221684 0.0221684 0.0221684
0.00273712 0.00273712 0.00273712
0.000271463 0.000271463 0.000271463
2.24887e-5 2.24887e-5 2.24887e-5
1.59922e-6 1.59922e-6 1.59922e-6
9.96062e-8 9.96062e-8 9.96055e-8
5.51839e-9 5.51838e-9 5.51834e-9
2.75295e-10 2.75295e-10 2.74665e-10
1.24898e-11 1.24848e-11 1.24487e-11
5.20199e-13 5.18696e-13 0.0
1.99499e-14 0.0 0.0
上で紹介した微分操作を区間演算と組み合わせることで、Chebyshev補間の微分を精度保証付き数値計算することができる。
using IntervalArithmetic,LinearAlgebra
function chebdiff(ia::Vector{Interval{T}}; I=[-1,1])where T # Input is Two-sided (inverval)
M = length(a)
ib = zeros(Interval{T},M+1)
for r = M-1:-1:1
ib[r] = ib[r+2] + 2*r*ia[r+1]
end
ib[1] /= 2.0
return ib[1:end-2]*(2/(interval(I[2])-I[1])) # Output is Two-sided (interval)
end
chebdiff (generic function with 2 methods)
f(x) = 1/(1+1000*(x+.5)^2)+1/sqrt(1+1000*(x-.5)^2)
# f(x) = exp(x)
# using SpecialFunctions
# f(x) = exp(erf(x^2)+x.^5).*sin(3*pi*x) + x
a = cheb(f)
maximum(radius,chebdiff(map(Interval,a)))
4.3298697960381105e-15
Clenshawのアルゴリズムで観測されたように補間の項数が多くなると区間演算を実行する回数が多くなり、区間の幅が増大する予想が考えられるが、今回は予想に反して微分操作で混入する誤差が小さくこの実装で十分である。
Chebyshev補間の微分は第2種Chebyshev多項式を用いると綺麗に表記できる。すなわち
$$ \frac{d}{dx}T_{n}(x) = nU_{n-1}(x) $$という関係式から（$U_n(x)$: 第2種Chebyshev多項式）
$$ p'(x)=\sum_{n=0}^{M-2}(n+1)a_{n+1} U_n(x) $$という関数の表現ができる。第1種Chebyshev多項式に拘らなければ、こちらの方が表現が簡潔で良い。
function chebdiff_secondkind(a; I=[-1,1]) # Input is Two-sided
M = length(a)
b = zeros(M-1)
for n = 0:M-2
b[n+1] = (n+1)*a[n+2]
end
return b*(2/(I[2]-I[1])) # Output is second kind (Two-sided)
end
chebdiff_secondkind (generic function with 1 method)
# f(x) = exp(x)
using SpecialFunctions
f(x) = exp(erf(x^2)+x.^5).*sin(3*pi*x) + x; dom = [-1.5,1]
a = cheb(f,dom)
b = chebdiff_secondkind(a,I=dom)
67-element Vector{Float64}:
1.5362508753531374
0.8726471810871159
0.6452024469103549
0.24253352085466193
-0.8857802833980725
-2.981737199638248
-2.9954225795004685
-3.2410567792303526
-3.5122672555296575
2.3806940321802212
4.151953577006191
-2.323259803879892
-1.9342326735501818
⋮
4.581952860245287e-11
-2.2754392996681254e-10
1.1738751397130197e-10
-3.831797897222334e-12
-3.024494088013351e-11
1.824550726092612e-11
-2.0409332598740264e-12
-3.818318795473811e-12
2.7827822964694075e-12
-5.183780436411e-13
-5.151540885415318e-13
3.9328540424321546e-13
係数だけを見ても分からないので、Clenshawのアルゴリズムを使って描画をしてみる。
function clenshaw_secondkind(a,x) # Clenshaw's algorithm
# a: (Two-sided) Chbyshev coefficients
# x: evaluating points in [-1,1]
n = length(a)-1
bk0 = 0
bk1 = 0
for r = (n+1):-1:1
tmp = 2x.*bk0 .- bk1 .+ a[r]
bk1 = bk0
bk0 = tmp
end
return bk0 #.+ bk1*(-x)
end
clenshaw_secondkind (generic function with 1 method)
一方は、chebdiff関数を使って描画した図、もう一方は、第2種Chebyshev多項式の係数を描画したもの。
plot_cheb(chebdiff(a,I=dom),I=dom,label="chebdiff")
x = range(dom[1],stop=dom[2],length=5000)
ξ = 2*(x.-dom[1])/(dom[2]-dom[1]) .- 1
fx = clenshaw_secondkind(b,ξ)
plot!(x,fx,label="second kind")
見た目はほぼ一致している。さらに、各点（5000点）での関数値を比較して差をとるとほぼ同じであることがわかる。
norm(eval_cheb(chebdiff(a,I=dom),x,I=dom)-fx,Inf)
8.526512829121202e-14
本資料は以下のような文献・Web ページ等を参考に書いています。