标题

Size: px

Start display at page:

Download "标题"

偃例满
5 years ago
Views:

1 第 53 卷第期厦门大学学报 ( 自然科学版 ) Vl.53 N. 04 年 3 月 JunlfXmenUnvesty (NtulSene) M.04 d:0.6043/j.ssn 反常次扩散问题的有限元逼近吴春红, 刘青霞 (. 厦门理工学院应用数学学院, 福建厦门 3604;. 厦门大学数学科学学院, 福建厦门 36005) 摘要 : 讨论一类反常次扩散问题, 进行了有限元数值模拟, 分别给出了其时间半离散时间空间全离散形式, 并且讨论了两种形式的稳定性收敛性. 最后给出数值例子显示所提出的数值方法的有效性. 关键词 : 反常次扩散问题 ; 有限元方法 ; 稳定性 ; 收敛性中图分类号 : O4.8 文献标志码 : A 文章编号 : (04) 预备知识自然界与工程应用中存在很多扩散现象, 如地下水溶质运移污染物在土壤中的迁移核废料的扩散药物在身体中的扩散石油渗流湍流等, 这些扩散现象一般不满足经典的 k 梯度扩散律, 通常称为反常扩散. 反常扩散现象很难用整数阶微分方程表征, 近年来研究者们广泛采用分数阶微分方程来描述反常扩散过程, 并求助数值方法进行模拟. 文献 [] 最早采用谱方法处理时间分数阶扩散方程, 随后进行推 [] 广, 将谱方法应用到空间时间分数阶扩散方程. 文献 [3] 采用有限元方法考虑空间分数阶对流 - 扩散方程. 更多的文献均采用有限差分方法处理分数阶次扩散 [4-6] 问题. 文献 [7] 采用 Admn 分裂方法给出时间分数阶 Klen-Gdn 方程的解析近似解. 本文尝试有限元方法处理分数阶导数施加于右侧扩散项的次扩散问题. 本文讨论下面带有初边值条件 u(x,0)=ω(x),x Ω=[0,L], () u(0,t)=u(l,t)=0,t [0,T] () 的反常次扩散问题 (ASE): u(x,t) t = u(x,t) D t -α,x Ω,t [0,T], x (3) 其中 0<α, >0,D t -α 表示关于 t 的 -α 阶 Remmn-Luvle 分数阶导数, 其定义为 : Dt -α u(x,t)= Γ(α) t t 0 其中 Γ( ) 是 Gmm 函数. u(x,τ) (t-τ) -α dτ, 本文采用有限元方法求解上述 ASE 问题, 分别给出了其时间半离散时间空间全离散形式, 讨论了两种形式的稳定性收敛性, 并且给出数值例子与有 [5] 限差分方法进行比较, 得到了较好的结果. ASE 半离散格式定义函数空间 : H (Ω)={v L (Ω), v x L (Ω)}, H 0(Ω)={v H (Ω),v Ω =0}, H m (Ω)={v L (Ω), k v x k L (Ω),0<k m}, (4) 其中 L (Ω) 表示平方可积的函数空间, 两个函数的内积定义为 (u,v)= Ωuvdx, 范数和半范定义为 u 0 = (u, 收稿日期 : 基金项目 : 国家自然科学基金 (0344,3094); 福建省自然科学基金 (03J00) * 通信作者 :wu@xmut.edu.n u), u k = k u k x 0,u k k = ( u 0 + = u ), 标准的 - 范数为 u =(u 0+ u ). 在文中定义 u '= u, 以及 - 范数 u' = (u 0 + u ), 其中 = 将在下面离散格式中给 Γ(α+)

