二阶偏微分方程matlab解答,二阶椭圆偏微分方程实例求解(附matlab代码).docx

二阶椭圆偏微分方程实例求解(附matlab代码).docx

《微分方程数值解法》期中作业实验报告二阶椭圆偏微分方程第一边值问题姓名：学号：班级：2013年11月19日1二阶椭圆偏微分方程第一边值问题摘要对于解二阶椭圆偏微分方程第一边值问题，课本上已经给出了相应的差分方程。而留给我的难题就是把差分方程组表示成系数矩阵的形式，以及对系数进行赋值。解决完这个问题之后，我在利用matlab解线性方程组时，又出现“outofmemory”的问题。因为99*99阶的矩阵太大，超出了分配给matlab的使用内存。退而求其次，当n=10，h=1/10或n=70，h=1/70时，我都得出了很好的计算结果。然而在解线性方程组时，无论是LU分解法或高斯消去法，还是gauseidel迭代法，都能达到很高的精度。关键字：二阶椭圆偏微分方程差分方程outofmemoryLU分解高斯消去法gauseidel迭代法一、题目重述解微分方程：22(,)(,)()(,)(,(,)1yxxyyxyyeueuuuxye已知边界:(0,),(),(0),(1)xe==求数值解,把区域分成，n=100[]G´2/,/0h注：老师你给的题F好像写错了，应该把改成。xye2yxe二、问题分析与模型建立2.1微分方程上的符号说明𝐴(𝑥,𝑦)=𝑒𝑦𝐵(𝑥,𝑦)=𝑒𝑥𝐶(𝑥,𝑦)=𝑥+𝑦𝐷(𝑥,𝑦)=𝑥‒𝑦𝐸(𝑥,𝑦)=1𝐹(𝑥,𝑦)=221yxyxyee2.2课本上差分方程的缺陷课本上的差分方程为：𝑎𝑖𝑗𝑢𝑖𝑗‒(𝑎𝑖‒1,𝑗𝑢𝑖‒1,𝑗+𝑎𝑖,𝑗‒1𝑢𝑖,𝑗‒1+𝑎𝑖+1,𝑗𝑢𝑖+1,𝑗+𝑎𝑖,𝑗+1𝑢𝑖,𝑗+1)=𝐹𝑖𝑗2{𝑎𝑖‒1,𝑗=ℎ‒2(𝐴𝑖‒1/2,𝑗+ℎ𝐶𝑖𝑗2)𝑎𝑖,𝑗‒1=ℎ‒2(𝐵𝑖,𝑗‒1/2+ℎ𝐷𝑖𝑗2)𝑎𝑖+1,𝑗=ℎ‒2(𝐴𝑖+1/2,𝑗‒ℎ𝐶𝑖𝑗2)𝑎𝑖,𝑗+1=ℎ‒2(𝐵𝑖,𝑗+1/2‒ℎ𝐷𝑖𝑗2)𝑎𝑖𝑗=ℎ‒2(𝐴𝑖+1/2,𝑗+𝐴𝑖‒1/2,𝑗+𝐵𝑖,𝑗‒1/2+𝐵𝑖,𝑗+1/2)+𝐸𝑖𝑗 举一个例子：当i=2,j=3时，；当i=3,j=3时，。但𝑎𝑖𝑗=𝑎23𝑎𝑖‒1，𝑗=𝑎23是，显然这两个不是同一个数，其大小也不相等。𝑎232.3差分方程的重新定义因此，为了避免2.2中赋值重复而产生的错误，我在利用matlab编程时，对这些系数变量进行了重新定义：𝑏𝑖𝑗=𝑎𝑖𝑗，𝑐𝑖𝑗=𝑎𝑖,𝑗+1，𝑑𝑖𝑗=𝑎𝑖,𝑗‒1，𝑔𝑖𝑗=𝑎𝑖+1,𝑗，𝑘𝑖𝑗=𝑎𝑖‒1,𝑗.2.4模型建立这里的模型建立就是把差分方程组改写成系数矩阵的形式。经过研究，我觉得写成如下的系数矩阵不仅看起来简单明了，而且在matlab编程时比较方便。系数矩阵为：Pu=f其中P是阶方阵，具体如下：(𝑛‒1)2(𝑏11𝑐110𝑑12𝑏12⋱0⋱⋱00⋱𝑐1,𝑛‒2𝑑1,𝑛‒1𝑏1,𝑛‒1𝑔11𝑔12𝑔13⋱𝑔1,𝑛‒1𝑘21𝑘22𝑘2,3⋱𝑘2,𝑛‒1𝑏21𝑐210𝑑22𝑏22⋱0⋱⋱00⋱𝑐2,𝑛‒2𝑑1,𝑛‒1𝑏2,𝑛‒100⋱𝑔𝑛‒2,1𝑔𝑛‒2,2𝑔𝑛‒23⋱𝑔𝑛‒2,,𝑛‒1𝑘𝑛‒1,1𝑘𝑛‒1,2𝑘𝑛‒1,3⋱𝑘𝑛‒1,𝑛‒1𝑏𝑛‒1,1𝑐𝑛‒1,10𝑑𝑛‒1,2𝑏𝑛‒1,2⋱0⋱⋱00⋱𝑐𝑛‒1,𝑛‒2𝑑𝑛‒1,𝑛‒1𝑏𝑛‒1,𝑛‒1)3而u是维的列向量，具体如下：(𝑛‒1)2u=(𝑢11𝑢12⋮𝑢1,𝑛‒1𝑢21⋮⋮𝑢𝑛‒1,𝑛‒1)而f是维的列向量，具体如下：(𝑛‒1)2𝑓=(𝑓11𝑓12⋮𝑓1,𝑛‒1𝑓21⋮⋮𝑓𝑛‒1,𝑛‒1)三、求解过程3.1对系数矩阵的分析对上述模型的求解就是对线性方程组的求解。通过观察，我发现P是一个对角占优的矩阵，这不仅确定了解的唯一性，还保证了迭代法的收敛性。此外，还可以确定进行LU分解，若使用高斯消去法还可以省去选主元的工作。3.2matlab编程因为矩阵阶数过大，所以此题的编程难点为构造系数矩阵，即对线性方程组的赋值。我采用的方法是分块赋值。对于P的赋值，过程如下：第一步：𝑏𝑐𝑑𝑖=(𝑏𝑖1‒𝑐𝑖10‒𝑑𝑖2𝑏𝑖2⋱0⋱⋱00⋱‒𝑐𝑖,𝑛‒2‒𝑑𝑖,𝑛‒1𝑏𝑖,𝑛‒1),𝑔𝑖=[𝑔𝑖1𝑔𝑖2⋮𝑔𝑖,𝑛‒1],𝑘𝑖=[𝑘𝑖1𝑘𝑖2⋮𝑘𝑖,𝑛‒1]4第二步：𝐵𝐶𝐷=(𝑏𝑐𝑑1𝑏𝑐𝑑200⋱𝑏𝑐𝑑𝑖),𝐺=[𝑔1𝑔2⋮𝑔𝑛‒2],𝐾=[𝑘2𝑘3⋮𝑘𝑛‒1]第三步：P=BCD-diag(G,99)-diag(K,99).P和f的具体构造见附录6.1主代码3.3编程求解过程中的问题3.3.1问题产生当按照老师要求，n=100，h=1/100时，运行编好的matlab程序时，会出现如图3.1的错误提示。图3.13.3.2问题分析在matlab的命令窗口输入“memory”，出现如图3.2的内存使用情况，可以得出：MemoryusedbyMATLAB:454MB(4.760e+008bytes)。，若不用稀疏矩阵定义P，经过粗略计算，我发现矩阵P就要占800MB左右的内存，加上其他数据，内存消耗至少在1G以上。但是我电脑上分配给matlab的内存只有：454MB，即使在关闭杀毒软件等大部分应用程序后，分配给matlab的内存也刚够1G。图3.23.3.3问题解决经过上网查找资料后，我找到了如下几个解决方法。1)尽量早的分配大的矩阵变量2)尽量避免产生大的瞬时变量，当它们不用的时候应该及时clear。3)尽量的重复使用变量(跟不用的clear掉一个意思)4)将矩阵转化成稀疏形式5)使用pack命令56)如果可行的话，将一个大的矩阵划分为几个小的矩阵，这样每一次使用的内存减少。7)增大内存针对本题，我觉得比较理想的解决方法是采用稀疏矩阵的方式定义P。这样可以有效的减小P的内存消耗。但是考虑到老师的这次期中作业主要是考察我们对二阶椭圆偏微分方程的理解与实例操作，而不是旨在考察我们的matlab编程能力。因此我在此，略作偷懒，把n改成10或70(75以上内存就不够用了)，适当的降低精度来得到结果。四、计算结果4.1当n=10,h=1/10时的结果取n=10,h=1/10时，matlab运行的部分结果如图4.1。表4.2为LU分解法和高斯消去法的部分结果(这两个直接法结果完全一样)，表4.3为迭代法的部分结果。图4.1i,j数值解真实值误差1，11.0100501453351.0100501670840.0000000217491，21.0202012644381.0202013400270.0000000