Pure Mathematics
Vol. 09  No. 10 ( 2019 ), Article ID: 33331 , 9 pages
10.12677/PM.2019.910140

The Quadratic Finite Volume Element Method for a Parameter Identification Problem of Parabolic Differential Equation

Juan Ma, Ling Li, Zhiguang Xiong

School of Mathematics and Computational Science, Hunan University of Science and Technology, Xiangtan Hunan

Received: Nov. 16th, 2019; accepted: Nov. 29th, 2019; published: Dec. 6th, 2019

ABSTRACT

In order to solve the problem of parameter identification for a class of parabolic differential equations with additional conditions, the quadratic finite volume element method is proposed for semi-discretization in space and the quadratic continuous finite element method for complete discretization in time. The numerical solution of unknown function and control parameter stability is derived, and the error analysis of the corresponding format is given. Finally, a numerical example is given to verify the stability and effectiveness of the proposed scheme.

Keywords:Parameter Identification Problem, Quadratic Finite Volume Element, Quadratic Continuous Finite Element, Totally Discrete Scheme, Stability

一类抛物微分方程参数识别问题的二次有限体积元方法

马娟,李玲,熊之光

湖南科技大学,数学与计算科学学院,湖南 湘潭

收稿日期:2019年11月16日;录用日期:2019年11月29日;发布日期:2019年12月6日

摘 要

针对一个具有附加条件的一类抛物微分方程的参数识别问题,首先提出了在空间上采用二次有限体积元法进行半离散,在时间方向上进行二次连续有限元全离散,导出了未知函数和控制参数稳定的数值解法,并给出了相应格式的误差分析,最后给出一个数值例子验证了所研究计算格式的稳定性和有效性。

关键词 :参数识别问题,二次有限体积元,二次连续有限元,全离散格式,稳定性

Copyright © 2019 by author(s) and Hans Publishers Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY).

http://creativecommons.org/licenses/by/4.0/

1. 前言

抛物微分方程参数识别问题是微分方程反问题中的一种,反问题是当微分方程中除了未知解以外,还有某个系数或初始条件或边界条件也是未知的。

本文将研究如下偏微分方程的参数识别问题

