利用中心差分格式数值求解导数
差分格式
§1. 差分 1．一阶导数的差分近似（差商）导数的定义： ()()()0000limx x f x f x f x x x ®-¢=-导数的近似：()()()10010f x f x f x x x -¢»- （当 1x 与 0x 足够接近时）这样的表达式称为差商，它可作为导数的近似，称为导数的差分近似。
误差分析 － 泰勒展开：将 ()1f x 在 0x 处做泰勒展开，有()()()()()()21001001012f x f x f x x x f x x x ⅱ?=+-+-+L于是()()()()1001010f x f x f x x x x x -¢-=--各种差分近似： 取 0h >（称为步长），则可以有 向前差分近似（相当于取 100x x h x =+>）()()()000f x h f x f x h+-¢»向后差分近似（相当于取 100x x h x =-<）()()()000f x f x h f x h--¢»中心差分近似 （前差近似与后差近似的算术平均）()()()0002f x h f x h f x h+--¢»2． 差分近似的一般形式差分近似的一般形式可写成()()()()()()()()022********* m m n n f x c f x c f x c f x hc f x c f x c f x c f x ------é¢?++êë+ù++++úûL L或简写为()()01nj j j mf x c f x h =-¢»å 称为一阶导数 ()0f x ¢ 的一个 1m n ++ 点差分近似。
这里0 ( , , 2 , 1 , 0 , 1 , 2 , , )j x x jh j m n =+=---L L差分近似的精度 ： 阶 定义：若()()()01n pj jj mf x c f x h h =-¢-=å 则称表达式 ()1nj j j mc f x h =-å 是一阶导数 ()0f x ¢ 的 p 阶差分近似。
求解三维空间分数阶对流扩散方程的Douglas-Gunn格式
!第"#卷第#期郑州大学学报!理学版"$%&’"#(%’#!)*#+年,月-./012340%56278.!(9:.;<7.=>."?9@.)*#+收稿日期!)*#AB*#B *C 基金项目!国家自然科学基金项目!##ED#)C)".作者简介!聂玉峰!#+CA %"#男#陕西西安人#教授#主要从事高性能数值计算方法’计算材料学’并行计算研究#=B I 97&&M \271N 2[L5.1>5.<2$通信作者&胡嘉卉!#+A*%"#女#河南郑州人#博士研究生#主要从事偏微分方程数值解研究#=B I 97&&05R 0NI 97&.2[L5.1>5.<2.求解三维空间分数阶对流扩散方程的!9D )?4>=U D ::格式聂玉峰#!!胡嘉卉#!)!!王俊刚#!#.西北工业大学计算科学研究中心!陕西西安D#*#)+$).河南工业大学理学院!河南郑州E"***#"摘要!由于分数阶导数的非局部性特征#在模拟反常扩散现象时使用分数阶偏微分方程具有更好的效果#但是分数阶导数的非局部性也给数值分析和计算带来了很大困难#尤其在多维空间情形下.通过对经典d%53&9S B V 522格式的推广#提出一种求解三维空间分数阶对流扩散方程!SL9<1\@9<:7%29&9>81<:7%2>7\\5S 7%21O59:7%2#;F G d ="的交替方向隐!9&:1@29:723>7@1<:7%27I L&7<7:#G d ^"差分格式#并用矩阵法证明了其稳定性和收敛性.用数值算例进一步验证了该格式在空间和时间方向均具有较高的二阶收敛精度#可以高效地求解三维;FG d =.关键词!三维;F G d =$G d ^格式$Y @92JB (7<%&S %2格式$d %53&9S B V 522格式$稳定性$收敛性中图分类号!e)E#文献标志码!G 文章编号!#CD#B CAE#!)*#+"*#B **EEB *D !"#!#*’#,D*"Q R .7S S 2.#CD#B CAE#’)*#A*##$%引言近年来#自然界中的反常扩散现象受到科研人员的广泛关注#为研究其独特的物理过程#常常用分数阶偏微分方程建立相应的数学模型.其中#在包含对流和超扩散两个物理过程的散布现象中#粒子束的传播与经典的布朗运动模型不再一致#此时把经典对流扩散方程中的空间二阶导数替换成分数阶导数构建的空间分数阶对流扩散方程!;FG d ="能更准确地模拟这一现象.分数阶导数或积分具有非局部性#这给相应方程的求解带来了很大困难#在大多数情况下很难得到解析解#因此研究可靠而有效的数值方法就显得尤为重要.目前常用的数值解法包括有限差分法(#X ))’有限元法(,X E )’有限体积法(")’配点法(C )以及谱方法(D )等.由于三维模型在科学研究中有广泛应用#本文考虑有限区域上带有零d 7@7<0&1:边界条件的三维;F G d =的数值求解.分数阶导数是一个非局部算子#这就使得离散;FG d =得到的线性系统的刚度矩阵不再是稀疏矩阵#导致计算工作量和存储量都非常大#尤其在多维空间情形下.目前求解三维;FG d =的数值方法还比较少.文献(A )采用了一种交替方向稳定法!9&:1@29:723>7@1<:7%27I L&7<7:I 1:0%>#G d ^"差分格式求解三维;F G d =#并提高了精度.文献(+)提出了一种求解三维分数阶扩散方程的G d ^差分格式.文献(#*)研究了一种求解三维空间分数阶扩散方程的快速迭代G d ^有限差分方法.在本文中#我们将提出一种求解三维;F G d =的有效的G d ^有限差分格式#这种方法是将经典的d %53&9S B V 522格式中的二阶中心差分算子推广为包含左’右分数阶导数离散算子及一阶中心差分算子在内的复杂算子得到的#同时给出了该格式的稳定性和收敛性的必要证明.最后用数值实验验证了理论分析的结果.&%三维F ’7!L 及其!9D )?4>=U D ::格式本文考虑三维;FG d =及它的初边值条件为-@!I #C #J #""-"-:#I Q R "I @!I #C #J #""#:)I R "I +@!I #C #J #""#P #C Q R !C@!I #C #J #""# Copyright©博看网 . All Rights Reserved.!第#期聂玉峰#等$求解三维空间分数阶对流扩散方程的d%53&9S B V 522格式P )C R !C +@!I #C #J #""#2#J QR %J @!I #C #J #""#2)J R %J +@!I #C #J #""#X #-@!I#C #J #""-I#X )-@!I #C #J #""-C #X ,-@!I#C #J #""-J#!!I #C #J #""#!I #C #J #""(*/!*#Z )#!#"@!I #C #J #""-*#!!I #C #J #""(-*/(*#Z )#!)"@!I #C #J #*"-@*!I #C #J "#!!I #C #J "(*#!,"其中&#K "#!#%K )$*-!I Q #I +"/!C Q #C +"/!J Q #J +"4S ,$:##:)#P ##P )#2##2)!*是,个空间方向的左’右扩散系数$X #’X )’X ,分别是,个空间方向的对流系数$!!I #C #J #""是源项.方程!#"中的分数阶导数是_71I 922B ]7%587&&1型的#即函数A!I "的"!#K "K )"阶_71I 922B ]7%587&&1左导数和右导数分别定义为I Q R "IA !I "-#’!).""-)-I ),I I Q!I.4"#."A !4">4和I R "I +A !I "-#’!).""-)-I ),I +I!4.I "#."A !4">41&,&%S (<@4::=V (9D Q(??<分数阶导数的离散设*#’*)’*,和,为正整数#3#-!I +.I Q "L *##3)-!C +.C Q "L *)#3,-!J +.J Q "L *,#.-Z L ,分别是一致的空间步长和时间步长#由此定义的空间和时间的剖分为I ;-I Q #;3#;-*###-#*##C G -C Q #G 3#G -*###-#*)#J 9-J Q #93#9-*###-#*,#"7-7.#7-*###-#,1设@7;#G #9表示@!I ;#C G #J 9#"7"的近似值#!7;#G #9-!!I ;#C G #J 9#"7"1采用文献(##)中的方法离散方程!#"中的分数阶导数1以I方向为例#有IQR "I@!I #C G #J 9#"7"I -I ;-#’!E .""3"##;##5-*?"5@!I ;.5###C G #J 9#"7"#_!3)#"#IR "I +@!I #C G #J 9#"7"I -I ;-#’!E .""3"##*#.;##5-*?"5@!I ;#5.##C G #J 9#"7"#_!3)#"#其中系数为"5-##5-*#.E #),."#5-##C .)"."#,,."#5-)#!5##",.".E 5,."#C !5.#",.".E !5.)",."#!5.,",."#5!,1记)#"#I @7;#G #9-#’!E .""3"##;##5-*?"5@7;.5###G #9#!E ")."#I @7;#G #9-#’!E .""3"##*#.;##5-*?"5@7;#5.##G #91!""同时#用中心差分近似对流项的一阶导数#记R "#I @7;#G #9-!@7;###G #9.@7;.##G #9"L )3#1!C "为了便于表示#进一步引入记号R ^"#I @7;#G #9-X #R "#I @7;#G #9#)^#"#I @7;#G #9-:#)#"#I @7;#G #9#)^."#I @7;#G #9-:))."#I @7;#G #91C 和J 方向的记号可以类似表示.&,+%三维F ’7!L 的有限差分近似用公式!E "i !C "离散空间导数#时间方向采用Y@92JB (7<%&S %2格式.记)"#I o -)^#"#I #)^."#I #R ^"#I #)!#C o -)^#!#C #)^.!#C #R ^!#C #)%#J o -)^#%#J #)^.%#J #R ^%#J #然后方程!#"就可以表示为!#..))"#I ..))!#C ..))%#J "@!I ;#C G #J 9#"7##"-!##.))"#I #.))!#C #.))%#J "@!I ;#C G #J 9#"7"#.!!I ;#C G #J 9#"7##L )"#+7##;#G #91!D "存在正常数$^#使得+7##;#G #9’$^.!.)#3)##3))#3),"1接下来#在方程!D "中用近似值@7;#G #9代替函数值@!I ;#C G #J 9#"7"#并去掉高阶项#得到方程!#"的全离散格式#"E Copyright©博看网 . All Rights Reserved.郑州大学学报!理学版"第"#卷!#..))"#I ..))!#C ..))%#J "@7##;#G #9-!##.))"#I #.))!#C #.))%#J "@7;#G #9#.!7##L );#G #91!A "!!为便于计算#下面构造G d ^差分格式.将方程!A "左边加上高阶项!.)E )"#I )!#C #.)E )"#I )%#J #.)E )!#C )%#J ",!@7##;#G #9.@7;#G #9"..,A)"#I )!#C )%#J !@7##;#G #9#@7;#G #9"#再把适当的部分移项到方程的右边并分解因式#得!#..))"#I "!#..))!#C "!#..))%#J "@7##;#G #9-!##.))"#I "!##.))!#C "!##.))%#J "@7;#G #9#.!7##L );#G #91!+"采用经典的d %53&9S B V 522格式分解式!+"得到G d ^格式#即!#..))"#I "!@+-!.)"#I #.)!#C #.)%#J "@7;#G #9#.!7##L );#G #9$!#*"!#..))!#C"!@++-!@+$!##"!#..))%#J "!@-!@++$!#)"!@-@7##;#G #9.@7;#G #91!#,"与经典的d %53&9S B V 522格式不同#,个方向的二阶中心差分算子在此处分别被替换为)"#I ’)!#C ’)%#J #它们是包含了左’右分数阶导数离散算子等在内的复杂算子#可以认为是经典d %53&9S B V 522格式在求解分数阶方程中的推广.接下来#我们将给出收敛性和稳定性的必要证明.+%收敛性和稳定性分析显然#如果在格式!#*"i !#,"中消去中间解变量#则得到格式!+"#即格式!#*"i !#,"和!+"是等价的.下面用矩阵法证明格式!+"是无条件稳定和收敛的.首先把方程!+"表示成矩阵形式#令"7-(@7######@7)#####-#@7*#.######@7##)###@7)#)###-#@7*#.##)###-#@7##*).####@7)#*).####-#@7*#.##*).####@7####)#@7)###)#-#@7*#.####)#-#@7##*).##)#@7)#*).##)#-#@7*#.##*).##)#-#@7##*).##*,.##@7)#*).##*,.##-#@7*#.##*).##*,.#)P #!#E "为了书写简单起见#我们用记号"7-!!(@7;#G #9);-##)#-#*#.#"G -##)#-#*).#"9-##)#-#*,.#表示式!#E "#类似的表示还有87-!!(!7;#G #9);-##)#-#*#.#"G -##)#-#*).#"9-##)#-#*,.#1记0^I -:#.)’!E .""3"#71710"#:).)’!E .""3"#71710P"#X #.E 3#71711#!#""0^C -P #.)’!E .!"3!)710!17#P ).)’!E .!"3!)710P!17#X ).E 3)71117#!#C "0^J -2#.)’!E .%"3%,0%1717#2).)’!E .%"3%,0P%1717#X ,.E 3,11717#!#D "其中0"#0!#0%是P %1L&7:4矩阵#0"和1分别表示为0"-"#?"**-**")?"#?"**-*",?")?"#?"*-*222222"*#.)---?"#?"*"*#.#?"*#.)?"*#.,-?")?"##1-*#*-**.#*#*-**.#*#-*222222***-*#***-.#*#其中&7是单位矩阵$符号1表示b @%21<J1@积(#)).0!’0%与0"类似.利用上述记号#式!+"可以写为!7.0^I "!7.0^C "!7.0^J""7##-!7#0^I "!7#0^C "!7#0^J""7#.87##L )1!#A "!!为了证明式!#A "的稳定性和收敛性#下面列出一些相关的引理和定理.CE Copyright©博看网 . All Rights Reserved.!第#期聂玉峰#等$求解三维空间分数阶对流扩散方程的d%53&9S B V 522格式引理&(#,)!一个7阶实矩阵0是正定的#当且仅当矩阵,-!0#0P"L )是正定的$,是正定的#当且仅当,的特征值都是正的.引理+(#,)!设0是一个7阶复矩阵#0T 表示H 的共轭转置#记,-!0#0T "L)#则对于0的任意特征值(#它的实部满足不等式(I 72!,"’+!(!0""’(I 9‘!,"#这里(I 72!,"和(I 9‘!,"分别表示,的最小和最大特征值.定理&(A )!设0"是式!#""中的P %1L&7:4矩阵#则对于0"的任意特征值(#有+!(!0"""K *#并且0"是负定矩阵1同时#+!(!=#0"#=)0P """K *#=##=)!*#=)##=))8*1引理-!设0(S 9/7#1(S 8/B #2(S T /[#则!011"12-01!112"1证明!此结论可以由b @%21<J1@积的定义直接得到.引理.!设0#1(S 9/7#2(S B /"#则有!0#1"12-012#112#21!0#1"-210#2111证明!此结论可以由b @%21<J1@积的定义直接得到.引理/(#))!设0(S 9/7#1(S 8/B #2(S 7/T #3(S B /"#则!011"!213"-02113!(S 98/T ""1引理0(#))!对于任意的矩阵0和1#有!011"P -0P 11P .引理1(#))!设矩阵0(S 7/7有特征值0(;17;-##矩阵1(S 9/9有特征值01G 19G -#1则矩阵011的97个特征值为(#1##-#(#19#()1##-#()19#-#(71##-#(7191为了叙述并证明下述引理和定理#记%,%表示矩阵的)B 范数.引理2%设0(S 7/7是正定矩阵#则对任意的+(S 且+]*#有%!7#+0".#%K #1证明!由矩阵)B 范数的定义#有%!7#+0".#%)-I 9‘I 8*!!7#+0".#I #!7#+0".#I "!I#I "1设C -!7#+0".#I #则有%!7#+0".#%)-I9‘C 8*!C #C "!!7#+0"C #!7#+0"C "-#I 72C 8*(##)+!0C #C "!C #C "#+)!0C #0C "!C#C ")K #1引理6%设0(S 7/7是正定矩阵#则对任意的+(S 且+]*#有%!7#+0".#!7.+0"%K #1证明!由矩阵)B 范数的定义并记C-!7#+0".#I #可得%!7.+0"!7#+0".#%)-I 9‘I 8*!!7.+0"!7#+0".#I #!7.+0"!7#+0".#I "!I #I "-I 9‘C 8*!!7.+0"C #!7.+0"C "!!7#+0"C #!7#+0"C "-I 9‘C 8*!C #C ".)+!0C #C "#+)!0C #0C "!C #C "#)+!0C #C "#+)!0C #0C "K #1!!引理&$%设0#1#7(S 7/7#0和1乘积可交换#且!7.0".#’!7.1".#存在#则!7#0"与!7.1".##!7.0".#与!7.1".#也是乘积可交换的.证明!首先#由01-10#不难验证!7a 0"!7.1"-!7.1"!7a 0"1所以有!7#0"!7.1".#-!7.1".#!7.1"!7#0"!7.1".#-!7.1".#!7#0"!7.1"!7.1".#-!7.1".#!7#0"#!7.0".#!7.1".#-!!7.1"!7.0"".#-!!7.0"!7.1"".#-!7.1".#!7.0".#1定理+%由式!#""i !#D "定义的矩阵0^I ’0^C ’0^J 是负定的.证明!记0I -:#.)’!E .""3"#0"#:).)’!E .""3"#0P"#X #.E 3#1#则有0PI-:#.)’!E .""3"#0P"#:).)’!E .""3"#0"#X #.E 3#1P #0I #0P I )-!:##:)".)’!E .""3"#0"#0P")#由定理##可得+!(!0I #0P I )""K *#也就是(!0I #0P I )"K *#而且由引理,#E 和C #有0^I #0^PI)-7171!0I #0PI)"#再根据引理D #可以得到(!0^I #0^PI)"K *1最后由引理)和引理#可得+!(!0^I""K *#且0^I 是负定的1类似地#可以证明0^C 和0^J 是负定的.定理-%差分格式!+"是无条件稳定的.DE Copyright©博看网 . All Rights Reserved.郑州大学学报!理学版"第"#卷证明!我们利用差分格式!+"的矩阵形式!#A "证明.设"7和#7分别是格式!#A "对应于初值"*和#*的解#记’7o -"7.#7-!!(’7;#G #9);-##)#-#*#.#"G -##)#-#*).#"9-##)#-#*,.#1由式!#A "可得’7满足方程!7.0^I "!7.0^C "!7.0^J"’7##-!7#0^I "!7#0^C "!7#0^J"’71!#+"!!矩阵0^I ’0^C ’0^J 是乘积可交换的.事实上#只需验证!#""i !#D "中b @%21<J1@积形式的矩阵都是乘积可交换的.由引理,和"#可得!71710""!710!17"-!71!710"""!71!0!17""-71!!710""!0!17""-71!0!10""#同时!710!17"!71710""-!71!0!17""!71!710"""-71!!0!17"!710"""-71!0!10""#等等1由0^I ’0^C ’0^J 乘积可交换并利用引理#*#方程!#+"可以写为’7-!!7.0^J".#!7#0^J ""7!!7.0^C".#!7#0^C ""7!!7.0^I".#!7#0^I""7’*1由定理)和引理+可知%!7.0^I ".#!7#0^I "%K ##所以!7.0^I ".#!7#0^I "的谱半径小于##因此当7)q 时#!!7.0^)".#!7#0^I ""7收敛到零矩阵#也就是对任意的7!*#!!7.0^I ".#!7#0^I ""7有界1同理可证!!7.0^C ".#!7#0^C ""7和!!7.0^J".#!7#0^J""7对任意的7!*有界.这就证明了差分格式!+"是无条件稳定的.定理.%设@!I ;#C G #J 9#"7"是问题!#"^!,"的解#@7;#G #9是差分格式!+"的解1记<7;#G #9-@!I ;#C G #J 9#"7".@7;#G #9#<7-!!(<7;#G #9);-##)#-#*#.#"G -##)#-#*).#"9-##)#-#*,.##那么存在一个正常数$#使得%<7%K$!.)#3)##3))#3),"1证明!记*7-!!(+7;#G #9);-##)#-#*#.#"G -##)#-#*).#"9-##)#-#*,.##那么方程!D "减去!+"的矩阵形式为!7.0^I "!7.0^C "!7.0^J "<7##-!7#0^I "!7#0^C "!7#0^J"<7#*7##1!)*"因为0^I ’0^C ’0^J 乘积可交换#根据引理#*#方程!)*"可以写为<7##-!7.0^J ".#!7#0^J "!7.0^C ".#!7#0^C ",!7.0^I ".#!7#0^I "<7#!7.0^J ".#!7.0^C ".#!7.0^I".#*7##1用矩阵的)B 范数作用上式两端#并利用引理A 和+#可得%<7%K %<*%##7.#X -*%*X ##%1因<*-$#所以存在正常数$#使得%<7%K$!.)#3)##3))#3),"1-%数值结果下面#我们通过两个数值算例验证本文所提出的数值格式的稳定性和收敛阶#也就是说格式是有效的#并在时间和空间方向都具有较高的二阶收敛精度.设"和"3分别表示问题!#"i !,"的解析解和采用格式!#*"i !#,"得到的数值解#用离散的Q q 和Q )范数计算全局截断误差#即%"3."%Q q o -I 9‘#’9’*,.##’G ’*).##’;’*#.#@,;#G #9.@!I ;#C G #J 9#Z"#%"3."%Q )o -!#*,.#9-##*).#G -##*#.#;-#@,;#G #9.@!I ;#C G #J 9#Z ")3#3)3,"#L )1!!算例&%在问题!#"i !,"中#取*-!*##"/!*##"/!*##"#Z-#$对流和扩散系数分别为X #-X )-X ,-.##:#-:)j P #-P )-2#-2)-#$初值取@*!I #C #J "-I ,!#.I ",C ,!#.C ",J ,!#.J ",1已知的解析解为@!I#C #J #""-1."I ,!#.I ",C ,!#.C ",J ,!#.J ",#由以上条件容易算出!!I #C #J #""1对常系数算例#取优化的步长比例,-*#-*)-*,进行测试.表#列出的数值结果表明#用格式!#*"i !#,"计算常系数问题!#"i !,"时#算法是无条件稳定的#而且在时间及空间方向都是二阶收敛的#这和理论分析的结果一致.算例+%在问题!#"i !,"中#取*-!*##"/!*##"/!*##"#Z -#$对流和扩散系数分别为X #-*’)"I #X )-*’)"C #X ,-*’)"J #:#-I "#:)-!#.I ""#P #-C !#P )-!#.C "!#2#-J %#2)-!#.J "%$初值取@*!I #C #J "-I ,!#.I ",C ,!#.C ",J ,!#.J ",1已知的解析解为@!I#C #J #""-1."I ,!#.I ",C ,!#.C ",J ,,!#.J ",#由以上条件!!I #C #J #""容易算出.表)列出了变系数算例)的数值结果#这里也取优化的步长比例,-*#-*)-*,进行测试#数值结果表明用格式!#*"i !#,"计算变系数问题!#"i !,"时#算法是无条件稳定的#而且在时间及空间方向也都具AE Copyright©博看网 . All Rights Reserved.!第#期聂玉峰#等$求解三维空间分数阶对流扩散方程的d%53&9S B V522格式!!表&!算例#在时刻"-##取,-*#-*)-*,的数值误差和收敛阶345*&!P011@@%@S92><%281@312<1%@>1@S\%@1‘9I L&1#9:"j#[7:0,-*#-*)-*,.%"3."%Q q+:"<%"3."%Q)+:"<"-#’) !-#’) %-#’)#Q#*#Q)*#Q E*#Q A*#Q#C*#’#*)D1X**D)’CD*#1X**AC’C,+"1X**+#’CCC,1X**+E’#D*"1X*#*X)’*EC#)’**D D#’++E E#’++A E#’+,C E1X**AE’D)E)1X**+#’#DA E1X**+)’+"))1X*#*D’,+C E1X*##X)’*,"))’**,)#’++D*#’++C+"-#’E !-#’" %-#’C #Q#*#Q)*#Q E*#Q A*#Q#C*#’)",E1X**D,’*#C+1X**AD’"*,)1X**+#’AAC E1X**+E’DE)"1X*#*X)’*"E D)’**D"#’++#+#’++#+)’*A+*1X**A"’*D+,1X**+#’)D*D1X**+,’)**C1X*#*A’*CE+1X*##X)’*E*##’+++*#’+A+)#’+AA C"-#’+ !-#’+ %-#’+#Q#*#Q)*#Q E*#Q A*#Q#C*)’,*"A1X**DE’D,+A1X**A#’*A""1X**A)’C)A"1X**+C’"#,D1X*#*X)’)A)E)’#)C")’*EC*)’*#)D,’"*"E1X**AD’,E,A1X**+#’D*,D1X**+E’#"A E1X*#*#’*,E"1X*#*X)’)""*)’#*D+)’*,E C)’**D#表+!算例)在时刻"-##取,-*#-*)-*,的数值误差和收敛阶345*+!P011@@%@S92><%281@312<1%@>1@S\%@1‘9I L&1)9:"j#[7:0,-*#-*)-*, .%"3."%Q q+:"<%"3."%Q)+:"<"-#’) !-#’) %-#’)#Q#*#Q)*#Q E*#Q A*#Q#C*#’*+)+1X**D)’C#E"1X**AC’EAE A1X**+#’C)#)1X**+E’*"A,1X*#*X)’*C,C)’*##E)’****#’++A#)’#*A C1X**AE’+C*C1X**+#’)),D1X**+,’*")*1X*#*D’C,**1X*##X)’*AD D)’*#+,)’**,E)’****"-#’E !-#’" %-#’C #Q#*#Q)*#Q E*#Q A*#Q#C*#’*A#D1X**D)’C#"E1X**AC’"#E+1X**+#’C,"D1X**+E’##*+1X*#*X)’*EA))’**")#’++,A#’++)E)’*,"A1X**AE’AE*+1X**+#’)**"1X**+,’**A D1X*#*D’"",C1X*##X)’*D)))’*##C#’++C E#’++,+"-#’+ !-#’+ %-#’+#Q#*#Q)*#Q E*#Q A*#Q#C*+’AE,"1X**A)’,C#,1X**A"’ADC E1X**+#’EA*D1X**+,’DEC#1X*#*X)’*"+C)’**C C#’+AA D#’+A)A#’C+*C1X**AE’*C#"1X**+#’*#,)1X**+)’""A,1X*#*C’EA,"1X*##X)’*"D")’**,##’+A"D#’+A*,有二阶收敛率#这和理论分析的结果是非常吻合的..%结论本文将求解三维整数阶抛物方程的经典d%53&9S B V522格式推广到分数阶#提出了一种求解三维;F G d=的有效的数值方法#并证明该格式具有无条件稳定性和较高的二阶收敛精度#必要而充足的数值实验验证了理论结果.最后#由于分数阶导数是非局部算子#对于多维空间问题的求解需要耗费较大的计算工作量和空间存储量#在今后的工作中#我们将考虑开展适当的快速算法#以减少计算花费和加快计算速度.参考文献!(#)!;e6;G=.G21‘L&7<7:0730%@>1@I1:0%>\%@\@9<:7%29&9>81<:7%2>7\\5S7%21O59:7%2S(-).-%5@29&%\<%I L5:9:7%29&L0M S7<S# )*#E#)DA&)"D X)DE.())!/T G(VT#]^6F#/T6G(V U#1:9&.(5I1@7<9&929&M S7S%\921[S L9<1B:7I189@79Z&1\@9<:7%29&%@>1@9>81<:7%2B>7S L1@S7%2 1O59:7%2(-).G LL&71>I9:01I9:7<S92><%I L5:9:7%2#)*#E#)E)&"E#X""*.+E Copyright©博看网 . All Rights Reserved.*"郑州大学学报!理学版"第"#卷(,)!=_$^($-#T=6=_(#_e e U-U.(5I1@7<9&9LL@%‘7I9:7%2%\9:7I1>1L12>12:#2%2&7219@#S L9<1B\@9<:7%29&>7\\5S7%21O59:7%2 (-).;^G?R%5@29&%225I1@7<9&929&M S7S#)**D#E"!)"&"D)X"+#.(E)!/T=(Vff#]^YU#/T G e/V.G2%:1%2:01\727:11&1I12:I1:0%>\%@:01S L9<1\@9<:7%29&9>81<:7%2>7\\5S7%21O59:7%2(-).Y%I L5:1@S92>I9:01I9:7<S[7:09LL&7<9:7%2S#)*#*#"+&#D#A X#D)C.(")!T=-G/^T#?e_e(=fP#]^6F.;:9Z7&7:M92><%281@312<1%\9\727:18%&5I1I1:0%>\%@:01S L9<1\@9<:7%29&9>81<:7%2B>7S L1@B S7%21O59:7%2(-).-%5@29&%\<%I L5:9:7%29&92>9LL&71>I9:01I9:7<S#)*#E#)""&CAE X C+D.(C)!虎晓燕#韩惠丽.重心插值配点法求解分数阶F@1>0%&I积分方程(-).郑州大学学报!理学版"#)*#D#E+!#"&#D X),.(D)!/T=(V?#]^6F#G(T$#1:9&.G0730B%@>1@S L1<:@9&I1:0%>\%@:01I5&:7B:1@I:7I1B\@9<:7%29&>7\\5S7%21O59:7%2S(-).G LBL&71>I9:01I9:7<9&I%>1&&723#)*#C#E*!D Q A"&E+D*X E+A".(A)!d=(Va T#Y T=(?T.=\\7<712:25I1@7<9&9&3%@7:0I S\%@:0@11B>7I12S7%29&\@9<:7%29&L9@:79&>7\\1@12:79&1O59:7%2S(-).-%5@29&%\<%I L5:9:7%29&I9:01I9:7<S#)*#E#,)!E"&,D#X,+#.(+)!Y T=(-#]^6F#]^6h#1:9&.(5I1@7<9&S7I5&9:7%2\%@:01:0@11B>7I12S7%2\@9<:7%29&S5ZB>7\\5S7%21O59:7%2(-).G LL&71> I9:01I9:7<9&I%>1&&723#)*#E#,A!#""&,C+"X,D*".(#*)aG(VT#d6(.F9S:9&:1@29:723B>7@1<:7%2\727:1>7\\1@12<1I1:0%>S\%@:0@11B>7I12S7%29&S L9<1B\@9<:7%29&>7\\5S7%21O59:7%2S (-).-%5@29&%\<%I L5:9:7%29&L0M S7<S#)*#E#)"A&,*"X,#A.(##)Y T=(?T#d=(Va T.GS1<%2>B%@>1@25I1@7<9&I1:0%>\%@:[%B>7I12S7%29&:[%B S7>1>S L9<1\@9<:7%29&<%281<:7%2>7\\5S7%2 1O59:7%2(-).G LL&71>I9:01I9:7<9&I%>1&&723#)*#E#,A!#,"&,)EE X,)"+.(#))]G6HG-.?9:@7‘929&M S7S\%@S<712:7S:S92>1237211@S(?).U07&9>1&L079&;^G?#)**".(#,)h6G_P=_e(^G#;G Y Y e_#;G]=_^F.(5I1@7<9&I9:01I9:7<S(?).H1@&72&;L@7231@B$1@&93#)**D.!9D)?4>=U D::’(:(;<!(H H<B<:E<F E N<@<H9B3N B<<=C(@<:>(9:4?F T4E<’B4E;(9:4?7C Q<E;(9:!(H H D>(9:L X D4;(9:(^=f5\123##T6-79057##)#aG(V-523923#!#’+<B<:823$<7"<8!48$49T@":";47:5)2;<72<#*48"3E<B"<87(45C"<237;2:5&7;A<8B;"C#M;N:7D#*#)+#$3;7:$ )’$455<?<4!)2;<72<#\<7:7&7;A<8B;"C4!Z<237454?C#O3<7?J34@E"***##$3;7:"75>;B4E;&d51:%:012%2B&%<9&7:M%\\@9<:7%29&>1@789:781S#\@9<:7%29&L9@:79&>7\\1@12:79&1O59:7%2S[1@1Z1::1@:%>1S<@7Z192%I9&%5S>7\\5S7%2L012%I129:092%:01@I1:0%>S.T%[181@#[07&112R%M723:01<%2B 812712<1\@%II9:01I9:7<9&I%>1&723#7:9&S%<95S1>&%:S%\:@%5Z&11S L1<79&&M72S%&8723I5&:7>7I12S7%29& <9S1S.G21\\7<712:25I1@7<9&9&3%@7:0I[9S L@%L%S1>\%@S%&8723:01:0@11B>7I12S7%29&S L9<1\@9<:7%29&9>81<:7%2>7\\5S7%21O59:7%2!;F G d="ZM3121@9&74723:01d%53&9S B V522S<01I1.;:9Z7&7:M92><%281@B312<1%\:01I1:0%>[1@1L@%81>ZM:01I9:@7‘I1:0%>.P01>1@781>9&:1@29:723>7@1<:7%27I L&7<7:!G d^"\727:1>7\\1@12<1S<01I109>:01S1<%2>%@>1@9<<5@9<M72Z%:0:7I192>S L9<1>7@1<:7%2S#@1S L1<:781&M.P011\\7<712<M92><%281@312<1%@>1@S[1@1\729&&M>1I%2S:@9:1>ZM S%I125I1@7<9&1‘9I L&1S.J<G K9B C>&:0@11B>7I12S7%29&;F G d=$G d^S<01I1$Y@92JB(7<%&S%2S<01I1$d%53&9S B V522S<01I1$ S:9Z7&7:M$<%281@312<1!责任编辑&方惠敏"Copyright©博看网 . All Rights Reserved.。
一阶导数边界处的高阶差分格式
一阶导数边界处的高阶差分格式在数学和科学领域中，一阶导数边界处的高阶差分格式是一个重要且复杂的概念。
它在数值计算和科学工程中有着广泛的应用。
本文将探讨一阶导数边界处的高阶差分格式，讨论其基本原理、应用场景以及个人理解。
1. 基本原理一阶导数边界处的高阶差分格式是用差分来逼近微分的方法。
在边界处求解一阶导数时，我们需要使用高阶的差分格式来更精确地逼近导数的值。
常见的高阶差分格式包括中心差分、向前差分和向后差分。
在边界处，特别是当函数在该点不可导或者导数变化剧烈时，高阶差分格式能够提高逼近的准确性。
2. 应用场景一阶导数边界处的高阶差分格式广泛应用于数值计算、科学工程和实际问题的求解中。
在物理学中，当需要计算边界处的导数时，高阶差分格式能够提供更精确的数值结果。
在工程领域中，处理边界条件时，高阶差分格式也能够有效地提高数值计算的准确性。
在金融领域和生物医学领域，一阶导数边界处的高阶差分格式也有着重要的应用价值。
3. 个人理解对于我个人而言，一阶导数边界处的高阶差分格式是一个具有挑战性但又十分重要的概念。
通过学习和应用高阶差分格式，我意识到数值计算中的边界条件处理非常关键，而高阶差分格式能够帮助我们更准确地处理这些边界条件。
在我的实际工作中，我也经常需要使用高阶差分格式来解决复杂的数值计算问题，因此对其原理和应用有着更深入的理解和实际经验。
总结一阶导数边界处的高阶差分格式是数值计算和科学工程中不可或缺的重要概念。
通过本文的探讨，我们对其基本原理和应用场景有了更深入的了解，并且探讨了个人对该概念的理解和应用经验。
在今后的工作和研究中，我将继续深入学习和探索一阶导数边界处的高阶差分格式，以提高自己在数值计算和科学工程领域的能力。
通过以上内容，我希望本文能够对一阶导数边界处的高阶差分格式有一个更全面、深刻和灵活的理解。
一阶导数边界处的高阶差分格式在数学和科学领域中扮演着重要的角色。
它是一种用差分逼近微分的方法，特别在边界处求解一阶导数时，能够提供更精确的数值结果。
二阶混合偏导差分格式的详解
二阶混合偏导差分格式的详解文章标题：深度解析二阶混合偏导差分格式正文：一、引言在数学和计算机科学领域中，二阶混合偏导差分格式是一种重要的数值计算方法，它在求解偏微分方程和数值模拟中发挥着关键作用。
本文将针对二阶混合偏导差分格式进行全面深入的解析，从简单的定义和原理出发，逐步深入探讨其数学推导和应用场景，以帮助读者更好地理解和应用这一重要的数值计算方法。
二、基本概念二阶混合偏导差分格式是一种数值计算方法，用于求解偏微分方程。
它通过将偏微分方程中的偏导数用差分表示，然后构建差分方程并进行数值求解，从而得到偏微分方程的近似解。
在二阶混合偏导差分格式中，一阶导数使用中心差分逼近，二阶导数则结合前向和后向差分来逼近，从而实现更高的数值精度和稳定性。
三、数学推导接下来，我们将通过数学推导来深入探讨二阶混合偏导差分格式的具体计算过程。
假设我们要求解的偏微分方程为一维扩散方程，即∂u/∂t = α∂²u/∂x²我们可以使用二阶混合偏导差分格式来离散化这个方程，首先对时间t 进行离散化，然后对空间x进行离散化，最终得到差分方程，并通过迭代求解得到数值解。
具体的推导过程将在以下几个步骤中详细展开。
第一步，对时间进行离散化，使用显式欧拉方法进行离散化求解。
第二步，对空间进行离散化，使用中心差分逼近一阶导数，结合前向和后向差分逼近二阶导数。
第三步，构建差分方程，将离散化的时间和空间方程组合在一起。
第四步，通过迭代求解差分方程，得到偏微分方程的数值解。
通过以上数学推导，我们可以清晰地了解二阶混合偏导差分格式的具体计算过程，以及其在求解偏微分方程中的应用。
四、应用场景二阶混合偏导差分格式主要适用于求解具有时间和空间耦合的偏微分方程，例如扩散方程、波动方程和热传导方程等。
在工程、物理、生物和环境科学等领域，这些偏微分方程广泛应用于模拟和预测各种复杂现象，而二阶混合偏导差分格式则成为了求解这些偏微分方程的重要数值方法之一。
数值模拟偏微分方程的三种方法：FDM、FEM及FVM
数值模拟偏微分方程的三种方法：FDM、FEM及FVM偏微分方程数值模拟常用的方法主要有三种：有限差分方法（FDM）、有限元方法（FEM）、有限体积方法（FVM），本文将对这三种方法进行简单的介绍和比较。
有限差分方法有限差分方法(Finite Difference Methods)是数值模拟偏微分方程最早采用的方法，至今仍被广泛运用。
该方法包括区域剖分和差商代替导数两个过程。
具体地，首先将求解区域划分为差分网格，用有限个网格节点代替连续的求解区域。
其次，利用Taylor级数展开等方法将偏微分方程中的导数项在网格节点上用函数值的差商代替来进行离散，从而建立以网格节点上的值为未知量的代数方程组。
该方法是一种直接将微分问题变为代数问题的近似数值解法，数学概念直观，表达简单，是发展较早且比较成熟的数值方法。
差商代替导数后的格式称为有限差分格式，从格式的精度来考虑，有一阶格式、二阶格式和高阶格式。
从差分的空间离散形式来考虑，有中心格式和迎风格式。
对于瞬态方程，考虑时间方向的离散，有显格式、隐格式、交替显隐格式等。
目前常见的差分格式，主要是以上几种格式的组合，不同的组合构成不同的差分格式。
差分方法主要适用于结构网格，网格的步长一般根据问题模型和Courant稳定条件来决定。
请输入标题有限元方法(Finite Element Methods)的基础是变分原理和分片多项式插值。
该方法的构造过程包括以下三个步骤。
首先，利用变分原理得到偏微分方程的弱形式(利用泛函分析的知识将求解空间扩大)。
其次，将计算区域划分为有限个互不重叠的单元(三角形、四边形、四面体、六面体等)。
再次，在每个单元内选择合适的节点作为求解函数的插值点，将偏微分方程中的变量改写成由各变量或其导数的节点值与所选用的分片插值基函数组成的线性表达式，得到微分方程的离散形式。
利用插值函数的局部支集性质及数值积分可以得到未知量的代数方程组。
有限元方法有较完善的理论基础，具有求解区域灵活(复杂区域)、单元类型灵活(适于结构网格和非结构网格)、程序代码通用(数值模拟软件多数基于有限元方法)等特点。
迎风差分法-概述说明以及解释
迎风差分法-概述说明以及解释1.引言1.1 概述迎风差分法是一种常用的数值计算方法，用于求解偏微分方程的数值解。
它以差分代替偏导数，将连续问题离散化为离散问题，通过逼近求解离散方程组来得到精确的数值解。
该方法的基本原理是将求解区域等分为小网格，利用节点上的函数值和它相邻节点上的函数值之间的差异来逼近导数的值。
根据所采用的向前差分、向后差分或中心差分的方式，可以得到迎风差分法的不同形式。
迎风差分法的应用场景非常广泛。
它可以用于求解一维或多维的空间问题，例如热传导方程、波动方程、扩散方程等。
此外，该方法还可以应用于流体力学、计算机模拟等领域，用于模拟复杂的物理现象和工程问题。
迎风差分法具有一些优点和缺点。
其优点之一是简单易懂，容易实现。
通过简单的差分运算，即可得到数值解。
此外，该方法的计算效率较高，可应对大规模离散问题。
然而，迎风差分法也存在一些缺点，例如对于非线性问题的处理相对困难，可能会出现数值耗散和数值耗散等问题。
综上所述，迎风差分法是一种重要的数值计算方法，能够有效地求解偏微分方程的数值解。
它在科学研究和工程实践中具有广泛的应用前景。
本文将深入探讨迎风差分法的基本原理、应用场景以及优缺点，并对其未来的发展做出展望。
文章结构部分的内容可以按照以下方式编写：1.2 文章结构本文将按照以下结构进行论述：第一部分是引言部分，包括概述、文章结构和目的。
引言部分将对迎风差分法进行简要介绍，说明文章的结构和目的，为读者提供文章的整体框架。
第二部分是正文部分，主要包括迎风差分法的基本原理、应用场景和优缺点。
在这一部分，将对迎风差分法的基本原理进行详细阐述，包括其基本概念、基本步骤和数学原理等。
随后，将介绍迎风差分法在各个领域的应用场景，包括物理学、工程学和计算机科学等。
同时，还将探讨迎风差分法的优点和缺点，对其进行客观评价，为读者提供全面的了解。
第三部分是结论部分，主要包括总结迎风差分法的重要性、对其展望和最终结论。
四点差分格式
四点差分格式
四点差分格式（Four Point Difference Formulation）是一种数值计算方法，用于求解偏微分方程。
该方法主要应用于有限差分法中，通过在求解域上离散化偏微分方程，将其转化为代数方程组进行求解。
四点差分格式通常采用四个网格点的值来逼近偏微分方程中的导数项，通过这种方式可以将偏微分方程离散化为代数方程组。
具体来说，对于一个偏微分方程中的导数项，我们可以用四个网格点的值来表示该导数项的逼近值。
这四个网格点通常包括当前点、前一点、后一点和下一点，形成一个矩形区域。
通过选取适当的权值，我们可以将这四个网格点的值组合起来，得到导数项的逼近值。
四点差分格式的优点是简单易懂，易于编程实现。
同时，由于该方法只涉及到四个网格点的值，因此计算量相对较小，适用于大规模的数值计算问题。
但是，四点差分格式也存在一些缺点，例如对于复杂边界条件和不规则求解域的处理不够灵活，可能会引入较大的数值误差等。
在实际应用中，需要根据具体的问题和要求选择合适的差分格式和离散化方法。
同时，还需要对离散化后的代数方程组进行适当的处理和优化，以提高计算效率和精度。
中国科学院大学2020年计算传热学期末考试题目
1. 数值求解不可压缩Navier-Stokes 方程时，我们讲了MAC 方法、Projection 方法和压力校正法，写出它们的求解步骤，并着重分析如何保证计算所得的速度场满足无散度约束条件。
求解可压缩流体力学方程，我们介绍了时间推进法，请简要介绍该方法，并解释为什么可以用时间推进法？（12分）2. 不可压流体的计算中，我们介绍了同位网格和交错网格。
请详细说明这两种网格布置的差异，并对比分析其优缺点。
同时解释可压缩流体计算中不需要引入交错网格的原因。
（8分）3. 数值解法的特性有相容性，稳定性，收敛性，守恒性等，请分别解释其含义并叙述它们之间的联系；解释对流稳定性和数值稳定性的区别。
（10分）4. 在专题报告中，我们介绍了浸没边界法和重叠网格法，请以二维圆柱绕流为例，简述贴体网格法、浸没边界法和重叠网格法的网格特点（可以绘制简略图表示网格分布特点），并简要说明浸没边界法和重叠网格法对比贴体网格法的优势。
（6分）5. 对于一阶线性对流方程0u u a t x ∂∂+=∂∂如果假设对流项的速度a 为常数。
（1）请写出其特征线方程和相容关系。
说明如果对流项分别采用中心差分格式和一阶迎风格式会各自带来什么影响？（2）请根据特征性理论，插值推导Lax-Friedrichs 格式；分析该格式的相容性并基于Von Neumann 分析法给出其稳定性条件。
（3）分析中心差分格式和Lax-Friedrichs 格式的耗散特性，并解释Lax-Friedrichs 格式相比中心差分格式，其稳定性条件得到改善的原因。
（24分）6. 考虑二维不可压缩无旋流动，试用张量公式证明速度u 可以用速度势函数表示，即u ϕ=∇, 将该形式代入连续性方程，可以得到一个关于速度势ϕ的拉普拉斯方程，离散该方程，并写出用Gauss-Seidel 迭代法辅之以亚松弛方法求解该离散方程的表达式。
（12分） 7. 考虑非定常Navier-Stokes 方程21Re i i i j j j j i u u u p u t x x x x ∂∂∂∂+=−∂∂∂∂∂基于SIMPLEC 的思想，给出两位直角坐标系均匀交错网格下的求解步骤及具体离散表达式。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
利用中心差分格式数值求解导数目录一、问题描述 (2)二、格式离散 (2)二阶导数中心差格式离散 (2)追赶法求解线性方程组简述 (3)计算流程图 (5)三、程序中主要符号和数组意义 (5)四、计算结果与讨论 (6)五、源程序 (9)一、问题描述利用中心差分格式近似导数22/dx y d ，数值求解 ()x dx y d 2sin 22= ()10≤≤x1/,0/10====x x y y 步长分别取 0001.0,001.0,01.0,05.0=∆x二、格式离散将x 轴上[0,1]之间的线段按上述步长，等步长的离散为n 个小段，包括端点，共n+1个网格节点，示意图如下：线段上边的数字表示x 轴上的坐标值，线段下边的数字表示节点编号，从0到n 编号。
二阶导数中心差格式离散211222)2sin(x y y y dx y d x i i i ∆+-==+- 整理为线性方程形式)2sin(2211x x y y y i i i ∆=+-+-其中，x ∆ 为空间离散步长；i=1,2，……，n-1包括边界条件的线性方程组如下：边界条件边界条件0)*)1(*2sin(2...................)**2sin(2..................)*1*2sin(2021221122100=∆-∆=+-∆∆=+-∆∆=+-=--+-n n n n i i i y x n x y y y x i x y y y x x y y y y 改写成矩阵形式：f Ay = 其中，⎪⎪⎪⎪⎪⎭⎪⎪⎪⎪⎪⎬⎫⎪⎪⎪⎪⎪⎩⎪⎪⎪⎪⎪⎨⎧----=1012112112112101 A ，⎪⎪⎪⎪⎭⎪⎪⎪⎪⎬⎫⎪⎪⎪⎪⎩⎪⎪⎪⎪⎨⎧=-n n i y y y y y y 110 ，⎪⎪⎪⎪⎭⎪⎪⎪⎪⎬⎫⎪⎪⎪⎪⎩⎪⎪⎪⎪⎨⎧=-n n i f f f f f f 110 系数矩阵A 中仅三对角线上的数值不全为0，其余位置上的数值全为0，是典型的对角占优的三对角矩阵，列向量f 中，)2sin(2x i x f i ∆∆=，且10==n f f ，作为边界条件。
追赶法求解线性方程组简述⎪⎪⎪⎪⎭⎪⎪⎪⎪⎬⎫⎪⎪⎪⎪⎩⎪⎪⎪⎪⎨⎧=⎪⎪⎪⎪⎪⎭⎪⎪⎪⎪⎪⎬⎫⎪⎪⎪⎪⎪⎩⎪⎪⎪⎪⎪⎨⎧----=---n n n n n i i i b a c b a c b a c b a c b A 111111001012112112112101对A 做Crout 分解，即LU A =,⎪⎪⎪⎪⎭⎪⎪⎪⎪⎬⎫⎪⎪⎪⎪⎩⎪⎪⎪⎪⎨⎧=--n n n n s r s r s r s r s L 1122110 ，⎪⎪⎪⎪⎭⎪⎪⎪⎪⎬⎫⎪⎪⎪⎪⎩⎪⎪⎪⎪⎨⎧=-111111210n t t t t U 其中i i i t r s ,,为待定常数，由矩阵乘法可以得到下面的式子：1,,3,2,1,/,,3,2,1,,,,3,2,1,100000-===-=====-n i s c t n i t a b s b c t b s ni a r i i i i i i i i i将对角占优三对角矩阵线性方程组f Ay =等价为如下两个方程组f Lg =，g Uy =求解对角占优三对角矩阵线性方程组的追赶法步骤：①输入数据i i i c b a ,,②计算i i t s ,③求解方程组f Lg =n i s g a f g b f g i i i i i ,,3,2,1,/)(/1000 =-==-④求解方程组g Uy = 0,1,,2,1,1 --=-==+n n i y t g y g y i i i i nn⑤输出T n y y y y ),,,(10 =计算流程图三、程序中主要符号和数组意义符号或数组意义A、B、C、D、filename 用于自动更改dat文件名的字符串变量h 离散步长n 离散网格数，共n+1个网格节点p 辅助变量，暂时记录网格节点上的y值数组x，y 离散节点的x，y坐标子程序数组a,b,c 记录系数矩阵占优对角线上的值子程序数组f 记录线性方程组常数向量子程序数组s,r,t,g 追赶法求解线性方程组需要用到的中间辅助向量四、计算结果与讨论不同步长的数值结果函数曲线与精确解的对比从对比图中可以看到，在所取的四个计算步长下，数值计算结果与精确解都符合得相当好Step Accuracy（Max-error）0.05 8.152404966677018E-0050.01 3.264604683472783E-0060.001 3.817221055912867E-0080.0001 9.862256566961491E-009不同网格步长的计算精度由相应步长下所有网格节点上数值解与精确解的误差的最大值来度量，即上表中的Max-error来度量。
从上表中可以看出，随着网格节点的加密，Max-error的数量级在降低，即数值解的精度提高。
上述四个步长中，将线段离散成一万个网格时，数值结果的精度最高。
五、源程序program mainimplicit nonecharacter(13) filename !定义了五个字符串变量，用于按输出数据的需要自动更改dat文件的文件名character(3) Acharacter(6) Bcharacter(4) Ccharacter(3) Dinteger :: n,ireal(8) :: h,error,preal(8),allocatable :: y(:),x(:) !声明可变长度数组，x、y轴坐标定义为可变数组，数组长度按需要自动更改write(*,*)"输入步长："read (*,*)h !读入空间步长n=NINT(1.0/h) !nint命令，取与浮点数最接近的整数allocate (y(0:n),x(0:n)) !指定可变数组的长度A="po-" !po 代表problemopen(unit=11,file="help.txt") !打开一个txt文件，用于帮助更改dat文件的文件名write(11,"(f6.4)")hrewind(11) ! 将文件的读写位置移回到文件的最前面read(11,"(A6)")BC=".dat"filename=A//B//Ccall subsolution(y,n,h) !调用追赶法求解线性方程组的子程序open(unit=10,file=filename)do i=0,nx(i)=real(i)*hend dowrite(10,*)'TITLE = "X - Y CURVE"' !写入到dat文件中的一段字符，便于用tecplot软件后处理计算数据write(10,*)'VARIABLES = "X", "Y"'write(10,"('ZONE T=""Problem1-',f8.5,'"", I=',I6,', F=POINT')")h,n+1do i=0,nwrite(*,"(F10.8,10x,F10.8)")x(i),y(i)write(10,"(F10.8,10x,F10.8)")x(i),y(i) !将数值解数据写入dat文件end doy(0)=0y(n)=1error=0.0do i=1,n-1p=y(i)y(i)=(1+0.25*sin(2.0))*(i*h)-0.25*sin(2.0*i*h) !求解节点精确值error=max(error,abs(p-y(i))) !寻找各个节点误差的最大值end dowrite(10,"('ZONE T=""Problem1-',f8.5,'-exact"", I=',I6,', F=POINT')")h,n+1do i=0,nwrite(*,"(F10.8,10X,F10.8)")x(i),y(i)write(10,"(F10.8,10X,F10.8)")x(i),y(i) !将精确解数据写入dat文件end doD="er-" !er 代表error，这里指精确值和数值计算值之间的差别filename=D//B//Copen(unit=12,file=filename)write(12,*)"max-error=",error !将误差最大值写入dat文件write(*,*)n,errorwrite(*,*)filenamestopend!主程序结束subroutine subsolution(y,n,h) !子程序implicit noneinteger ::n,ireal(8) :: hreal(8) :: y(0:n) !数组和变量的声明不能同时进行real(8),allocatable :: a(:),b(:),c(:),s(:),t(:),f(:),g(:) !声明可变长度数组allocate(a(1:n),b(0:n),c(0:n-1),s(0:n),t(0:n-1),f(0:n),g(0:n)) !指定可变数组的长度a=1 !对数组a，b，c赋值b=-2c=1a(n)=0b(0)=1b(n)=1c(0)=0f(0)=0 !对数组f 赋值f(n)=1do i=1,n-1f(i)=h*h*sin(2.0*real(i)*h)end dos(0)=b(0) !计算s，tt(0)=c(0)/b(0)do i=1,n-1s(i)=b(i)-a(i)*t(i-1)t(i)=c(i)/s(i)end dos(n)=b(n)-a(n)*t(n-1)g(0)=f(0)/b(0) !追过程do i=1,ng(i)=(f(i)-a(i)*g(i-1))/s(i) end doy(n)=g(n) !赶过程do i=n-1,1,-1y(i)=g(i)-t(i)*y(i+1)end doy(0)=0y(n)=1returnend subroutine subsolution。