微分方程数值解法课程设计---抛物型方程问题的差分格式.doc

目 录 一、问题的描述5二、算法设计及流程图52.1 算法设计52.2 流程图6三、算法的理论依据及其推导63.1 截断误差分析63.2 稳定性分析7四、数值结果及分析7 五、总结9六、附件（源代码）10 抛物型方程问题的差分格式一、问题的描述有限差分方法就是一种数值解法，它的基本思想是先把问题的定义域进行网格剖分，然后在网格点上，按适当的数值微分公式把定解问题中的微商换成差商，从而把原问题离散化为差分格式，进而求出数值解。此外，还要研究差分格式的解的存在性和唯一性、解的求法、解法的数值稳定性、差分格式的解与原定解问题的真解的误差估计、差分格式的解当网格大小趋于零时是否趋于真解(即收敛性)，等等。偏微分方程边值问题的差分法是物理上的定常问题，其定解问题为各种边值问题， 即要求解在某个区域内满足微分方程，在边界上满足给定的边界条件。常系数扩散方程的差分解法可归结为选取合理的差分网格，建立差分格式求解。常系数扩散问题的有限差分格式求常系数扩散问题 （1.1）的近似解，其初始条件为 二、算法设计及流程图2.1 算法设计运用加权隐式格式求解常系数扩散问题（1.1），（1.6）步骤1 输入初始值，确定加权隐式格式的参数；步骤2 定义向量A，把初边值条件离散，得到，j=0,1,J的值存入向量A步骤3 利用加权隐式差分格式由第n层计算第n+1层，建立相应线性方程组，求解并且存入向量A;步骤4 计算到t=1，输出u2.2 流程图开始输入初始值求h=(maxx-minx)/(n-1)求u0(j) = PrIniU(minx+(j-1)*h)利用加权隐式格式求u输出u;结束三、算法的理论依据及其推导3.1 截断误差分析常系数扩散问题（1.1）的加权隐式格式如下： ，（1.6）其中，把（1.6）改写为便于计算的形式= （1.7）其中。下面来求出差分格式（1.6）的截断误差，设u(x,t)是方程（1.1）的充分光滑的解，在处进行Taylor级数展开并经化简有.由此可以看出，当时，截断误差为。特别引起注意的是的情况，此时差分格式的截断误差是，即差分格式是二阶精度的。3.2稳定性分析我们用Fourier方法分析差分格式（1.6）式的稳定性。容易求出（1.6）式的增长因子为根据von Neumann 条件是稳定性的充要条件，因此只验证，右端一定成立，只需要考虑左端的不等式，即因，因此只需要，这就是差分格式（1.6）的稳定性限制。此条件也可以写得更明确些，即加权隐式格式稳定的条件是，当；无限制， 当 四、数值结果及分析采用加权隐式差分格式求解其精确解为。取轴方向的网格步长分别为轴方向的网格步长为。计算在的近似解。（1）取0.5在 Matlab 命令窗口键入：u = peParabWegImp(1,0.5,0.01,11,0,1,0,0,10) 输出结果为：u = Columns 1 through 9 0 0.1203 0.2248 0.3076 0.3607 0.3790 0.3607 0.3076 0.2248 Columns 10 through 11 0.1203 0(2) 取0在 Matlab 命令窗口键入：u = peParabWegImp(1,0,0.01,11,0,1,0,0,10)输出结果为：u = Columns 1 through 10 0 -0.8819 5.5939 -7.9892 9.9736 -9.5804 9.9736 -7.9892 5.5939 -0.8819 Column 11 0(3) 取1在 Matlab 命令窗口键入：u = peParabWegImp(1,1,0.01,11,0,1,0,0,10)输出结果为：u = Columns 1 through 10 0 0.1215 0.2310 0.3180 0.3738 0.3930 0.3738 0.3180 0.2310 0.1215 Column 11 0（4）取0.2在 Matlab 命令窗口键入：u = peParabWegImp(1,0.2,0.01,11,0,1,0,0,10)输出结果为：u = Columns 1 through 9 0 0.1208 0.2228 0.3017 0.3549 0.3708 0.3549 0.3017 0.2228 Columns 10 through 11 0.1208 0实验结果分析：所得答案u为u在的近似解，第一到第十一列分别对应在第一到第十一个网格点定解问题的近似解。取0,0.2时，所得结果不稳定，取0.5,1时，所得结果稳定。五、总结通过两周的课程设计，收获颇丰。在设计和修改的过程中，我不仅巩固了以前所学过的知识，而且知道了很多在书本上所没有学到过的知识。通过这次课程设计使我懂得了理论与实际相结合是很重要的，只有理论知识是远远不够的，只有把所学到的理论知识与实践相结合起来，从理论中得出结论，将结论辅助与理论，才能真正学到知识并应用。 这次课程设计是对matlab功能的进一步熟悉和了解，尤其在函数的运用方面提高了很多。在这期间遇到了很多问题，譬如定义的函数名和文件保存的名称不同，及子文件保存的目录与主函数不一致或者在主窗口运行时，路径不相同，导致程序不可运行等等问题，但是在老师的指导和同学的帮助下，这些问题一一得到了解决，在此要感谢老师和同学的帮助。六、附件（源代码）function u = peParabWegImp(c,sita,dt,n,minx,maxx,lbu,rbu,M)format long;h = (maxx-minx)/(n-1);u0(1) = lbu;u0(n) = rbu;for j=2:n-1 u0(j) = PrIniU(minx+(j-1)*h);end u1 = u0;for k=1:M A = zeros(n-2,n-2); cb = zeros(n-2,1); cb(1) = -u0(2) -(1 - sita)*(u0(3)-2*u0(2)+u0(1)*dt*c/h/h/2 . - sita*dt*c*lbu/h/h; cb(n-2)=-u0(n-1)-(1-sita)*(u0(n)-2*u0(n-1)+u0(n-2)*dt*c/h/h/2 . - sita*dt*c*rbu/h/h; for i=2:n-3 cb(i) = -u0(i+1) -(1 - sita)*(u0(i+2)-2*u0(i+1)+u0(i)*dt*c/h/h; end A(1,1) = -2*sita*dt*c/h/h -1; A(1,2) = sita*dt*c/h/h ; for i=2:n-3 A(i,i-1) = sita*dt*c/h/h ; A(i,i) = - 2*sita*dt*c/h/h -1 ; A(i,i+1) = sita*dt*c/h/h ; end A(n-2,n-2) = - 2*sita*dt*c/h/h -1; A(n-2,n-3) = sita*dt*c/h/h; u1(2:(n-1) = Acb; u0 = u1