{ u ( x , t ) t 2 u ( x , t ) x 2 = p ( t ) u ( x , t ) , x ( a , b ) , t > 0 , u ( a , t ) = 0 , u ( b , t ) = 0 , t 0 , u ( x , 0 ) = f ( x ) , x [ a , b ] , u ( x * , t ) = E ( t ) , a < x * < b , t 0. (1)

其中 f ( x ) , E ( t ) 为充分光滑已知函数, p ( t ) , u ( x , t ) 为未知函数。该问题可以描述热源的扩散、粒子振动、控制理论等许多物理现象,文 [1] [2] 利用Adomain分解法研究了该方程的解,文 [3] 中用几种有限差分法求解该方程,讨论并比较了它们的准确性和稳定性,文 [4] 中使用了二次拟插值法解决抛物型方程的参数识别问题,文 [5] [6] 对抛物型参数识别问题给出有限体积元方法的求解格式,本文将用二次有限体积元法求解该问题。本文将作如下安排:第2节给出空间上的二次有限体积元的半离散格式。第3节给出时间上的二次连续有限元的全离散格式。第4节给出全离散格式的误差分析,讨论其稳定性和收敛性。第5节给出数值例子来说明该方法的有效性。

2. 有限体积元格式

首先对问题(1)作如下函数变换

r ( t ) = e 0 t p ( s ) d s , w ( x , t ) = r ( t ) u ( x , t ) , (2)

则问题(1)变换如下:

{ w ( x , t ) t 2 w ( x , t ) x 2 = 0 , x ( a , b ) , t > 0 , w ( a , t ) = 0 , w ( b , t ) = 0 , t 0 , w ( x , 0 ) = f ( x ) , x [ a , b ] , r ( t ) = w ( x * , t ) / E ( t ) , a < x * < b , t 0. (3)

该定解问题中方程不含参数函数,比原问题要简单,而且前三式直接为关于 w ( x , t ) r ( t ) 的正问题,继而就可以求出参数 p ( t )

下面给出空间上二次有限体积元离散的计算格式,对空间区间 I = [ a , b ] 作如下剖分:

设单元 I j = [ x j 1 , x j ] ,步长 h j = x j x j 1 , h = max 1 j n h j ,单元中点为 x j 1 / 2 = ( x j + x j 1 ) / 2 ,单元内四分点为 x j ± 1 / 4 = x j ± h j / 4 ,并要求该剖分 J h 是拟一致剖分,即存在常数 C > 0 使得 h C h j ,此外再设

I 0 = [ x 0 , x 1 / 4 ] , I 2 i + 1 = [ x i + 1 / 4 , x i + 3 / 4 ] ( i = 0 , 1 , 2 , , n 1 ) , I 2 i = [ x i 1 / 4 , x i + 1 / 4 ] ( i = 1 , 2 , , n 1 ) , I 2 n = [ x n 1 / 4 , x n ] ,

那么所有的 I i ( i = 0 , 1 , , 2 n ) 构成 J h 的一个对偶剖分 J h * ,每个 I i ( i = 0 , 1 , , 2 n ) 为体积控制元。设试探函数空间为 U h = { u | u C ( I ) , u I j } ,对每个t, w ( x , t ) 在单元 I j 上的基函数为:

{ φ j 1 ( x ) = 2 h j 2 ( x j x ) ( x j 1 / 2 x ) , φ j 1 / 2 ( x ) = 4 h j 2 ( x j x ) ( x x j 1 ) , φ j ( x ) = 2 h j 2 ( x x j 1 / 2 ) ( x x j 1 ) ,

则定义插值算子 Π 2 : C ( I ) U h 使得对任意的 w ( x , t ) U h 表示为

Π 2 w h = j = 0 n w j ( x ) φ j ( x ) + j = 1 n w j 1 / 2 ( x ) φ j 1 / 2 ( x ) .

再设函数空间 V h 为分片常数空间,对每个t, w ( x , t ) 在基于节点 x j , x j + 1 / 2 上的对偶单元 I i * 的基函数为:

ψ j ( x ) = { 1 , x I 2 j , 0 , , j = 0 , 1 , 2 , , n , ψ j + 1 / 2 ( x ) = { 1 , x I 2 j + 1 , 0 , , j = 0 , 1 , 2 , , n 1 ,

则定义插值算子 Π 2 * : C ( I ) V h ,使得对任意的 w ( x , t ) V h 表示为

Π 2 * w h = j = 0 j = n w j ( t ) ψ j ( x ) + j = 1 j = n w j + 1 / 2 ( t ) ψ j + 1 / 2 ( x ) ,

即对任意 w h U h ,在单元 I j 可以表示为

w h = w j 1 ( t ) φ j 1 ( x ) + w j 1 / 2 ( t ) φ j 1 / 2 ( x ) + w j ( t ) φ j ( x ) ,

其中 w j ( t ) = w h ( x j , t ) , w j 1 / 2 ( t ) = w h ( x j 1 / 2 , t ) ,且 w 0 ( t ) = w h ( a , t ) , w n ( t ) = w h ( b , t )

则(3)的二次有限体积元半离散为:固定变量t,寻求 w h ( t ) U h ,使得,

{ I i * w h ( x , t ) t d x I i * 2 w h ( x , t ) x 2 d x = 0 , i = 1 , , 2 n 1 , w 0 ( t ) = w n ( t ) = 0 , w h ( 0 ) = Π 2 f ( x ) U h ,

{ x i + 1 / 4 x i + 3 / 4 w h ( x , t ) t d x x i + 1 / 4 x i + 3 / 4 2 w h ( x , t ) x 2 d x = 0 , i = 0 , 1 , , n 1 , x i 1 / 4 x i + 1 / 4 w h ( x , t ) t d x x i 1 / 4 x i + 1 / 4 2 w h ( x , t ) x 2 d x = 0 , i = 1 , , n 1 , w 0 ( t ) = w n ( t ) = 0 , w h ( 0 ) = Π 2 f ( x ) U h ,

引入如下记号:

w ( t ) = ( w 1 / 2 ( t ) w 1 ( t ) w 3 / 2 ( t ) w 2 ( t ) w n 1 ( t ) w n 1 / 2 ( t ) ) T ,

α 0 = ( f ( x 1 / 2 ) f ( x 1 ) f ( x 3 / 2 ) f ( x 2 ) f ( x n 1 ) f ( x n 1 / 2 ) ) T ,

则可以得到二次体积元半离散的形式,用矩阵表示为如下满足初始条件的一阶线性常微分方程组

{ A w ˙ ( t ) = B w ( t ) , w ( 0 ) = α 0 . (4)

其中

A = h 48 [ 22 1 0 0 0 0 0 0 0 0 5 16 5 1 0 0 0 0 0 0 0 1 22 1 0 0 0 0 0 0 0 1 5 16 5 1 0 0 0 0 0 0 0 1 22 1 0 0 0 0 0 0 0 1 5 16 0 0 0 0 0 0 0 0 0 0 1 22 1 0 0 0 0 0 0 0 1 5 16 5 0 0 0 0 0 0 0 0 1 22 ] ( 2 n 1 ) × ( 2 n 1 ) ,

B = 1 h [ 4 2 0 0 0 0 0 0 0 0 2 4 2 0 0 0 0 0 0 0 0 2 4 2 0 0 0 0 0 0 0 0 2 4 2 0 0 0 0 0 0 0 0 2 4 2 0 0 0 0 0 0 0 0 2 4 0 0 0 0 0 0 0 0 0 0 2 4 2 0 0 0 0 0 0 0 0 2 4 2 0 0 0 0 0 0 0 0 2 4 ] ( 2 n 1 ) × ( 2 n 1 ) .

3. 二次连续有限元全离散格式

将时间部分有界区间 J = [ 0 , T ] 作拟一致剖分: 0 = t 0 < t 1 < < t j < < t n 1 = T ,时间单元 J m = [ t m 1 , t m ] 的中点为 t m 1 / 2 = ( t m + t m 1 ) / 2 ,步长为 l m = t m t m 1 ,且记 l = max 1 m n 1 { l m / 2 } 。对任意 U = ( u 1 , u 2 , , u 2 n 1 ) T , V = ( v 1 , v 2 , , v 2 n 1 ) T ,为了以下研究方便,可以定义一个向量运算“ ”: U V = ( u 1 v 1 , u 2 v 2 , u 3 v 3 , , u 2 n 1 v 2 n 1 ) T ,其它运算同普通的矩阵和向量运算。

因为矩阵A是严格对角占优的,所以A是可逆的,记 A 1 B = G ,则(4)中 A w ˙ ( t ) = B w ( t ) 可化为:

w ˙ ( t ) = G w (t)

将上式两边与 Z ( t ) = ( z 1 , z 2 , z 3 , , z 2 n 1 ) T ( P 1 ( t ) ) 2 n 1 作“ ”运算,并在单元 J m 上积分,有

J m w ˙ ( t ) Z ( t ) d t = J m G w ( t ) Z ( t ) d t . (5)

定义 w ( t ) 在t方向上的二次连续有限元 D ( t ) 使:

J m D ˙ ( t ) Z ( t ) d t = J m G D ( t ) Z ( t ) d t . (6)

D ( t ) = B 0 m + ( t t m 1 ) B 1 m + ( t t m 1 ) 2 B 2 m , t [ t m 1 , t m ] ,其中 B 0 m , B 1 m , B 2 m R N 为待定常数向量,当 t = t m 1 时, B 0 m = D ( t m 1 ) = D m 1 为已知。我们可以分别取 Z ( t ) = ( 1 , 1 , 1 , , 1 ) 2 n 1 T ( P 1 ( t ) ) 2 n 1 Z ( t ) = ( t t m 1 , t t m 1 , t t m 1 , , t t m 1 ) 2 n 1 T ( P 1 ( t ) ) 2 n 1 ,即可得到关于 B 1 m B 2 m 的广义方程组:

{ l m B 1 m + l m 2 B 2 m = G [ l m D m 1 + 1 2 l m 2 B 1 m + 1 3 l m 3 B 2 m ] , 1 2 l m 2 B 1 m + 2 3 l m 3 B 2 m = G [ 1 2 l m 2 D m 1 + 1 3 l m 3 B 1 m + 1 4 l m 4 B 2 m ] ,

既得

{ B 1 m = [ I N 1 2 l m G + 1 12 ( l m G ) 2 ] - 1 [ I N 1 2 l m G ] G D m 1 , B 2 m = 1 2 [ I N 1 2 l m G + 1 12 ( l m G ) 2 ] - 1 G 2 D m 1 .

其中 I N N = ( 2 n 1 ) 阶单位矩阵。

将上式求得的 B 1 m , B 2 m 代入 D ( t ) = B 0 + B 1 ( t t m 1 ) + B 2 ( t t m 1 ) 2 中,并取 t = t m 得到:

D m = D m 1 + l m B 1 m + l m 2 B 2 m = [ I N 1 2 l m G + 1 12 ( l m G ) 2 ] - 1 [ I N + 1 2 l m G + 1 12 ( l m G ) 2 ] D m 1 . (7)

由于 D 0 = w ( 0 ) = α 0 ,所以用此递推的方法就可以求出每个时间层的全离散数值解 D m 。此方法就是二次连续有限元法的全离散格式,是一种显示单步格式。

最后可以得到数值解 r m , u j m , p m 的计算格式,由(2)和(3)式可知:

r ( t ) = w ( x * , t ) E ( t ) , u ( x , t ) = w ( x , t ) r ( t ) , p ( t ) = r ( t ) r ( t ) .

我们已由式(7)得到 w ( x * , t ) 的数值解 w k m ,由此可得以下数值解 r m , u j m , p m 的计算格式:

r m = w k m E ( t m ) , u j m = w j m r m , p m = r m + 1 r m 1 2 r m l m , (8)

其中 r m , u j m , w j m , p m 分别对应为 r ( t m ) , u ( x j , t m ) , w ( x j , t m ) , p ( t m ) 的数值解。

4. 误差分析

下面对空间上的二次有限体积元半离散解和全离散格式进行误差分析,根据文 [7] [8] 我们有如下引理:

引理4.1 设 w ( x , t ) 为式(3)的精确解, w h ( x , t ) 是式(4)的二次有限体积元半离散解,则有如下模误差估计:

w ( x j , t ) w h ( x j , t ) C h 3 .

引理4.2 w h ( x j , t ) 是式(4)的二次有限体积元半离散解, w j m 是式(6)的二次有限元全离散解,则有如下误差估计:

w h ( x j , t k ) w j m C l 4 .

引理4.3 设 w ( x , t ) 为式(3)的精确解, w j m 是式(6)的二次有限元全离散解则有如下误差估计:

w ( x j , t k ) w j m C h 3 + C l 4 .

定理4.1设 u ( x , t ) 为式(1)的精确解, u j m 是由式(8)得到的二次有限元全离散解, p ( t m ) , p m 分别为参数 p ( t ) 的精确解与数值解, r ( t m ) , r m 分别为参数 r ( t ) 的精确解与数值解,则有如下误差估计:

r ( t m ) r m C 1 h 3 + C 1 l 4 ,

u ( x j , t m ) u j m C 2 h 3 + C 2 l 4 ,

p ( t m ) p m c 1 h 3 l + c 2 l 2 .

证明 r ( t ) 的误差估计:

r ( t m ) r m = 1 E ( t m ) ( w ( x * , t m ) w k j ) 1 E ( t m ) ( w ( x * , t m ) w h ( x * , t m ) + w h ( x * , t m ) w k j ) M ( C h 3 + C l 4 ) C 1 h 3 + C 1 l 4 .

u ( x , t ) 的误差估计:由于 r ( t m ) r m C 1 h 3 + C 1 l 4 < r ( t m ) / 2

r ( t ) = e 0 t p ( s ) d s > 0 ,

则有

r m = r m r ( t m ) + r ( t m ) r ( t m ) r ( t m ) r m 1 2 r ( t m ) .

从而有

u ( x j , t m ) u j m = w ( x j , t m ) r ( t m ) w k m r m = w ( x j , t m ) ( r m r ( t m ) ) + r ( t m ) ( w ( x j , t m ) w k m ) r ( t m ) r m w ( x j , t m ) r ( t m ) r m ( C 1 h 3 + C 1 l 4 ) + 1 r m ( C h 3 + C l 4 ) ( M 2 w ( x j , t m ) r 2 ( t m ) + 2 r ( t m ) ) ( C h 3 + C l 4 ) C 2 h 3 + C 2 l 4

p ( t ) 的误差分析:由Taylor展开式得

r ( t m ) = r ( t m + 1 ) r ( t m 1 ) 2 l m + o ( l m 2 )

从而有:

p ( t m ) p m = r ( t m + 1 ) r ( t m 1 ) 2 r ( t m ) l m o ( l m 2 ) r ( t m ) + r m + 1 r m 1 2 r m l m = 1 2 l m r ( t m ) ( r ( t m + 1 ) r m + 1 ) + ( r ( t m + 1 ) r ( t m 1 ) ) ( r ( t m ) r m ) + r ( t m ) ( r ( t m 1 ) r m 1 ) r ( t m ) r m + o ( l m 2 ) r ( t m ) 1 l m ( 2 r ( t m ) ( C 1 h 3 + C 1 l 4 ) + r ( t m + 1 ) r ( t m 1 ) r 2 ( t m ) ( C 1 h 3 + C 1 l 4 ) ) + o ( l m 2 ) r ( t m ) c 1 h 3 l + c 2 l 2

5. 数值例子

本小节求解如下参数识别问题的数值例子,通过误差数据来说明计算格式的有效性。

{ u ( x , t ) t 2 u ( x , t ) x 2 = p ( t ) u ( x , t ) , x [ 0 , 1 ] , t > 0 , u ( 0 , t ) = 0 , u ( 1 , t ) = 0 , t 0 , u ( x , 0 ) = sin ( π x ) , x [ 0 , 1 ] , u ( 0.25 , t ) = 2 2 e t 2 , t 0.

已知问题的精确解为: u ( x , t ) = e t 2 sin ( π x ) , p ( t ) = π 2 2 t

先作代换消去 p ( t ) 后,在空间上用二次有限体积元法,在时间上用二次连续有限元全离散进行数值计算,取均匀空间步长为h和均匀时间步长 Δ t 表1给出了当 t = 1 时,网格比为 Δ t / h = 1 的节点在 N = 100 , 200 , 400 , 800 的绝对误差值。表2给出了当 t = 1 时,网格比为 Δ t / h 2 = 10 2 的节点在 N = 100 , 200 , 400 , 800 的绝对误差值。

Table 1. Absolute errors of u ( x , 1 ) , p ( t ) when the mesh ratio is Δ t / h = 1

表1. 网格比为 Δ t / h = 1 u ( x , 1 ) , p ( t ) 的绝对误差

Table 2. Absolute errors of u ( x , 1 ) , p ( t ) when the mesh ratio is Δ t / h = 10 2

表2. 网格比为 Δ t / h = 10 2 u ( x , 1 ) , p ( t ) 的绝对误差

表1表2可以看出用本章研究的数值方法进行抛物方程参数识别是有效而可靠的。它的精度不仅与剖分密度有关,还与时间和空间的相对密度有关,也就是说,当剖分密度越大,时间和空间的相对密度越大,则利用该方法的精度越高。

基金项目

本文为国家自然科学基金项目(11571102)资助课题。

文章引用

马 娟,李 玲,熊之光. 一类抛物微分方程参数识别问题的二次有限体积元方法
The Quadratic Finite Volume Element Method for a Parameter Identification Problem of Parabolic Differential Equation[J]. 理论数学, 2019, 09(10): 1139-1147. https://doi.org/10.12677/PM.2019.910140

参考文献

  1. 1. Tatari, M. and Dehghan, M. (2006) Identifying a Control Function in Parabolic Partial Differential Equations from Overspecified Boundary Data. Computers and Mathematics with Applications, 53, 1933-1942. https://doi.org/10.1016/j.camwa.2006.01.018

  2. 2. Tatari, M., Dehghan, M. and Razzaghi, M. (2006) Determina-tion of a Time-Dependent Parameter in a One-Dimensional Quasi-Linear Parabolic Equation with Temperature Overspecification. International Journal of Computer Mathematics, 83, 905-913. https://doi.org/10.1080/00207160601117420

  3. 3. Dehghan, M. (2004) Parameter Determination in a Partial Differential Equation from the Overspecified Data. Mathematical and Computer Modelling, 41, 196-213. https://doi.org/10.1016/j.mcm.2004.07.010

  4. 4. Ma, L.-M. and Wu, Z.-M. (2010) Identifying the Temperature Distribution in a Parabolic Equation with Overspecified Data Using a Multiquadric Quasi-Interpolation Method. Chinese Physics B, 19, 1-6. https://doi.org/10.1088/1674-1056/19/1/010201

  5. 5. Xiong, Z.G., Deng, K. and Liu, Z.S. (2015) The Finite Volume Element Method for a Parameter Identification Problem. Journal of Ambient Intelligence and Humanized Computing, 6, 533-539. https://doi.org/10.1007/s12652-014-0238-7

  6. 6. 陈国荣, 王雪玲, 熊之光. 一类参数识别问题的有限体积元计算[J]. 衡阳师范学院学报, 2011, 32(3): 32-35.

  7. 7. 熊之光, 刘晓奇, 邓康. 抛物方程初边值问题连续有限元的超收敛性[J]. 数学的实践与认识, 2007, 37(11): 141-147.

  8. 8. Xiong, Z.G. and Deng, K. (2017) A Quadratic Triangular Finite Volume Element Method for a Semilinear Elliptic Equation. Advances in Applied Mathematics and Mechanics, 9, 186-204. https://doi.org/10.4208/aamm.2014.m63

期刊菜单