圣维南方程的3次B样条拟插值数值解_钱江.pdf

h t t p:/ww w.j s j k x.c o mD O I:1 0.1 1 8 9 6/j s j k x.2 2 0 8 0 0 1 1 8到稿日期:2 0 2 2-0 8-1 2 返修日期:2 0 2 2-1 1-0 9通信作者:张鼎(2 0 1 3 1 2 0 1 0 0 0 9h h u.e d u.c n)圣维南方程的3次B样条拟插值数值解钱 江张 鼎河海大学理学院 南京2 1 1 1 0 0(q i a n j i a n g j o b h h u.e d u.c n)摘 要 首先针对不同阶的连续可导函数,对3次B样条拟插值算子进行相应的误差估计。其次将3次B样条拟插值方法用于求解圣维南方程,利用3次B样条拟插值的一阶导数近似圣维南方程的空间导数,同时使用向前差分近似其一阶时间导数,求出其数值解。最后将所得结果与4阶龙格库塔法和蛙跳格式所得数值解进行对比分析,结果表明3次B样条拟插值方法具有一定的优越性。关键词:B样条;样条拟插值;圣维南方程;偏微分方程数值解中图法分类号 O 2 4 1.8 2 N u m e r i c a l S o l u t i o no fS a i n t-V e n a n tE q u a t i o nb yC u b i cB-s p l i n eQ u a s i-i n t e r p o l a t i o nQ I ANJ i a n ga n dZ HANGD i n gC o l l e g eo fS c i e n c e,H o h a iU n i v e r s i t y,N a n j i n g2 1 1 1 0 0,C h i n a A b s t r a c t F i r s t l y,t h ee r r o re s t i m a t e so fc u b i cs p l i n eq u a s i-i n t e p o l a t i n go p e r a t o r sa r ed e r i v e df o rc o n t i n u o u sd i f f e r e n t i a l f u n c t i o nw i t hd i f f e r e n to r d e r s.S e c o n d l y,c u b i cB-s p l i n eq u a s i-i n t e r p o l a t i o n i su s e dt og e t t h en u m e r i c a l s o l u t i o no fS a i n t-V e n a n t e q u a t i o n.S p e c i f i c a l l y,t h ed e r i v a t i v e so f t h eq u a s i-i n t e r p o l a t i o na r eu s e dt oa p p r o x i m a t e t h es p a t i a l d e r i v a t i v eo f t h ed e p e n d e n t v a r i a b l ea n df o r w a r dd i f f e r e n c em e t h o d i su s e dt oa p p r o x i m a t et h et i m ed e r i v a t i v eo f t h ed e p e n d e n tv a r i a b l e.F i n a l l y,t h en u m e r i c a l s o l u t i o n sa r ec o m p a r e dw i t ht h e s o l u t i o no b t a i n e db y t h e f o u r t ho r d e rR u n g e-K u t t am e t h o da n d t h e l e a p f r o gs c h e m e.T h e nn u m e r i c a l e x a m-p l e ss h o wt h a t c u b i cs p l i n eq u a s i-i n t e p o l a t i n gm e t h o dh a ss o m ea d v a n t a g e s.K e y w o r d s Bs p l i n e,S p l i n eq u a s i-i n t e p o l a t i o n,S a i n t-V e n a n t e q u a t i o n,N u m e r i c a l s o l u t i o n so fp a r t i a l d i f f e r e n t i a l e q u a t i o n s 1 引言样条函数是保持一定连续阶的分段或分片的多项式函数1,由于其具有良好的性质,如保凸性、光滑性等,因此在计算机几何设计、微分方程数值解等诸多领域有着较好的应用。拟插值方法作为函数逼近中的重要方法,由于其无须求解方程组的优点而得到了广泛的研究与应用。W a n g在文献2 中对多元样条函数及其应用进行了论述,而后又构建了非均匀2型三角剖分上的二元2次拟插值算子3-4。文献5 给出了一类非均匀二元3次样条空间S1,23(2)m n)的样条拟插值,并在后续给出了其逼近误差6。文献7 对有界区间上的离散拟插值进行了进一步的研究,验证了B样条拟插值方法对某类函数一阶导数的逼近效果优于有限差分法。在此基础上,诸多学者尝试使用一元3次B样条拟插值方法求解方 程 的 数 值 解8-1 4,如D e g a s p e r i s-P r o c e s i方 程、B u r g e r s-H u x l e y方程、波动 方 程 等,且 都 取 得 了 较 好 的 数 值 结 果。文献1 5 利用光滑余因子方法给出了4次B样 条基 函数及4次B样条拟插值,并将其用于求解对流扩散方程。文献1 6 将3次B样条拟插值方法用于求 解 二维 的非 定常平流扩散方程。圣维南方程1 7作为水动力学中的基本方程,具有十分重要的意义。圣维南方程所描述的物理现象在现实生活中广泛存在,且其在洪水预报、水路规划、环境保护等方面都起着一定的作用。但除了在极少数特殊情况下,圣维南方程一般难以求出解析解。因此,诸多学者针对不同情形下的圣维南方程,求解其数值解1 8-2 0。例如,文献2 1 提出了一种计算效率高的混合有限体积/有限差分法,用于求解一维明渠流中的圣维南方程;文献2 2 采用P r e d i c t o r-C o r r e c t o r方法进行求解;文献2 3 分别使用有限差分法以及特征线法求解圣维南方程,结果发现有限差分法的结果准确度高于特征线法,且随着时间以及距离的增加,两者的误差也逐渐增大;文献2 4 将圣维南方程组转化为常微分方程,使用4阶龙格库塔法进行求解;文献2 5 采用G o d u n o v格式的有限体积法求解双曲守恒型圣维南方程;文献2 6 在P r e i s s m a n n差分格式研究的基础上,提出了4种优化算法并进行了仿真验证。据笔者所知,至今鲜有文章研究3次B样条拟插值的方法在圣维南方程组中的应用,故本文考虑将3次B样条拟插值方法用于求解圣维南方程的数值解。本文回顾3次B样条拟插值,针对不同阶的连续可导函数进行相应的误差分析,并利用3次B样条拟插值求解圣维南方程。具体步骤为:利用3次B样条拟插值的一阶导数来近似微分方程组的空间导数,利用欧拉方法近似微分方程组的时间导数。第2节引入了单变量3次B样条基函数及3次B样条拟插值的构造格式。第3节针对不同光滑度的函数,对3次B样条拟插值算子进行误差估计。第4节使用3次B样条拟插值方法求解圣维南方程的数值解,并与4阶龙格库塔法及蛙跳格式进行对比。2 3次B样条拟插值本文的主要内容涉及3次B样条基函数与样条拟插值,故在本节对二者进行简单的介绍。给定参数轴x轴的一个剖分XN=xj=a+j h,0jN,可由光滑余因子协调法1或d eB o o r-C o x公式1递推得到3次B样条基函数Bj,3。Bi,0(x)=1,xxi,xi+1)0,xxi,xi+1)Bj(x)=(x-xj)3(xj+1-xj)(xj+2-xj)(xj+3-xj),xxj,xj+1)(x-xj)2(xj+2-x)(xj+2-xj)(xj+2-xj+1)(xj+3-xj)+(x-xj)(xj+3-x)(x-xj+1)(xj+3-xj)(xj+3-xj+1)(xj+2-xj+1)+(xj+4-x)(x-xj+1)2(xj+2-xj+1)(xj+3-xj+1)(xj+4-xj+1),xxj+1,xj+2)(x-xj)(xj+3-x)2(xj+3-xj)(xj+3-xj+1)(xj+3-xj+2)+(x-xj+1)(xj+3-x)(xj+4-x)(xj+3-xj+2)(xj+3-xj+1)(xj+4-xj+1)+(xj+4-x)2(x-xj+2)2(xj+4-xj+1)(xj+4-xj+2)(xj+3-xj+2),xxj+2,xj+3)(xj+4-x)3(xj+4-xj+1)(xj+4-xj+2)(xj+4-xj+3),xxj+3,xj+4)(1)用B(4),B(3),B(2)分别表示4重、3重及2重节点情形的3次B样条基函数。3次B样条基函数图像(包含重节点)如图1所示。图2为非均匀情形下的B样条基函数,具体的节点为0,1,3,4,7。图3为节点为0,4,6,7,1 0情形下的B样条基函数图像。对于区间I=a,b,给定均匀划分XN=xj=a+j h,0jn,h=b-an。用Sd(Xn)表示d次一元样条空间。Bj的支撑集记为s u p p(Bj)=xj-d-1,xj。在两侧端点处增加重节点x-d=x-d+1=x-d+2=x-1=x0=a,b=xn=xn+1=xn+d。d次一元B样条拟插值的一般形式为:Qdf=n+dj=1jBj(x)f(2)其中,j(f)与x无关,其为函数f在s u p p(Bj)区域中某些节点处离散值fi的线性组合。本文考虑使用3次B样条拟插值求解圣维南方程,设fi=f(xi),i=0,1,2,n,x-3=x-2=x-1=x0 x1x2xn=xn+1=xn+2=xn+3,空间步长h=xi-xi-1,i=1,2,n,可得:定理17 一元3次B样条基函数的拟插值算子为:Qnf=n+3j=1j(f)Bj,3(x)(3)使其在S3(Xn)上精确成立,即对f=1,x,x2,x3成立,可得其中各项系数为:1(f)=f0,2(f)=11 8(7f0+1 8f1-9f2+2f3),j(f)=16(-fj-3+8fj-2-fj-1),j=3,n+1n+2(f)=11 8(2fn-3-9fn-2+1 8fn-1+7fn)n+3(f)=fn 利用3次B样条拟插值的一阶导数近似因变量的空间导数,即Qnf=n+3j=1jB j(x),可得:Qnf(xi)=1h-1 16f0+3f1-32f2+13f3(),i=01h-13f0-12f1+f2-16f3(),i=11h-11 2fi-2-23fi-1+23fi+1-11 2fj+2(),2in-21h13fn+12fn-1-fn-2+16fn-3(),i=n-11h1 16fn-3fn-1+32fn-2-13fn-3(),i=n图1 3次B样条基函数F i g.1 C u b i cB-s p l i n eb a s i s f u n c t i o n s图2 非均匀节点为0,1,3,4,7时的三次B样条基函数F i g.2 C u b i cB-s p l i n eb a s i s f u n c t i o n sw h e nt h en o n-u n i f o r mn o d e sa r e0,1,3,4,7621C o m p u t e rS c i e n c e计算机科学V o l.5 0,N o.4,A p r.2 0 2 3图3 非均匀节点为0,4,6,7,1 0时的三次B样条基函数F i g.3 C u b i cB-s p l i n eb a s i s f u n c t i o n sw h e nt h en o n-u n i f o r mn o d e sa r e0,4,6,7,1 03 样条拟插值误差分析考虑到文献7 给出了样条拟插值误差分析的粗略结论,在本节中,我们将利用连续模与极大范数,计算出3次B样条拟插值算子误差分析的具体表达式。对于区间I=x0,xn,x0 x1xn,步长h=xi+1-xi,i=0,1,2,n-1。定义连续模和极大范数为:I(f,r)=s u p|f(x)-f(u)|:x,uI,|x-u|r(4)fI=fL(I)=e s s s u pfI|f|(5)定理2 若f(x)C(I),则 对于 足够 大的 正 整数n,Qn(f)的误差估计满足:(1)当i=0,1,n-2,n-1时,f-QnfIi2I(f,4h)(6)(2)当i=2,3,n-3时,f-QnfIi53I(f,4h)(7)证明:当xI0,由B样 条 基 函 数 的 单 位 分 解 性 以 及式(2)可得:f-Qnf=f(x)(B(3)-3(I0)+B(2)-2(I0)+B(1)-1(I0)+B0(I0)-Qnf则f-QnfI0=I(f,4h)(B(3)-3(I0)+2B(2)-2(I0)+53B(1)-1(I0)+53B0(I0)2I(f,4h)当xI1,In-2,In-1时,与上述情形类似。当xIi,i=2,3,n-3时,f-QnfI0=I(f,4h)53Bi-3(Ii)+53Bi-2(Ii)+(53Bi-1(Ii)+53Bi(Ii)=53I(f,4h)定理3 若f(x)C1(I),则对于足够大的正整数n,Qn(f)的误差估计满足:(1)当i=0,1,n-2,n-1时,f-QnfIi32h If,h2()(8)(2)当i=2,3,n-3时,f-QnfIi43h If,h2()(9)证明:取函数f(x)C1(I),xI0=x0,x1,将f(x)在区间中点x1/2处进行一阶泰勒展开,此时我们可以得到:f(x)=T1,0+(f(0)-f(x1/2)(x-x1/2)即:f(x)-T1,0=(f(0)-f(x1/2)(x-x1/2)其中,0I(x,x1/2)x0,x1,T1,0=f(x1/2)+f(x1/2)(x-x1/2)+f(0)(x-x1/2)。注意到QnfI0=s u pf0|Qn(f)|f|=2,故有:f-QnfI0=f-T1,0-Qn(f)+Qn(T1,0)I03f-T1,0I03|x-x1/2|If,h2()3h2If,h2()当xI1,In-2,In-1时,同理可证。当xIi,i=2,3,n-3时,f(x)=T1,i+(f(i)-f(xi+1/2)(x-xi+1/2),可得:f(x)-T1,i=(f(i)-f(xi+1/2)(x-xi+1/2)其中T1,i=f(xi+1/2)+f(xi+1/2)(x-xi+1/2),iI(xi,xi+1/2)xi,xi+1,QnfIi=s u pf0|Qn(f)|f|=53。此 时可得:f(x)-Qn(f)=f-T1,i-Qn(f)+Qn(T1,i)83f-T1,iIi83|x-xi+1/2|If,h2()4h3If,h2()定理4 当f(x)C2(I),对于任意大的正整数n:(1)当i=0,1,n-2,n-1时,f-Qn(f)Ii3h28If,h2()(1 0)(2)当i=2,3,n-3时,f-Qn(f)Iih33If,h2()(1 1)证明:同理,对于f(x)C2(x),利用2阶泰勒多项式展开,可以得到:T2,i(x)=f(xi+1/2)+f(xi+1/2)(x-xi+1/2)+f(xi+1/2)(x-xi+1/2)2f(x)=T2,i(x)+12(f(i)-f(xi+1/2)(x-xi+1/2)2其中,iI(x,xi+1/2),xIi,i=0,1,n-1,此时可得:f-Qn(f)Ii83f-T2,iIi86|x-xi+1/2|2If,h2()h23If,h2(),xIi,i=2,3,n-3f-Qn(f)Ii3f-T2,iIi32|x-x1/2|2If,h2()3h28If,h2(),xIi,i=0,1,n-2,n-1定理5 当f(x)C3(I),对于任意大的正整数n:721钱 江,等:圣维南方程的3次B样条拟插值数值解(1)当i=0,1,n-2,n-1时,f-Qn(f)Iih31 6If,h2()(1 2)(2)当i=,2,3,n-3时f-Qn(f)Iih31 8If,h2()(1 3)证明:当f(x)C3(I),我们在区间中点处进行三阶泰勒多项式展开。T3,i(x)=f(xi+1/2)+f(xi+1/2)(x-xi+1/2)+12f(xi+1/2)(x-xi+1/2)2+16f(xi+1/2)(x-xi+1/2)3f(x)=T3,i(x)+16(f(i)-f(xi+1/2)(x-xi+1/2)3其中,iI(x,xi+1/2)Ii,xIi,i=0,1,n-1,此时可得:f-Qn(f)Ii83f-T3,iIi81 8|x-xi+1/2|3If,h2()h31 8If,h2(),xIi,i=2,3,n-3f-Qn(f)Ii3f-T3,iIi36|x-x1/2|3If,h2()h31 6If,h2(),xIi,i=0,1,n-2,n-1定理6 当f(x)C4(I),对于任意大的正整数n:(1)当i=0,1,n-2,n-1时,f-Qn(f)Iih41 2 8If(4)I(1 4)(2)当i=,2,3,n-3时,f-Qn(f)Iih41 4 4If(4)I(1 5)证明:此时,我们使用f(x)的4阶泰勒公式来进行误差估计:f(x)=T4,i+14!f(4)(i)(x-xi+1/2)4其中iI(x,xi+1/2),可得:f-Qn(f)Ii3f-T4,i 34!|x-xi+1/2|4f(4)Ih41 2 8f(4)I,xIi,i=0,1,n-2,n-1f-Qn(f)Ii83f-T4,i 87 2|x-xi+1/2|4f(4)Ih41 4 4f(4)I,xIi,i=2,3,n-34 圣维南方程数值求解本节我们将利用3次样条拟插值对圣维南方程进行数值模拟,圣维南方程的基本形式为1 7:At+Qx=0(1 6)Qt+xQ2A()+g Ahx+sf-s0()=0(1 7)其中,式(1 6)为连续方程,式(1 7)为动量方程。由于现实中实际问题的多样性,圣维南方程产生了多种形态,但其实质并未改变。我们考虑如下形式的圣维南方程2 3:t+(h u)x=0(h u)t+xh u2+12gc o s()h2()=gs i n()h-Cfu|u|h其中,为地表高度,h为水面高度,u是水流的水平速度,是河床倾角,d0为水深,g是重力加速度,h=d0+。同时Cf是底部摩擦参数,其值取决于常数的选择2 3。我们首先对圣维南方程组进行离散,使用向前差分近似因变量的时间导数,利用3次B样条拟插值的一阶导数近似因变量的空间导数。此处我们模拟波浪在浅水区域的传播,因此,需取=0,Cf=0,此时圣维南方程组的离散形式为:hk+1j=hkj-t(hkj(ux)kj+ukj(hx)kj)(1 8)uk+1j=ukj-tukj(ux)kj-tg(hx)kj)(1 9)给定初值条件h(x,0)=Atc o s(ktx-t),其 中,At=0.2 5,=0.5,kt=0.4 5 5 6,u(x,0)=0,g=9.8。取时间步长 t=0.0 1,空间步长 x=0.5,x0,1 0,将数值代入计算可以得到如下的波传播图像(见图4),同时使用4阶龙格库塔法与蛙跳格式差分进行求解,得到的图像如图5和图6所示。图4 t=0.0 1,x=0.5时三次B样条拟插值法所得的高度h随时间t与空间x的变化曲线F i g.4 C h a n g ec u r v eo f t h eh e i g h tho b t a i n e db yt h ec u b i cB-s p l i n eq u a s i-i n t e r p o l a t i o nm e t h o dw i t ht i m eta n ds p a c exw h e nt=0.0 1,x=0.5图5 t=0.0 1,x=0.5时四阶龙格库塔法所得的高度h随时间t与空间x的变化曲线F i g.5 C h a n g ec u r v eo f t h eh e i g h tho b t a i n e db yt h e f o u r t h-o r d e rR u n g e-K u t t am e t h o dw i t ht i m eta n ds p a c exw h e n t=0.0 1,x=0.5821C o m p u t e rS c i e n c e计算机科学V o l.5 0,N o.4,A p r.2 0 2 3图6 t=0.0 1,x=0.5时蛙跳格式所得的高度h随时间t与空间x的变化曲线F i g.6 C h a n g ec u r v eo f t h eh e i g h tho b t a i n e db yt h e l e a p f r o gs c h e m ew i t ht i m eta n ds p a c exw h e n t=0.0 1,x=0.5将3次B样条拟插值方法及4阶龙格库塔法所得的数值结果进行对比,可以得到差值图像,如图7和图8所示。图7 t=0.0 1,x=0.5两种方法所得h的差值F i g.7 D i f f e r e n c eo fhw h e n t=0.0 1,x=0.5图8 t=0.5时3种方法对应的曲线F i g.8 C o r r e s p o n d i n gc u r v e so f t h r e em e t h o d sw h e nt=0.5各方法在时刻t=0.5,x分别位于不同位置时的具体数值如表1所列。表1 3种方法在不同x处的对应值T a b l e1 C o r r e s p o n d i n gv a l u e so f t h r e em e t h o d sa td i f f e r e n txx3次B样条拟插值方法4阶龙格库塔法蛙跳格式0.00.2 3 5 1 1 2 9 3 80.3 8 5 3 5 6 2 9 90.3 1 4 5 6 6 4 3 51.00.2 1 5 2 8 2 2 8 90.2 3 1 2 6 3 4 6 30.2 3 6 2 9 5 3 3 61.50.1 9 0 1 7 8 2 4 90.1 9 8 8 5 1 0 2 50.2 0 1 6 1 6 7 3 63.00.0 6 4 0 5 4 0 7 50.0 9 0 4 1 9 0 3 30.0 6 0 4 4 2 4 0 45.0-0.1 5 8 6 7 8 5 7 8-0.1 2 9 8 4 5 0 3-0.1 3 5 2 6 1 0 3 96.5-0.2 6 1 8 1 1 6 0 8-0.2 6 3 1 3 9 6 8-0.2 2 1 2 5 2 8 9 37.0-0.2 6 6 7 2 7 3 3 4-0.2 6 9 2 1 9 1 8-0.2 2 9 9 2 8 6 8 9取时间步长 t=0.0 0 1,空间步长 x=0.5,使用3次B样条拟插值方法、4阶龙格库塔法及蛙跳格式,得到的图像如图9-图1 1所示。图1 2和图1 3分别为时间t=0.0 1,t=0.1时刻3种方法所对应的曲线。3次B样条拟插值方法与4阶龙格库塔法所得结果的差值图像如图1 4所示。图9 t=0.0 0 1,x=0.5时三次B样条插值法所得高度h随时间t与空间x的变化曲线F i g.9 C h a n g ec u r v eo f t h eh e i g h tho b t a i n e db yt h ec u b i cB-s p l i n eq u a s i-i n t e r p o l a t i o nm e t h o dw i t ht i m eta n ds p a c exw h e n t=0.0 0 1,x=0.5图1 0 t=0.0 0 1,x=0.5时四阶龙格库塔法所得的高度h随时间t与空间x的变化曲线F i g.1 0 C h a n g ec u r v eo f t h eh e i g h tho b t a i n e db yt h e f o u r t h-o r d e rR u n g e-K u t t am e t h o dw i t ht i m eta n ds p a c exw h e n t=0.0 0 1,x=0.5图1 1 t=0.0 0 1,x=0.5时蛙跳格式所得的高度h随时间t与空间x的变化曲线F i g.1 1 C h a n g ec u r v eo f t h eh e i g h tho b t a i n e db yt h e l e a p f r o gs c h e m ew i t ht i m eta n ds p a c exw h e nw h e n t=0.0 0 1,x=0.5图1 2 t=0.0 1时,3种方法的对应曲线F i g.1 2 C o r r e s p o n d i n gc u r v e so f t h r e em e t h o d sw h e nt=0.0 1图1 3 t=0.1时,3种方法的对应曲线F i g.1 3 C o r r e s p o n d i n gc u r v e so f t h r e em e t h o d sw h e nt=0.1921钱 江,等:圣维南方程的3次B样条拟插值数值解图1 4 t=0.0 0 1,x=0.5两种方法所得h的差值F i g.1 4 D i f f e r e n c eo fhw h e n t=0.0 0 1,x=0.1各方法在时刻t=0.1,x分别位于不同位置时的具体数值如表2所列。表2 3种方法在不同x处的对应值T a b l e2 C o r r e s p o n d i n gv a l u e so f t h r e em e t h o d sw i t hd i f f e r e n txx3次B样条拟插值方法4阶龙格库塔法蛙跳格式0.00.2 4 9 3 5 0 2 9 00.2 7 8 3 2 1 0 6 20.2 5 2 6 0 8 3 2 01.00.2 2 4 1 1 6 2 8 90.2 2 5 3 3 7 0 3 30.2 2 9 4 0 6 5 4 31.50.1 9 3 7 2 8 8 9 90.1 9 6 1 3 6 3 0 50.1 9 7 6 8 9 5 5 83.00.0 5 1 2 2 2 3 3 20.0 5 7 3 9 2 2 9 50.0 5 1 5 0 7 4 5 85.0-0.1 6 2 3 2 7 6 5 8-0.1 5 8 4 1 7 6 8 8-0.1 5 9 7 1 3 1 3 26.5-0.2 4 6 5 4 3 0 8 0-0.2 4 6 3 4 8 2 4 9-0.2 4 0 4 7 5 0 8 17.0-0.2 5 0 3 4 5 4 0 9-0.2 5 0 3 3 5 8 9 3-0.2 4 4 0 7 6 3 0 7更换时间步长,取时间步长 t=11 0-6,空间步长 x=0.5,t0,51 0-5,x0,1 0,其余条件不变,得到如下图像。其中图1 5-图1 7分别为使用3次B样条拟插值、4阶龙格库塔法及蛙跳格式所得数值结果;图1 8为三者在时刻t=21 0-5的对应曲线。图1 5 t=11 0-6,x=0.5时三次B样条插值法所得的高度h随时间t与空间x的变化曲线F i g.1 5 C h a n g ec u r v eo f t h eh e i g h tho b t a i n e db yt h ec u b i cB-s p l i n eq u a s i-i n t e r p o l a t i o nm e t h o dw i t ht i m eta n ds p a c exw h e n t=11 0-6,x=0.5图1 6 t=1 1 0-6,x=0.5时四阶龙格库塔法所得的高度h随时间t与空间x的变化曲线F i g.1 6 C h a n g ec u r v eo f t h eh e i g h tho b t a i n e db yt h e f o u r t h-o r d e rR u n g e-K u t t am e t h o dw i t ht i m eta n ds p a c exw h e n t=11 0-6,x=0.5图1 7 t=11 0-6,x=0.5时蛙跳格式所得的高度h随时间t与空间x的变化曲线F i g.1 7 C h a n g ec u r v eo f t h eh e i g h tho b t a i n e db yt h e l e a p f r o gs c h e m ew i t ht i m eta n ds p a c exw h e n t=11 0-6,x=0.5图1 8 t=21 0-5时3种方法的对应曲线F i g.1 8 C o r r e s p o n d i n gc u r v e so f t h r e em e t h o d sw h e nt=21 0-5特别地,取时间步长 t=0.2,空间步长 x=0.5,3次B样条拟插值方法,4阶龙格库塔法及蛙跳格式所得图像分别如图1 9-图2 1所示。图1 9 t=0.2,x=0.5时三次B样条插值法所得的高度h随时间t与空间x的变化曲线F i g.1 9 C h a n g ec u r v eo f t h eh e i g h tho b t a i n e db yt h ec u b i cB-s p l i n eq u a s i-i n t e r p o l a t i o nm e t h o dw i t ht i m eta n ds p a c exw h e n t=0.2,x=0.5图2 0 t=0.2,x=0.5时四阶龙格库塔法所得的高度h随时间t与空间x的变化曲线F i g.2 0 C h a n g ec u r v eo f t h eh e i g h tho b t a i n e db yt h e f o u r t h-o r d e rR u n g e-K u t t am e t h o dw i t ht i m eta n ds p a c exw h e n t=0.2,x=0.5031C o m p u t e rS c i e n c e计算机科学V o l.5 0,N o.4,A p r.2 0 2 3图2 1 t=0.2,x=0.5时蛙跳格式所得的高度h随时间x与空间t的变化曲线F i g.2 1 C h a n g ec u r v eo f t h eh e i g h tho b t a i n e db yt h e l e a p f r o gs c h e m ew i t ht i m eta n ds p a c exw h e n t=0.2,x=0.5对比相同步长时,3种方法所得图像及表1、表2可以发现,3次B样条拟插值方法所得数值解和4阶龙格库塔法、蛙跳格式所得数值解误差较小,且此时3次B样条拟插值方法所得数值解更稳定,得到的图像更光滑。通过图8中的对比可知,相比4阶龙格库塔法,3次B样条拟插值方法与蛙跳格式所得数值解更为接近,且当x2.5,4 时,两者的差值较小。相比图2 0和图2 1,图2 1拐点出现的时间更晚,且拐点的波动幅度相对较小。图2 2和图2 3中不光滑点的数量明显多于图1 9。同时,在数值模拟过程中,我们发现当x0,7,x=0.5,t=0.0 0 1时,3次B样条拟插值方法相比在其他条件时更加稳定。结束语 本文回顾了3次B样条拟插值方法,并利用不同阶的连续可导 函数,对拟 插 值算 子进 行了 误差 估 计。应用3次B样条拟插值方法求解圣维南方程组,获得其数值解,与4阶龙格 库 塔法 和蛙 跳 格式 相比,所得 的数 值结果误差较小。当3种方法所得的数值结果都保持稳定时,3种方法计 算 所 得 结 果 也 并 无 较 大 误 差。当 利 用 向 前 差分近似时间导数以及利用3次B样条拟插值方法求解时,相比4阶龙格库塔法和蛙跳格式差分法,减少了所需要的边界条件。且在上述实验中,我们发现当时间x0,7,x=0.5,t0.0 0 1时,3次B样条拟插值方法具有更好的稳定性。本文尚未对3次B样条 拟插 值 方 法 求 解 圣 维 南 方 程的稳定条件进行理论分析,后续可以考虑更换近似时间导数的方法,如求解 圣 维南 方程 组 常用 的蛙 跳格 式等,结合3次B样条拟插值方法进行求解,以期得到更为精确且稳定的数值解。参 考 文 献1WAN GR H,L ICJ,Z HUCG.C o m p u t a t i o n a lg e o m e t r yT u t o-r i a lM.B e i j i n g:S c i e n c eP r e s s,2 0 0 8:6-2 2 1.2WAN G R H.M u l t i v a r i a t es p l i n ef u n c t i o na n di t sa p p l i c a t i o nM.B e i j i n g:S c i e n c eP r e s s,1 9 9 4:1-4 5 5.3WAN GRH,L UY.Q u a s i-i n t e r p o l a t i n go p e r a t o r sa n dt h e i ra p-p l i c a t i o n s i nh y p e r s i n g u l a ri n t e g r a l sJ.J o u r n a lo fC o m p u t a-t i o n a lM a t h e m a t i c s,1 9 9 8,1 6(4):3 3 7-3 4 4.4WAN GR H,L U Y.Q u a s i-i n t e r p o l a t i n go p e r a t o r si nS12(2)m n)o nn o n u n i f o r mt y p e-2t r i a n g u l a t i o n sJ.N u m e r i c a lM a t h m a t i e sAJ o u r n a l o fC h i n e s eU n i v e r s i t i e s,1 9 9 9,(2):9 7-1 0 3.5Q I ANJ,WAN GR H,Z HUCG,e t a l.O ns p l i n eq u a s i-i n t e r p o-l a t i o n i nc u b i c s p l i n e s p a c eS1,23(2)m n)J.S c i e n t i aS i n i c a(M a t h e-m a t i c a),2 0 1 4,4 4(7):7 6 9-7 7 8.6Q I ANJ,WAN GF.O nt