一维热传导方程的差分格式

1、偏微分方程数值解课程论文微分方程数值解课程论文偏微分方程数值解课程论文偏微分方程数值解课程论文学生姓名1:许慧卿学生姓名2：向裕学生姓名3：邱文林学生姓名4：高俊学生姓名5：赵禹恒学生姓名6：刘志刚学院：学号：20144329学号：20144327学号：20144349学号：20144305学号：20144359学号：20144346理学院偏微分方程数值解课程论文偏微分方程数值解课程论文专业：14级信息与计算科学指导教师：陈红斌2017年6月25日一维热传导方程的差分格式论文一、微分方程数值解课程论文的格式1)引言：介绍研究问题的意义和现状2)格式：给出数值格式3)截断误差：给出数值格式的截断

2、误差4)数值例子：按所给数值格式给出数值例子5)参考文献：论文所涉及的文献和教材二、微分方程数值解课程论文的评分标准1)文献综述：10分；2)课题研究方案可行性：10分；3)数值格式：20分；4)数值格式的算法、流程图：10分；5)数值格式的程序：10分；6)论文撰写的条理性和完整性：10分；7)论文工作量的大小及课题的难度：10分；8)课程设计态度：10分；9)独立性和创新性：10分。评阅人:一维热传导方程的差分格式一维热传导方程的差分格式 -一维热传导方程的差分格式1引言考虑如下一维非齐次热传导方程Dirichlet初边值问题=+/U,0，cxdydtof0tT,(1.1)(圮0)=於)，

3、cxd.(1.2)u(c,t)=a(t)yu(d,t)=p(t0tT(1.3)的有限差分方法，其中a为正常数久x),&,0(/)为已知常数，久c)=tz(O),久d)=0(O).称(1.2)为初值条件，(1.3)为边值条件.本文将给出(1.1)(1.3)的向前民0格式，向后Euler格式和Crank-Nicolson格式，并给出其截断误差和数值例子.经对比发现，Crank-Nicolson格式误差最小，向前民0格式次之，向后民0格式误差最人.2差分格式的建立2.1向前EW刃格式将区间c,d作M等分，将0,T作N等分，并记h=(d-c)/M,t=T/N,x=c+jhf0jM,tk=kr,0N.分

