二阶常微分方程的差分格式

中第并哮书去七方《微分方程数值解》

课程论文学生姓名1:许慧卿学生姓名2：向裕学生姓名3：邱文林学生姓名4：高俊学生姓名5：赵禹恒学生姓名6：刘志刚学院：学号：20144329学号：20144327学号：20144349学号：20144305学号：20144359学号：20144346理学院专业： 14级信息与计算科学 指导教师： 陈红斌 2017年6月25日《偏微分方程数值解》课程论文《二阶常微分方程的差分格式》论文一、《微分方程数值解》课程论文的格式引言：介绍研究问题的意义和现状格式：给出数值格式截断误差：给出数值格式的截断误差数值例子：按所给数值格式给出数值例子参考文献：论文所涉及的文献和教材二、《微分方程数值解》课程论文的评分标准文献综述：10分；课题研究方案可行性：10分；数值格式：20分；数值格式的算法、流程图：10分；数值格式的程序：10分；论文撰写的条理性和完整性：10分；论文工作量的大小及课题的难度：10分;课程设计态度：10分；独立性和创新性：10分。评阅人:-2-

二阶常微分方程的差分格式二阶常微分方程的差分格式1引言考虑如下二阶常微分方程两点边值问题d2u(x)d2u(x)+u(x)q(x)=f(x)9a<x<b,=a,u(b)=p,的有限差分方法，其中q(x),f(x)为［。，1上的连续函数，(/(X)>0;a,p为给定常数.本文将对(1),(2)差分格式及紧差分格式进行求解，并给出截断误差与数值例子.2差分格式及紧差分格式的建立2.1差分格式将区间［a,。］作M等分，记力=@—a)/M,x^a+ih,0<i<M,/?为网格步长，为为网格结点，％为网格，称定义在鼻上的函数为网格函数.在网格结点占处考虑方程(1),有-u\xi)+qQJ〃(％)=/(%), l<z<M-1, (3)〃(大o)=a, 〃(・%)=△任取一个结点苍，(i=l,2,3,…，A7-1),将〃(/J,〃(％_］)在占处运用泰勒级数展开/ ，/>“"GM?ir(x.)h5w(4)(x.)/?4 .”(Z+1)=〃(苍)+ugh+---+—--+—--+o(h),"%)="用+/(0(叫+ ；)+ 3： +I+。()利用上面两式得-3-

