文章目录
引言
起因是数值分析的一次大作业,在学习了分段三次一元 H e r m i t e Hermite Hermite插值并编程实践后,老师布置了如下两个题目:
要求:
自己构造相应的插值函数是困难的,而且不好验证,于是先转向文献调研,再根据调研结果进行分析和编程
最初想法与文献调研
由于题目给出了
f
(
x
,
y
)
,
∂
f
∂
x
(
x
,
y
)
,
∂
f
∂
y
(
x
,
y
)
,
∂
2
f
∂
x
∂
y
(
x
,
y
)
f(x,y),\dfrac{\partial f}{\partial x}(x,y),\dfrac{\partial f}{\partial y}(x,y),\dfrac{\partial^2 f}{\partial x\partial y}(x,y)
f(x,y),∂x∂f(x,y),∂y∂f(x,y),∂x∂y∂2f(x,y),并且需要考虑分块插值,一个很自然的想法是考虑
H
e
r
m
i
t
e
Hermite
Hermite插值,当然,
H
e
r
m
i
t
e
Hermite
Hermite插值可以看作是
L
a
g
r
a
n
g
e
Lagrange
Lagrange和
N
e
w
t
o
n
Newton
Newton插值的极限形式,这两种插值也可以考虑,而在文献查阅过程中,考虑
L
a
g
r
a
n
g
e
Lagrange
Lagrange插值一般不考虑导数;而大部分文献似乎并没有将
N
e
w
t
o
n
Newton
Newton插值作为单独的一类插值,而是“
N
e
w
t
o
n
Newton
Newton形式的插值”[1]:
p
(
x
,
y
)
=
∑
i
=
0
m
∑
j
=
0
n
f
[
x
0
.
⋯
,
x
i
;
y
0
,
⋯
,
y
j
]
∏
h
=
0
i
−
1
(
x
−
x
h
)
∏
k
=
0
j
−
1
(
x
−
x
k
)
p(x,y)=\sum_{i=0}^m\sum_{j=0}^n f[x_0.\cdots,x_i;y_0,\cdots,y_j]\prod_{h=0}^{i-1}(x-x_h)\prod_{k=0}^{j-1}(x-x_k)
p(x,y)=i=0∑mj=0∑nf[x0.⋯,xi;y0,⋯,yj]h=0∏i−1(x−xh)k=0∏j−1(x−xk)
但是[1]中提供的关于
N
e
w
t
o
n
Newton
Newton插值的文献大部分是不能访问的,又或者是没有基于导数的。
当然,其实最暴力的方法是线性方程组求导[2],尽管题目不允许,但这种方法可以帮助我们判断插值方程是否唯一(对于第二题);样条插值也是一种可以考虑的方法,但是如果 H e r m i t e Hermite Hermite插值可行且线性方程组的解唯一,就没必要构造样条插值了。
针对题目,仅考虑二元 H e r m i t e Hermite Hermite插值,而事实上,在实际应用中,一般也考虑二元插值,最多考虑三元,二元插值也是研究的最多的[3]。
二元 H e r m i t e Hermite Hermite插值的两种“范式”
回顾一元
H
e
r
m
i
t
e
Hermite
Hermite插值的构造过程(方便起见,只考虑一阶导),对于已知
f
(
x
0
)
,
f
(
x
1
)
,
⋯
f
(
x
n
)
;
f
′
(
x
0
)
,
f
′
(
x
1
)
,
⋯
f
′
(
x
n
)
f(x_0),f(x_1),\cdots f(x_n);f'(x_0),f'(x_1),\cdots f'(x_n)
f(x0),f(x1),⋯f(xn);f′(x0),f′(x1),⋯f′(xn)的情形,构造插值函数
h
(
x
)
=
∑
i
=
0
n
A
i
(
x
)
f
(
x
i
)
+
∑
i
=
0
n
B
i
(
x
)
f
′
(
x
i
)
h(x)=\sum_{i=0}^nA_i(x)f(x_i)+\sum_{i=0}^nB_i(x)f'(x_i)
h(x)=i=0∑nAi(x)f(xi)+i=0∑nBi(x)f′(xi)
使其满足
h
(
x
i
)
=
f
(
x
i
)
,
h
′
(
x
i
)
=
f
′
(
x
i
)
,
i
=
0.1
,
⋯
,
n
h(x_i)=f(x_i),h'(x_i)=f'(x_i),i=0.1,\cdots,n
h(xi)=f(xi),h′(xi)=f′(xi),i=0.1,⋯,n,可以构造
A
i
(
x
)
=
[
1
−
2
l
i
′
(
x
i
)
(
x
−
x
i
)
]
l
i
2
(
x
)
B
i
(
x
)
=
(
x
−
x
i
)
l
i
2
(
x
)
\begin{align} &A_i(x)=[1-2l_i'(x_i)(x-x_i)]l_i^2(x)\\ &B_i(x)=(x-x_i)l_i^2(x) \end{align}
Ai(x)=[1−2li′(xi)(x−xi)]li2(x)Bi(x)=(x−xi)li2(x)
其中,
l
i
(
x
)
l_i(x)
li(x)为
L
a
g
r
a
n
g
e
Lagrange
Lagrange插值基函数
以下介绍两种二元 H e r m i t e Hermite Hermite插值“范式”(作者自称),方便起见,按作者起名为 A h l i n Ahlin Ahlin范式[4]和 C h a w l a − J a y a r a j a n Chawla-Jayarajan Chawla−Jayarajan范式[5]。前者的推导思想与一元情形类似,而后者是基于前者的改进
在构造插值函数的过程中,两者都参考了 S p i t z b a r t Spitzbart Spitzbart[6]在已知任意阶导数的情况下一元 H e r m i t e Hermite Hermite插值的构造
S p i t z b a r t Spitzbart Spitzbart的构造
假设已知
f
(
k
)
(
x
j
)
,
j
=
0
,
1
,
⋯
,
n
;
k
=
0
,
1
,
⋯
,
r
j
f^{(k)}(x_j),j=0,1,\cdots ,n;k=0,1,\cdots,r_j
f(k)(xj),j=0,1,⋯,n;k=0,1,⋯,rj,设
p
j
(
x
)
=
(
x
−
x
0
)
r
0
+
1
⋯
(
x
−
x
n
)
r
n
+
1
g
j
(
x
)
=
[
p
j
(
x
)
]
−
1
p_j(x)=(x-x_0)^{r_0+1}\cdots(x-x_n)^{r_n+1}\\ g_j(x)=[p_j(x)]^{-1}
pj(x)=(x−x0)r0+1⋯(x−xn)rn+1gj(x)=[pj(x)]−1
可以构造
n
+
∑
j
=
0
n
r
j
n+\sum\limits_{j=0}^nr_j
n+j=0∑nrj次多项式函数
H
(
x
)
H(x)
H(x),满足
H
(
k
)
(
x
j
)
=
f
(
k
)
(
x
j
)
H^{(k)}(x_j)=f^{(k)}(x_j)
H(k)(xj)=f(k)(xj)
H
(
x
)
H(x)
H(x)如下:
H
(
x
)
=
∑
j
=
0
n
∑
k
=
0
r
j
A
j
k
f
(
k
)
(
x
j
)
H(x)=\sum_{j=0}^n\sum_{k=0}^{r_j}A_{jk}f^{(k)}(x_j)
H(x)=j=0∑nk=0∑rjAjkf(k)(xj)
其中
A
j
k
(
x
)
=
p
j
(
x
)
(
x
−
x
j
)
k
k
!
∑
t
=
0
r
j
−
k
1
t
!
g
i
(
t
)
(
x
j
)
(
x
−
x
j
)
t
A_{jk}(x)=p_j(x)\dfrac{(x-x_j)^k}{k!}\sum_{t=0}^{r_j-k}\dfrac{1}{t!}g_i^{(t)}(x_j)(x-x_j)^t
Ajk(x)=pj(x)k!(x−xj)kt=0∑rj−kt!1gi(t)(xj)(x−xj)t
容易验证:
A
j
1
=
(
x
−
x
j
)
p
j
(
x
)
p
j
(
x
j
)
=
(
x
−
x
j
)
l
j
2
(
x
)
A_{j1}=(x-x_j)\dfrac{p_j(x)}{p_j(x_j)}=(x-x_j)l_j^2(x)\\
Aj1=(x−xj)pj(xj)pj(x)=(x−xj)lj2(x)
A j 0 = p j ( x ) [ g j ( x j ) + g j ′ ( x j ) ( x − x j ) ] = p j ( x ) [ 1 p j ( x j ) + p j ′ ( x j ) [ p j ( x j ) ] 2 ( x − x j ) ] = [ 1 − p j ′ ( x j ) p j ( x j ) ( x − x j ) ] l j 2 ( x ) = [ 1 − 2 l i ′ ( x i ) ( x − x i ) ] l i 2 ( x ) \begin{align} A_{j0}&=p_j(x)[g_j(x_j)+g'_j(x_j)(x-x_j)]\\ &=p_j(x)\left[\dfrac{1}{p_j(x_j)}+\dfrac{p'_j(x_j)}{[p_j(x_j)]^2}(x-x_j)\right]\\ &=\left[1-\dfrac{p'_j(x_j)}{p_j(x_j)}(x-x_j)\right]l^2_j(x)\\ &=[1-2l_i'(x_i)(x-x_i)]l_i^2(x) \end{align} Aj0=pj(x)[gj(xj)+gj′(xj)(x−xj)]=pj(x)[pj(xj)1+[pj(xj)]2pj′(xj)(x−xj)]=[1−pj(xj)pj′(xj)(x−xj)]lj2(x)=[1−2li′(xi)(x−xi)]li2(x)
这与我们的构造结果是一致的
A h l i n Ahlin Ahlin范式
已知
f
(
p
,
q
)
(
x
r
,
y
i
)
,
p
,
q
=
0
,
⋯
,
k
−
1
;
r
,
i
=
1
,
⋯
,
n
f^{(p,q)}(x_r,y_i),p,q=0,\cdots,k-1;r,i=1,\cdots,n
f(p,q)(xr,yi),p,q=0,⋯,k−1;r,i=1,⋯,n,假设插值函数唯一,记为
p
(
x
,
y
)
p(x,y)
p(x,y),则
∂
p
+
q
∂
x
p
∂
y
q
p
(
x
r
,
y
i
)
=
f
(
p
,
q
)
(
x
r
,
y
i
)
\dfrac{\partial^{p+q}}{\partial x^p\partial y^q}p(x_r,y_i)=f^{(p,q)}(x_r,y_i)
∂xp∂yq∂p+qp(xr,yi)=f(p,q)(xr,yi)
具体的推导过程比较复杂,可以参照原文(因为我没完全看懂),作者先是在复平面上考虑插值余项,进而构造插值函数,最后考虑
x
,
y
x,y
x,y一阶偏导和混合偏导的特殊情形,从一个非常复杂的式子得到一个漂亮的结果:
f
(
x
,
y
)
≈
p
(
x
,
y
)
=
∑
i
=
1
n
∑
r
=
1
n
h
r
(
x
)
g
i
(
y
)
f
(
x
r
,
y
i
)
+
∑
i
=
1
n
∑
r
=
1
n
h
r
(
x
)
g
~
i
(
y
)
∂
∂
y
f
(
x
r
,
y
i
)
+
∑
i
=
1
n
∑
r
=
1
n
h
~
r
(
x
)
g
i
(
y
)
∂
∂
x
f
(
x
r
,
y
i
)
+
∑
i
=
1
n
∑
r
=
1
n
h
~
r
(
x
)
g
~
i
(
y
)
∂
2
∂
x
∂
y
f
(
x
r
,
y
i
)
\begin{align} f(x,y)\approx p(x,y)=&\sum_{i=1}^n\sum_{r=1}^n h_r(x)g_i(y)f(x_r,y_i)+\sum_{i=1}^n\sum_{r=1}^n h_r(x)\widetilde{g}_i(y)\dfrac{\partial}{\partial y}f(x_r,y_i)\\ &+\sum_{i=1}^n\sum_{r=1}^n \widetilde{h}_r(x)g_i(y)\dfrac{\partial}{\partial x}f(x_r,y_i)+\sum_{i=1}^n\sum_{r=1}^n \widetilde{h}_r(x)\widetilde{g}_i(y)\dfrac{\partial^2}{\partial x\partial y}f(x_r,y_i) \end{align}
f(x,y)≈p(x,y)=i=1∑nr=1∑nhr(x)gi(y)f(xr,yi)+i=1∑nr=1∑nhr(x)g
i(y)∂y∂f(xr,yi)+i=1∑nr=1∑nh
r(x)gi(y)∂x∂f(xr,yi)+i=1∑nr=1∑nh
r(x)g
i(y)∂x∂y∂2f(xr,yi)
其中,
h
i
(
x
)
=
[
1
−
2
l
r
′
(
x
r
)
(
x
−
x
r
)
]
l
r
2
(
x
)
h
~
i
(
x
)
=
(
x
−
x
r
)
l
r
2
(
x
)
g
i
(
y
)
=
[
1
−
2
m
r
′
(
y
i
)
(
y
−
y
i
)
]
m
i
2
(
y
)
h
~
i
(
y
)
=
(
y
−
y
i
)
m
i
2
(
y
)
\begin{align} &h_i(x)=[1-2l_r'(x_r)(x-x_r)]l_r^2(x)\\ &\widetilde{h}_i(x)=(x-x_r)l_r^2(x)\\ &g_i(y)=[1-2m_r'(y_i)(y-y_i)]m_i^2(y)\\ &\widetilde{h}_i(y)=(y-y_i)m_i^2(y)\\ \end{align}
hi(x)=[1−2lr′(xr)(x−xr)]lr2(x)h
i(x)=(x−xr)lr2(x)gi(y)=[1−2mr′(yi)(y−yi)]mi2(y)h
i(y)=(y−yi)mi2(y)
l
r
(
x
)
,
m
i
(
y
)
l_r(x),m_i(y)
lr(x),mi(y)分别为关于
x
,
y
x,y
x,y的
L
a
g
r
a
n
g
e
Lagrange
Lagrange插值基函数
上式与一元的情形是相似的,且容易验证,满足 ∂ p + q ∂ x p ∂ y q p ( x r , y i ) = f ( p , q ) ( x r , y i ) , p , q = 0 , 1 \dfrac{\partial^{p+q}}{\partial x^p\partial y^q}p(x_r,y_i)=f^{(p,q)}(x_r,y_i),p,q=0,1 ∂xp∂yq∂p+qp(xr,yi)=f(p,q)(xr,yi),p,q=0,1
C h a w l a − J a y a r a j a n Chawla-Jayarajan Chawla−Jayarajan范式
C h a w l a − J a y a r a j a n Chawla-Jayarajan Chawla−Jayarajan范式的推导过程与一元情形类似,对于 A h l i n Ahlin Ahlin范式改进的地方在于,考虑了对于 x , y x,y x,y已知条件不对称的情形
已知
f
(
k
,
l
)
(
x
i
,
y
j
)
,
i
=
0
,
⋯
,
m
,
j
=
0
,
⋯
,
n
;
k
=
1
,
⋯
,
r
i
,
l
=
1
,
⋯
,
r
j
f^{(k,l)}(x_i,y_j),i=0,\cdots,m,j=0,\cdots,n;k=1,\cdots,r_i,l=1,\cdots,r_j
f(k,l)(xi,yj),i=0,⋯,m,j=0,⋯,n;k=1,⋯,ri,l=1,⋯,rj,假设插值函数唯一,记为
H
M
,
N
(
x
,
y
)
H_{M,N}(x,y)
HM,N(x,y),其中
M
=
m
+
∑
i
=
0
m
r
i
,
N
=
n
+
∑
j
=
0
n
r
j
M=m+\sum\limits_{i=0}^m r_i,N=n+\sum\limits_{j=0}^n r_j
M=m+i=0∑mri,N=n+j=0∑nrj,构造
H
m
,
n
(
x
,
y
)
=
∑
i
=
0
m
∑
j
=
0
n
∑
k
=
0
r
i
∑
l
=
0
s
j
A
i
k
(
x
)
B
j
l
(
y
)
f
(
k
,
l
)
(
x
i
,
y
j
)
H_{m,n}(x,y)=\sum_{i=0}^m\sum_{j=0}^n\sum_{k=0}^{r_i}\sum_{l=0}^{s_j}A_{ik}(x)B_{jl}(y)f^{(k,l)}(x_i,y_j)
Hm,n(x,y)=i=0∑mj=0∑nk=0∑ril=0∑sjAik(x)Bjl(y)f(k,l)(xi,yj)
其中
A
i
k
(
x
)
,
B
j
l
(
y
)
A_{ik}(x),B_{jl}(y)
Aik(x),Bjl(y)作者参照了
S
p
i
t
z
b
a
r
t
Spitzbart
Spitzbart的研究,这里不再赘述
对于已知两对
x
,
y
x,y
x,y的函数值、偏导和混合偏导的特殊情形,即
∂
k
+
l
∂
x
k
∂
y
l
H
2
m
−
1
,
2
n
−
1
(
x
i
,
y
j
)
=
f
(
k
,
l
)
(
x
i
,
y
j
)
,
i
,
j
=
0
,
1
;
k
=
0
,
1
,
⋯
,
m
−
1
,
l
=
0
,
1
,
⋯
,
n
−
1
\dfrac{\partial^{k+l}}{\partial x^k\partial y^l}H_{2m-1,2n-1}(x_i,y_j)=f^{(k,l)}(x_i,y_j),i,j=0,1;k=0,1,\cdots,m-1,l=0,1,\cdots,n-1
∂xk∂yl∂k+lH2m−1,2n−1(xi,yj)=f(k,l)(xi,yj),i,j=0,1;k=0,1,⋯,m−1,l=0,1,⋯,n−1,作者给出了以下公式:
H
2
m
−
1
,
2
n
−
1
=
(
x
−
x
0
)
m
(
y
−
y
0
)
n
∑
k
=
0
m
−
1
∑
l
=
0
n
−
1
A
k
l
(
x
−
x
1
)
k
k
!
(
y
−
y
1
)
l
l
!
+
(
x
−
x
0
)
m
(
y
−
y
1
)
n
∑
k
=
0
m
−
1
∑
l
=
0
n
−
1
B
k
l
(
x
−
x
1
)
k
k
!
(
y
−
y
0
)
l
l
!
+
(
x
−
x
1
)
m
(
y
−
y
0
)
n
∑
k
=
0
m
−
1
∑
l
=
0
n
−
1
C
k
l
(
x
−
x
0
)
k
k
!
(
y
−
y
1
)
l
l
!
+
(
x
−
x
1
)
m
(
y
−
y
1
)
n
∑
k
=
0
m
−
1
∑
l
=
0
n
−
1
D
k
l
(
x
−
x
0
)
k
k
!
(
y
−
y
0
)
l
l
!
\begin{align} H_{2m-1,2n-1}&=(x-x_0)^m(y-y_0)^n\sum_{k=0}^{m-1}\sum_{l=0}^{n-1}A_{kl}\dfrac{(x-x_1)^k}{k!}\dfrac{(y-y_1)^l}{l!}\\ &+(x-x_0)^m(y-y_1)^n\sum_{k=0}^{m-1}\sum_{l=0}^{n-1}B_{kl}\dfrac{(x-x_1)^k}{k!}\dfrac{(y-y_0)^l}{l!}\\\\ &+(x-x_1)^m(y-y_0)^n\sum_{k=0}^{m-1}\sum_{l=0}^{n-1}C_{kl}\dfrac{(x-x_0)^k}{k!}\dfrac{(y-y_1)^l}{l!}\\ &+(x-x_1)^m(y-y_1)^n\sum_{k=0}^{m-1}\sum_{l=0}^{n-1}D_{kl}\dfrac{(x-x_0)^k}{k!}\dfrac{(y-y_0)^l}{l!}\\ \end{align}
H2m−1,2n−1=(x−x0)m(y−y0)nk=0∑m−1l=0∑n−1Aklk!(x−x1)kl!(y−y1)l+(x−x0)m(y−y1)nk=0∑m−1l=0∑n−1Bklk!(x−x1)kl!(y−y0)l+(x−x1)m(y−y0)nk=0∑m−1l=0∑n−1Cklk!(x−x0)kl!(y−y1)l+(x−x1)m(y−y1)nk=0∑m−1l=0∑n−1Dklk!(x−x0)kl!(y−y0)l
其中
A
k
l
=
[
∂
k
+
l
∂
x
k
∂
y
l
f
(
x
,
y
)
(
x
−
x
0
)
m
(
y
−
y
0
)
n
]
x
−
x
1
,
y
=
y
1
A_{kl}=\left[\dfrac{\partial^{k+l}}{\partial x^k\partial y^l} \dfrac{f(x,y)}{(x-x_0)^m}(y-y_0)^n\right]_{x-x_1,y=y_1}
Akl=[∂xk∂yl∂k+l(x−x0)mf(x,y)(y−y0)n]x−x1,y=y1
B
k
l
,
C
k
l
,
D
k
l
B_{kl},C_{kl},D_{kl}
Bkl,Ckl,Dkl类似
解题与编程
范式的选择
对于上述两种范式,需要根据需求选择一个。从编程的难易程度而言, A h l i n Ahlin Ahlin范式是优于 C h a w l a − J a y a r a j a n Chawla-Jayarajan Chawla−Jayarajan范式的,因为 A h l i n Ahlin Ahlin范式与一元的 H e r m i t e Hermite Hermite插值公式相似,并且方便定义函数和构造循环体;而对于 C h a w l a − J a y a r a j a n Chawla-Jayarajan Chawla−Jayarajan范式,比较头疼的是 A k l A_{kl} Akl的计算,尽管对于特定的 k , l k,l k,l, A k l , B k l , C k l , D k l A_{kl},B_{kl},C_{kl},D_{kl} Akl,Bkl,Ckl,Dkl可以一起实现,但其中偏导的部分在程序中是难以实现的,对于每阶偏导需要单独手算然后根据手算的结果代入数据,就非常繁琐。而且,当插值函数唯一时,这两种范式的表达式理论上是统一的。
于是我们在解题时选择 A h l i n Ahlin Ahlin范式。
分块的实现
回顾一元分段三次
H
e
r
m
i
t
e
Hermite
Hermite插值,每个子区间
[
x
i
,
x
i
+
1
]
[x_i,x_{i+1}]
[xi,xi+1]上的
p
3
(
x
)
p_3(x)
p3(x)可以写成
p
3
(
x
)
=
(
1
−
2
x
−
x
i
x
i
−
x
i
+
1
)
(
x
−
x
i
+
1
x
i
−
x
i
+
1
)
2
f
(
x
i
)
+
(
1
−
2
x
i
+
1
−
x
x
i
−
x
i
+
1
)
(
x
−
x
i
x
i
−
x
i
+
1
)
2
f
(
x
i
+
1
)
+
(
x
−
x
i
)
(
x
−
x
i
x
i
−
x
i
+
1
)
2
f
′
(
x
i
)
+
(
x
−
x
i
+
1
)
(
x
−
x
i
x
i
+
1
−
x
i
)
2
f
′
(
x
i
+
1
)
\begin{align} p_3(x)=&\left(1-2\dfrac{x-x_i}{x_i-x_{i+1}}\right)\left(\dfrac{x-x_{i+1}}{x_i-x_{i+1}}\right)^2f(x_i) +\left(1-2\dfrac{x_{i+1}-x}{x_i-x_{i+1}}\right)\left(\dfrac{x-x_i}{x_i-x_{i+1}}\right)^2f(x_{i+1})\\ &+(x-x_i)\left(\dfrac{x-x_i}{x_i-x_{i+1}}\right)^2f'(x_i)+ (x-x_{i+1})\left(\dfrac{x-x_i}{x_{i+1}-x_{i}}\right)^2f'(x_{i+1}) \end{align}
p3(x)=(1−2xi−xi+1x−xi)(xi−xi+1x−xi+1)2f(xi)+(1−2xi−xi+1xi+1−x)(xi−xi+1x−xi)2f(xi+1)+(x−xi)(xi−xi+1x−xi)2f′(xi)+(x−xi+1)(xi+1−xix−xi)2f′(xi+1)
于是,类似地,可以构造二元分块
3
×
3
3\times 3
3×3次
H
e
r
m
i
t
e
Hermite
Hermite插值,为方便起见,对于
x
x
x引入插值基函数:
ϕ
0
(
x
)
=
{
(
1
−
2
x
−
x
0
x
0
−
x
1
)
(
x
−
x
1
x
0
−
x
1
)
2
,
x
∈
[
x
0
,
x
1
]
0
,
e
l
s
e
ϕ
i
(
x
)
=
{
(
1
−
2
x
−
x
i
x
i
−
x
i
+
1
)
(
x
−
x
i
+
1
x
i
−
x
i
+
1
)
2
,
x
∈
[
x
i
−
1
,
x
i
)
(
1
−
2
x
i
+
1
−
x
x
i
−
x
i
+
1
)
(
x
−
x
i
x
i
−
x
i
+
1
)
2
,
x
∈
[
x
i
,
x
i
+
1
]
0
,
e
l
s
e
ϕ
n
(
x
)
=
{
(
1
−
2
x
n
−
x
x
n
−
1
−
x
n
)
(
x
−
x
n
−
1
x
n
−
1
−
x
n
)
2
,
x
∈
[
x
n
−
1
,
x
n
]
0
,
e
l
s
e
ϕ
0
~
(
x
)
=
{
(
x
−
x
0
)
(
x
−
x
0
x
i
−
x
1
)
2
,
x
∈
[
x
0
,
x
1
]
0
,
e
l
s
e
ϕ
i
~
(
x
)
=
{
(
x
−
x
i
)
(
x
−
x
i
x
i
−
x
i
+
1
)
2
,
x
∈
[
x
i
−
1
,
x
i
)
(
x
−
x
i
+
1
)
(
x
−
x
i
x
i
+
1
−
x
i
)
2
2
,
x
∈
[
x
i
,
x
i
+
1
]
0
,
e
l
s
e
ϕ
n
~
(
x
)
=
{
(
x
−
x
n
)
(
x
−
x
n
−
1
x
n
−
x
n
−
1
)
2
,
x
∈
[
x
n
−
1
,
x
n
]
0
,
e
l
s
e
\phi_0(x)= \begin{cases} \left(1-2\dfrac{x-x_0}{x_0-x_{1}}\right)\left(\dfrac{x-x_{1}}{x_0-x_{1}}\right)^2,&x\in[x_0,x_1]\\ 0,&else \end{cases}\\ \phi_i(x)= \begin{cases} \left(1-2\dfrac{x-x_i}{x_i-x_{i+1}}\right)\left(\dfrac{x-x_{i+1}}{x_i-x_{i+1}}\right)^2,&x\in[x_{i-1},x_i)\\ \left(1-2\dfrac{x_{i+1}-x}{x_i-x_{i+1}}\right)\left(\dfrac{x-x_i}{x_i-x_{i+1}}\right)^2,&x\in[x_i,x_{i+1}]\\ 0,&else \end{cases}\\ \phi_n(x)= \begin{cases} \left(1-2\dfrac{x_{n}-x}{x_{n-1}-x_{n}}\right)\left(\dfrac{x-x_{n-1}}{x_{n-1}-x_{n}}\right)^2,&x\in[x_{n-1},x_n]\\ 0,&else \end{cases}\\ \widetilde{\phi_0}(x)= \begin{cases} (x-x_0)\left(\dfrac{x-x_0}{x_i-x_{1}}\right)^2,&x\in[x_0,x_1]\\ 0,&else \end{cases}\\ \widetilde{\phi_i}(x)= \begin{cases} (x-x_i)\left(\dfrac{x-x_i}{x_i-x_{i+1}}\right)^2,&x\in[x_{i-1},x_i)\\ (x-x_{i+1})\left(\dfrac{x-x_i}{x_{i+1}-x_{i}}\right)^22,&x\in[x_i,x_{i+1}]\\ 0,&else \end{cases}\\ \widetilde{\phi_n}(x)= \begin{cases} (x-x_{n})\left(\dfrac{x-x_{n-1}}{x_n-x_{n-1}}\right)^2,&x\in[x_{n-1},x_n]\\ 0,&else \end{cases}\\
ϕ0(x)=⎩
⎨
⎧(1−2x0−x1x−x0)(x0−x1x−x1)2,0,x∈[x0,x1]elseϕi(x)=⎩
⎨
⎧(1−2xi−xi+1x−xi)(xi−xi+1x−xi+1)2,(1−2xi−xi+1xi+1−x)(xi−xi+1x−xi)2,0,x∈[xi−1,xi)x∈[xi,xi+1]elseϕn(x)=⎩
⎨
⎧(1−2xn−1−xnxn−x)(xn−1−xnx−xn−1)2,0,x∈[xn−1,xn]elseϕ0
(x)=⎩
⎨
⎧(x−x0)(xi−x1x−x0)2,0,x∈[x0,x1]elseϕi
(x)=⎩
⎨
⎧(x−xi)(xi−xi+1x−xi)2,(x−xi+1)(xi+1−xix−xi)22,0,x∈[xi−1,xi)x∈[xi,xi+1]elseϕn
(x)=⎩
⎨
⎧(x−xn)(xn−xn−1x−xn−1)2,0,x∈[xn−1,xn]else
对于
y
y
y:
ψ
0
(
y
)
=
{
(
1
−
2
y
−
y
0
y
0
−
y
1
)
(
y
−
y
1
y
0
−
y
1
)
2
,
y
∈
[
y
0
,
y
1
]
0
,
e
l
s
e
ψ
i
(
y
)
=
{
(
1
−
2
y
−
y
i
y
i
−
y
i
+
1
)
(
y
−
y
i
+
1
y
i
−
y
i
+
1
)
2
,
y
∈
[
y
i
−
1
,
y
i
)
(
1
−
2
y
i
+
1
−
y
y
i
−
y
i
+
1
)
(
y
−
y
i
y
i
−
y
i
+
1
)
2
,
y
∈
[
y
i
,
y
i
+
1
]
0
,
e
l
s
e
ψ
n
(
y
)
=
{
(
1
−
2
y
n
−
y
y
n
−
1
−
y
n
)
(
y
−
y
n
−
1
y
n
−
1
−
y
n
)
2
,
y
∈
[
y
n
−
1
,
y
n
]
0
,
e
l
s
e
ψ
0
~
(
y
)
=
{
(
y
−
y
0
)
(
y
−
y
0
y
i
−
y
1
)
2
,
y
∈
[
y
0
,
y
1
]
0
,
e
l
s
e
ψ
i
~
(
y
)
=
{
(
y
−
y
i
)
(
y
−
y
i
y
i
−
y
i
+
1
)
2
,
y
∈
[
y
i
−
1
,
y
i
)
(
y
−
y
i
+
1
)
(
y
−
y
i
y
i
+
1
−
y
i
)
2
2
,
y
∈
[
y
i
,
y
i
+
1
]
0
,
e
l
s
e
ψ
n
~
(
y
)
=
{
(
y
−
y
n
)
(
y
−
y
n
−
1
y
n
−
y
n
−
1
)
2
,
y
∈
[
y
n
−
1
,
y
n
]
0
,
e
l
s
e
\psi_0(y)= \begin{cases} \left(1-2\dfrac{y-y_0}{y_0-y_{1}}\right)\left(\dfrac{y-y_{1}}{y_0-y_{1}}\right)^2,&y\in[y_0,y_1]\\ 0,&else \end{cases}\\ \psi_i(y)= \begin{cases} \left(1-2\dfrac{y-y_i}{y_i-y_{i+1}}\right)\left(\dfrac{y-y_{i+1}}{y_i-y_{i+1}}\right)^2,&y\in[y_{i-1},y_i)\\ \left(1-2\dfrac{y_{i+1}-y}{y_i-y_{i+1}}\right)\left(\dfrac{y-y_i}{y_i-y_{i+1}}\right)^2,&y\in[y_i,y_{i+1}]\\ 0,&else \end{cases}\\ \psi_n(y)= \begin{cases} \left(1-2\dfrac{y_{n}-y}{y_{n-1}-y_{n}}\right)\left(\dfrac{y-y_{n-1}}{y_{n-1}-y_{n}}\right)^2,&y\in[y_{n-1},y_n]\\ 0,&else \end{cases}\\ \widetilde{\psi_0}(y)= \begin{cases} (y-y_0)\left(\dfrac{y-y_0}{y_i-y_{1}}\right)^2,&y\in[y_0,y_1]\\ 0,&else \end{cases}\\ \widetilde{\psi_i}(y)= \begin{cases} (y-y_i)\left(\dfrac{y-y_i}{y_i-y_{i+1}}\right)^2,&y\in[y_{i-1},y_i)\\ (y-y_{i+1})\left(\dfrac{y-y_i}{y_{i+1}-y_{i}}\right)^22,&y\in[y_i,y_{i+1}]\\ 0,&else \end{cases}\\ \widetilde{\psi_n}(y)= \begin{cases} (y-y_{n})\left(\dfrac{y-y_{n-1}}{y_n-y_{n-1}}\right)^2,&y\in[y_{n-1},y_n]\\ 0,&else \end{cases}\\
ψ0(y)=⎩
⎨
⎧(1−2y0−y1y−y0)(y0−y1y−y1)2,0,y∈[y0,y1]elseψi(y)=⎩
⎨
⎧(1−2yi−yi+1y−yi)(yi−yi+1y−yi+1)2,(1−2yi−yi+1yi+1−y)(yi−yi+1y−yi)2,0,y∈[yi−1,yi)y∈[yi,yi+1]elseψn(y)=⎩
⎨
⎧(1−2yn−1−ynyn−y)(yn−1−yny−yn−1)2,0,y∈[yn−1,yn]elseψ0
(y)=⎩
⎨
⎧(y−y0)(yi−y1y−y0)2,0,y∈[y0,y1]elseψi
(y)=⎩
⎨
⎧(y−yi)(yi−yi+1y−yi)2,(y−yi+1)(yi+1−yiy−yi)22,0,y∈[yi−1,yi)y∈[yi,yi+1]elseψn
(y)=⎩
⎨
⎧(y−yn)(yn−yn−1y−yn−1)2,0,y∈[yn−1,yn]else
于是,每个子块
[
x
i
,
x
i
+
1
]
×
[
y
j
,
y
j
+
1
]
[x_i,x_{i+1}]\times [y_j,y_{j+1}]
[xi,xi+1]×[yj,yj+1]上的
p
3
,
3
(
x
,
y
)
p_{3,3}(x,y)
p3,3(x,y)可以写成
p
3.3
(
x
,
y
)
=
∑
j
=
1
n
∑
i
=
1
n
ϕ
i
(
x
)
g
j
(
y
)
f
(
x
i
,
y
j
)
+
∑
j
=
1
n
∑
i
=
1
n
ϕ
i
(
x
)
ψ
~
j
(
y
)
∂
∂
y
f
(
x
i
,
y
j
)
+
∑
j
=
1
n
∑
i
=
1
n
ϕ
~
i
(
x
)
g
j
(
y
)
∂
∂
x
f
(
x
i
,
y
j
)
+
∑
j
=
1
n
∑
i
=
1
n
ϕ
~
i
(
x
)
ψ
~
j
(
y
)
∂
2
∂
x
∂
y
f
(
x
i
,
y
j
)
⋯
(
∗
)
\begin{align} p_{3.3}(x,y)=&\sum_{j=1}^n\sum_{i=1}^n \phi_i(x)g_j(y)f(x_i,y_j)+\sum_{j=1}^n\sum_{i=1}^n \phi_i(x)\widetilde{\psi}_j(y)\dfrac{\partial}{\partial y}f(x_i,y_j)\\ &+\sum_{j=1}^n\sum_{i=1}^n \widetilde{\phi}_i(x)g_j(y)\dfrac{\partial}{\partial x}f(x_i,y_j)+\sum_{j=1}^n\sum_{i=1}^n \widetilde{\phi}_i(x)\widetilde{\psi}_j(y)\dfrac{\partial^2}{\partial x\partial y}f(x_i,y_j) \end{align} \cdots (*)
p3.3(x,y)=j=1∑ni=1∑nϕi(x)gj(y)f(xi,yj)+j=1∑ni=1∑nϕi(x)ψ
j(y)∂y∂f(xi,yj)+j=1∑ni=1∑nϕ
i(x)gj(y)∂x∂f(xi,yj)+j=1∑ni=1∑nϕ
i(x)ψ
j(y)∂x∂y∂2f(xi,yj)⋯(∗)
问题1的解决
根据上述分析,算法如下:
用$C/C++$进行编程实现,顺便对于区间的选择和分块插值次数做了泛化,代码如下:
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
void data(double *x, double *y, double *f, double *dx, double *dy, double *dxy, int n) {
// 生成sin(x^2y+1)的dx,dy,dxy
int i, j;
for (i = 0; i <= n; i++) {
for (j = 0; j <= n; j++) {
f[i * n + j] = sin(x[i] * x[i] * y[j] + 1);
dx[i * n + j] = 2 * x[i] * y[j] * cos(x[i] * x[i] * y[j] + 1);
dy[i * n + j] = x[i] * x[i] * cos(x[i] * x[i] * y[j] + 1);
dxy[i * n + j] = 2 * x[i] * cos(x[i] * x[i] * y[j] + 1) - 2 * pow(x[i], 3) * y[i] * sin(x[i] * x[i] * y[j] + 1);
}
}
}
double lagrange(double *z, double zs, int n, int i) {
// 求lagrange基函数
int j;
double t = 1;
for (j = 0; j < n; j++) {
if (i != j) {
t *= ((zs - z[j]) / (z[i] - z[j]));
}
}
return t;
}
double lagrange_(double *z, int i, int n) {
// 求lagrange基函数导数
int j;
double r = 0;
for (j = 0; j < n; j++) {
if (i != j) {
r += (1 / (z[i] - z[j]));
}
}
return r;
}
double Hermite(double *x, double *y, double *f, double *dx, double *dy, double *dxy, int m, double xs, double ys, int s) {
int i, j, k, t;
double r = 0, h, h_, g, g_;
double *x_, *y_;
x_ = (double *) malloc(sizeof(double) * (s - 1));
y_ = (double *) malloc(sizeof(double) * (s - 1));
for (i = 0; i < m; i++) {
for (j = 0; j < m; j++) {
if (x[i] <= xs && xs <= x[i + s - 2] && y[j] <= ys && ys <= y[j + s - 2] && i % (s - 2) == 0 && j % (s - 2) == 0 && i + s - 2 <= m && j + s - 2 <= m) {
for (k = 0; k < s - 1; k++) {
x_[k] = x[i + k];
y_[k] = y[j + k];
}
for (k = 0; k <= s - 2; k++) {
for (t = 0; t <= s - 2; t++) {
h = (1 - 2 * lagrange_(x_, k, s - 1) * (xs - x_[k])) * lagrange(x_, xs, s - 1, k) * lagrange(x_, xs, s - 1, k);
h_ = (xs - x_[k]) * lagrange(x_, xs, s - 1, k) * lagrange(x_, xs, s - 1, k);
g = (1 - 2 * lagrange_(y_, t, s - 1) * (ys - y_[t])) * lagrange(y_, ys, s - 1, t) * lagrange(y_, ys, s - 1, t);
g_ = (ys - y_[t]) * lagrange(y_, ys, s - 1, t) * lagrange(y_, ys, s - 1, t);
r += (h * g * f[(i + k) * m + (j + t)] + h * g_ * dy[(i + k) * m + (j + t)] + h_ * g * dx[(i + k) * m + (j + t)] + h_ * g_ * dxy[(i + k) * m + (j + t)]);
}
}
}
}
}
free(x_);
free(y_);
return r;
}
int main() {
int n = 5; // 区间n等分
int s = 3; // s次插值
double *x, *y, *dx, *dy, *dxy, *f;
double xs, ys;
double left = 0, right = 1;
int i, j;
double value, result;
xs = 1.0 / 3.0;
ys = 2.0 / 3.0;
x = (double *) malloc(sizeof(double) * (n + 1));
dx = (double *) malloc(sizeof(double) * (n + 1) * (n + 1));
y = (double *) malloc(sizeof(double) * (n + 1));
dy = (double *) malloc(sizeof(double) * (n + 1) * (n + 1));
dxy = (double *) malloc(sizeof(double) * (n + 1) * (n + 1));
f = (double *) malloc(sizeof(double) * (n + 1) * (n + 1));
// 等分区间
for (i = 0; i <= n; i++) {
x[i] = ((right - left) / n) * i + left;
y[i] = ((right - left) / n) * i + left;
}
data(x, y, f, dx, dy, dxy, n);
value = sin(xs * xs * ys + 1);
printf("Actual Value: %.15lf\n", value);
result = Hermite(x, y, f, dx, dy, dxy, n, xs, ys, s);
printf("Hermite Value: %.15lf\n", result);
printf("error: %.15lf", abs(result - value));
free(x);
free(dx);
free(y);
free(dy);
free(f);
free(dxy);
return 0;
}
结果为:
问题2的解决
计算插值多项式在特定点的值
如果要满足题目所需条件,即满足
p
(
x
,
y
)
=
f
(
x
,
y
)
,
∂
∂
x
p
(
x
,
y
)
=
∂
∂
x
f
(
x
,
y
)
,
∂
∂
y
p
(
x
,
y
)
=
∂
∂
y
f
(
x
,
y
)
p(x,y)=f(x,y),\dfrac{\partial}{\partial x}p(x,y)=\dfrac{\partial}{\partial x}f(x,y),\dfrac{\partial}{\partial y}p(x,y)=\dfrac{\partial}{\partial y}f(x,y)
p(x,y)=f(x,y),∂x∂p(x,y)=∂x∂f(x,y),∂y∂p(x,y)=∂y∂f(x,y),显然,
(
∗
)
(*)
(∗)式是满足条件的,因此只需将上述代码中Hermite
函数中的
r += (h * g * f[(i + k) * m + (j + t)] + h * g_ * dy[(i + k) * m + (j + t)] + h_ * g * dx[(i + k) * m + (j + t)] + h_ * g_ * dxy[(i + k) * m + (j + t)]);
改成
r += (h * g * f[(i + k) * m + (j + t)] + h * g_ * dy[(i + k) * m + (j + t)] + h_ * g * dx[(i + k) * m + (j + t)]);
即可
结果为:
连续性和可微性的讨论
对于连续性,在 ( x i , x i + 1 ) , i = 0 , 1 , ⋯ , n − 1 (x_i,x_{i+1}),i=0,1,\cdots,n-1 (xi,xi+1),i=0,1,⋯,n−1,由于 p 3 ( x ) p_3(x) p3(x)为多项式函数,显然连续;而在点 x = x i , i = 0 , 1 , ⋯ , n x=x_i,i=0,1,\cdots,n x=xi,i=0,1,⋯,n处, p 3 ( x ) = 0 p_3(x)=0 p3(x)=0,显然连续。
综上, p 3 ( x ) p_3(x) p3(x)在 [ x 0 , x n ] [x_0,x_n] [x0,xn]上连续。
对于可微性,由于 ∂ 2 ∂ x ∂ y f ( x , y ) \dfrac{\partial^2}{\partial x\partial y}f(x,y) ∂x∂y∂2f(x,y)没有数据(实际上是按0处理了),则 ∀ x ∈ [ x 0 , x n ] \forall x\in[x_0,x_n] ∀x∈[x0,xn], ∂ 2 ∂ x ∂ y p 3 ( x , y ) ≡ 0 \dfrac{\partial^2}{\partial x\partial y}p_3(x,y)\equiv 0 ∂x∂y∂2p3(x,y)≡0;由于 p 3 ( x ) p_3(x) p3(x)为多项式函数,显然其关于 x , y x,y x,y的偏导数存在,故关于 x , y x,y x,y可微。
综上, p 3 ( x ) p_3(x) p3(x)在 [ x 0 , x n ] [x_0,x_n] [x0,xn]上可微。
对于问题2的疑惑
从上述对于问题2的讨论中可以发现,由于
∂
2
∂
x
∂
y
f
(
x
,
y
)
\dfrac{\partial^2}{\partial x\partial y}f(x,y)
∂x∂y∂2f(x,y)没有数据,我们的处理是直接将
∂
2
∂
x
∂
y
f
(
x
,
y
)
\dfrac{\partial^2}{\partial x\partial y}f(x,y)
∂x∂y∂2f(x,y)设置成0,这就导致了我们的误差在问题1结果的基础上又加上了
∑
j
=
1
n
∑
i
=
1
n
ϕ
~
i
(
x
)
ψ
~
j
(
y
)
∂
2
∂
x
∂
y
f
(
x
i
,
y
j
)
\sum_{j=1}^n\sum_{i=1}^n \widetilde{\phi}_i(x)\widetilde{\psi}_j(y)\dfrac{\partial^2}{\partial x\partial y}f(x_i,y_j)
j=1∑ni=1∑nϕ
i(x)ψ
j(y)∂x∂y∂2f(xi,yj)
这直接导致误差大了十倍甚至九倍,是否合理呢?
首先,我们的已知条件由16个缩减到了12个,于是我们只能构造12个方程,因此我们无法构造 p ( x , y ) = ∑ i = 0 3 ∑ j = 0 3 c i j x i y j p(x,y)=\sum\limits_{i=0}^3\sum\limits_{j=0}^3c_{ij}x^iy^j p(x,y)=i=0∑3j=0∑3cijxiyj,因为有16个未知数。我们只能构造 p ( x ) = ∑ i = 0 2 ∑ j = 0 2 c i j x i y j + q ( x , y ) p(x)=\sum\limits_{i=0}^2\sum\limits_{j=0}^2c_{ij}x^iy^j+q(x,y) p(x)=i=0∑2j=0∑2cijxiyj+q(x,y),其中 q ( x , y ) q(x,y) q(x,y)包含 x k y 3 − k , k = 0 , 1 , 2 , 3 x^ky^{3-k},k=0,1,2,3 xky3−k,k=0,1,2,3中的三项并乘一个常数,因此插值函数是不唯一的
而在Tomas Sauer和Yuan Xu的研究中[7],对于多元 H e r m i t e Hermite Hermite插值的不同插值条件,即已知不同方向导数的情形,做了进一步的讨论,其引入了适定性(posedness)的概念(事实上,现代对于插值的大部分研究都会引入并考虑这一概念),对于类似问题2的情形做了进一步的阐释
总之,由于数据的缺失,误差是无法避免的
误差分析
回顾一元分段三次
H
e
r
m
i
t
e
Hermite
Hermite插值的误差分析,其余项为:
R
(
x
)
=
f
(
4
)
(
ξ
i
)
4
!
(
x
−
x
i
)
2
(
x
−
x
i
+
1
)
2
,
x
i
<
ξ
i
<
x
i
+
1
R(x)=\dfrac{f^{(4)}(\xi_i)}{4!}(x-x_i)^2(x-x_{i+1})^2,x_i<\xi_i<x_{i+1}
R(x)=4!f(4)(ξi)(x−xi)2(x−xi+1)2,xi<ξi<xi+1
则其最大误差限为:
∣
R
(
x
)
∣
=
∣
f
(
4
)
(
ξ
i
)
4
!
(
x
−
x
i
)
2
(
x
−
x
i
+
1
)
2
∣
≤
h
4
384
max
x
i
<
x
<
x
i
+
1
∣
f
(
4
)
(
x
)
∣
,
x
i
<
ξ
i
<
x
i
+
1
|R(x)|=\left|\dfrac{f^{(4)}(\xi_i)}{4!}(x-x_i)^2(x-x_{i+1})^2\right|\le \dfrac{h^4}{384}\max\limits_{x_i<x<x_{i+1}}|f^{(4)}(x)|,x_i<\xi_i<x_{i+1}
∣R(x)∣=∣
∣4!f(4)(ξi)(x−xi)2(x−xi+1)2∣
∣≤384h4xi<x<xi+1max∣f(4)(x)∣,xi<ξi<xi+1
对于
A
h
l
i
n
Ahlin
Ahlin范式和
C
h
a
w
l
a
−
J
a
y
a
r
a
j
a
n
Chawla-Jayarajan
Chawla−Jayarajan范式,两个范式的作者分别给出了各自的误差分析
A
h
l
i
n
Ahlin
Ahlin是在复平面中分析的,考虑包含
x
,
x
0
,
x
1
,
⋯
,
x
n
−
1
x,x_0,x_1,\cdots,x_{n-1}
x,x0,x1,⋯,xn−1的
J
o
r
d
a
n
Jordan
Jordan曲线
C
1
C_1
C1和
y
,
y
0
,
y
1
,
⋯
,
y
n
−
1
y,y_0,y_1,\cdots,y_{n-1}
y,y0,y1,⋯,yn−1的
J
o
r
d
a
n
Jordan
Jordan曲线,且
k
−
1
k-1
k−1为已知的最高次偏导,则插值余项为:
f
(
x
,
y
)
−
p
(
x
,
y
)
=
[
λ
(
x
)
]
k
2
π
i
∫
C
1
f
(
s
,
y
)
[
λ
(
s
)
]
k
(
s
−
x
)
d
s
+
[
μ
(
y
)
]
k
2
π
i
∫
C
2
f
(
x
,
t
)
[
μ
(
s
)
]
k
(
t
−
y
)
d
t
−
[
λ
(
x
)
μ
(
y
)
]
k
2
π
i
∫
C
1
∫
C
2
f
(
s
,
t
)
[
λ
(
s
)
μ
(
t
)
]
k
(
s
−
x
)
(
t
−
y
)
d
s
d
t
\begin{align} f(x,y)-p(x,y)=&\dfrac{[\lambda(x)]^k}{2\pi i}\int_{C_1}\dfrac{f(s,y)}{[\lambda(s)]^k(s-x)}ds +\dfrac{[\mu(y)]^k}{2\pi i}\int_{C_2}\dfrac{f(x,t)}{[\mu(s)]^k(t-y)}dt\\ &-\dfrac{[\lambda(x)\mu(y)]^k}{2\pi i}\int_{C_1}\int_{C_2}\dfrac{f(s,t)}{[\lambda(s)\mu(t)]^k(s-x)(t-y)}dsdt \end{align}
f(x,y)−p(x,y)=2πi[λ(x)]k∫C1[λ(s)]k(s−x)f(s,y)ds+2πi[μ(y)]k∫C2[μ(s)]k(t−y)f(x,t)dt−2πi[λ(x)μ(y)]k∫C1∫C2[λ(s)μ(t)]k(s−x)(t−y)f(s,t)dsdt
其中
λ
(
x
)
=
∏
i
=
0
n
−
1
(
x
−
x
i
)
μ
(
y
)
=
∏
i
=
0
n
−
1
(
y
−
y
i
)
\lambda(x)=\prod_{i=0}^{n-1}(x-x_i)\\ \mu(y)=\prod_{i=0}^{n-1}(y-y_i)
λ(x)=i=0∏n−1(x−xi)μ(y)=i=0∏n−1(y−yi)
C
h
a
w
l
a
−
J
a
y
a
r
a
j
a
n
Chawla-Jayarajan
Chawla−Jayarajan范式的余项为:
R
M
,
N
(
x
,
y
)
=
α
(
x
)
(
M
+
1
)
!
∂
M
+
1
∂
x
M
+
1
f
(
ξ
,
y
)
+
β
(
y
)
(
N
+
1
)
!
∂
N
+
1
∂
y
N
+
1
f
(
x
,
η
)
−
α
(
x
)
β
(
y
)
(
M
+
1
)
!
(
N
+
1
)
!
∂
M
+
N
+
2
∂
x
M
+
1
∂
y
N
+
1
f
(
ξ
,
η
)
,
x
i
<
ξ
<
x
i
+
1
,
y
j
<
η
<
y
j
+
1
\begin{align} R_{M,N}(x,y)=&\dfrac{\alpha(x)}{(M+1)!}\dfrac{\partial^{M+1}}{\partial x^{M+1}}f(\xi,y)+ \dfrac{\beta(y)}{(N+1)!}\dfrac{\partial^{N+1}}{\partial y^{N+1}}f(x,\eta)\\ &-\dfrac{\alpha(x)\beta(y)}{(M+1)!(N+1)!}\dfrac{\partial^{M+N+2}}{\partial x^{M+1}\partial y^{N+1}}f(\xi,\eta),\\ &x_i<\xi<x_{i+1},y_j<\eta<y_{j+1} \end{align}
RM,N(x,y)=(M+1)!α(x)∂xM+1∂M+1f(ξ,y)+(N+1)!β(y)∂yN+1∂N+1f(x,η)−(M+1)!(N+1)!α(x)β(y)∂xM+1∂yN+1∂M+N+2f(ξ,η),xi<ξ<xi+1,yj<η<yj+1
其中
α
(
x
)
=
(
x
−
x
0
)
r
0
+
1
⋯
(
x
−
x
m
)
r
m
+
1
β
(
y
)
=
(
y
−
y
0
)
s
0
+
1
⋯
(
y
−
y
m
)
s
n
+
1
\alpha(x)=(x-x_0)^{r_0+1}\cdots(x-x_m)^{r_m+1}\\ \beta(y)=(y-y_0)^{s_0+1}\cdots(y-y_m)^{s_n+1}
α(x)=(x−x0)r0+1⋯(x−xm)rm+1β(y)=(y−y0)s0+1⋯(y−ym)sn+1
由于 A h l i n Ahlin Ahlin范式的误差分析是在复平面上分析,需要用到留数定理进行误差分析;而 C h a w l a − J a y a r a j a n Chawla-Jayarajan Chawla−Jayarajan范式的余项相对简洁的多,而且和一元 H e r m i t e Hermite Hermite插值的余项更加相似。
我们只考虑插值公式唯一的情形(问题1),对于不唯一的情形(问题2)可以在插值公式唯一的情形加上或减去缺失的项。
由于插值公式唯一,所以插值余项理论上也是唯一的,于是我们选择更为简洁的 C h a w l a − J a y a r a j a n Chawla-Jayarajan Chawla−Jayarajan范式的余项进行误差分析。
则问题1的余项为:
R
3
,
3
(
x
,
y
)
=
(
x
−
0.2
)
2
(
x
−
0.4
)
2
4
!
∂
4
∂
x
4
f
(
ξ
,
y
)
+
(
y
−
0.6
)
2
(
y
−
0.8
)
2
4
!
∂
4
∂
y
4
f
(
x
,
η
)
−
(
x
−
0.2
)
2
(
x
−
0.4
)
2
(
y
−
0.6
)
2
(
y
−
0.8
)
2
4
!
4
!
∂
8
∂
x
4
∂
y
4
f
(
ξ
,
η
)
\begin{align} R_{3,3}(x,y)=&\dfrac{(x-0.2)^2(x-0.4)^2}{4!}\dfrac{\partial^{4}}{\partial x^{4}}f(\xi,y)+ \dfrac{(y-0.6)^2(y-0.8)^2}{4!}\dfrac{\partial^4}{\partial y^{4}}f(x,\eta)\\ &-\dfrac{(x-0.2)^2(x-0.4)^2(y-0.6)^2(y-0.8)^2}{4!4!}\dfrac{\partial^{8}}{\partial x^{4}\partial y^{4}}f(\xi,\eta) \end{align}
R3,3(x,y)=4!(x−0.2)2(x−0.4)2∂x4∂4f(ξ,y)+4!(y−0.6)2(y−0.8)2∂y4∂4f(x,η)−4!4!(x−0.2)2(x−0.4)2(y−0.6)2(y−0.8)2∂x4∂y4∂8f(ξ,η)
则最大误差限为:
∣
R
3
,
3
(
x
,
y
)
∣
≤
∣
0.
2
4
384
∂
4
∂
x
4
f
(
ξ
,
y
)
∣
+
∣
0.
2
4
384
∂
4
∂
x
4
f
(
x
,
η
)
∣
+
∣
0.
2
8
38
4
2
∂
8
∂
x
4
∂
y
4
f
(
ξ
,
η
)
∣
\begin{align} |R_{3,3}(x,y)|&\le \left|\dfrac{0.2^4}{384}\dfrac{\partial^{4}}{\partial x^{4}}f(\xi,y) \right|+\left|\dfrac{0.2^4}{384}\dfrac{\partial^{4}}{\partial x^{4}}f(x,\eta) \right|+\left|\dfrac{0.2^8}{384^2}\dfrac{\partial^{8}}{\partial x^{4}\partial y^{4}}f(\xi,\eta) \right|\\ \end{align}
∣R3,3(x,y)∣≤∣
∣3840.24∂x4∂4f(ξ,y)∣
∣+∣
∣3840.24∂x4∂4f(x,η)∣
∣+∣
∣38420.28∂x4∂y4∂8f(ξ,η)∣
∣
对于
f
(
x
,
y
)
=
s
i
n
(
x
2
y
+
1
)
,
0.2
<
x
<
0.4
,
0.6
<
y
<
0.8
f(x,y)=sin(x^2y+1),0.2<x<0.4,0.6<y<0.8
f(x,y)=sin(x2y+1),0.2<x<0.4,0.6<y<0.8,利用
S
a
g
e
m
a
t
h
Sagemath
Sagemath求导可得:
sage: x, y = var('x,y')
sage: f = sin(x^2*y+1)
sage: f.diff(x, 4)
16*x^4*y^4*sin(x^2*y + 1) - 48*x^2*y^3*cos(x^2*y + 1) - 12*y^2*sin(x^2*y + 1)
sage: f.diff(y, 4)
x^8*sin(x^2*y + 1)
sage: f.diff(x, 4, y, 4)
16*x^12*y^4*sin(x^2*y + 1) - 304*x^10*y^3*cos(x^2*y + 1) - 1740*x^8*y^2*sin(x^2*y + 1) + 3360*x^6*y*cos(x^2*y + 1) + 1680*x^4*sin(x^2*y + 1)
代入得:
∣
R
3
,
3
(
x
,
y
)
∣
≤
∣
0.
2
4
384
[
16
ξ
4
y
4
s
i
n
(
ξ
2
y
+
1
)
−
48
ξ
2
y
3
c
o
s
(
ξ
2
y
+
1
)
−
12
y
2
s
i
n
(
ξ
2
y
+
1
)
]
∣
+
∣
0.
2
4
384
x
8
s
i
n
(
x
2
η
+
1
)
∣
+
∣
0.
2
8
38
4
2
[
16
ξ
12
η
4
s
i
n
(
ξ
2
η
+
1
)
−
304
ξ
10
η
3
c
o
s
(
ξ
2
η
+
1
)
−
1740
ξ
8
η
2
s
i
n
(
ξ
2
η
+
1
)
+
3360
ξ
6
η
c
o
s
(
ξ
2
η
+
1
)
+
1680
ξ
4
s
i
n
(
ξ
2
η
+
1
)
]
∣
≤
3.806762660674352
×
1
0
−
5
\begin{align} |R_{3,3}(x,y)|&\le \left|\dfrac{0.2^4}{384}[16\xi^4y^4sin(\xi^2y + 1) - 48\xi^2y^3cos(\xi^2y + 1) - 12y^2sin(\xi^2y + 1)]\right|\\ &+\left|\dfrac{0.2^4}{384}x^8sin(x^2\eta + 1) \right|+\Bigg|\dfrac{0.2^8}{384^2}[16\xi^{12}\eta^4sin(\xi^2\eta + 1) - 304\xi^{10}\eta^3cos(\xi^2\eta + 1)\\ &- 1740\xi^8\eta^2sin(\xi^2\eta + 1) + 3360\xi^6\eta cos(\xi^2\eta + 1) + 1680\xi^4sin(\xi^2\eta + 1)]\Bigg|\\ &\le 3.806762660674352\times 10^{-5} \end{align}
∣R3,3(x,y)∣≤∣
∣3840.24[16ξ4y4sin(ξ2y+1)−48ξ2y3cos(ξ2y+1)−12y2sin(ξ2y+1)]∣
∣+∣
∣3840.24x8sin(x2η+1)∣
∣+∣
∣38420.28[16ξ12η4sin(ξ2η+1)−304ξ10η3cos(ξ2η+1)−1740ξ8η2sin(ξ2η+1)+3360ξ6ηcos(ξ2η+1)+1680ξ4sin(ξ2η+1)]∣
∣≤3.806762660674352×10−5
参考文献
[1] M. Gasca and T. Sauer, “On the history of multivariate polynomial interpolation,” J. Comput. Appl. Math., vol. 122, no. 1–2, pp. 23–35, 2000.
[2] Wikipedia contributors. (2022, September 21). Bicubic interpolation. In Wikipedia, The Free Encyclopedia. Retrieved 08:11, October 9, 2022, from https://en.wikipedia.org/w/index.php?title=Bicubic_interpolation&oldid=1111555946.
[3] Wikipedia contributors. (2022, August 14). Multivariate interpolation. In Wikipedia, The Free Encyclopedia. Retrieved 08:20, October 9, 2022, from https://en.wikipedia.org/w/index.php?title=Multivariate_interpolation&oldid=1104415950.
[4] A. C. Ahlin, ”A bivariate generalization of Hermite’s interpolation formula”, Math. Comp.
18 (1964), 264-273.
[5] M. M. Chawla and N. Jayarajan, “A generalization of Hermite’s interpolation formula in two variables”, J. Austral. Math. Soc., vol. 18, no. 4, pp. 402-410, Dec. 1974.
[6] A. Spitzbart, “A generalization of Hermite’s interpolation formula”, Amer. Math. Monthly 67 (1960), 42-6.
[7] T. Sauer and Yuan Xu, “On multivariate Lagrange interpolation”, Math. Comp. 64 (1995) 1147–1170.