偏微分方程数值解法

精品文档偏微分方程数值解法课 程 设 计题 目： 六点对称差分格式解热传导方程的初边值问题姓 名： 王晓霜 学 院： 理学院 专 业： 信息与计算科学 班 级： 0911012 学 号： 091101218 指导老师：翟方曼2012年12月14日一、题目用六点对称差分格式计算如下热传导方程的初边值问题已知其精确解为二、理论1考虑的问题考虑一维模型热传导方程(1.1) ，其中为常数。是给定的连续函数。（1.1）的定解问题分两类： 第一，初值问题（Cauchy 问题）：求足够光滑的函数，满足方程（1.1）和初始条件：（1.2） , 第二，初边值问题（也称混合问题）：求足够光滑的函数，满足方程（1.1）和初始条件： , 及边值条件 ， 假定和在相应的区域光滑，并且于，两点满足相容条件,则上述问题有唯一的充分光滑的解。现在考虑边值问题（1.1），（1.3）的差分逼近取 为空间步长，为时间步长，其中，是自然数，, ； , 将矩形域分割成矩形网格。其中 表示网格节点；表示网格内点（位于开矩形中的网格节点）的集合；表示位于闭矩形中的网格节点的集合；表示-网格边界点的集合。表示定义在网点处的待求近似解，。注意到在节点处的微商和差商之间的下列关系（）：2区域网格剖分取空间步长和时间步长，其中都是正整数。用两族平行直线和将矩形域分割成矩形网格，网格节点为。以表示网格内点集合，即位于开矩形的网点集合；表示所有位于闭矩形的网点集合；是网格界点集合。其次，用表示定义在网点的函数，3建立相应差分格式数值分析中，Crank-Nicolson方法是有限差分方法中的一种，用于数值求解热方程以及形式类似的偏微分方程。它在时间方向上是隐式的二阶方法，数值稳定。该方法诞生于20世纪，由John Crank与Phyllis Nicolson发展。向前差分格式 ， =0 向后差分格式 ， =0将向前差分格式和向后差分格式做算术平均，得到的差分格式称之为六点对称格式，也称为Grank-Nicholson格式： ， =0进一步， +=+按层计算：首先，取，则利用初值和边值=0，来确定出第一层的，即求解方程组：+=+，=0。求出，在由,取，可利用，解出，。如此下去，即可逐层算出所有，。若记三、截断误差=。注意：=又两式相加而+故有。四、稳定性与收敛性 抛物方程的两层差分格式可以统一写成向量形式:(1.7) 其中，和是阶矩阵。我们假定可逆，即（1.7）是唯一可解的。对于显格式，等于单位矩阵。三层格式可以通过引入新变量化成两层格式。假设差分解的初始值（其实可以是任一层的值）有误差，以后各层计算没有误差，让我们来考察初始误差对以后各层的影响。令和分别是以和为初始值由差分格式（1.7）得到的两组差分解，则满足（1.8） 因此，按初值稳定应该意味着。这就导致如下定义： 假设，我们称差分格式（2.1）按初值稳定，如果存在正常数和，使得以下不等式成立：（1.8） ， 这里是上的某一个范数，例如 类似地，假设，我们称差分格式（2.1）按右端稳定，如果存在正常数和，使得以下不等式成立：（1.8） ， 可以证明，差分格式若按初值稳定，则一定按右端稳定。因此，这时我们简单地称差分格式稳定。前面讨论的向前差分格式当网比时稳定，当时不稳定。这就意味着给定空间步长以后，时间步长必须足够小，才能保证稳定。而向后差分格式和Grank-Nicholson格式（1.6）则对任何网比都是稳定的，时间步长可以取得大一些，从而提高运算效率。如果某个差分格式的截断误差当和趋于0时随之趋于0，则称这个差分格式是相容的。可以证明：若差分格式是相容的和稳定的，则它是收敛的，并且差分解与微分解之间误差的阶等于截断误差的阶。因此，对任何网比，向后差分格式（1.6）有收敛阶。五、结论对于扩散方程（包括许多其他方程），可以证明Crank-Nicolson方法无条件稳定。但是，如果时间步长与空间步长平方的比值过大（一般地，大于1/2），近似解中将存在虚假的振荡或衰减。基于这个原因，当要求大时间步或高空间分辨率的时候，往往会采用数值精确较差的后向欧拉方法进行计算，这样即可以保证稳定，又避免了解的伪振荡。由本题可以总结出，抛物型方程的六点对称差分法所得的数值解能够较好地逼近方程的精确 解，且区域剖分得越细，即步长越小，数值解与精确解的误差就越小，数值解越逼近精确解。六、附录取，则，满足稳定性条件取，则，满足稳定性条件另取，则，亦满足稳定性条件程序代码format longa=2;l=1;T=1;N=5;M=100;h=l/N;to=T/M;r=(a*to)/h2;for j=1:N+1 x(j)=(j-1)*h; for k=1:M+1 t(k)=(k-1)*to; u(j,k)=exp(x(j)+2*t(k); endendu %求解精确解for j=1:N+1 x(j)=(j-1)*h; us(j,1)=exp(x(j);endfor k=1:M+1 t(k)=(k-1)*to; us(1,k)=exp(2*t(k); us(N+1,k)=exp(1+2*t(k);endfor k=2:M+1 for j=2:N us(j,k)=r*us(j-1,k-1)+(1-2*r)*u(j,k-1)+r*us(j+1,k-1); endendus %求解数值解for k=1:M+1 for j=1:N+1 R(j,k)=abs(u(j,k)-us(j,k); endendR %计算误差Rmax=max(max(R) %求误差的最大值精确解与数值解的比较：x=0:0.1:1;hold onplot(x,u(:,M+1),b);plot(x,us(:,M+1),y);title(t=1，h=1/10，=1/400时精确解和数值解的比较)text(0.05,21,蓝：精确解);text(0.05,20,黄：数值解);hold off取不同步长时的误差比较：x=0:1/5:1;y=0:1/10:1;z=0:1/20:1;hold onplot(x,R(:,M+1),b);hold offN分别取5,