2 次元 Padua 点上の多項式補間2010SE058

稲葉三奈

指導教員：杉浦 洋

はじめに1L.Bos ら [1] は、正方形領域上の２次元多項式補間に用いる標本点集合;Padua 点を提案した. 本研究では、Padua点上の Chebyshev 補間と，その係数である Padua 補間を求める高速アルゴリズムを学び, その特性を数値実験で明らかにすることを目標とする.

Padua 点2

［定理 3.1］整数 k ，l ≥ 0，k + l ≤ 2n について，Qn (Tkl ) = Q(Tkl )((k ，l) ̸= (0，2n))，
Qn (T0,2n ) = 1 ̸= 0 = Q(T0,2n )が成立する．Padua 則は 2n − 1 次則である．

4

Padua 点は正方形領域 [−1，1] ⊂ R 上の有限点集合で22

次のように定義される．

An :={ (rnπln(n + 1))

3Padua 則であるので，式 (6) により，{
akl = dk1dl Qn (Tkl f )， (k ，l) ̸= (0，n)，a0n = 2d01dn Qn (T0n f )・正方形領域 Ω = [−1, 1]2 上の積分

Q(f ) =

1π2∫∫
f (x，y) √Ωとなる．//

11
√dxdy
2
1−x
1 − y25N ′′∑fj = f (γn (πj/N ))(A は頂点)，(A は境界点)(A は内点)(3)

は初項と末項に 1/2 を掛けた総和と定義する．すなn ′′
∑
k=0

ak =

1
a0 +2n−1∑
k=1

1
ak + an2(4)n(n+1)∑1πlQn (f ) =))． (5)f (γn (n(n + 1)
n(n + 1)l=0(8)

N ′′2∑
πjk
bk =fj cos
N j=0N(0 ≤ k ≤ N )(9)とする．bk は高速 cosine 変換で高速に計算できる．これを用いて，k, l ≥ 0, k + l ≤ n において，(k ，l ≤ 0［命題 3.2］次の公式が成り立つ．
′′(0 ≤ j ≤ N )と書き，その台形則 cosine 変換をで定義する．これを Padua 則と呼ぶ．
わち，(7)の値をwA f (A)，

1/2
1× 1
wA =n(n + 1) 2∑′′(0 ≤ k ≤ N )で定義される．N = n(n + 1) とし，Padua 点上の関数 fA∈An・2∑πkj
ak =fj cos
N j=0N(2)・Padua 点 An 上の積分則をQn (f ) =高速アルゴリズム

高速 Padua 補間の具体的なアルゴリズムにつて述べる．
N + 1 個のデータ fj (0 ≤ j ≤ N ) の台形則 cosine 変換
ak は，を考える．また，2 変数 Chebyshev 多項式を

Tkl (x，y) = Tk (x)Tl (x)(k ≥ 0，l ≥ 0) で定義する．
［命題 3.1］整数 k ，l ≥ 0 について，{
1 (k = l = 0)，Q(Tkl ) = δk0 δl0 =0 otherwise．//

1Q(Tkl f )dk dl{
2 k=0dk =1 k ̸= 0(1)(0 ≤ l ≤ n(n + 1)) には重複がある．)

πln(n+1)

Padua 点 An 上の n 次 Chebyshev 補間を Padua 補間

akl =要素数は，n(n) + 1)/2．
((γnPadua 補間という．補間係数}0 ≦ l ≦ n(n + 1)) ，γn (t) = (cos nt，cos(n + 1)t)(6)

，k + l ≤ n) のとき，

( ( ))N
1 ∑πjakl =fj Tkl γn
N dk d′l j=0N′′=)1 (bnk+(n+1)l + b|nk−(n+1)l|
4dk dl(10)

1

表 1 計算時間 (秒)
0.001

nN
通常計算
130.016
260.
4
150.
8
450.10916153
0.967
32561
12.262642145213.066
1288385—
256 33153—
512 131841—
1024 525825
—

高速計算

10 -6

0.
0.
0.
0.
0.016
0.062
0.218
0.89
4.196
24.71187.00110 -9

10 -12100

200300400

図 1 f (x, y) = e−a(x +y ) ，a = 1 のとき
22

ÈErÈ
1
0.001
10-6
10-9となる．ここで，d′l =10-12{

2, l = 0, n,1, otherwise

50

100

150

200N

図 2 f (x, y) = sin a(x + y)，a = 1 のときである．
ある．破線は C 補間，実線は超補間の結果である．項数を数値実験

6揃えて比較すると，Padua 補間は C 補間より精度が良い
５節で示したアルゴリズムを Mathematica 上に実現し， ことがわかる．また，項数が多くなるにつれてその差は広実験を行った．この実験は FUJITSU のノートパソコン がる．

FMV-S8390，CPU は intel Core2 Duo P8700，2.53GHzである．OS は Windows 7 で Mathematica ver8.0.1.0 上

7.2振動関数でのプログラムを作成した．

6.1

計算時間

f (x, y) = sin a(x + y),

Padua 補間係数を求める今回の高速アルゴリズムと，通常のアルゴリズムの計算時間を比べる．

n = 2m (0 ≤ m ≤ 10) として係数計算に要した CPU 時間を測定した．結果を表 1 に示す．今回のアルゴリズムは通常のアルゴリズムと比べて圧倒的に速い．n = 64 では計算時間は約 1/1000 となった．通常の方法では計算時間が長すぎて測定を断念した n = 128〜1024 でも，効率的に2 に示す．横軸は項数 N ，縦軸は最大絶対誤差である．破線は C 補間，実線は Padua 補間の結果である．Padua 補間のが精度が良いことがわかる．しかし，丸め誤差の影響により Padua 補間のグラフが水平になっているので項数が多くなると C 補間のがやや精度が良くなる．C 補間の方がやや丸め誤差の影響を受けにくいと思われる．

8補間精度

a = 1, 2, 4の補間を行い，両者を比較した．a = 1 のときの結果を図計算できた．

7振動関数の補間おわりに
本研究では，Mathematica を用いて Padua 補間係数を直積型 Chebyshev 補間 (C 補間) と Padua 超補間 (超補間) の精度を比較した．C 補間の次数 n は 1 ≤ n ≤ 90 で補間多項式の項数は N = (n + 1)2 ，Padua 補間の次数 nは 1 ≤ n ≤ 90 で補間多項式の項数 N = (n + 1)(n + 2)/2求める高速アルゴリズムを実現した．直積型 Chebyshev
補間と精度を比較し，その特徴を調べた．テスト関数としてピーク関数と振動関数を取り上げた．同じ項数で比較すると，Padua 補間が優れていた．

である．近似多項式の計算コストは項数に比例する．被近参考文献

似関数としては，ピーク関数 f (x，y) = e−a(x +y ) と振動9関数 f (x，y) = sin a(x + y) を選んだ．

[1] L.Bos , M.Caliari , S.De Marchi , M.Vianello , Y.Xu :Bivariate Lagrange interpolation at the Padua points: The generating curve approach. Journal of Approximation Theory,Vol.143, pp.15-25(2006).

27.12ピーク関数の補間

ピーク関数

f (x, y) = e−a(x +y ) ,22

a = 1, 8, 16, 25の補間を行い，両者の精度を比較した．a = 1 のときの結果を図 1 に示す．横軸は項数 N ，縦軸は最大絶対誤差で[2] 杉浦洋：数値計算の基礎と応用 [新訂版]，サイエンス社，2009，pp.14-17