二阶常微分方程的差分格式〃(％])+ =2〃区)+〃”(为)始+--T—+。(/?4).JL乙整理得〃―丫…-J=/⑺+飞％。々). ⑷由⑷式得〃«)=出上竿2士四墟_[*£+。断)]. (5)将式⑸代入式⑶得」品)-2%)+，，院)+依)(6)〃",(丫)人2其中凡二——d^+o(/)为差分方程的截断误差.12用/代替〃(%),并舍去截断误差得如下差分方程TOC\o"1-5"\h\zw..-2m.+u... /、 ~、 /、-――K——+q(%)%=/(Xf.). (7)方程(7)结合边值条件，得如下差分格式-――2"1 +q(xj〃,=/(%)， 1<Z<M-l. (8)lr〃o=a, 11m=B•2.2紧差分格式现在定义平均算子—+1。/+匕+1), 1<z<Af-1.JL乙(A"，=<, i=0,M.■•在结点工处考虑方程(1),有(10)一〃”(菁)+q(xj〃(z)=/(%),0<i<M.(10)用算子A作用于(10)式，整理可得二阶常微分方程的差分格式112dx2112dx2।小「心)।d"Q+J

dx2dx2+'夕(％_])〃(4])+与q(xDJL乙 JL乙(14)+：q(z+J〃(x,+J=][/(苍.J+i°/(z)+(14)JL乙 JL乙根据引理团,当函数g(x)eC'[c—/"+可，有1 1 /74—[g\c-h)+l0g\c)+g\c+h)]=—[g(c+h)-2g(c)+g(c-h)]+—g{6>^.),(15)结合(14),(15)两式，整理有-p-[〃(％)—2“(七)+H(X/+1)]+，■国(Xi)〃(x，T)+10qQ,)〃(x,)+式4”(£*)]TOC\o"1-5"\h\z1 /74=—[/(Vi)+10/(Xi)+/C\+1)]- "⑹(A)• (16)JL乙 乙\J其中4g(Zt，/[1<z<M-1j/(a0)=a,〃(％)=p.将(16)代入(10)中，有"+1-[q(xi_l)u(<xi_l)+10贝王)%(％)+q(xM)u(x^)]n 12=今[”“1)+10""‘)+"'+用+K'其中为=-h其中为=-h4240〃⑼(。)为差分方程的截断误差.(17)(17)(18)(19)舍去截断误差，并用/代替”(%),可得如下差分方程u—24.+u. 1一--一■十——十B国QtMt+iOg(XM+夕(九北一』=II【/(4J+1°/(w)+/(刈方程(17)结合边值条件，可得如下.紧差分格式-+也+'团(4】),*+1°q(x,M+夕(七+1)/+J

=:[/(%)+10/(x,)+/(菁+1]，1</<M-1,JL乙・5・

二阶常微分方程的差分格式3差分格式及紧差分格式的求解1差分格式的求解将差分格式(8),(9)中(8)整理为(20)可将式(20)表示为4〃二下,其中’2+/q(xj-1-12+/*(4)-1尸=(/("+a，2+/?一夕(九) -1-1 2+h2q(x2)-1(20)可将式(20)表示为4〃二下,其中’2+/q(xj-1-12+/*(4)-1尸=(/("+a，2+/?一夕(九) -1-1 2+h2q(x2)-1-12+lrq(xM_2)-1一1 2+/r(7(xv.1)>、U./(七)=/r-12+/厂式/_?)-1-12+该系数矩阵为三对角矩阵，用追赶法进行求解较好.3.2紧差分格式的求解将差分格式(18),(19)中(18)整理为11/、］ 「210/J「11/1丁五式巧G—h一夕区)"=5九+岑力+*力-1, \<j<M-1.JL乙 JL乙 JL乙将上式整理为=尸,・6・二阶常微分方程的差分格式其中210 -11后+逐式"京+豆虱&）二十工虱1,）马+山式&）人=h212 -h212 3T1, 、210,、庐+》(％7)庐+小⑸川……11、o■01012112112101211211210h，712人/".I>1、.乜）JL乙0

*07T/（%）\1Z/-1布+丘夕(％-2)210庐+正式j-J1121I21121I210fM-212/\/w-l>/（X。）12

00■

.

*0/（L）124数值算例算例1应用差分格式(8),(9)计算如下两点边值问题-7-二阶常微分方程的差分格式TOC\o"1-5"\h\z一〃"(x)+w(x)=e'(sinx-2cosx),0<x<^, (21)“(0)=0,”(乃)=0. (22)该定解问题的精确解为〃(x)=e'smx.将区间［0,%］作M等分，并记〃二//M,%=法，差分格式为+w=eXi(smx-2COSXJ,1<z<M-I,ltQ=°， =5卜面将给出取不同步长时在结点处所得的最大模误差及阶1)最大模误差Em=max|h,-w(xz)|, (23)i2)阶E

rate=log.——. (25)■E2M从下面表1可以看出，步长缩小到原来的1/2时，最大模误差缩小到原来的1/4.下面图1给出了取不同步长时所得精确解与数值解的曲线.表1算例1取不同等分数目时最大模误差及其对应的阶MEIErate81.532e-l4.00042.0023163.824e-24.00182.0006329.555e-34.00052.0017642.388e-34.00012.00001285.971e-44.00002.0000二阶常微分方程的差分格式X XX X图1算例1取不同等分数目时数值解与精确解的对比曲线由图1可知，当等分数目M逐渐增大时，数值解与精确解的曲线拟合程度愈佳，当M增大到一定程度时，两者几乎重合。算例2应用差分格式(18),(19)计算算例1所给的两点边值问题(21),(22).将区间［0,4］作M等分，并记/?二4/用，玉=，力，OKx<4.差分格式为_2%+/+1)+：(〃,T_2%+/+1)+：(〃,TJL乙+10%+%+J-2cosx/t)+1Oe'(sin毛一2cosxf.)+e%(sinxiU-2cosx/+1)〃o=°，“M=。•下面表2给出了取不同步长时在结点处所得模误差及其阶.可以看出，当步长缩小到原来的1/2时，最大模误差约缩小到原来的1/16.下图2给出了算例2取不同步长时所得精确解与数值解的曲线.-9-

二阶常微分方程的差分格式表2算例2取不同等分数目时最大模误差及其对应的阶ME/Erate81.874e-315.53123.9571161.207e-415.91093.9919327.583e-615.96093.9965644.751e-715.99783.99981282.970e-815.99473.9995图2算例2取不同等分数目时的数值解与精确解的对比曲线♦参考文献[1]孙志忠.偏微分方程数值解法[M].第二版.北京：科学出版社，2005.6-7.[2]李荣华.偏微分方程数值解法[M].第二版.北京：高等教育出版社，2005.[3]刘卫国.MATLAB程序设计教程[M].第二版.北京：中国水利水电出版社，2010.-10-二阶常微分方程的差分格式♦程序流程图总流程图及雅可比流程图如下:程序代码附录A差分格式代码建立函数文件Dirichlet,m,如下：functionU=Dirichlet(M,h,qx,fx,uO,uM)断10,uM为边值，M为区间［aO,b0_」,.的等分数目，h为步长，其中qx,fx为(a,b)上的连续函数A=zeros(M-l); 为构造矩阵A为(M-l)X(M-l)的零矩阵for 先给矩阵A主对角线上赋值二阶常微分方程的差分格式A(i,i)=2+h*h*qx(i);endforj=l:M-2 %给矩阵A对角线序号为1的位置赋值A(j,j+l)=-l；endfork=2:M-l 先给矩阵A对角线序号为-1的位置赋值A(kfk-1)=-1；endB=zeros(MT,1); 际构造向量B为(MT)X1的零矩阵fory=2:M-2B(y,1)=h*h.*fx(y); 先给向量B从2到\1-2行位置上赋值endB(lzl)=h*h.*fx(l)+uO; %给向量B第一个元素赋值l)=h*h.*fx(M-l)+UM;%给向量B最后一个元素赋值disp(,数值解ux为：')； %AU=B,即U=inv(A)*BU=(A\B)1; ，输出数值解的值命令行窗II：»disp('输入定义区间［a,b］：')输入定义区间［a,b］：>a=inputCa=，);a=0>b=input('b=')；b=pi»M二input('输入等分数目M：');输入等分数目M：8»h=(b-a)/M;x=a+h：h：b-h;qx=x./x;ux=exp(x).*sin(x);-12-二阶常微分方程的差分格式fx=exp(x),*(sin(x)-2.*cos(x));>uO=inputCx=a处的函数值u(a)：')x二a处的函数值u(a)：0uO=0>uM=input('x=b处的函数值u(b)：')x二b处的函数值u(b)：0uM=0>Dirichlet(M,h,qx,fx,uO,uM)数值解ux为：Ans=0.5303 1.4769 2.8905 4.6704 6.4287 7.3225 5.8940»disp('精确解ux为：');精确解ux为：»uxux=0.5667 1.5509 3.0009 4.8105 6.5819 7.4605 5.9796紧差分格式代码建立函数文件DirichletJ.m,如下：functionU=DirichletJ(hzMzqxzfxzu0zuM,q0zqmzf0zfm) 为边值/M为区间［a0,b0］上的等分数目，h为步长，其中qx,fx为(a,b)上的连续函数A=zeros(M-1); *构造矩阵A为(MT)X(M-1)的零矩阵fori=l:M-l 为给矩阵A主对角线上赋值A(i,i)=2./hA2+5*qx(i)./6;endforj=l:M-2 为给矩阵A对角线序号为1的位置赋值A(j,j+1)=T•/(h*h)+qx(j+1)./12;end-13-二阶常微分方程的差分格式fork=2:M-1fork=2:M-1为给矩阵A对角线序号为-1的位置赋值A(k,k-1)=-1./(h*h)+qx(k-l)./12;A(k,k-1)=-1./(h*h)+qx(k-l)./12;endB=zeros(MT,MT);%构造向量B为(M-l)X(M-l)的零矩阵endB=zeros(MT,MT);%构造向量B为(M-l)X(M-l)的零矩阵C=zeros(MT,1);学构造向量C为(M-l)XI的零矩阵forc=l:M-1电给矩阵B主对角线上赋值endC(c,1)=fx(c);为为给矩阵B主对角线上赋值endforjl=l:M-2为endforjl=l:M-2为给矩阵B对角线序号为1的位置赋值end学学给矩阵B对角线序号为T的位置赋值endendD=B*C；为生成右边矩阵DD=B*C；D(lzl)=D(lrl)+fO./12-uO./(h*h)-qO*uO./12;留给向量D第一个元素赋值D(M-l,1)=D(M-l, /12-uM./(h*h)-qm*uM./12;留给向量D最后一个元素赋值disp(,数值解disp(,数值解ux为：1;%AU=DZKPU=inv(A)*DU=(A\D)1U=(A\D)1；命令行窗II：»disp('输入定义区间')输入定义区间[a,b]:>>a=input('a=')；a=0>>b=inputCb=，);b=pi-14-

二阶常微分方程的差分格式»M二i