向前差分格式

偏微分方程数值解法作业 题        目：  抛物型方程最简差分格式        学　  院：      理  学  院          专业：  2008级信息与计算科学二班  姓        名：        姜智                学  号：        20083718              指导老师：        陈宏斌                2011年6月 一 考虑问题及条件 考虑一维热传导方程： …………（1.1） 其中 是正常数， 是给定的连续函数。 现在考虑第二类初边值问题的差分逼近： 初始条件： …………（1.2） 边值条件： ， ， ………（1.3） 假设 和 在相应区域光滑，并且在 满足相容条件，使上述问题有惟一充分光滑的解。 用向前差分格式计算如下热传导方程的初边值问题 已知其精确解为 u(x,t)=1-x. 二 网格剖分与差分格式 （1）.区格剖分 取空间步长 和时间步长 ，其中 都是正整数。用两族平行直线 和 将矩形域 分割成矩形网格，网格节点为 。以 表示网格内点集合，即位于开矩形的网点集合； 表示所有位于闭矩形 的网点集合； 是网格界点集合。 其次，用 表示定义在网点 的函数， （2）.微分方程的离散，建立相应差分格式 将方程在节点 离散化， ，   …………（1.4） 对充分光滑的解 ，由Taylor展式： …………（1.5） …………（1.6） …………（1.7） （1.5）移项得： …………（1.8） （1.6）（1.7）相加得： …………（1.9） 将（1.8）（1.9）代入（1.4）得： …………（1.10） 其中， 舍去 ，得到逼近（1.1）的向前格式差分方程： ，   ……（1.11） 其中， ， 其中j=1,2,…,N-1, k=0,1,2,…,M-1.以 表示网比,则可将（1.11）式化为： 矩阵形式： . 3、截断误差 将（1.8）（1.9）代入（1.4）得：    …………                                                                            （1.10） 其中， 而得到的向前逼近差分方程是舍去 的结果，固 即为向前差分方程的局部截断误差。 4、稳定性 如果在应用差分格式（1.1）时，计算右端函数 有误差 ，计算初值 时有误差 ，令 ，则计算可得 ， ，由定理可知，当网比 时，有 . 上式说明当 和 很小时，误差 也很小，即当网比 时，差分格式是稳定的。 5、数值例子 实现结果： 6、参考文献 [1] 李荣华 偏微分方程数值解法  高等教育出版社 2005。 [2] 孙志忠 《偏微分方程数值解法》科学出版社，  2005。 [3] 汤怀民、胡建伟，《微分方程数值解法》，南开大学出版社，1990。 [4] 李庆扬 《数值分析》 清华大学出版社 ，2008 。 [5] 余德浩，汤华中编著：《微分方程数值解法》（第三版）科学出版社2004.6。 7、附件 1.流程图 程序实现流程图 2、详细 1.建立函数文件及执行文件 文件1：fun_ux0.m      用于处理初值条件； 文件2：fun_u0t.m        用于处理边值条件； 文件3：fun_u1t.m        用于处理边值条件； 文件4：fun_u.m          计算真值； 文件5：fun_f.m          计算f(x) 文件 6：fun.m 2.程序实现代码 2.1 输入 输入时间上限T，时间方向的节点数M；输入空间上限l，空间方向的节点数N；输入大于0的常数a.T,l大于0；M，N为正整数. 2.2 设置初值 h=l/N;v=T/M; 节点： x=0+h:h:1;t=v:v:T; 网比：r=a*v/(h*h); u0=ones(1,N-1);% u0为初值,f是一个1*（N-1）的矩阵 for k=1:N-1 u0(k)=fun_ux0(x(k))*u0(k); f(k)=v*fun_f(x(k)); end f(1)=f(1)+r*fun_u0t(0);%边值条件 f(N-1)=f(N-1)+r*fun_u1t(0); A=diag(r*ones(1,N-2),-1)+diag((1-2*r)*ones(1,N-1))+diag(r*ones(1,N-2),1); 2.3 计算节点的值 u1=A*u0'+f'; u0=u1 u(k,:)=u1; f(1)=f(1)-r*fun_u0t(0); f(N-1)=f(N-1)-r*fun_u1t(0); for k=2:M f(1)=f(1)+r*fun_u0t(t(k-1)); f(N-1)=f(N-1)+r*fun_u1t(k-1); u1=A*u0+f';              u0=u1; u(k,:)=u1; f(1)=f(1)-r*fun_u0t(t(k-1)); f(N-1)=f(N-1)-r*fun_u1t(k-1); end