利用中心差分格式数值求解导数
差分格式
§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. FDM 1.1概念有限差分方法(FDM)是计算机数值模拟最早采用的方法,至今仍被广泛运用。
该方法将求解域划分为差分网格,用有限个网格节点代替连续的求解域。
有限差分法以Taylor级数展开等方法,把控制方程中的导数用网格节点上的函数值的差商代替进行离散,从而建立以网格节点上的值为未知数的代数方程组。
该方法是一种直接将微分问题变为代数问题的近似数值解法,数学概念直观,表达简单,是发展较早且比较成熟的数值方法。
1.2差分格式(1)从格式的精度来划分,有一阶格式、二阶格式和高阶格式。
(2)从差分的空间形式来考虑,可分为中心格式和逆风格式。
(3)考虑时间因子的影响,差分格式还可以分为显格式、隐格式、显隐交替格式等。
目前常见的差分格式,主要是上述几种形式的组合,不同的组合构成不同的差分格式。
差分方法主要适用于有结构网格,网格的步长一般根据实际地形的情况和柯朗稳定条件来决定。
1.3构造差分的方法构造差分的方法有多种形式,目前主要采用的是泰勒级数展开方法。
其基本的差分表达式主要有三种形式:一阶向前差分、一阶向后差分、一阶中心差分和二阶中心差分等,其中前两种格式为一阶计算精度,后两种格式为二阶计算精度。
通过对时间和空间这几种不同差分格式的组合,可以组合成不同的差分计算格式。
2. FEM 2.1概述有限元方法的基础是变分原理和加权余量法,其基本求解思想是把计算域划分为有限个互不重叠的单元,在每个单元内,选择一些合适的节点作为求解函数的插值点,将微分方程中的变量改写成由各变量或其导数的节点值与所选用的插值函数组成的线性表达式,借助于变分原理或加权余量法,将微分方程离散求解。
采用不同的权函数和插值函数形式,便构成不同的有限元方法。
2.2原理有限元方法最早应用于结构力学,后来随着计算机的发展慢慢用于流体力学、土力学的数值模拟。
在有限元方法中,把计算域离散剖分为有限个互不重叠且相互连接的单元,在每个单元内选择基函数,用单元基函数的线形组合来逼近单元中的真解,整个计算域上总体的基函数可以看为由每个单元基函数组成的,则整个计算域内的解可以看作是由所有单元上的近似解构成。
(完整版)有限差分方法概述
有限差分法(Finite Difference Method,简称FDM)是数值方法中最经典的方法,也是计算机数值模拟最早采用的方法,至今仍被广泛运用。
该方法将求解域划分为差分网格,用有限个网格节点代替连续的求解域。
有限差分法以Taylor级数展开等方法,把控制方程中的导数用网格节点上的函数值的差商代替进行离散,从而建立以网格节点上的值为未知数的代数方程组。
该方法是一种直接将微分问题变为代数问题的近似数值解法,数学概念直观,表达简单,是发展较早且比较成熟的数值方法。
对于有限差分格式,从格式的精度来划分,有一阶格式、二阶格式和高阶格式。
从差分的空间形式来考虑,可分为中心格式和逆风格式。
考虑时间因子的影响,差分格式还可以分为显格式、隐格式、显隐交替格式等。
目前常见的差分格式,主要是上述几种形式的组合,不同的组合构成不同的差分格式。
差分方法主要适用于有结构网格,网格的步长一般根据实际地形的情况和柯朗稳定条件来决定。
构造差分的方法有多种形式,目前主要采用的是泰勒级数展开方法。
其基本的差分表达式主要有三种形式:一阶向前差分、一阶向后差分、一阶中心差分和二阶中心差分等,其中前两种格式为一阶计算精度,后两种格式为二阶计算精度。
通过对时间和空间这几种不同差分格式的组合,可以组合成不同的差分计算格式。
下面我们从有限差分方法的基本思想、技术要点、应用步骤三个方面来深入了解一下有限差分方法。
1.基本思想有限差分算法的基本思想是把连续的定解区域用有限个离散点构成的网格来代替,这些离散点称作网格的节点;把连续定解区域上的连续变量的函数用在网格上定义的离散变量函数来近似;把原方程和定解条件中的微商用差商来近似,积分用积分和来近似,于是原微分方程和定解条件就近似地代之以代数方程组,即有限差分方程组,解此方程组就可以得到原问题在离散点上的近似解。
然后再利用插值方法便可以从离散解得到定解问题在整个区域上的近似解。
在采用数值计算方法求解偏微分方程时,再将每一处导数由有限差分近似公式替代,从而把求解偏微分方程的问题转换成求解代数方程的问题,即所谓的有限差分法。
求解广义Black-Scholes方程的指数型差分法
1 引 言
近年来 . 资者 对创 新产 品 的需求 1 投 3益增 加 , 证券 交 易所 之 间在 产 品开 发上 的竞 争 1趋激 烈 , 国 3 各 ( 区) 地 竞相推 出种类 繁多 的期 权产 品 。 权是 指买方 在支付 了期权 费后取得 的在 约定 的到期 1或 到期 1 期 3 3 之前 以约定 的价格买 人或卖 出一定 数量某种 标 的资产 的权 利 。 只能在选 择权 到期 日当天履 约的称为 欧式期 权 ; 而能 够在选 择权 到期 日当天或 到期前 任何一 天履 约 的, 称之为美 式期权 。 基于单 一标 的资产价 格 的欧式 期权 , 价格 满足一个 关 于时间和 标 的资产价格 的二 其 阶抛物方 程 . 称其 为 Ba k S h l 方程 。对于标 准 的欧 式期权 , l k和 S h ls早 已给 出了解析形 式 的 lc — c oe s Ba c c oet 】 定 价公式 。然而还有 很多 复杂 的期 权类 型 , 它们涉 及到各 种动态 随机 变化 的参 数 、 种标 的物 、 多 复杂 的收 益 函数 以及 提前 执行可能 等 , 很难 得到 它们 的解 析解 。 因此 , 发展 各种计 算期权 价格 的数值计 算方法 具有
域 ( ,∞) 0 + 将变换 成 ( ∞, ∞) 一 + 。因此 , 当求解变换 后 的问题时 , 数值求 解必须 引人 两个额外 的参数 来将无
穷 区域截 断成有 限 区间。 在数值计 算 中 , 截断 右边无穷 区域成 有限 区间是不 可避免 的 , 是截断 左边无穷 但 区域来 人为 的除去奇性 , 这将会 引起 数值误差 。 另外 , 在变换 后 的定 义 区域上 使用 等距 网格将导 致原来 的
(.) 26 的解 ( ,) 足如下估 计 t满
二维热传导问题的时间中心差分法求解
二维热传导问题的时间中心差分法求解下载提示:该文档是本店铺精心编制而成的,希望大家下载后,能够帮助大家解决实际问题。
文档下载后可定制修改,请根据实际需要进行调整和使用,谢谢!本店铺为大家提供各种类型的实用资料,如教育随笔、日记赏析、句子摘抄、古诗大全、经典美文、话题作文、工作总结、词语解析、文案摘录、其他资料等等,想了解不同资料格式和写法,敬请关注!Download tips: This document is carefully compiled by this editor. I hope that after you download it, it can help you solve practical problems. The document can be customized and modified after downloading, please adjust and use it according to actual needs, thank you! In addition, this shop provides you with various types of practical materials, such as educational essays, diary appreciation, sentence excerpts, ancient poems, classic articles, topic composition, work summary, word parsing, copy excerpts, other materials and so on, want to know different data formats and writing methods, please pay attention!二维热传导问题的时间中心差分法求解1. 引言热传导问题是工程领域中常见的研究课题,涉及到如何有效地计算和模拟热量在空间中的传播和分布。
数值传热学部分习题问题详解2
习题4-2一维稳态导热问题的控制方程:022=+∂∂S xTλ 依据本题给定条件,对节点2采用二阶精度的中心差分格式,节点3采用第三类边界条件具有二阶精度的差分格式,最后得到各节点的离散方程: 节点1: 1001=T节点2: 1505105321-=+-T T T 节点3:75432=+-T T求解结果:852=T ,403=T对整个控制容积作能量平衡,有:02150)4020(15)(3=⨯--⨯=∆+-=∆+x S T T h x S q f f B即:计算区域总体守恒要求满足习题4-5在4-2习题中,如果25.03)(10f T T h -⨯=,则各节点离散方程如下:节点1: 1001=T节点2: 1505105321-=+-T T T节点3:25.03325.032)20(4015])20(21[-⨯+=-⨯++-T T T T对于节点3中的相关项作局部线性化处理,然后迭代计算; 求解结果:818.822=T ,635.353=T (迭代精度为10-4)迭代计算的Matlab 程序如下:x=30; x1=20;while abs(x1-x)>0.0001a=[1 0 0;5 -10 5;0 -1 1+2*(x-20)^(0.25)]; b=[100;-150; 15+40*(x-20)^(0.25)]; t=a^(-1)*b; x1=x; x=t(3,1); end tcal=t习题4-14充分发展区的温度控制方程如下:)(1rTr r r x T uc p ∂∂∂∂=∂∂λρ 对于三种无量纲定义w b w T T T T --=Θ、∞∞--=ΘT T T T w 、w w T T T T --=Θ∞进行分析如下1)由wb wT T T T --=Θ得:w w b T T T T +Θ-=)(由T 可得:x T x T x T T T x T w b w w b ∂∂Θ-+∂∂Θ=∂+Θ-∂=∂∂)1(])[(rT r T T r T T T r T w w b w w b ∂∂Θ-+∂Θ∂-=∂+Θ-∂=∂∂)1()(])[( 由b T 与r 无关、Θ与x 无关以及x T ∂∂、rT∂∂的表达式可知,除了w T 均匀的情况外,该无量纲温度定义在一般情况下是不能用分离变量法的; 2)由∞∞--=ΘT T T T w 得:∞∞+Θ-=T T T T w )(由T 可得:xT x T T T x T w w ∂∂Θ=∂+Θ-∂=∂∂∞∞])[(rT r T T r T T T r T w w w ∂∂Θ+∂Θ∂-=∂+Θ-∂=∂∂∞∞∞)(])[( 由b T 与r 无关、Θ与x 无关以及x T ∂∂、rT∂∂的表达式可知,在常见的四种边界条件中除了轴向及周向均匀热流const q w =的情况外,有0=∂∂rT w离变量法的;3)由wwT T T T --=Θ∞得:w w T T T T +Θ-=∞)(由T 可得: xT x T T T x T w w w ∂∂Θ-=∂+Θ-∂=∂∂∞)1(])[(rT r T T r T T T r T w w w w ∂∂Θ-+∂Θ∂-=∂+Θ-∂=∂∂∞∞)1()(])[( 同2)分析可知,除了轴向及周向均匀热流const q w =的情况外,有0=∂∂rT w,该无量纲温度定义是可以用分离变量法的;习题4-181)采用柱坐标分析,写出统一的稳态柱坐标形式动量方程:S r r r r r r x x w r v r r r u x +∂∂∂∂+∂∂∂∂+∂∂∂∂=∂∂+∂∂+∂∂)(1)(1)()(1)(1)(θφλθφλφλφρθφρφρ x 、r 和θ分别是圆柱坐标的3个坐标轴,u 、v 和w 分别是其对应的速度分量,其中x 是管内的流动方向;对于管内的层流充分发展有:0=v 、0=w ,0=∂∂xu;并且x 方向的源项:x p S ∂∂-= r 方向的源项:r pS ∂∂-=θ方向的源项:θ∂∂-=pr S 1 由以上分析可得到圆柱坐标下的动量方程: x 方向: 0)(1)(1=∂∂-∂∂∂∂+∂∂∂∂x pu r r r u r r r θλθλ r 方向:0=∂∂r pθ方向:0=∂∂θp边界条件: R r =,0=u0=r ,0=∂∂r u ;对称线上,0=∂∂θu 不考虑液体的轴向导热,并简化分析可以得到充分发展的能量方程为:)(1)(1θλθλρ∂∂∂∂+∂∂∂∂=∂∂Tr r r T r r r x T uc p 边界条件: R r =,w q r T =∂∂λ;0=r ,0=∂∂rTπθ/0=,0=∂∂-θλT2)定义无量纲流速:dxdp R uU 2-=λ并定义无量纲半径:R r /=η;将无量纲流速和无量纲半径代入x 方向的动量方程得:0))1((1))1((122=∂∂-∂-∂∂∂+∂-∂∂∂xp U dx dp R R R R U dx dp R RR R θληλθηηλληηη 上式化简得:01)1(1)(1=+∂∂∂∂+∂∂∂∂θηθηηηηηU U 边界条件:1=η,0=U0=η,0=∂∂ηU ;对称线上,0=∂∂θU定义无量纲温度:λ/0R q T T b-=Θ其中,0q 是折算到管壁表面上的平均热流密度,即:Rq q wπ=0; 由无量纲温度定义可得: b T Rq T +Θ=λ将T 表达式和无量纲半径η代入能量方程得:)(1)(100θληλθηηλληηηρ∂Θ∂∂∂+∂Θ∂∂∂=∂∂R q R R R R q R R R x T uc b p 化简得:)1(1)(10θηθηηηηηρ∂Θ∂∂∂+∂Θ∂∂∂=∂∂x T u c q R b p (1)由热平衡条件关系可以得:mm m b m p b p p RU U q R u u R q A u u dx dT A u c x T u c x T uc 020221221)(===∂∂=∂∂ππρρρ 将上式代入式(1)可得:)1(1)(12θηθηηηηη∂Θ∂∂∂+∂Θ∂∂∂=m U U 边界条件: 0=η,0=∂Θ∂η;1=η,R q q w πη10==∂Θ∂0=θ,0=∂Θ∂θ;πθ=,0=∂Θ∂θ单值条件:由定义可知:0/0=-=ΘλR q T T b b b 且: ⎰⎰Θ=ΘAAb UdAUdA即得单值性条件:0=Θ⎰⎰AA UdAUdA 3)由阻力系数f 及Re 定义有:228)(21/Re ⎪⎭⎫ ⎝⎛=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡-=D D U D u u dx dp D f e m e m me νρ 且:m W b m W b m W R q T T D T T q Nu ,0,,0~2)/(2Θ=-=-=λλ5-21.一维稳态无源项的对流-扩散方程如下所示:xx u 22∂∂Γ=∂∂φφρ (取常物性)边界条件如下:L L x x φφφφ====,;,00上述方程的精确解如下:11)/(00--=--⋅PeL x Pe L e e φφφφ Γ=/uL Pe ρ 2.将L 分成20等份,所以有:∆=Pe 20图示如下:1 2 3 4 5 6 ………… …………… 17 18 19 20 21 对于中心差分、一阶迎风、混合格式和QUICK 格式分别分析如下: 1) 中心差分中间节点: 2)5.01()5.01(11-∆+∆++-=i i i P P φφφ 20,2Λ=i2) 一阶迎风中间节点: ∆-∆++++=P P i i i 2)1(11φφφ 20,2Λ=i3) 混合格式当1=∆P 时,中间节点:2)5.01()5.01(11-∆+∆++-=i i i P P φφφ20,2Λ=i当10,5=∆P 时,中间节点: 1-=i i φφ 20,2Λ=i 4) QUICK 格式*12111)35(8122121⎥⎦⎤⎢⎣⎡---++++++=+--∆∆-∆∆+∆i i i i i i i P P P P P φφφφφφφ 2≠i *1111)336(8122121⎥⎦⎤⎢⎣⎡--++++++=+-∆∆-∆∆+∆i i i i i i P P P P P φφφφφφ 2=i 5-3乘方格式:⎪⎪⎩⎪⎪⎨⎧<-≤≤--+≤≤->=∆∆∆∆∆∆∆∆10,010,)1.01(100,)1.01(10,055P P P P P P P P D a e E当1.0=∆P 时有:951.0)1.01.01()1.01(55=⨯-=-=∆P D a eE因为:301.0/3)()()()()()(===Γ=Γ=∆eee e e e e e e P u x u u x D ρδρρδ 所以:5297.2830951.0951.0=⨯==e E D a由系数关系式∆=-P D a D a eEw W 可得: 53.3130)951.01.0()(=⨯+=⨯+=∆w eEW D D a P a 且: 205.01.010=⨯=∆∆=txa P p ρ 当采用隐式时1=f ,因此可得:0597.62253.315297.280=++=++=P W E P a fa fa a同理可得当10=∆P 时有:0=E a ,3=W a ,5=P a5-5二维稳态无源项的对流-扩散问题的控制方程:)()()()(yy x x y v x u ∂∂Γ∂∂+∂∂Γ∂∂=∂∂+∂∂φφφρφρφφ 对于一阶迎风、混合、乘方格式的通用离散方程:S S N N W W E E P P a a a a a φφφφφ+++=其中:[]0,)(e e e E F P A D a -+=∆ []0,)(w w w W F P A D a +=∆ []0,)(n n n N F P A D a -+=∆ []0,)(s s s S F P A D a +=∆5-71)QUICK 格式的界面值定义如下:⎪⎪⎩⎪⎪⎨⎧-+=-+=)36(81)36(81WW P W w W E P e φφφφφφφφ0>u 对(5-1)式dxdx d d dx u d )()(φφρΓ=积分可得: w e w e dxd dx d u u )()()()(φφφρφρΓ-Γ=-对流项采用QUICK 格式的界面插值,扩散项采用线性界面插值,对于0>u 及均分网格有:)]()([]))(36())(36[(81x x u u W P w P E e w WW P W e W E P ∆-Γ-∆-Γ=-+--+φφφφρφφφρφφφ 整理得:WW w W w e w E e e P w e w e u u u x u x x x u u φρφρρφρφρρ)(81])(43)(81[])(83[)]()(83)(43[-++∆Γ+-∆Γ=∆Γ+∆Γ+-上式即为QUICK 格式离散得到的离散方程;2)要分析QUICK 格式的稳定性,则应考虑非稳平流方程:xut ∂∂-=∂∂φφ 在t ∆时间间隔内对控制容积作积分:⎰⎰⎰⎰∆+∆+∂∂-=∂∂t t t e w e w tt tdxdt x u dtdx xφφ得:dt u dx tt tw e e wttt ⎰⎰∆+∆+--=-)()(φφφφφ随时间变化采用阶梯显式,随空间变化采用QUICK 格式得:t u x WW P W W E P tP t t P ∆+---+-=∆-∆+)]3636(81[)(φφφφφφφφ整理得:xu t ni n i n i n i ni n i ∆+-+-∆---++87332111φφφφφφ对于初始均匀零场,假设在),(n i 点有一个扰动n i ε; 对1+i 点写出QUICK 格式的离散方程:xu tni n i n i n i n i n i ∆+-+-∆--+++++8733121111φφφφφφ可得:ni n i xt u εφ∆∆=++8711 对1-i 点分析可得:ni n i xt u εφ∆∆-=+-8311 由于扩散对扰动的传递恒为正,其值为ni x t ερ2∆Γ∆,所以根据符号不变原则有: 0)/)83(2≥∆Γ∆+∆∆-ni n i n i xt x t u εερε 整理得到QUICK 格式的稳定性条件为:38≤∆P 5-91)三阶迎风格式采用上游两个节点和下游一个节点的值来构造函数界面插值形式,所以定义如下:⎩⎨⎧<++=>++=00u c b a u c b a EEE P e W P E e φφφφφφφφ根据上述定义,在0>u 时对控制容积内的对流项作积分平均可得:])()([1)(11WW W P E e w w e c b c a b a xxdx x x φφφφφφφ--+-+∆=-∆=∂∂∆⎰由表2-1式可知三阶迎风格式的差分格式:xxni n i n i n i ni ∆+-+=∂∂--+1221264211,φφφφφ 由控制容积积分法得到的对流项离散格式应与Taylor 离散展开得到的离散格式具有相同的形式和精度,所以比较可得:61,65,31-===c b a所以三阶迎风格式的函数插值定义为:⎪⎪⎩⎪⎪⎨⎧<-+=>-+=06165310616531u u EE E P e W P E e φφφφφφφφ2)由上述分析可知,得到的三阶迎风格式的插值定义与给出节点上导数表达式的定义在形式上显然是一致的;6-1二维直角坐标中不可压缩流体的连续方程及动量方程如下:⎪⎪⎪⎪⎪⎩⎪⎪⎪⎪⎪⎨⎧+∂∂∂∂+∂∂∂∂+∂∂-=∂∂+∂∂+∂∂+∂∂∂∂+∂∂∂∂+∂∂-=∂∂+∂∂+∂∂=∂∂+∂∂)3()()()()()()2()()()()()()1(0vu S y y v x x v y p y vv x vu t v S y y u x x u x p y uv x uu tu y v x u ηηρρρηηρρρ假设常粘性,则0==v u S S ;对公式(2)及(3)分别对y x ,求偏导得:⎪⎪⎩⎪⎪⎨⎧∂∂+⎪⎪⎭⎫ ⎝⎛∂∂∂∂+⎪⎪⎭⎫ ⎝⎛∂∂∂∂-=⎪⎪⎭⎫ ⎝⎛∂∂∂∂+⎪⎭⎫ ⎝⎛∂∂∂∂+⎪⎭⎫ ⎝⎛∂∂∂∂⎪⎪⎭⎫⎝⎛∂∂∂∂+∂∂+⎪⎭⎫ ⎝⎛∂∂∂∂-=⎪⎪⎭⎫ ⎝⎛∂∂∂∂+⎪⎭⎫ ⎝⎛∂∂∂∂+⎪⎭⎫ ⎝⎛∂∂∂∂33222233)()()()()()(y v x v y y p y y vv y x vu y t v y y u x x u x p x y uv x x uu x t u x ηηρρρηηρρρ 两式相加得并变换积分顺序有:⎥⎦⎤⎢⎣⎡⎪⎪⎭⎫ ⎝⎛∂∂+∂∂∂∂+⎥⎦⎤⎢⎣⎡⎪⎪⎭⎫ ⎝⎛∂∂+∂∂∂∂+⎪⎪⎭⎫⎝⎛∂∂+∂∂-=⎥⎦⎤⎢⎣⎡⎪⎪⎭⎫ ⎝⎛∂∂+∂∂+∂∂∂∂+⎪⎪⎭⎫ ⎝⎛∂∂+∂∂+∂∂∂∂+⎪⎪⎭⎫ ⎝⎛∂∂+∂∂∂∂y v x u yy v x u xy p x p x v u x u v y v v y y v u y u v x u u x y v x u t 2222222222ηηρρ利用连续方程有:⎪⎪⎭⎫ ⎝⎛∂∂+∂∂-=⎥⎦⎤⎢⎣⎡⎪⎪⎭⎫ ⎝⎛∂∂+∂∂∂∂+⎪⎪⎭⎫ ⎝⎛∂∂+∂∂∂∂2222y p x p x v u y v v y y u v x u u x ρ ⎪⎪⎭⎫ ⎝⎛∂∂+∂∂-=⎥⎦⎤⎢⎣⎡∂∂∂∂-∂∂∂∂+∂∂+∂∂+∂∂∂∂22222222222y p xp y v x u y v x u y v x u x v y u ρ 最后即得:⎥⎦⎤⎢⎣⎡∂∂∂∂-∂∂∂∂=⎪⎪⎭⎫ ⎝⎛∂∂+∂∂x v y u y v x u y p x p ρ222226-4假设5*=P p ,则有:5105*-=-=e u 5.3)05(7.0*=-⨯=n v由连续性条件有:s w n e v u v u +=+按SIMPLE 算法有:'''*5)(P E P e e e p p p d u u +-=-+= '''*7.05.3)(P n P n n n p p p d v v +=-+=将上两式代入连续性方程中有:20507.05.35''+=+++-P P p p计算得:06.42'=P p所以:06.4706.425'*=+=+=P P P p p p06.371006.47=-=-=E P e p p u 94.32)006.47(7.0)(7.0=-⨯=-=N P n p p v6-5假设250*3=p ,150*6=p ,所以各点的流量为:⎪⎪⎪⎩⎪⎪⎪⎨⎧-=-⨯==-⨯=-=-⨯=-=-⨯==-⨯=11)15040(1.020)150250(2.024)25010(1.04)270250(2.010)250275(4.0*****E DC B A Q Q Q Q Q 上述流量满足动量方程,但并不满足连续性方程,所以对流量修正:⎪⎪⎪⎩⎪⎪⎪⎨⎧-⨯+-=-⨯+=-⨯+-=-⨯+-=-⨯+=)(1.011)(2.020)(1.024)(2.04)(4.010'6'5'6'3'3'4'2'3'3'1p p Q p p Q p p Q p p Q p p Q ED C B A 对节点3作质量守恒有:B DC A Q Q Q Q +=+即得:)(2.04)(2.020)(1.024)(4.010'2'3'6'3'3'4'3'1p p p p p p p p -⨯+--⨯+=-⨯+--⨯+对节点3作质量守恒有:F E D Q Q Q =+即得:20)(1.011)(2.020'6'5'6'3=-⨯+--⨯+p p p p联立求解上两式有:70.48'3-=p ,13.69'6-=p修正后的压力为:3.20170.48250'3*33=-=+=p p p 87.8013.69150'6*66=-=+=p p p修正后的流量为:⎪⎪⎪⎩⎪⎪⎪⎨⎧-=-⨯==-⨯=-=-⨯=-=-⨯==-⨯=09.4)87.8040(1.009.24)87.803.201(2.013.19)3.20110(1.074.13)2703.201(2.048.29)3.201275(4.0ED C B A Q Q Q Q Q由)(76p p C Q F F -=。
有限差分求解1维流体方程
有限差分求解1维流体方程全文共四篇示例,供读者参考第一篇示例:有限差分法是求解偏微分方程的一种常用数值方法,它将连续的微分方程离散化为有限个点上的代数方程组,通过对这些代数方程进行求解得到数值解。
在流体力学中,有限差分法被广泛应用于求解流体方程,其中涉及了流体的运动和力学性质。
本文将着重介绍有限差分法求解一维流体方程的方法和步骤。
一维流体方程是描述流体在一维空间中运动的数学模型,通常可以描述为一维流体方程组。
有限差分法的基本思想是将一维流体方程中的时间和空间进行离散化,将连续的一维空间划分为有限个网格点,时间也进行离散化为有限个时间步长,通过有限差分近似代替微分算子,并在每个网格点上建立代数方程组,最终通过求解这些代数方程组得到数值解。
具体来说,有限差分法求解一维流体方程的步骤如下:1. 确定求解区域和边界条件:首先需要确定求解区域的大小和边界条件,包括流体的初始状态和边界条件。
这些信息将决定了求解的范围和边界条件的设定。
2. 离散化:将一维空间和时间进行离散化,将空间和时间分别划分为有限个网格点和时间步长。
这一步是有限差分法的核心思想,通过离散化将连续的微分方程转化为离散的代数方程。
3. 近似微分算子:在每个网格点上近似代替微分算子,例如将一维流体方程中的导数项用差分近似代替,通常采用中心差分、前向差分或后向差分等方式。
通过这种方式可以将微分方程转化为代数方程。
4. 建立代数方程组:在每个网格点上建立代数方程组,将近似的微分算子代入原始方程中,然后利用相邻网格点之间的关系建立代数方程,得到形如Ax=b的线性方程组。
5. 求解线性方程组:通过求解线性方程组得到数值解,这可以采用各种求解线性方程组的方法,如直接法、迭代法或者共轭梯度法等。
6. 可视化和分析:最后通过对数值解进行可视化和分析,可以得到流体在一维空间中的运动情况和物理性质,例如流速、压强等。
第二篇示例:有限差分法是一种常用的数值计算方法,常用于求解偏微分方程。
差分法基本原理范文
差分法基本原理范文差分法是一种数值求解微分方程的方法,它将微分方程转化为离散形式,通过离散点上的函数值之差来近似计算导数,从而得到微分方程的数值解。
差分法的基本原理可以总结为以下几个步骤:1.网格划分:将求解区域划分为若干有限个小区域,每个小区域称为一个网格单元。
通常采用均匀网格划分,将区域划分为有限个等距的小区域。
2.网格节点:在每个网格单元的边界上选择一个或多个节点,节点是离散区域内的点。
节点数量取决于所选择的差分格式。
通常要求节点密度足够高,以确保数值解的精度。
3.差分逼近:使用差分公式对微分方程中的导数进行近似。
差分公式的选择是差分法的核心。
常见的差分公式包括:中心差分、向前差分、向后差分等等。
不同的差分公式对应着不同的差分格式,如前向差分格式、后向差分格式、中心差分格式等等。
4.离散化方程:根据差分逼近的结果,将微分方程中的导数用离散点上的函数值之差来代替,得到离散的差分方程。
离散化的过程将微分方程转化为代数方程组,可以通过求解代数方程组来得到数值解。
5.边界条件:确定边界条件,在差分方程中将其加入到方程组中。
边界条件通常是指在求解区域边界上的已知函数值或导数值。
6.求解代数方程组:根据离散化方程和边界条件,得到一个代数方程组。
通过数值方法,如高斯消元法、迭代法等,求解得到方程组的数值解。
7.结果输出和误差分析:根据求解得到的数值解,可以进行结果输出、误差分析和收敛性研究。
通常需要对数值解进行采样,与解析解进行比较,以评估差分法的精度和稳定性。
差分法的优点包括:简单易用,计算效率高,适用于各种类型的微分方程,比如常微分方程、偏微分方程及椭圆、抛物、双曲型方程等等。
然而,差分法也存在一些限制,如数值方法的稳定性与精度受节点密度、步长选择和差分格式的影响。
因此,在实际应用中,需要根据具体问题的特点和求解要求选择合适的差分格式和参数,以获得满足要求的数值解。
一种椭圆形方程的差分格式及迭代法求解
一种椭圆形方程的差分格式及迭代法求解作者:于红张宏来源:《新课程学习·下》2015年第03期摘要:有限差分法是解偏微分方程的一个重要数值方法。
对正方形域上的Laplace方程的第一边值问题用差分法建立了其差分格式,并用Jacobi迭代法、Gauss-Seidel迭代法和超松弛迭代法(SOR法)对该差分格式进行求解。
对三种迭代法进行编程并上机实践,求得相应数值解,通过表格对运行结果进行了比较。
关键词:差分格式;Jacobi迭代法;Gauss-seidel迭代法;超松弛迭代法一、问题及其差分格式考虑正方形域上的Laplace方程的第一边值问题:uxx+uyy=0,(0容易验证方程(1)的精确解为:u(x,y)=■sinπy.在xoy平面上作两组平行直线:x=x0+ih,y=y0+jh(i,j=0,±1,±2…).(x0,y0)是xoy平面上的任意一点,取(x0,y0)为坐标原点(0,0),步长h=■,这样,整个平面就被这两组平行直线构成的正方形网格所覆盖,两组平行直线的交点称为网格结点,只考虑属于正方形区域[0,1;0,1]的结点。
若一个结点的四个相邻接点都属于[0,1;0,1],则称此结点为内部结点;若一结点的四个相邻结点至少有一个不属于[0,1;0,1],则称此结点为边界结点。
在每一个内部结点上,用二阶中心差代替问题中的二阶导数:(uxx)■=(■)■≈■=■(uyy)■=(■)■≈■=■则有:(■)■+(■)■≈■由此得到(1)的差分格式为:uij=■(ui+1, j+ui-1, j+ui, j+1+ui, j-1).(i,j=1,2,…n-1),其中u(0,y)=u(x,0)=u(x,1)=0,u(1,y)=sinπy.二、三种迭代方法及收敛性比较对线性方程组Ax=b,系数矩阵A=(aij)n×n非奇异,且A的主对角元aij≠0(i=1,2…n).1.Jacobi迭代法将A分裂成A=D-(D-A),其中D=diag(a11,a22…ann),于是方程Ax=b可以写成Dx=(D-A)x+b或x=(E-D-1A)x+D-1b (2)令B=E-D-1A,g=D-1b,则(2)可写成:x=Bx+g这样得到了迭代公式:xk=Bxk-1+g(k=1,2…)(3)2.Gauss-Seidel迭代法将A分裂成:A=D(E-L)-DU其中:D=diag(a11,a22,…ann),L=0 0 … 0 0-■ 0 … 0 0-■ -■ … 0 0 … … …-■ -■ … -■ 0U=0 -■ … -■ -■0 0 … -■ -■ … … …0 0 … 0 -■0 0 … 0 0对比(2)和(3)则得出:B=L+U于是方程Ax=b可以写成:D(I-L)x=DUx+b (4)显然D和I-L都是非奇异的,因此可以用(I-L)-1D-1左乘上式的两端,得出:x=(I-L)-1+Ux+(I-L)-1D-1b由此得Gauss-Seidel迭代法的迭代公式:■由此可见Gauss-Seidel迭代法是Jacobi迭代法的修正。
解四阶抛物型方程的两层显式差分格式
解四阶抛物型方程的两层显式差分格式四阶抛物型方程是指具有四个导数项的抛物型偏微分方程,可以写成如下形式:u_t = a*u_xx + b*u_xxx + c*u_xxxx + f(x,t)其中,u表示未知函数,t表示时间,x表示空间,a、b、c为系数,f(x,t)为已知的源项函数。
为了数值求解这类方程,我们可以使用显式差分格式。
显式差分格式是指通过将方程中的导数项用差分运算进行离散化,将连续的偏微分方程转化为离散的差分方程。
两层显式差分格式指使用两个时间层次的差分方程进行迭代求解。
下面我们将介绍两种常用的两层显式差分格式:双边五点差分格式和五点中心差分格式。
1.双边五点差分格式(BDF5)双边五点差分格式采用五点差分近似导数,其中时间层次的差分使用五阶向前差分,空间层次的差分使用五阶中心差分,可以得到如下差分方程:(u_i^(n+1)-u_i^n)/Δt=a*(u_{i-2}^n-4u_{i-1}^n+6u_i^n-4u_{i+1}^n+u_{i+2}^n)/(Δx^2)+b*(u_{i-2}^n-2u_{i-1}^n+2u_{i+1}^n-u_{i+2}^n)/(2Δx^2)+c*(u_{i-2}^n-u_{i-1}^n-u_{i+1}^n+u_{i+2}^n)/(Δx^2)+f_i^n其中,i表示空间格点的索引,n表示时间层次的索引,Δt和Δx 分别表示时间和空间的步长,u_i^n表示在第n个时间层次上的第i个空间点的解,f_i^n表示在第n个时间层次上的第i个空间点的源项。
2.五点中心差分格式(CD5)五点中心差分格式采用五点差分近似导数,其中时间层次的差分使用五阶前后向差分,空间层次的差分使用五阶中心差分,可以得到如下差分方程:(u_i^(n+1)-u_i^(n-1))/(2Δt)=a*(u_{i-2}^n-4u_{i-1}^n+6u_i^n-4u_{i+1}^n+u_{i+2}^n)/(Δx^2)+b*(u_{i-2}^n-2u_{i-1}^n+2u_{i+1}^n-u_{i+2}^n)/(2Δx^2)+c*(u_{i-2}^n-u_{i-1}^n-u_{i+1}^n+u_{i+2}^n)/(Δx^2)+f_i^n这两个差分方程可以通过逐步迭代求解,用现有的时间层次上的解来计算下一个时间层次上的解。
一维热传导方程的差分法
一维热传导方程的差分法1. 引言1.1 介绍一维热传导方程的差分法一维热传导方程是描述物体内部温度分布随时间变化的数学模型。
差分法是一种常用的数值解法,通过将时间和空间进行离散化,将偏微分方程转化为差分方程,从而可以通过计算机进行数值求解。
在一维热传导方程的差分法中,我们通常将时间和空间分别进行离散化,将连续的温度变化转化为离散的温度值。
通过迭代计算,可以得到物体内各个离散点的温度随时间的变化情况。
差分法的优点在于可以较好地模拟物体内部温度分布的变化,同时可以较快地得到数值解,对于复杂的边界条件和非线性问题也有较好的适用性。
通过研究一维热传导方程的差分法,可以更好地理解物体内部温度分布的变化规律,为工程实践提供有效的数值模拟手段。
同时也可以探讨数值解法的稳定性和收敛性,为进一步的数值模拟研究提供参考。
通过不断改进差分法的算法和技术,可以更准确地预测物体内部温度变化,为工程设计和科学研究提供有力支持。
1.2 研究背景一维热传导方程是描述热量在一维空间内传递和分布的数学模型,广泛应用于工程领域和物理学中。
研究热传导方程的差分法是为了解决实际问题中复杂边界条件和非线性情况下的热传导问题,以及对传热过程进行数值模拟和分析。
在工程实践中,热传导问题经常出现在各种材料的传热过程中,例如石油钻井中地下油层的温度分布、金属材料的焊接过程中的温度控制等。
研究热传导方程的差分法可以帮助工程师们更好地理解热传导过程,优化工程设计,提高生产效率。
研究热传导方程的差分法还可以为其他科学领域提供理论支持和数值计算方法。
在地质学中用于模拟地热传导过程、在气象学中用于模拟大气环流等。
深入研究一维热传导方程的差分法对于推动科学研究和解决实际问题具有重要意义。
1.3 研究目的研究目的是通过对一维热传导方程的差分法进行深入分析和研究,探索其在实际工程和科学问题中的应用潜力。
具体来说,我们的研究目的包括以下几个方面:我们希望能够建立一种有效的数学模型,用以描述和解决一维热传导问题,为实际问题的数值模拟提供理论基础。
matlab差分法解偏微分方程
Matlab 差分法解偏微分方程1.引言解偏微分方程是数学和工程领域中的一项重要课题,它在科学研究和工程实践中具有广泛的应用。
而 Matlab 差分法是一种常用的数值方法,用于求解偏微分方程。
本文将介绍 Matlab 差分法在解偏微分方程中的应用,包括原理、步骤和实例。
2. Matlab 差分法原理差分法是一种离散化求解微分方程的方法,通过近似替代微分项来求解微分方程的数值解。
在 Matlab 中,差分法可以通过有限差分法或者差分格式来实现。
有限差分法将微分方程中的导数用有限差分替代,而差分格式指的是使用不同的差分格式来近似微分方程中的各个项,通常包括前向差分、后向差分和中心差分等。
3. Matlab 差分法步骤使用 Matlab 差分法解偏微分方程一般包括以下步骤:(1)建立离散化的区域:将求解区域离散化为网格点或节点,并确定网格间距。
(2)建立离散化的时间步长:对于时间相关的偏微分方程,需要建立离散化的时间步长。
(3)建立离散化的微分方程:使用差分法将偏微分方程中的微分项转化为离散形式。
(4)建立迭代方程:根据离散化的微分方程建立迭代方程,求解数值解。
(5)编写 Matlab 代码:根据建立的迭代方程编写 Matlab 代码求解数值解。
(6)求解并分析结果:使用 Matlab 对建立的代码进行求解,并对结果进行分析和后处理。
4. Matlab 差分法解偏微分方程实例假设我们要使用 Matlab 差分法解决以下一维热传导方程:∂u/∂t = α * ∂^2u/∂x^2其中 u(x, t) 是热传导方程的温度分布,α 是热扩散系数。
4.1. 离散化区域和时间步长我们将求解区域离散化为网格点,分别为 x_i,i=1,2,...,N。
时间步长为Δt。
4.2. 离散化的微分方程使用中心差分格式将偏微分方程中的导数项离散化得到:∂u/∂t ≈ (u_i(t+Δt) - u_i(t))/Δt∂^2u/∂x^2 ≈ (u_i-1(t) - 2u_i(t) + u_i+1(t))/(Δx)^2代入原偏微分方程可得离散化的微分方程:(u_i(t+Δt) - u_i(t))/Δt = α * (u_i-1(t) - 2u_i(t) + u_i+1(t))/(Δx)^24.3. 建立迭代方程根据离散化的微分方程建立迭代方程:u_i(t+Δt) = u_i(t) + α * Δt * (u_i-1(t) - 2u_i(t) + u_i+1(t))/(Δx)^24.4. 编写 Matlab 代码使用以上建立的迭代方程编写 Matlab 代码求解热传导方程。
二阶导数的四阶精度格式的推导
二阶导数的四阶精度格式的推导一、概述1. 二阶导数在数学和工程领域中广泛应用,而对于高精度计算来说,需要设计具有较高精度的数值格式来求解二阶导数。
本文将对二阶导数的四阶精度格式进行推导,以提供在实际计算中的参考。
二、推导过程2. 假设要计算函数f(x)在点x处的二阶导数,我们可以采用离散化的方式,即通过一定步长h对x进行逼近,然后利用差分格式得到二阶导数的近似值。
3. 我们首先利用中心差分格式求解一阶导数,得到:f'(x) ≈ (f(x+h) - f(x-h)) / (2h)4. 再对一阶导数利用中心差分格式求解二阶导数,得到:f''(x) ≈ (f'(x+h) - f'(x-h)) / (2h)≈ (f(x+2h) - 2f(x) + f(x-2h)) / (4h^2)5. 上述是常用的二阶导数的离散格式,其精度为二阶。
为了提高精度,我们可以利用更高阶的差分格式来逼近二阶导数。
6. 假设我们利用五点差分格式进行逼近,即利用前后各两个点进行逼近,得到:f''(x) ≈ (-f(x+2h) + 16f(x+h) - 30f(x) + 16f(x-h) - f(x-2h)) / (12h^2)7. 上述格式即为二阶导数的四阶精度格式。
通过将更多近邻点的信息纳入差分逼近中,可以提高逼近的精度。
三、数值实例8. 为了验证四阶精度格式的有效性,我们可以通过数值实例来进行对比。
假设我们要计算函数f(x)=sin(x)在x=1处的二阶导数。
9. 利用标准的中心差分格式进行计算,取步长h=0.1:f''(1) ≈ (sin(1+0.1) - 2sin(1) + sin(1-0.1)) / (0.1)^2≈ (0.453 - 2*0.841 + 1.411) / 0.01≈ -0.53610. 利用四阶精度格式进行计算,取步长h=0.1:f''(1) ≈(-sin(1+0.2) + 16sin(1+0.1) - 30sin(1) + 16sin(1-0.1) - sin(1-0.2)) / (12*0.1)^2≈ (-0.818 + 16*0.891 - 30*0.841 + 16*0.782 - 0.745) / 0.xxx≈ -0.84111. 通过对比可以看出,利用四阶精度格式得到的结果更加接近真实值sin''(1)=-sin(1)=-0.841的精度更高,验证了四阶精度格式的有效性。
- 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。