（应用数学专业论文）带有第三类边界条件的热传导方程逆时问题数值解.pdf

摘要 在很多工程应用领域里，人们常常遇到热传导方程的逆时问题，即由某一时 刻的温度场去确定该时刻以前的温度分布由于这类问题通常是不适定的，因此 求解此问题需要利用各种正则化方法。处理该类问题的经典方法是把原问题转化 为求解第一类积分方程的问题，再利用正则化方法求解。众所周知，正则化方法 中的难点之一是正则化参数的选取由于正向热传导问题温度场随时间是指数衰 减的，因此逆向热传导问题中对正则化参数的选取的要求更为严格 本论文研究一维空间的具有第三类边界条件的热传导方程 甓= 鑫( 女( z ) 舞) ， ( 。，t ) ( n ，6 ) ( o ，卅 u x ( a ，t ) 一 ( n ) u ( o ，t ) = 0 u 。( b ，t ) + h ( b ) u ( b ，t ) = 0 u ( z ，0 ) = g ( 。) ，z ( 。，6 ) 的逆时问题数值求解。本文分别研究了两类问题：一是由某一时刻的温度u ( 。，t ) 来求出初始温度9 ( z ) ，将求解第一类积分方程的最小模解的方法用于确定正则化 参数，并据此对逆时热传导问题进行了数值试验，取得了令人满意的数值结果 我们的数值结果表明，借助于正则化参数的适当选取，即使初始温度场的振荡性 较强或光滑性较弱，也能获得相对而言比较好的重建结果另一类是由某一时刻 的温度u ( z ，t ) 来求出初始温度u ( x ，t o ) ，t o ( o ，t ) ，利用对数凸性估计建立解的条 件稳定性，在g ( z ) 的先验条件下给出正则化参数的选取标准及正则化解的收敛 速度 本文的结构如下第一章介绍了数学物理反问题的背景。第二章通过特征函 数和特征值构造积分算子方程a g ( x ) = u ( x ，t ) ，并且说明该方程是不适定的。第 三章乖j 用t i k h o n o v 正则化方法求解积分方程，并讨论正则化参数的选取问题 第四章是数值实验，三个不同的例子验证了提出的方法的数值有效性 关键词：热传导，逆时问题，t i k h o n o v 正则化，正则化参数，r o b i n 边界 a b s t r a c t t h eb a c k w a r dh e a tc o n d u c t i o np r o b l e mi sat y p i c a li n v e r s ep r o b l e mw h i c ha i m st o d e t e r m i n et h ei n i t i a lt e m p e r a t u r ed i s t r i b u t i o nf r o mt h ek n o w nt e m p e r a t u r eg i v e na t s o m e t i m e t 0 t h i s p r o b l 一 ，, e m i ss e v e r ei l lp o s e dd u et ot h e e x p o n e n t i a l l yd ecreasingof t e m p e r a t u r e t h e r e f o r et h er e g u l a r i z i n gp a r a m e t e rs h o u l db ec h o s e nc a r e f u l l y t h e p u r p o s eo ft h i st h e s i si st os t u d yt h e1 - db a c k w a r dh e a tc o n d u c t i o np r o b l e m w i t hr o b i nb o u n d a r yc o n d i t i o nf r o mt h ef o l l o w i n gm o d e lp r o b l e m o o ( z ，t ) ( o ，6 ) ( 0 ，司 x ( a ，b ) w ec o n s i d e rt h ed e t e r m i n a t i o no ft h et e m p e r a t u r ed i s t r i b u t i o nu ，) ，0st 0 t w o m e t h o df o rt h ec h o i c eo ft h er e g u l a r i z i n gp a r a m e t e ra r eg i v e n t h ef i r s tm e t h o dd e - p e n d so nt h em o r o z o vd i s c r e p a n c yp r i n c i p l ef o rw h i c h w e g i v ea ni t e r a t i v em e t h o dt o d e t e r m i n et h er e g u l a r i z i n gp a r a m e t e r ，w h e nw es o l v et h ei n i t i a lv a l u eu ( x ，0 ) = 9 ( z ) f r o mu ( o ，t ) t h es e c o n do n ei st od e t e r m i n et h er e g u l a r i z i n gp a r a m e t e rf r o mt h e c o n d i t i o n a ls t a b i f i t yo fs o l u t i o nw h e nw er e c o v e ru ( z ，t o ) f r o mu ( z ，t ) f o r0 o o ( t t h e c o n v e r g e n c er a t eo ft h er e g u l a r i z i n gs o l u t i o ni s a l s oa n a l y z e du n d e rs o m en p r i o r i a s s u m p t i o no n9 ( z ) t h eo u t l i n eo ft h i st h e s i si sa sf o l l o w s w ef i r s t l yg i v et h eb a c k g r o u n d so ni n - v e r s ep r o b l e m s ，e s p e c i a l l yo nt h eb a c k w a r dc o n d u c t i o np r o b l e mi nc h a p t e r1 t h e n a ni l l p o s e d i n t e g r a le q u a t i o na 9 ( o ) = n ( z ，t ) f o rt h eb a c k w a r dh e a tc o n d u c t i o ni s c o n s t r u c t e di nc h a p t e r2 ，i nc h a p e r3 ，w ec o n s i d e rt h em e t h o d so fd e t e r m i n i n gr e g u - l a r i z a t i o np a r a m e t e rb yc o n d i t i o n a ls t a b i l i t ya n dd e s c r e p a n c yp r i n c i p l e ，w h e nw es o l v e 乱( z ，t o ) f o r0 0 时的温度场u ( 。，t ) 是已知的，要确定在0 到t 之间的温度场“( z ，t ) 。这 就是所说的热传导方程逆时问题由于这类问题是不适定的，因此求解要比正问题复杂 的多。 l 恤h孔嚣m 叫 i忡删艄舢旷卅归 2 逆时热传导问题的难点( 问题的强不适定性) 来源于正向热传导问题温度场随时间的 指数衰减性质时间越长，初始温度的特性就越难体现，从而由测量数据反演初始温度 难度就越大。该特性可以从简单的模型问题中看出在上述长杆温度的数学模型中，取 ( ，b ) = ( 0 ，7 r ) ，( o ) = 1 ，h ( o ) = 0 ，h ( ) = 1 ，初始温度g ( z ) = c o s ( 5 z ) 正问题温度u ( x ，t ) 在 不同时刻的分布见图1 1 从图中可以看出初始温度的较强振荡性，经过一秒后已很难 体现这说明对于的热传导方程逆时问题，使用的观测时刻相对说来不能取的太大 1 2 逆时热传导的已有工作 近年来，这方面的文章有很多在f 1 1 中，r c h a p k o 利用正则n e w t o n 法讨论了二 维热传导方程边界的重构问题h ，h a n 2 1 利用最小能量技巧改进了边界元方法，解决了 一维热传导方程逆时问题，得到较为满意的结果。c h a p k o 和h a n 均考虑第一类边界条 件的热传导方程，而在一般情况下考虑阻尼边界条件则更为实际1 9 9 9 年，m u n i z 3 1 则 分别用三种方法求解一维具有第一类边界条件的热传导方程逆时问题，第一种是先把原 问题转化为第一类积分方程，再直接求解积分方程，效果很差；第二种是用向后欧拉差 分隐格式，在输入数据无扰动时效果不错，但输入数据有小扰动时则很不理想；第三种 是用t i k h o n o v 正则化方法，得到较满意的结果刘继军教授【4 ，5 ，1 5 】同样用t i k h o n o v 3 正则化方法解决了二维具有第一类边界条件的热传导方程求初始温度的问题，并给出解 的条件稳定性和收敛速度。由于使用t i k h o n o v 正则化方法时，会涉及正则化参数的选取 问题，【3 和【5 】均未详细讨论此问题，m u n i z 只有取不同的正则化参数进行数值实验， 比较结果，找出较好的正则化参数，这样很费时费力对于正则化参数的选取，近年来 也有不少文献出版。程晋教授6 1 在2 0 0 0 年提出的先建立解的某种条件稳定性，再由条 件稳定性确定正则化参数的方法得到了广泛的应用正则化解收敛到精确解的速度依赖 于该条件稳定性。j i a n l ix i e 在【8 中利用双参数的策略决定正则化参数，并且对椭圆型 方程进行数值实验，结果令人满意。肖庭延在【7 】中阐述了决定正则化参数的高阶收敛算 法，本文的数值方法利用其中提出的高阶迭代格式进行数值计算，得到了满意的结果。 1 3 本文的工作 本文探讨了一维具有第三类边界条件的热传导方程求初值的问题。我们的目标是通 过在某个时刻t 的已知的温度场“( 蜀t ) 来确定 f 镑= 鑫( ( 。) 爱) ( z ，t ) ( 。，6 ) ( o ，t 】 “z ( 。，) 一九( 。) “( 。t ) 。o o l ( 。) o 1 ) lu 。( 6 ，t ) + h ( b ) u ( b ，t ) = 00 t s zh ( b ) 0 、。 【u ( z ，0 ) = 9 ( z )o t o tz ( a ，b ) 中的温度场“( z ，t ) ，0 t t 。 文章主要内容安排如下第二章首先利用( 1 3 1 ) 的特征函数和特征值来构造积分算 子a ：x y ： ( a 9 ) ( 。) = k ( x ，y ) g ( y ) d y = u ( z ，t ) 这样把( 1 3 。1 ) 的逆时问题转化为由算予方程求9 ( 。) x 的阅题由于a 的核函数k ( 。，y ) 非常光滑，该问题是不适定的。该分析同时也揭示了由u ( z ，t ) 确定u ( 。，t o ) ，t o ( 0 ，t ) 的 不适定性第三章利用t i k h o n o v 正则化方法由u ( z ，t ) 求解u ( z ，t ) ，t 【o ，t ) ，在第三章 第二节里，我们由u ( z ，t ) 来反演初值9 ( z ) ，并对正则化参数n 采用后验取法，即利用 m o r o z o v 偏差原理决定a ，引入了计算a 的迭代格式并证明它是至少三阶收敛的。在第 三章第三节里，我们由u ( x ，t ) 来反演“( z ，t o ) ，0 0 ，。( a ，b ) 则( 1 3 1 ) 可转化为下方程： f 嚣= 一是( 女( z ) 罄) f ( n ，t ) 一h ( n ) w ( 0 1 t ) = 0 i 毗( 6 ，t 7 ) + h ( b ) w ( b ，t 7 ) = 0 【l r = 。= “( 。，t ) 这样原问题由u 功求( 1 3 1 ) 中的初值9 ( z ) ，转化为由( 2 2 2 ) 来求 ( z ，t ) = ”( 。，0 ) 9 ( 。) 。( 2 2 1 ) 可写为： a 一1 ”扛，0 ) = 0 ，t ) 下面的结果揭示了反问题的不稳定性为表达式的方便，记 m = 2 卜) 坼) 以则) 州啪( 6 ) 以邶) + z 6m ) ( o w 。( x 。o ) ) 2 出| 定理2 2 1 ：( 2 2 2 ) 的解”( z ，t ) 满足下面的估计： ，哟慨蚰刊艮h 圹印( 丽) 怛 冈 吣 c 。 “ 罗鼍0吣吣d 6 址：令 日( t ，) - m ，删。) 2 上，t 引2 d z ( 2 。2 3 ) 利用分部积分和方程( 2 2 ，2 ) 可得 掣= 。z 6 出 = 2 陆( b ) 砌) 以”，) 十凇a ) 舻( 叫，) 十z 6 蛐) ( 警) 2 如】 m 4 ) 箸笋= a 蚋( b , t r m 哟m “叫，) + f m 蚓盎 同样对上式右端最后一项分部积分得； a 2 h 矿( t ) o t i = 4z 6 罅如 ( 2 删 2 j d 1 、1 这样我们有： 愀妁卅f ) ) 2 = a f 上6 a 。小如一( r w w t , d x ) h ( t 1 1 j 2w wd 2 卜0 ) 日“( t ) 一( ( t ，) ) 2 = 4 i 出f ”如一( t。i ， f oj o oa l 最后一个不等式由s c h w a r z s 不等式得到。 该估计表明( 1 0 9 h ( t ) ) ”0 ，即l o g ( h ( t ，) ) 是凸函数，从而有下式成立： 螋学兰陋 t = 0 ，即抡即) e x p 【需2 t 十7 一l、j + 、一 、7 l【u ) i 由日的表达式( 2 2 3 ) ，即得所要证明的不等式，证毕。 令t ，：t ，并且注意到 ( z ，0 ) = u ( z ，t ) ，u ( z ，t ) = 9 ( z ) ，上述不等式可化为： 怕( ) 悒孙( t 川2 e x p ( 斋) ( 2 2 - 6 ) 丽：) k ( 。) 。( 啦) + h ( 6 ) 女( b ) u z ( b ，丁) + 上6k ( z ) ( o u 瓦( x t ) ) 2 d z 从( 2 2 6 ) 式看出，当i i u ( ，t ) n f ：一0 时，怕( ) 伲z o 。成立例如，我们取u ( z ，习2 寺s 饥( z ) ，自0 ) = l ，( n ，6 ) = ( o ，7 r ) ，则由( 2 ，2 6 ) 得： m ) 慨咿) 赤e x p 【c n 2 】 随着j 的增大， ( ，t ) u c z = 丢杀一o ，1 9 ( ) i i i ：一o 。 这说明当精确解u ( z ，t ) 的近似解“6 ( 。，t ) 满足： i i u 6 ( ，t 1 ) 一u ( ，t ) l l l 。6 ， 7 这就是说由a 一1 ( z ，t ) 来确定矿( z ) 是不稳定的 在实际应用中，我们只能通过测量来得到u ( x ，t ) 的近似值“6 ( z ，t ) ，这样不可避免 的会有测量误差于是一个无法避免的问题就是如何通过u 6 ( z ，t ) 来求g ( z ) 的近似值 矿( 。) ，并且使得当6 0 时，有矿( z ) 一g ( z ) 本文即利用t i k h o n o v 正则化方法求此近 似解 第三章正则化方法及正则化参数的选取 3 1t i k h o n o v 正则化方法 设a ：x y 是紧算子，x ，y 均为h i l b e r t 空间。在本文中x ，y 均为l 2 ( n ，b ) 。考 虑方程 a g ( x ) = ( 。，t )( 3 1 ，1 ) 对有误差的右端项( z ，t ) ，假定扰动数据( 。，丁) 满足j j “6 ( ，t ) 一u ( ，t ) ij 6 下面我们由( z ，t ) 来求9 ( z ) 的近似值矿( z ) 。从上一章的分析我们知道矿( z ) 不能 由a - 1 u 6 ( x ，t ) 直接确定。 为了保证矿( $ ) 对( 。，t ) 的连续依赖性，必须构造a 。的正则化解算子吼：y x 。有关正则化解算子的定义见9 1 所谓的t i k h o n o v 正则化方法就是求业( z ) ，使得下面t i k h o n o v 泛函值极小： 碟( 9 ) = i i a g u 6 1 1 2 。( 。，b ) 十a l l 9 i i 色( 。，6 ) 其中的n 称为正则化参数。 关于t i k h o n o v 泛函有如下标准结果( 【9 t h e o r e m2 1 1 ) 引理3 1 1 ：设x ，】，是h i l b e r t 空间，a ：x y 为有界线性算子，则 ( 1 ) 砭( 9 ) 在x 上存在唯一的极小元9 ：； ( 2 ) 业x 满足 叼：+ a + a 9 ：= a + 其中n o ，a + 为a 的伴随算子。 在前面我们已知a 是自伴算子，从而对我们的问题，( 3 1 3 ) 可写为 a 菇十a 2 9 ：= a 这样，我们不再对( 2 1 2 ) 求解，而是转向( 3 1 4 ) 。 在该方法中，正则化解算子和正则化解谚分别为 r d = ( o ，+ a 2 ) 一1 a ，g ：= r 。u 5 = ( a j + a 2 ) 一1 a u 6 由正则化解和精确解的经典误差分析( 业一9 1 1 = 1 i r 。u 6 9 | l i i r 。u 6 i l r 。l 札6 一叫i + l l r 。a 9 8 一咚u l + i i r 。u g l l 9 9 曼6 1 1 r a l l + i i 凡a g 9 1 1 ( 3 1 2 ) ( 3 1 3 ) ( 3 14 ) 9 可知，正则化解g ：与精确解g 的误差分为两项，第一项是由输入数据的误差6 0 而产 生的解的误差，它被正则化解算子的模i l r 。| | 放大了，对任何给定的6 0 ，当n 一0 时， 训r 。i i o 。；第二项表示正则化解算子逼近不适定算子a _ 1 在精确右端数据u ( z ，t ) 处产生的误差i i ( 如一a - 1 ) u 忆当a 一0 时，它趋于0 。一方面，近似解对输入数据的误 差d 0 的稳定性要求需要| | 凰| | 很小( 由正则化解算子的性质知：o 不能太小) ； 另一方面，正则化解算子凰对不适定算子a 。的逼近性要求i i ( r 。一a 。) u 1 1 很小( 即。 越小越好) 。这样就产生了一个问题：正则化参数a 如何选取，使得 刚r 圳十j i 风a g g | | 极小。 这样在利用正则化方法求解过程中，确定合适的正则化参数成了核心问题。其基本 要求是：对于待求解施加定性的或定量的信息，以使正则化参数与原始输入数据的误差 水平相匹配。 确定正则化参数a 通常有先验和后验两种策略先验策略是在求解( 3 1 4 ) 前，先取 a 的值。从正则化解算子如的定义可知，d 不会太大，应在区间( 0 ，1 ) 内在求解( 3 1 4 ) 前，在区间( 0 ，1 ) 中随意取值作为正则化参数a ，通过( 3 1 4 ) 求出正则化解g ：，如果9 ： 与精确解g 的估计范围之外，则另取a 值。这样不停的试，以期能求得g 的近似值。先 验策略便于理论分析，但一般说来关于精确解的假定条件难以验证，故数值求缌上采用 后验的方法居多。后验策略就是在解( 3 1 4 ) 过程中确定正则化参数o 本文使用两种方 法：( 1 ) 在由u ( 札t ) 反演初值9 ( 。) 时，利用m o r o z o v 偏差原理决定q ；( 2 ) 在由u ( 。，t ) 反演初值u ( x ，t o ) ，0 0 时，对任 意的y y ，方程a g = “在偏差6 下存在唯一的最小模解，并且最小模解连续地依赖于 右端项。我们也可以得到最小模解收敛速度的估计( 【9 1t h e o r e m2 1 7 ) 。 在实际计算中，对( 3 1 4 ) 求最小模解时，首先要通过( 3 , 2 2 ) 式确定n 的取值，在得 到。的确定值后，代回( 3 1 4 ) 式求解，这就是所谓的正则化参数。的后验策略。 一般情况下，总认为原始输入数据的大小比误差水平大得多，即1 1 i l 6 ，否则原 始输入数据就会完全被噪音覆盖，因而没有意义在这种情况下，可通过n e w t o n 方法来 数值求解( 3 2 2 ) 的根 本文利用一个迭代格式( 【7 】) 来确定n 的取值。通过简单计算可得： g 弘) ：一z ( 警，以) ，g ) ：一z ( 警，蘸) 一。an 警旷+ 辔) 1 ， 其中以，譬，搭分别由下列三式确定： ( 。a 2 ) 醒= a u d ，( a i + a 2 ) 警= 魂( a i + a 2 ) 辔_ - 2 警 【7 l 中提出的迭代格式如下 g ( a k ) g ( o k ) 瓯“2 吼一研面( 瓦丽 并指出该格式是局部三阶收敛的，但未给出详细证明过程事实上此迭代格式是局部至 少三阶收敛的。关于收敛阶定义如下： 定义：设序列如) 收敛于矿，并记e k = 矿一z k ，( k = 0 ，1 ，2 ) 如果存在常数p 1 及非零常数c ，使得： l i mt e k + l ：c 一。 则称序列 。k ) 是p 阶收敛的 在证明收敛性之前，先给出p 阶收敛定理( 【1 3 】，定理2 4 ) 作为引理 引理3 2 1 ：设迭代格式为z k + 1 = 妒( 钆) ， = 1 ，2 ，3 ) 若妒( z ) 在矿附近的某个邻 域内有v ( p 1 ) 阶导数，且 _ p ( ) ( 。+ ) = 0 ，( k = 1 ，2 ，3 ，一，p 一1 ) ，妒( p ) ( z + ) 0 则迭代格式在z + 附近是p 阶局部收敛的 显然a u 6 0 ，否则( 3 1 4 ) 只有零解我们可以知道g ( n ) 在区间( 0 ，1 ) 中有单根， g 7 ( 。) 0 ( 见1 9 1t h e o r e m2 1 6 ) 设矿为g ( a ) 在( o ，1 ) 内的单根，即g ( o + ) = 0 为近似 n 4 所使用的迭代程序为 令 堡! 塑型堕2 口( 。k ) 2 一 g ( 口k ) g ”( o k ) g ( 女) 1 丽f 哥 m 阳一器嚼1 原迭代格式变为：。k + 1 = f ( 口k ) 再记 如) = 铲，啦) 则上式化为f ( a ) = o s ( ) t ( a ) ，且 如) = 器 f 。( a ) = 1 一一( a ) t ( a ) 一s ( a ) t ( o ) f ”( a ) = 一8 1 1 ( o ) t ( o ) 一2 s ( n ) r ( o ) 一s ( a ) t ”( a ) 通过简单计算得 s 协) = l - t 弘) - _ 搿+ g ( 坝群一搿弘) = 杀 由于g 7 ( n + ) = 0 ，所以： _ o ! 以。卜篱 删乩以_ 1 ，s 一搿 邓卜1 ，t 协弘扣卜嬲 f ( o + ) = 1 一t ( n ) s 7 ( o + ) = 1 1 = 0 脚卜川a 忡卜2 s 铀f ( 。卜篇一2 器= o 同样我们也可以得到f ”，( a + ) ：塾塑盟；嚣警二2 1 兰盟，若无附加条件，则无法判断f ( a + ) 是否为零。这样由引理( 3 2 1 ) 知，上述迭代格式是至少三阶收敛的 在第四章数值实验中，通过上述迭代格式，我们先得到正则化参数n 值，进而得到 最小模解口! 。 3 3 利用条件稳定性确定正则化解 本节我们讨论由u ( z ，t ) 反演u ( z ，t o ) 的逆时问题，其中0 e 2 + 1 ，取口= 6 2 时，( 3 l4 ) 的正则 化解g ：满足下式： f 嘉( 以) 诺6 2 ，l l a 9 ：一a g l ls ( c o + 1 ) 6 证明当n = d 2 时，显然有： f 嘉b ) = i l a g ( ) 一“6 ( ，t ) 1 1 2 + 6 2 1 1 9 ( ) 1 1 2 = i i “( ，t ) 一“6 ( ，t ) 1 1 2 十6 2 1 1 9 ( ) 1 1 2 兰6 2 + 硎9 ( ) 1 1 2 。( e 2 + 1 ) 6 2 由于9 ：是使得t i k h o n o v 泛函垛( 9 ) 在l 2 ( n ，6 ) 上达到最小值，所以有 f 嘉( 蘸) 冬f 嘉( 9 ) ( e 2 十1 ) 6 2 c 0 占2 此式意味着 i i a g ：( ) 一u 6 ( ，t ) i i c b 6 ，i i g ：( ) l lsc o 这样可得： i l a g ：( ) 一a g ( ) 1 1 i i a g ：( ) 一u 6 ( ，t ) i i 十| | u 6 ( ，t ) 一a k ) ( c o + 1 ) 6 1 3 证毕。 当a = 6 2 时，可由( 3 1 4 ) 式得到正则化解9 ：( z ) ，它是初值g ( z ) 的近似解。这里我 们没有得到初值的正则化解g ：( z ) 的收敛速度但是把g ：( z ) 作为初值代入( 1 3 1 ) 中，求 解此正问题可以得到温度场( z ，t o ) ，t o ( o ，t ) ，对于这样的( 。，t o ) ，如果我们对精确 的初值g ( z ) 加一些限制，则由已建立的条件稳定性和定理3 3 3 可以得到对( 。，t o ) 的下 面的收敛速度估计。 定理3 3 4 设g ( x ) p ( e ) ，对于上述( z ，t o ) ，有下式成立： i i ( ，t o ) 一u ( ，t o ) 1 1 “- o m r 其中0 为正常数 证明对= 6 2 ，由定理3 3 3 知： j | u 6 ( ，o ) 1 l = | | g ：( 训c o l i u ( ，0 ) l | = 1 1 9 ( ) i | e c o 即9 i ( z ) ，9 ( x ) p ( c o ) 。由定理3 3 2 和定理3 3 3 可得 1 1 u 6 ( ，t o ) 一u ( - ，o ) 1 1 ( 2 c o ) l l a g ：( ) 一a 9 ( ) i i 争 墨( 2 c o ) ( c o + 1 ) 争6 挈 0 ，6 0 ，e 0 和最大迭代次数。，令k ：0 ( 2 ) ：离散并求解下列方程： g ( 。) = i i 且9 ：“1 1 2 搿，g 协) = 也( 警，或) 昨卜。( 豢，：) 一缸憎1 2 十辔) ( 4 ) ：利用迭代格式 啉一t 一丽鹪 求下一迭代点q + 1 ( 5 ) ：若l o z k + 1 一n k f e 或迭代次数女= k 。，则停止；否则，令k = + 1 ，回到( 2 ) 这里迭代初值o o = o 0 1 ，e = 1 0 一，。= 5 0 ，m = 1 0 0 ，t = 0 0 0 5 当u 6 = u 十o 0 0 1 s i n ( 2 x 1 ) 时，迭代5 次停止，+ = 9 5 1 5 1 4 5 3 1 0 4 相对误差数 值结果见图4 ，4 。 当= u + o o l s i n ( 2 x 1 ) 时，迭代2 次停止。o + = 8 5 5 4 1 0 4 6 1 0 一、，相对误差数 值结果见图4 5 。 从这两张图和由我们迭代得到的。数值结果可知，利用3 1 节中的迭代格式确定正 则化参数。，进而得到的正则化解，虽然达不到最优，但仍然能得到比较满意的精度 4 2扩散系数较小时逆时问题数值实现 小扩散系数的热传导方程有广泛的物理背景例如，在考古学中人们利用碳元素的 衰变周期来确定文物的年代，由于碳元素的衰变周期很长，因此该问题就表现为小扩散 磁卜g 一 = 丝卿 a + a 擎 一 = 赋i a+让a = 蘸 a+甜 算 0 汁0 图4 4 ：6 = 00 0 1 时取不同的相对误差 图4 5 ：6 = 0 0 1 时取不同a 的相对误差 2 l 系数的热传导方程。考虑如下问题： f 象= k 2 貉 ( z ，t ) ( o ，n ) ( o ，t 】 “z ( o ，。) 3 o o t ( 4 2 1 1 i “。( ? r ，t ) + u ( 7 r ，t ) = 00 t t 。 【u ( x ，0 ) = c o s ( 5 z )o ( 0 ，7 r ) 其中= o 1 。由于需要某一个时刻的温度u ( z ，t ) 作为反演过程的输入数据，但是无法 直接写出该问题的解析解，所以我们先用有限差分方法解正问题。 把区间( o ，”) 等分2 0 0 份，则空间步长为h = 盎，时间步长为r = i 0 7 3 。在内点采 用四阶紧格式，边界采用二阶偏心差商，将( 4 2 1 ) 转为下面差分方程； 其中 击( 盈u ；l - 一 + 1 0 巩“：一 十盈u ) = 2 。2 “：一 1 t 1 9 9 = ! 生1 2 1 h 互二生= 0 当t n 二= 鼍t 舞出+ “1 0 0 = 0 u ? = c o s ( s h i ) 0 i 2 0 0 u ：。h ，l r ) ，u ：一 ：竺! 堂2 ， 瓯屯掣，磋。：t 照掣 得到的数值解精度为o ( r 2 + h 2 ) 给定t 后，可从差分方程( 4 2 2 ) 求出u ( z ，t ) 。下面来构造算子a 。 ( 4 , 2 1 ) 对应的特征值问题为： 鬻兰鬻意：。，嚣； 显然( 4 2 3 ) 的通解是： ( z ) = a o c o s 、+ b 0 5 饥、，恹z ，其中a o ，b o 均为常数。由( 4 2 4 ) 知：c 0 8 弧7 r = 瓜s 讯瓠7 r ，令1 = 弧7 r ，得： 咖( 7 ) = 1 ( 4 2 5 ) 用n e w t o n 迭代法解上式，令f ( z ) = c t g x 一。i x ，迭代格式如下： m 。x n 一烈 显然，方程( 4 2 5 ) 有无穷多个根，本文取前2 0 个正根 这样特征值和单位正交的特征函数分别为： 佤= 垒7 r ， 驴去一佤。 、f 。 记为7 n ，n = 1 ，2 ，2 0 。 n = 1 2 - 2 0 其中l n = 兰+ i 去s f n 2 瓶。由( 2 1 4 ) ，( 2 1 、5 ) 式得： k ( x ，g ) = c o s 、瓦zc o s 硒e x p ( 一h 2 t ) ， 由( z ) 2 o ( 刚) g ( 口j 由。z 荟瓦1 c o s 瓜c o sv 瓦ve x p ( 一a n 2 t ) 9 国) d r ， g ( 。) = k ( z ，y ) k ( y ，z ) g ( z ) d z d y 同样把区间( 0 ，”) 作m 等分，利用复化梯形公式离散积分，则( 3 1 4 ) 写为： ( z ) + 筹j 曼= 0 0 j ( z ，蜥) 匡n t ( 胁矧以) = 三mi 尹= 0 q k ( z 矧“( 执用，咖气z ) + 磊凸j ( 蜥l 蚤啦( 协崩) 以) j2 三f q k ( 矧“( 执用， 再对z 离散，可得类似( 4 1 2 ) 的矩阵方程 a 2 ：正则化参数的先验取法 首先取a = o 0 0 0 1 ，6 = 0 时，对于不同的丁，在反演过程中反演的初值u ( z ，t ) 加上 1 0 的随机扰动后反演结果和初始值比较如图4 6 从图中可以知道当t 较小时反演的 结果不错，能够很好的反映初始温度分布振荡的待点和振幅大小，但丁取较大的值时结 果不太好这是由于正向热传导问题温度场是随时间是指数衰减( 图1 1 ) ，时间越长，正 问题的解就越趋向于常数，反演的结果自然就越差。取t = 0 5 ，区间( 0 ，”) 等分2 0 0 份即 m = 2 0 0 ，设有扰动的右端项为：= “+ 5 s i n ( 2 x 一1 ) 数值结果如表4 2 1 ： 表4 2 1 ：不同扰动的数值结果 从表4 2 1 可以看出当扰动为零时，数值解与精确解的绝对误差大部分为o ( 1 0 - 3 ) 。由于 在利用差分方程( 4 2 2 ) 求解“( z ，0 5 ) 时误差为o ( h 2 十r 2 ) ，在本例中就是o ( 1 0 _ 4 ) ，也就 图4 6 ：由不同时刻的温度场反演的初始温度 是说反演过程中未加扰动的初始数据已有o ( 1 0 “) 的误差，所以这样的数值结果是不错 的当扰动增大时，误差也随之增大 表4 2 2 反映了对于不同的正则化参数n 的数值解和精确解的相对误差图4 7 则表 现了对于不同的正则化参数a 各个节点上的误差，此时5 取0 0 0 5 。表4 , 2 2 租图4 7 均 说明正则化参数不能太大或太小。 表4 2 2 ：不同和5 的相对误差 b 2 ：正则化参数的后验取法 我们仍然采用第一节中的迭代格式，这里迭代初值o e 0 = 0 0 1 ，e = 1 0 一，n z = 5 0 ，m = 2 0 0 ，t = 0 5 。 当u 6 ：u + o 0 0 5 s i n ( 2 x 一1 ) 时，迭代4 次停止矿= 3 7 5 4 2 0 8 6 1 0 ，相对误差数 值结果见图4 8 。 m l j 暮& e a l 寻苦u l 图4 7 ：6 = o0 0 5 时取不同q 的反演结果在各点的误差 2 4 当u 6 = u 十o o l s i n ( 2 x 一1 ) 时，迭代5 次停止o + = 7 0 6 8 6 5 8 8 1 0 ，相对误差数 值结果见图4 9 。 数值结果不是最优，但精度还是令人满意的 图4 8 ：6 = od 0 5 时取不同。的相对误差 图4 9 ：5 = o 0 1 时取不同d 的相对误差 4 3 初值有不可导点的逆时问题数值实现 本节取的初值是有不可导点的连续函数， f 象= k 2 貉 ( z ，t ) ( o ，”) ( o ，刁 i ( o ，t ) = 0 0 t t iu 。( 7 r ，t ) + ，缸( 7 r ，t ) = 0 0 t 曼t 【u ( z ，0 ) = 9 ( z ) z ( 0 ， ) 茁【0 ，吾l z ( 量，7 r 】 显然，9 ( 。) 在。= 处不可导类似4 2 节中求解正问题的方法，我们在内点用四阶紧 格式，边界用二阶偏心差商，时间步长取r = 1 0 ，空间步长为h = 盎，取定t ，可以解 出u ( z ，t ) 作为反演的初始数据，同样按4 2 节的方式离散方程( 3 1 4 ) ，得到类似( 4 1 2 ) 的矩阵方程。 a 3