2 66 厦门大学学报 ( 自然科学版 ) 04 年出. 本文中定义的 - 范数与标准 - 范数之间的关系为 : u' = u 0 + u = u 0 + Γ(α+) u Cu { } ). u (Cu =, Γ(α+) (5) 假设 ASE()~ (3) 的解满足 :u(,t) H 0 (Ω), u(x, ) C ([0,T]). x. ASE 半离散格式首先利用时间间隔 tn=nτn=0,,,,n,τ= T N ø 离散时间变量, 其中 τ 为时间步长. 对式 (3) 两端从 tn 到 tn+ 进行积分, 得到 : u(x,tn+)=u(x,tn)+ Γ(α) τ 0 tn 0 U(x ἠ ) dη+ x (tn+ -η )-α Γ(α) v(x ἠ ) x (tn -η )-α dη, ( 6) 其中 v(x,t)=u(x,t+τ)-u(x,t). 为了方便, 引入记号 : bs =(s+) α -s α,,,,,n, (7) 并对式 (6) 右端后两项分别进行如下逼近 : Γ(α) τ 0 u(x ἠ ) dη= x (tn+ -η )-α τ α bn u(x,τ) +R Γ(α+) n+ x, R n+ Cbnτ +α, Γ(α) tn v(x ἠ ) dη= 0 x (tn -η )-α n- τ Γ(α+) α bn-s- R n+, R n+ Cτ, v(x,ts+) x + 其中 C=M' / Γ(+α),C=M / α,m' 和 M 分别 u(x, ) 是对时间 t 的一阶导数和二阶导数的界. x 将上述两个逼近形式代入式 (6) 可得 : u(x,τ) u(x,tn+)=u(x,tn)+bn + x n- bn-s- u(x,ts+) u(x,ts+) - x x ø 其中 = Γ(α+), 截断误差为 : R n+ C * (τ +bnτ α+ )(C * = +R n+, (8) {C,C}). (9) 下面引入关于系数性质的引理 [5] : 引理式 (7) 定义的系数 bs 满足 : ()b0=,bn >0,n=0,,, ; ()bn >bn+,n=0,,, ; () 存在只与时间 T 和 α 有关的正系数 C', 使得 τ C'bnτ α,n=,,. 由引理, 截断误差式 (9) 满足 : R n+ C bnτ α+ (C =C * (C'+)), (0) 可知 C 与 τ 无关. 假设 un 为 u(x,tn) 的逼近解,ASE ()~(3) 的隐式差分格式为 : u n+ =u n + u n+ x n- (bs+ + -bs) u n-s x, () 再由边界条件 u n+ (0)=u n+ (L)=0, 则 ASE(8) 的弱解形式是 : 找 u n+ H 0(Ω), 使 (u n+,v)=(u n,v)- un+ bs+) un-s n- + (bs - x, v, v H 0(Ω). () x ø. ASE 半离散格式的稳定性定的. 定理对任意的 τ, 离散格式 () 是无条件稳 n- 证明定义 E n = u n 0 + bs u n-s.n=0 时,(u,v)=(u 0,v)- u. 令 v=u, 有 u 0 ( u 0 0+ u 0)- u, 即 E u 0 0. 当 n 时, 取 v=u n+, 有 u n+ 0 ( u n 0 + u n+ 0)- u n+ + n- (bs -bs+)( u n-s + u n+ ). n- 由系数 bs 的定义及 (bs -bs+)=-bn, 有 u n+ 0 n + bs u n+-s u n n- 0 + bs u n-s, 由引理 () 以及 = Γ(α+), 有 E n+ E n E u 0 0, 即 u n 0 u ASE 半离散格式的误差分析引理 (Pné 不等式 ) 假设 Ω R 为有界连通开区域, 则存在仅与区域测度有关的正常数 CΩ,

3 第期吴春红等 : 反常次扩散问题的有限元逼近 67 使得 v 0 CΩ v, v H 0(Ω). 记 u 为 ASE()~(3) 的精确解,{u n } N n=0 是带有初始条件 u 0 (x)=u(x,0) 的半离散问题式 () 的解, 并且令 η n (x)=u(x,tn)-u n, 由式 (8) () 得 : n+ ( η,v)=( ηn,v)- η n+ n- s 取 v=η n+, 有 : + (bs -bs+) η n-s x, v + (R n+,v), x ø η n+ 0 ( η n 0 + η n+ 0)- η n+ + n- (bs - bs+)( η n-s + η n+ )+ (R n+ ἠ n+ ). 记 Y n = η n n- 0 bs + η n-s, 应用 ( ϕ, ψ ) ϕ ψ 0 和引理得 : Y n+ Y n + C Ω R n+ 0. (3) bn 由截断误差 (0) 和 = Γ(α +), 式 (3) 变为 Y n+ Γ(α +)C Y n +^Cbnτ α+, 这里 ^C= ΩC σω, 与 n 和 τ 无关,σΩ 表示区域 Ω 的长度. 由于 Y 0 = η 0 0=0, 有 Y n+ n ^Cτ α+ bs C sτ (C s =^CT α ), 即 Cs 与 τ 无关. 综上得到半离散问题 () 的误差分析 : 定理假设 u 是 ASE()~ (3) 的精确解, {u n } N n=0 是带有初始条件 u 0 (x)=u(x,0) 的半离散问题 () 的解, 则有 u n -u(,tn) 0 Csτ. 记号 : 接下来考虑 ASE()~ (3) 的变分问题. 引入 A(u n+,v) =(u n+,v)+ un+, n- f =u n (bs+ + -bs) u n-s x, ( 4) 以及 (v)=(f,v), 弱形式 () 可写为 : A(u n+,v)=(v). (5) 容易证明双线性形式 A(, ) 在 H 0(Ω) H 0(Ω) 上是强制的连续的, 即有下面的引理 : 引理 3 A(v,v) v ', v H 0(Ω). 引理 4 A(u,v) C(τ)u ' v ', u,v H 0(Ω), 其中 C(τ) 为只与 τ 有关的正常数. 由 Lx-Mlm 定理可得下面定理 : 定理 3 对于固定的 τ, 离散问题 (5) 有唯一解 u n+ H 0(Ω). 3 ASE 的全离散问题假设 T 是传输区域 Ω 的剖分, 是网格参数, Pl(K) 是 K 上的 l 次多项式,X l = {v H (Ω) ( Ω ) v K Pl (K ), K T }.ASE(5) 的 Glekn 离散格式为 : 找 u n+ X=H 0 X, l 使 (u n+,v)+ un+ x, v =(u,v)+ n n- (bs -bs+) un-s 引理 5 x, vn, v X.(6) x ø ( 逼近性质 )u H l+ (Ω),T l u 是插值函数 [8], 则有 u-t l u CT l u l+,ci 为只与区域 Ω 有关的正常数. 引入正交投影算子 Π l :H X l, 则 (Π l ψ, v)+ Πl ψ x, v =( ψ,v)+ ψ x, v, v X. (7) 由式 (5) 和逼近性质可得 u-π l u ' CΠ l u l+, (8) 其中 CΠ = CuCT, 与 u 和无关. 定理 4 假设 {u n } N n=0 是带有初始条件 u 0 =Π l u 0 的全离散问题 (6) 的解,{u n } N n=0 是半离散问题 () 的解, 假设 u n+ H l+ (Ω), 则 : 得到 : u n -u n ' Cτ - l u s l+. (9) 0 s n 证明由正交投影算子 Π l 的定义 (7) 以及 () (Π l u n+,v)+ Πl u n+ x, v =(u n,v)+ n- (bs -bs+) un-s x, v, v X, x ø (0) 引入记号 :e n =u n -u n ἠ n =u n -Π u n,ē n =Π u n - u n, 显然 e n =η n +ē n. () 取 v=ē n+, 并利用不等式 ( ϕ, ψ ) ϕ 0+ 4 ψ 0, 式 (0) 与 (6) 相减得 n- ē n+ 0 + ē n+ e n 0 (bs + - bs+) e n-s -bn e n+, () 由 () 和 () 得 :

4 68 厦门大学学报 ( 自然科学版 ) 04 年 n- e n+ 0 + e n+ e n 0 (bs + - bs+) e n-s -bn ē n+ + η n+ 0 + η n+ n+ +( η,ē n+ )+ η n+ 令 E n = e n, ēn+. x x ø n- 0 + bs e n-s, 并利用不等式 ( ϕ, ψ ) ϕ ψ 0 和引理得 : E n+ E n + η n+ 0 + η n+ + C Ω bn η n+ 0 + η n+, bn 利用引理, 有 E n+ C τ - ' 0 s n+ ηs, C 为只与 T 有关的正常数. 于是有 e n+ ' Cτ - '. 0 s n+ ηs 利用引理 5, 可得结论 : u n+ -u n+ ' Cτ - l 0 s n+ u s l+. [8] 下面利用 Ntse 技巧得到 L(Ω) 中的最优估计. 假设 w 为下列变分问题的解 : A(v,w)=(v,u n+ -u n+ ), v H 0(Ω), 取 v=u n+ -u n+, 则 u n+ -u n+ 由式 () 和 (6) 可得 (3) 0 =A(u n+ -u n+,w), (4) A(u n+ -u n+,v)=0, v X, 再由式 (8) 定理 4 有 u n+ -u n+ Cτ - l 0 =A(u n+ -u n+,w -v) 0 s n+ u s l+ w l+, (5) 此处常数 C 与式 (8) 和定理 4 中的常数有关, 与 τ 和均无关. 下面来估计 w l+,w 满足式 (3), 则 w 是下列边值问题的解 : ì w - w ï =u n+ -u n+ í x ï îw(0)=w(l)=0,, 首先估计 w '. 式 (3) 中取 v=w, 再由引理 3 得 : 3 时, k w = k- w x k x, 故 k- w k C()u n+ -u n+ 0. 由此得 w l+ C()u n+ -u n+ 0, 则式 (5) 为 : u n+ -u n+ 0 C(τ) l+ u s l+,(6) 0 s n+ 其中常数 C (τ) 仅与 τ 有关. 由定理 3 以及式 (6) 可得下面定理 : 定理 5 设 {u n } N n=0 是带有初始条件 u 0 =Π l u 0 的全离散问题 (6) 的解,u(x,t) 是 ASE()~ (3) 的解, 且 u H ( H l+ (Ω) H0(Ω),[0,T] l ), 则有 u(x,tn)-u n 4 数值例子 0 Csτ+ C(τ) l+ u s l+. 0 s n 例利用有限元法近似下列问题 : u(x,t) t = u(x,t) D t -α +f(x,t), x x [0,π],t>0, u(0,t)=u(π,t)=0,t>0,u(x,0)=0, x [0,π], (7) 其中 =,f(x,t)=snx (+α)t α + Γ (+α) Γ(+3α) t3α ø 其精确解为 u(x,t)=t +α snx. 本例均计算 t= 时的误差, 并定义精确解与数值解的误差为 ε= u n -u n 0. 表给出了二阶有限差分方法和采用 P(K ) 的有限元方法的误差以及收敛率, 从中可以看出, 虽然有限差分方法可以得到二阶格式, 但是显然没有有限元方法的误差精度高. 从图 3 可以看出, 只要选取到合适的基函数, 有限元方法可以达到更高阶的精度. 图是误差与时间步长的关系图, 图 3 是误差与空间步表 τ= 000 ἀ=0.8, 有效元法与有限差分方法的误差收敛率 Tb. Teendnveeneffnteelementmetd ndfntedfeenemetdwenτ= 000 ἀ=0.8 有限元法有限差分方法误差收敛率 /% 误差收敛率 /%. w 0 w ' u n+ -u n+ un+ -u n+ 0, 则 w = w ' w 0. 由边值问题得 = x [ w -(u n+ - u n+ )], 故 w = w x [ w 0 + u n+ -u n+ 0], 则 w un+ -u n+ 0. 当 k

图5描述了式 ( 8)不同阶数α 的次扩散情形,可以看出随着分数阶导数阶数的增加,扩散速度在加快. π,采用 P ( 图固定 = K )的数值误差 00 π. Tenume le sbyus e fp ( K )wen= 00 图 4 α=0.7 时式( 8)的数值解.4 Tenume ls l u t n f( 8)wenα=0.

5 第期吴春红等:反常次扩散问题的有限元逼近长的关系图. 可以看出,针对不同的阶数 α,格式对时间步长具有一阶精度,对空间步长分别具有二阶和三阶精度. 69 u( x, 0)= 6x, 0 x 0.5,.5 x, -x +3x +, 0 u( 0, t)=u(, t)=0, t>0 的问题: { (,) u( x, t) -α u x t, 0 x, t >0. =Dt t x ( 8) 图4 描述了α=0.7 时式( 8)解的次扩散情况,可以看出,随着时间的增加,源在扩散. 图5描述了式 ( 8)不同阶数α 的次扩散情形,可以看出随着分数阶导数阶数的增加,扩散速度在加快. π,采用 P ( 图固定 = K )的数值误差 00 π. Tenume le sbyus e fp ( K )wen= 00 图 4 α=0.7 时式( 8)的数值解.4 Tenume ls l u t n f( 8)wenα=0.7,采用 P ( 图固定τ= K )的数值误差 000. Tenume le sbyus e fp ( K )wenτ= 000 图 5 α [ 0, ], t=.0 时式( 8)的数值解.5 Tenume ls l u t n f( 8)wenα [ 0, ] 5 结,采用 P ( 图 3 固定τ= K )的数值误差 Tenume le sbyus e fp ( K )wenτ= 000 例考虑带有初边值条件: 论本文针对分数阶导数施加于右侧扩散项的次扩散问题,采用有限元法分别给出了其时间半离散时间空间全离散形式,并且讨论了两种形式的稳定性收敛性,数值例子中与有限差分方法进行比较,得到

6 70 厦门大学学报 ( 自然科学版 ) 04 年了较好的结果. 通过本文的研究, 不仅探索了次扩散问题有限元方法的处理方式和技巧, 同时为求解其他类型的次扩散问题提供了思路. 参考文献 : [] Ln Yumn,Xu Cunju.nte dfeene/spetlppxmtnsftetme-ftnldfusn euqtn [J].JCmputPys,007,5(): [] LXnjun,XuCunju.A spe-tmespetlmetd ftetmeftnldfusnequtn[j].siamjnumeanl,009,47(3):08-3. [3] Zen Yunyn,LCnpn,ZZenn.A nten te fnte element metd f te spe-ftnl dvetndfusnequtn[j].cmputesnd MtemtswtAppltns,00,59(5): [4] Lu,Yn C,Bue K.Numel metd nd nlytl tenque f te mdfed nmlus subdfusnequtnwtnnlnesuetem[j].j CmputApplMt,009,3(): [5] Zun P,Lu, An V,et l.new slutn nd nlytltenquesftempltnumelmetdf tenmlussubdfusnequtn[j].siam JNume Anl,008,46(): [6] Yu Q,Lu,An V,etl.Slvnlnendnn-lne spe-tmeftnletn-dfusnequtnsbyte Admndempstn metd[j].intentnljunl f Numel Metdsn Ennen,008,74(): [7] 郭鹏, 王艺红, 陶春兴, 等. 时间分数阶 Klen-Gdn 型方程的解析近似解 [J]. 科技导报,009,7():47-5. [8] QutenA,VlA.Numelppxmtnfptl dfeentlequtns[m].new Yk:Spne,997. nteelementappxmtnfteanmlussub-dffusnpess WU Cun-n,LIU Qnḡx (.SlfAppled Mtemts,XmenUnvestyfTenly,Xmen3604,Cn;.SlfMtemtlSene,XmenUnvesty,Xmen36005,Cn) Abstt: Reentlyftnldfusnequtnsewdelyusedtdesbenmlusdfusnpesses,tenteesefdffusnpessesplysnmptntlenmnyfeldssusenneen,pyss,et.Intsppe,wensdesub-dfusnequtnbyfnteelementmetd.Tesem-dseteppxmtnndfuldseteppxmtneppsed.Andtestbltynd nveeneedsussed.nlysmenumelexmplesepesentedtdemnstteteefetvenessfteetlnlyss. Keywds: nmlussub-dfusnpess ;fnteelementmetd;stblty;nveene

标题

第 35 卷第期西南大学学报 ( 自然科学版 ) 3 年月 Vol.35 No. JouralofSouthwestUiversity (NaturalScieceEditio) Feb. 3 文章编号 :673 9868(3) 69 4 一类积分型 Meyer-KiḡZeler-Bzier 算子的点态逼近赵晓娣, 孙渭滨宁夏大学数学计算机学院, 银川 75 摘要 : 应用一阶 DitziaṉTotik