4、别称力和为空间步长和时间步长.用两组平行直线x=xj,0jM,t=tk,0kN将。分割成矩形网格.记Q/z=xy|OjM,r=tk0kN,Q,ir=Q.,xQr.称(，4)为结点.定义耳“上的网格函数Q=)|0jM,OkN,其中在结点(x.,rj处考虑方程(I.I),有尬心川）_dt_dx2将山+J以结点（山）为中心关于t运用泰勒级数展开，有)=“(,/)+/(Xj山)整理有）-du（xjytk）ro2w（xpr,）=4;+O（T）.rdt2dt2(2.1)(2.2)w(Vi，-)=u(山)+0(，4)3+”(，)()r(，L)(-府2!3!+“()(抑+耐),4!U(X田,L)=W(Xj,/

5、J+u(x.(,JFI)/?3uH2!3!/4)H4!皿)由上述两式可得（T，h）-2（山）+（饬人）_刊（山）.Ir巩山）+o（F）h2将(2.2),(2.3)两式代入(2.1)中,得(f+J-U(fJMU,.=Cl-Hdx212dx4(2.3)(S山)-2(wJ+M%,+/(&)+疋h2(2.4)再将（T，/J，（%山）分别以结点（卩山）为中心关于X运用泰勒级数展开，有一维热传导方程的差分格式一维热传导方程的差分格式 #-一维热传导方程的差分格式一维热传导方程的差分格式 -其中R：=_箋。竹严）*r_c（丁）*冋为方程.I）的截断误差.(2.5)3112dx42dr舍去截断误差，用0代替“

6、（,）,得到如下差分方程匸过=/爲_2与+；yI1,0kN-l.rlr结合初边值条件，可得如卜差分格式严1+小1JM-ilkN1,(2.6)W-=0(巧)，0JA/,(2.7)o=MJ,UM=0(4)，lkN.(2.8)2.2向后EWw差分格式仝叫丄心如），口g”(2.9)在结点（厂山+J处考虑方程（1.1）,有dtdx2将u（xjytk）以（，4+J为中心关于f运用泰勒级数展开，有2心沁）=2心川+J+讥巧山+J（r）+将上式整理得（如）-2心川）色心厂心）了刊匕）,、:=;+O（T）.rdt2or-再将山+J，（X田如）分别以（有,心）为中心关于兀运用泰勒级数展开，有匕山+J（-疔+（,）

7、（-府2!+3!w(V1、心)=u(Xj山+J+/(Xf丿,心）（-力）+-(2.10)一维热传导方程的差分格式一维热传导方程的差分格式 #-一维热传导方程的差分格式一维热传导方程的差分格式 #-d4!“（，心：H心如）+心，“+叫弓严2+叭；加由上述两式可得W（V1，hJ-2匕,L+J+（%=刊（f+J+夕（兀，）+2（2口）lr齐+12茨+久八将（2.10）,（2.11）两式代入（29）中，得一维热传导方程的差分格式一维热传导方程的差分格式 #一维热传导方程的差分格式一维热传导方程的差分格式 # r=a疋+心,心）+用其中R：=W）_宾。（号+J+0匸+冋为方程（2.9）的截断误差.2dt

8、212dx4舍去截断误差，用it代替,得到如下差分方程U广亿丄疋mt,N.T/广结合初边值条件，町得如卜差分格式i学-尤_屹;-2m*+1+ur=a示/t+i=MXik由上述两式得）-21心）+心+/）=员心山）+64（.%）+h2dx2+12dx4+同理，将Mt，“+J，u（+m+J分别以（山+J为中心关于*泰勒级数展开，(2.19)整理得“（S，心）-2（,心）+（6,心）_刊（，心）h2夕（卩儿丄）+0（肝）lr此时，分别将（f+J,UXj,tk）以U（X尸t刊餉心八gj.九（以+山丿Tdx21Tdx12dx4奸力）为中心关于f泰勒级数展开，有）r,1肌心2）fI2丿+dx2dx2dt2

9、2!dx2dt2空(Xj山)_空(Xj,S+HJ|空(Xj,如dx1dx2dx2dt,1夕山+1/2+2!dx2dt2(2.20)（山加（f|尸曲（f|3）（218）rdt245/3再将”（t，4）,（兮）分别以（山）为中心关于x运用泰勒级数展开，有2!+3!（八（八“亠八叫S）（如心）（-府一维热传导方程的差分格式一维热传导方程的差分格式- #- #一维热传导方程的差分格式一维热传导方程的差分格式- #- 矶心gj_1矶心M+J刊匕山）、1心山+”）dx22+0Fk丿2dx2dt2JJ利用上述两式得+。（尸）(2.21)利用(2.19),(2.20)两式，整理有dx1*dx202叭Xj山+J

10、。妆(Xj,tk)u(%山)-2w(Xj山)+“(Xy+I山)h212dx4(2.22)I“(T,心)一2(山+J+U(X和,心)h2夕(厂山)+乔12dx4结合(2.21),(2.22)两式，整理得一维热传导方程的差分格式一维热传导方程的差分格式- #- #2/?2(+l4+l)1心，心22dx2dt2(2丿20%(形山+”2)M(尸+J-2u(,g)+”tk)-2uxj,tk)+uxj+tk)一维热传导方程的差分格式一维热传导方程的差分格式- #- #一维热传导方程的差分格式一维热传导方程的差分格式- #- #1(h2巩，4)h2夕(212dx4V将(2.18),(2.23)式代入(2.1

11、7)得J如)/宀1+3+O/厂)12去4)(2.23)一维热传导方程的差分格式一维热传导方程的差分格式- #- #一维热传导方程的差分格式一维热传导方程的差分格式- #- #(山)2山)+(Vi，)2h2(2.24)Ms”心)-2M,“+j+M+/+j2/r+d+jXrtk)+/厂其中/r山)十/r夕(山+】)+4(卩山+”12dx2dt2212a?)fr24a?一维热传导方程的差分格式一维热传导方程的差分格式- #- #一维热传导方程的差分格式一维热传导方程的差分格式- #- #+0(尸+IF)为方程(2.17)的截断误差.舍去截断误差，用比代替(山)，可得如下差分方程I讥一idi占一2i(

12、k+U=cir2/?J川；-2喟+唸；，2Jj1jM-i0kN-l.(2.24)一维热传导方程的差分格式一维热传导方程的差分格式- - #一维热传导方程的差分格式一维热传导方程的差分格式- - 结合初边值条件，可得如卜差分格式喟屹厂2心+唏|屹;-2喟+唸;r2h22h21jM-i0kN-L,(2.25)(2.26)(2.27)町=倾吗)，o;A/,0=ag，屹=0(4)，.kN.3差分格式的求解3.1向前EWw厂格式记r=d%2,称厂为步长比.差分格式(2.6)(2.8)中(2.6)可改写为喟=r+(1-2r)w*+ru*+1+,(3.1)1j;W-1,0kN-l.将(3.1)写成如下形式U

13、k+l=A-Uk+r-Fk+r-Fl.其中叫彳SJ,，心J,0=彳,心心J,戸=/，/f，允JF严凶，0,0,，0,啲:计，TOC o 1-5 h z4-2rr、rl-2rrA=r.rr12rX/(Af1)3.2向后格式记r=%2,称厂为步长比.差分格式(2.14)(2.16)中(2.14)可改写为=r(l+2r)wj+1+r：；：+r-f：*,(3.2)1jM-l,UkVN.将(3.2)写成如下形式A-Uk+l=-Uk-T-Fk+l-r-F,.其中一维热传导方程的差分格式一维热传导方程的差分格式一维热传导方程的差分格式一维热传导方程的差分格式严+1广1t./+l.i+l.i+11M1心、上M

14、T1jM-l,OKN(3.3)=严,才，朋，坊=喟，0,0,0,隙打计，一(1+2广)r、r-(l+2r)rrr_(1+2厂力3.3Crank-Nicolson格式记/=%,称r为步长比.差分格式(2.25)(2.27)中(2.25)可改写为-r/2喑-(1+r)+1+厂/2唸;=2“二+(1-诃+厂/2临+Gf一维热传导方程的差分格式一维热传导方程的差分格式一维热传导方程的差分格式一维热传导方程的差分格式将(33)写成如下形式-A.-Uk1=AcUk+r-FkW2+r(Ft+FJ其中=咁1，W2+1，,喘J，U*=4-17，,0,0,o,町;_叩,代=仁0,0,0,啜爲川一维热传导方程的差分

15、格式一维热传导方程的差分格式一维热传导方程的差分格式一维热传导方程的差分格式严+1/2=厂+1/2yA+1/2厂+1/2丁l-rr/2r/2l-rr/2r/2r/2r/21广丿(M-1户(MT)nr/2r/2-(1+厂)r/2r/24数值例子io-在本文中，将应用向前Euler格式(2.6)-(2.8),向后Euler格式(2.14)-(2.16)和Crank-Nicolson格式(2.25)(2.27)分别来求解如下算例.算例现有如下定解问题红斗巴0*1,dtdx20r1,u(x,0)=e0 xl,(0,Z)=e2t,z/(l,r)=e1+2/,0f1.上述定解问题的精确解为u(xyt)=e

16、x+21.4.1向前Euler格式在表4.1.1中给出了取步长/?=1/10和=1/200(步长比7=1/2)时计算得到的部分数值结果.表4.1.2给出了取步长/?=1/10和:tl/lOO(步长比r=l)时计算得到的部分数值结果，随着计算层数的增加，误差越来越人，数值结果是发散的.经步长比，的多次取值发现，当厂1/2时，其数值结果是发散的.表4.1.3给出了厂=1/2时，取不同步长，数值解的最大误差E*(/?,T)=max|u(Xj,tk)-u-.lkN从表4.1.3可以看出当空间步长缩小到原来的1/2,时间步长缩小到原来的1/4时，最大误差约缩小到原来的1/4.图4.1.1给出了心1时的数

17、值解曲线(/?=1/10,r=1/200)与精确解曲线.图4.1.2给出了=1时不同步长数值解的误差曲线图(厂=1/2),由图4.1.2可知，当步长比厂保持不变，随着时间与空间方向上的加密，误差越来越小.图4.1.3给出了不同步长的误差曲面图.表4.1.1算例在部分结点处数值解、精确解和绝对误差(/?=l/10,r=l/200)kCM)数值解精确解丨精确解-数值解20(0.5,0.1)1.8028101.8039881.178e-3一维热传导方程的差分格式一维热传导方程的差分格式- - 40(0.5,0.2)2.2014382.2033961.959e-360(0.5,0.3)2.688651

18、2.6912342.583e-380(0.5,0.4)3.2838563.2870813.225e-3100(0.5,0.5)4.0108864.0148503.965e-3120(0.5,0.6)4.8988974.9037494.852e-3140(0.5,0.7)5.9835235.9894525.929e-3160(0.5,0.8)7.3082907.3155347.243e-3180(0.5,0.9)8.9263668.9352138.847e-3200(0.5,1.0)10.90268710.9134941.081e3表4.1.2算例在部分结点处数值解、精确解和绝对误差(/7=l/

19、10,r=l/100)kCM)数值解精确解1精确解-数值解1(0.5,0.01)1.5216741.5219622.879e-42(0.5,0.02)1.5521231.5527075.845e-43(0.5,0.03)1.5831841.5840748.901e-44(0.5,0.04)1.6148701.6160741.205e35(0.5,0.05)1.6473861.6487211.336e-36(0.5,0.06)1.6797851.6820282.242e-37(0.5,0.07)1.7160571.7160075.051e-58(0.5,0.08)1.7414461.750673

20、9.227e-39(0.5,0.09)1.8070891.7860382.105e210(0.5,0.10)1,7441431.8221197.798e-211(0.5,0.11)2.0925461.8589282.336e-112(0.5,0.12)1.1649231.8964817.316e-113(0.5,0.13)4.1483931.9347922.213e+014(0.5,0.14)-4.7094151.9738786.684e+015(0.5,0.15)21.9980932.0137531.998e+l16(0.5,0.16)-57.4223992.0544335.948e+l17

21、(0.5,0.17)178.2946232.0959361.762e+218(0.5,0.18)-518.1435032.1382765.203e+2表4.1.3算例取不同步长时数值解的最人误差(r=l/2)/?TEMC耳7,3)/E咖)1/101/2001.178e-2*1/201/8002.973e-33.9641/401/32007.431e-44.0001/801/128001.858e-44.000一维热传导方程的差分格式一维热传导方程的差分格式- #- #一维热传导方程的差分格式一维热传导方程的差分格式- #- #图4.1.1算例取/=1时的数值解曲线(/?=1/10,r=1/20

22、0)与精确解曲线0012h1/10,f1/200-.-.-h=1/2O,r=1.i800/、h=1/40.t=1J32O0/*/Z0.01一维热传导方程的差分格式一维热传导方程的差分格式- - 一维热传导方程的差分格式 图4.1.2算例取t=1时不同步长数值解的误差曲线(r=l/2):lh=1Q.T=200ih=20t=8D00.042=O01.图4.1.3算例取不同步长数值解的误差曲面(r=l/2)4.2向后EWer格式在表4.2.1和表4.2.2中分别给出a/10,r=l/200(步长比厂=1/2)和/2=1/10,r=l/100(步长比r=l)时计算得到的部分数值结果，发现两表的数值结果

23、都是收敛的，经步长比的多次取值发现，无论尸取任何值，其数值结果都是收敛的.表4.2.3给出了厂=1/2时，取不同步长，数值解的最大误差Eg(九r)=maxa|w(x.山)一町|lkN从表4.2.3可以看出当空间步长缩小到原来的1/2,时间步长缩小到原来的1/4时，最人误差缩小到原来的1/4.表4.2.4给出了厂=1时的类似结论.图4.2.1给出了心1时的数值解曲线(/?=l/10,r=l/200)与精确解曲线.图4.2.2给出了=1时，取不同步长数值解的误差曲线图(r=l/2),由图4.2.2可知，当步长比尸保持不变，随着时间与空间方向上的加密，误差越来越小图4.2.3和图4.2.4分别当给出

24、了步长比厂=1/2和厂=1下不同步长的误差曲面图.表4.2.1算例在部分结点处数值解、精确解和绝对误差(/?=1/10“=1/200)kW)数值解精确解1精确解-数值解20(0.5,0.1)1.8053431.8039881355e340(0.5,0.2)2.2056742.2033962.278e-360(0.5,0.3)2.6942582.6912343.023e-380(0.5,0.4)3.2908673.2870813.785e-3100(0.5,0.5)4.0195094.0148504.659e-3120(0.5,0.6)4.9094534.9037495.704e-3140(0.

25、5,0.7)5.9964255.9894526.973e-3160(0.5,0.8)7.3240527.3155348.518e3180(0.5,0.9)10.9262038.9352131.041e-2200(0.5,1.0)10.92620310.9134941.271e-2表4.2.2算例在部分结点处数值解、精确解和绝对误差(/?=l/10,r=l/100)kW)数值解精确解1精确解-数值解10(0.5,0.1)1.7885001.7860382.461e-320(0.5,0.2)2.1857422.1814724269e-330(0.5,0.3)2.6701712.6644565.71

26、5e-340(0.5,0.4)3.2615513.2543747.177e-350(0.5,0.5)3.9837453.9749028.843e-360(0.5,0.6)4.8657884.8549561.083e270(0.5,0.7)5.9430985.9298561.324e-280(0.5,0.8)7.2589217.2427431.618e-290(0.5,0.9)8.8660688.8463061.976e-2100(0.5,1.0)10.82904110.8049032.414e-2表4.2.3算例取不同步长时数值解的最人误差(r=l/2)hTEMC耳7,4“/$(/)1/101

27、/2001.386e2*1/201/8003.509e-33.9501/401/32008.779e-43.9971/801/128002.195e-43.999表4.2.4算例取不同步长时数值解的最大误差(r=l)hTEx(h,r)1/101/1002.659e-2*1/201/4006.744e-33.9431/401/16001.688e33.955一维热传导方程的差分格式一维热传导方程的差分格式- - #一维热传导方程的差分格式 #一维热传导方程的差分格式一维热传导方程的差分格式- #- #一维热传导方程的差分格式 #1/801/64004.222e-43.999一维热传导方程的差分格

28、式一维热传导方程的差分格式- #- #一维热传导方程的差分格式 #一维热传导方程的差分格式一维热传导方程的差分格式- #- #一维热传导方程的差分格式 #图4.2.1算例取21时的数值解曲线(力=1/10=1/200)与精确解曲线一维热传导方程的差分格式一维热传导方程的差分格式- #- #一维热传导方程的差分格式 #一维热传导方程的差分格式一维热传导方程的差分格式- #- #一维热传导方程的差分格式 #一维热传导方程的差分格式一维热传导方程的差分格式- #- #一维热传导方程的差分格式 #一维热传导方程的差分格式一维热传导方程的差分格式- #- #一维热传导方程的差分格式 #图4.2.2算例取

29、f=l时不同步长数值解的误差曲线(r=l/2)一维热传导方程的差分格式一维热传导方程的差分格式- #- #-一维热传导方程的差分格式一维热传导方程的差分格式一维热传导方程的差分格式- #- #-一维热传导方程的差分格式0.0140.012-h=1；10.r=1.200h=1J20,r=vaooh=1J40T=13200-Tx)nS001.0008.0.006r0004、0002、040.6040.6一维热传导方程的差分格式一维热传导方程的差分格式- #- #-一维热传导方程的差分格式一维热传导方程的差分格式一维热传导方程的差分格式- - -一维热传导方程的差分格式图4.2.3算例取不同步长数值

30、解的误差曲面卩=1/2)hvw.1nooh=V2Ot=1MOOh=”8.T=meoo0.0310025.0.02.图4.2.4算例取不同步长数值解的误差曲面(r=1)4.3Crank-Nicolson格式在表4.3.1中给出了取步长/7=l/10,r=l/10时计算得到的部分数值结果，表4.3.2给出了取步长/?=l/100,r=l/100时计算得到的部分数值结果，发现两表的数值结果都是收敛的，经步长比的多次取值发现，无论尸取任何值，其数值结果都是收敛的.表4.3.3给出了取不同步长时，所得数值解的最大误差Eg(九r)=maxa|w(x.|正左N从表可以看出当空间步长缩小到原来的1/2,时间步

31、长缩小到原来的1/2时，最人误差缩小到原来的1/4.图4.3.1给出了f=1时的数值解曲线(/?=l/10,r=l/10)与精确解曲线图4.3.2给出了=1时，取不同步长数值解的误差曲线图，由图4.3.2可知，当步长比广保持不变,随着时间与空间方向上的加密，误差越来越小.图4.3.3给出了取不同步长时所得数值解的误差曲面图.图4.3.4给出了向前Euler,向后和Crank-Nicolson三种方法数值解的误差对比图(/?=l/10,r=l/200),由图4.3.4分析可得，当时间方向f与空间方向x都取相同步长时，Crank-Nicolson方法误差最小，精度最高，向前EWw方法次之，向后EW

32、w方法所得误差最犬，精度最低.表4.3.1算例在部分结点处数值解、精确解和绝对误差(/?=l/10,r=l/10)kW)数值解精确解1精确解-数值解1(0.5,0.1)1.8231141.8221199.949e-42(0.5,0.2)2.2271962.2255411.656e33(0.5,0.3)2.7204142.7182822.132e-34(0.5,0.4)3.3227763.3201172.659e-35(0.5,0.5)4.0584604.0552003.260e-36(0.5,0.6)4.9570214.9530323.988e-37(0.5,0.7)6.0545206.049

33、6474.873e-38(0.5,0.8)7.3950087.3890565.952e-39(0.5,0.9)9.0322849.0250137.271e-310(0.5,1.0)11.03205611.0231768.880e-3表4.3.2算例在部分结点处数值解、精确解和绝对误差(/2=l/100,r=l/100)k(M)数值解精确解1精确解-数值解110(0.5,0.1)1.9937261.9937151.081e520(0.5,0.2)2.4351472.4351301.749e-530(0.5,0.3)2.9742972.9742972.296e-540(0.5,0.4)3.6328

34、153.6327862.864e-550(0.5,0.5)4.4371314.4370963.521e560(0.5,0.6)5.9195245.4194814.308e-570(0.5,0.7)6.6194216.6193695.266e-580(0.5,0.8)8.0849798.0849156.433e-590(0.5,0.9)9.8750169.8749387.857e-5100(0.5,1.0)12.06137212.0612769.597e-5一维热传导方程的差分格式一维热传导方程的差分格式- - #一维热传导方程的差分格式表4.3.3算例取不同步长时数值解得最人误差hTEx(2/

35、?,2r)/E,x(/2,r)1/101/109.588e-3*1/201/202.429e-33.94571/401/406.078e-43.99611/801/801.520e-43.99901/1601/1603.799e-53.99981/3201/3209.499e-63.99991/6401/6402.375e-640000图4.3.1算例取f=1时数值解曲线(/?=1/10=1/10)与精确解曲线一维热传导方程的差分格式一维热传导方程的差分格式2O2O一维热传导方程的差分格式一维热传导方程的差分格式一维热传导方程的差分格式2O2O一维热传导方程的差分格式0.010.0090.00

36、80.0070.0060.0050.0040.00300020.0010IVVVB1V、T=h=1/1O/、T=h=1/20/T=h=1/40/T=h=1/6O././一一一一一I/I、*、/00.10.20.30.40.50.60.70.80.91图4.3.2算例取f=l时不同步长数值解的误差曲线T=h=V10T=h=1i*20T=h=1M00.010.008、1|图4.3.3算例取不同步长数值解的误差曲面一维热传导方程的筮分格式一维热传导方程的筮分格式- - -一维热传导方程的差分格式图4.3.4算例取t=1时三种方法数值解的误差对比曲线(/?=l/10,r=l/200)5总结在本文中，给

37、出了向前Euler,向后Euler和Crank-Nicolson三种差分格式，分别对抛物型方程进行求解，其中向前格式是条件稳定的，后两者是无条件稳定的.经对比发现，Cw伙-TWco/son格式的精度最高，误差最小，向前总“少格式次之，向后总少格式精度最低，误差最大.参考文献孙志忠.偏微分方程数值解法M.第二版.北京：科学出版社，2005.李荣华.偏微分方程数值解法M.第二版.北京：高等教育出版社，2005.刘卫国.MATLAB程序设计教程M.第二版.北京：中国水利水电出版社，2010.程序流程图1向前氏少法一维热传导方程的差分格式一维热传导方程的差分格式 # 一维热传导方程的差分格式一维热传导

38、方程的差分格式一维热传导方程的差分格式 # #一维热传导方程的差分格式2向后Euler法一维热传导方程的差分格式一维热传导方程的差分格式 #一维热传导方程的差分格式一维热传导方程的差分格式一维热传导方程的差分格式 # #一维热传导方程的差分格式3Crank一Nicolson法一维热传导方程的差分格式一维热传导方程的差分格式 一维热传导方程的差分格式程序代码附录A(向前乩/少法)建立函数文件ForwardEuler.m,如下：%求解一维非齐次热传导方程ut-uxx=e(x+2t),0 xl,0t=l,%u(x,0)=e*x,0=x=l,%u(0,t)=e*(2t),u(l,t)=e(l+2t),

39、0t0)%其中M,N分別为x与t方向上的网格数。%uxy为梢确解，uij为数值解，eps为数值解与精确解的绝对误差.u=zeros(N+l,M+l);%初始化矩阵u,存储数值解的值.x=0:h2:l;%h2=(l-0)/Mt=0:hl:T;%hl=(T-0)/N%给边值条件赋值u(0,t)=e(2t),u(l,t)=e(l+2t).fork=l:N+lu(k,l)=exp(2.*t(k);u(k,M+l)=exp(1+2.*t(k);end%给初值条件赋值u(x,0)=exforj=l:M+lu(l,j)二exp(x(j);endf=zeros(N+l,M+l);%初始化源函数f(x,t).%

40、给源函数(SourceFunction)赋值，f(x,t)=e(x+2t)fork=l:N+lforj=l:M+lf(k,j)=exp(x(j)+2.*t(k);endend%r=(a*hl)/(h22)为步长比，这里取a=l,即r=hl/(h22);%当rWl/2时，数值解是越来越接近梢确解的，%当r1/2时，两者的误差会越来越大.r二hl/(h22);fork=l:Nforj=2：Mu(k+l,j)=r.*(u(k,j-l)+u(k,j+1)+(l-2*r).*u(k,j)+hl.*f(k,j);endenduxy二zeros(N+l,M+l);冬给梢确解赋初值.eps=zeros(N+l

41、,M+l);$给数值解与梢确解的绝对误差赋初值.fork=l:N+lforj=l:M+luxy(k,j)二exp(x(j)+2.*t(k);$给梢确解賦值.eps(k,j)=abs(u(k,j)-uxy(k,j);瓯数值解与精确解的绝对误差.endend一维热传导方程的差分格式一维热传导方程的差分格式- - #-一维热传导方程的差分格式枫B(向后法)建立函数文件BackwordEu1er.如下：%求解一维非齐次热传导方程ut-uxx=e(x+2t),0 xl,0t=l,%u(x,0)=e*x,0=x=l,%u(0,t)=e(2t),u(l,t)=e(l+2t),0t0)%其中M,N分别为x与t

42、方向上的网格数。%UX为梢确解，U为数值解，eps为数值解与精确解的绝对误差.一维热传导方程的差分格式一维热传导方程的差分格式- #- #-一维热传导方程的差分格式一维热传导方程的差分格式一维热传导方程的差分格式- #- -一维热传导方程的差分格式U=zeros(M+D*(N+l),1);r=c./h2;aal=zeros(1,M-2);fori=l:M-2aal(i)=r;endaa2=zeros(1,M-l);fori=l:M-laa2(i)=-l-2*r;endaa3=zeros(1,M-2);fori=l:M-2aa3(i)=r;enda=zeros(M-l,M-l);A=diag(a

43、al,-1)+diag(aa2)+diag(aa3,1);k=l;%记录行数count二1;%记录总行数fori=l:M+l%给口(x,0)赋值U(count,l)=exp(x(i);count=count+l;endU0=zeros(M-l,1);fori=l:M-l%U向虽:初始值U0(i,l)=exp(x(i+l);end玄给矩阵A对角线序号为-1的位置赋值玄给矩阵A主对角线上元素赋值.%给矩阵A对角线序号为1的位置赋值%构建矩阵A.一维热传导方程的差分格式一维热传导方程的差分格式- #- #-一维热传导方程的差分格式一维热传导方程的差分格式一维热传导方程的差分格式- #- #-一维热传

44、导方程的差分格式fl二zeros(M-l,1);一维热传导方程的差分格式一维热传导方程的差分格式 一维热传导方程的差分格式- -f二zeros(M-l,1);whilek=N+lfl(l,l)=exp(2*t(k+l);fl(M-l,l)=exp(l+2*t(k+l);Fl=r*fl;forj=l:M-1f(j,l)=exp(x(j+l)+2*t(k+1);endF=c*f;Ul=inv(A)*(-UO-F-Fl);U(count,1)二exp(2*t(k+1);%给最左边界赋值count=count+l;fori=l:M-lU(count,l)=Ul(i,1);count=count+l;endU(count,l)=exp(l+2*t(k+l);%给最右边界赋值count=count+l;k二k+1;%进去下一行UO二U1;endU=reshape(r,M+1,N+1);p二1;%求梢确解fori二1:N+1forj=l:M+lux(p,l)=exp(x(j)+2*t(i);p二p+1;endendux=reshape(ux,M+l,N+l);eps=abs(ux-U);%精确解与数值解的绝对误差附录C(Crank一Nicolso