偏微分方程数值解第一次大作业.docx

偏微分方程数值解——第一次大作业作业一：用有限差分法计算：要求：L随机选择光滑（随机），且统计L2-误差规模尽可能大根据题意：将线段0~L均匀分成M段，则网格步长，节点分布为共M+1个，，可得原方程的差分形式：其中，下标k表示在处取值，（）表示在与（与）中点处的取值。化简上式可得：写成矩阵形式：，得，计算时，取，，，，可以解得原方程的解析解（精确解）为由于A为对称三对角阵且主对角占优，可以用直接法求解此方程组。此外，此方程组还可以用经典的Jacobi迭代，共轭梯度（CG）法以及预处理CG法分别来求解，并比较它们的统计误差以及计算时间。但是，注意到当M很大时，由于我们选取的函数形式，有，所以对于Jacobi迭代的系数矩阵有（M越大，越趋近于1），也即此时Jacobi迭代的收敛速度很慢，所以此时没有必要进行Jacobi迭代的计算了。程序中只对直接法、共轭梯度（CG）法和预处理CG法的统计误差和计算时间进行了比较。由于系数矩阵A为对称三对角阵，所以有以下算法：记直接法：，，以上过程其实为消元过程，再用回代过程求出方程组的各变量：，共轭梯度法：令对直到收敛完成，结束。特别注意：由于A为对称三对角，为了减少计算过程的乘法次数以及计算机的存储空间，采用两个数组来表示矩阵A，则有其中，，。预处理共轭梯度法：采用的是Jacobi迭代预处理，取令对直到收敛完成，结束注意：仍用CG法中提到的方法简化计算，其实就是解线性方程组，为了提高计算速度采用直接法解此方程组。至此，所有方法的算法都已经给出。在程序的调试过程中发现，当M6000时，CG法的收敛速度已经很慢了，这是由于此时系数矩阵A的非常趋近于1。因此程序中对于CG算法取的M最大值为5000次。上述，分别为A的最大，最小特征值。结果如下所示：表一、计算时间（s）M*100012345678910直接法0.0020.0010.0010.0010.0010.0010.0010.0010.0020.002CG法2.3710.6726.2037.2067.12/////预处理CG0.0040.0140.0350.0950.2060.3730.7841.1963.8834.593表二、计算误差M*100012345678910直接法5.09E-071.80E-079.80E-086.43E-084.49E-083.54E-082.41E-082.63E-082.17E-081.25E-08CG法5.09E-071.80E-079.80E-086.43E-084.49E-08/////预处理CG5.09E-071.80E-079.80E-086.43E-084.50E-083.54E-082.42E-082.64E-082.17E-081.27E-08图1.M与计算误差的双对数由表一可知，用直接法来解此问题是最快捷的，其计算时间在M6000内基本不随M变化，CG法和预处理CG法随着规模的增大逐渐增大，同时预处理CG法比CG法快了许多。当然，由于这个矩阵为对称三对角矩阵，是一个大型稀疏矩阵，CG算法以及预处理CG算法的优势不能体现出来。但是，如果系数矩阵仅仅是一个大型矩阵，CG算法尤其是预处理CG算法的优势就非常明显了。至于Jacobi迭代方法，本人也进行了尝试，发现其在M=500时计算时间就达到了15s；至于M=1000时，计算运行了5min还是不能达到收敛，而且其计算误差也比前三者大了至少一个量级，因此对于大型系数矩阵，Jacobi方法不具有可行性。对于计算精度，从表一可以看出三种方法的误差基本差不多，同时从图1可以看出当M6000时，计算误差随着M的增大稳定的减小，但是在M6000以后减小幅度不再稳定，甚至误差还会增大，这是由于计算机的固有计算精度所致。再增大M，舍入误差的影响将会变的明显。再往后继续增加计算规模，计算误差只能维持10-8左右，这与之前做的第三次作业的结果是一致的，都是与计算机步长有关。作业二：同样考虑二维问题：要求：用Jacobi迭代用CG迭代用预处理CG迭代统计误差和计算时间最大规模不小于此类问题在实际中可以当作变导热系数、有内热源的二维稳态导热问题。为解决此问题，做矩形网格的差分格式。取定x轴和y轴方向的步长和，设的四个邻点为，，，，用表示的中点，，见下图。原方程等价于在阴影区域上对原方程进行积分，利用梯形公式有类似的有此外，有利用上面的四个积分近似式并用差商代替其中的偏导数，可以得到原方程的差分表达式将网格的顶点先按j由小到大，j相同的按i由小到大的次序排列，这种排列方式叫“自然顺序排列”，上述差分形式可以得到一个线性方程组，写成矩阵形式，则有，，其中且可记其中，，系数矩阵A有以下几个特点：A是块三对角矩阵，共有5条对角线上有非零元素；A是不可约对角占优的（由于）A为稀疏矩阵由于