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

二阶椭圆偏微分方程实例求解(附matlab代码)用的是五点差分法。

《微分方程数值解法》期中作业实验报告

二阶椭圆偏微分方程第一边值问题

姓名： 学号： 班级：

2013年11月19

日

二阶椭圆偏微分方程实例求解(附matlab代码)用的是五点差分法。

二阶椭圆偏微分方程第一边值问题

摘要

对于解二阶椭圆偏微分方程第一边值问题，课本上已经给出了相应的差分方程。而留给我的难题就是把差分方程组表示成系数矩阵的形式，以及对系数进行赋值。解决完这个问题之后，我在利用matlab解线性方程组时，又出现“out of memory”的问题。因为99*99阶的矩阵太大，超出了分配给matlab的使用内存。退而求其次，当n=10，h=1/10或n=70，h=1/70时，我都得出了很好的计算结果。然而在解线性方程组时，无论是LU分解法或高斯消去法，还是gauseidel迭代法，都能达到很高的精度。

关键字：二阶椭圆偏微分方程差分方程out of memory LU分解高斯消去法gauseidel迭代法

一、题目重述

解微分方程：

(eyux(x,y))x (exuy(x,y))y (x y)ux(x,y) (x y)uy(x,y) u(x,y) ye xe e y x 1 e

2y

2x

xy

2

2

xy

已知边界:u(0,y)=1,u(1,y)=ey,u(x,0)=1,u(x,1)=ex

求数值解, 把区域G=[0,1] [0,1]分成h1=1/100,h2=1/100，n=100 注：老师你给的题F好像写错了，应该把y2ex x2ey改成y2ey x2ex。

二、问题分析与模型建立

2.1微分方程上的符号说明

, = , = , = + , =

, =1 , = y2ey x2ex exy y2 x2 1 exy

2.2课本上差分方程的缺陷

课本上的差分方程为：

1, 1, + , 1 , 1+ +1, +1, + , +1 , +1 =

二阶椭圆偏微分方程实例求解(附matlab代码)用的是五点差分法。

1, = 2 1/2, + , 1= 2 , 1/2+ +1, = 2 +1/2,

2 2

2

, +1= 2 , +1/2 2

2

= +1/2, + 1/2, + , 1/2+ , +1/2 +

举一个例子：当i=2,j=3时， = 23；当i=3,j=3时， 1， = 23。但是，显然这两个 23不是同一个数，其大小也不相等。

2.3差分方程的重新定义

因此，为了避免2.2中赋值重复而产生的错误，我在利用matlab编程时，对这些系数变量进行了重新定义： = ， = , +1， = , 1， = +1, ， = 1, .

2.4模型建立

这里的模型建立就是把差分方程组改写成系数矩阵的形式。经过研究，我觉得写成如下的系数矩阵不仅看起来简单明了，而且在matlab编程时比较方便。

系数矩阵为：Pu=f

其中P是（n 1）阶方阵，具体如下：

11 11 110

120 12 12

13 0

1, 2

0

1, 11, 1 1, 1

21

21 210

220 22 22

2,30 2, 2

0 1, 1 2, 1 2, 1

23

2,, 1

0

1, 2 1, 1 1, 1

2

2,1

2,2

1,1

1,2

1,3

1, 1

1,1 1,20 1,1 1,2 0

而u是（n 1）维的列向量，具体如下：

2

二阶椭圆偏微分方程实例求解(附matlab代码)用的是五点差分法。

11 12

1, 1

= 21

1, 1

而f是（n 1）维的列向量，具体如下：

11 12

2

1, 1

=

21

1, 1

三、求解过程

3.1对系数矩阵的分析

对上述模型的求解就是对线性方程组的求解。通过观察，我发现P是一个对角占优的矩阵，这不仅确定了解的唯一性，还保证了迭代法的收敛性。此外，还可以确定进行LU分解，若使用高斯消去法还可以省去选主元的工作。

3.2matlab编程

因为矩阵阶数过大，所以此题的编程难点为构造系数矩阵，即对线性方程组的赋值。我采用的方法是分块赋值。对于P的赋值，过程如下：

第一步：

1 10 1 1

0 2 2 2 2

bcdi= 0, = , = , 2

, 1 , 10

, 1 , 1

第二步：

二阶椭圆偏微分方程实例求解(附matlab代码)用的是五点差分法。

bcd1

BCD=

第三步：

0bcd2

1 2

,G= 2 ,K= 3

2 1bcdi

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.1

3.3.2问题分析

在matlab的命令窗口输入“memory”，出现如图3.2的内存使用情况，可以得出：Memory used by MATLAB: 454 MB (4.760e+008 bytes)。，若不用稀疏矩阵定义P，经过粗略计算，我发现矩阵P就要占800MB左右的内存，加上其他数据，内存消耗至少在1G以上。但是我电脑上分配给matlab的内存只有：454 MB，即使在关闭杀毒软件等大部分应用程序后，分配给matlab的内存也刚够1G。

图3.2

3.3.3问题解决

经过上网查找资料后，我找到了如下几个解决方法。 1）尽量早的分配大的矩阵变量

2）尽量避免产生大的瞬时变量，当它们不用的时候应该及时clear。 3）尽量的重复使用变量（跟不用的clear掉一个意思） 4）将矩阵转化成稀疏形式 5）使用pack命令

6）如果可行的话，将一个大的矩阵划分为几个小的矩阵，这样每一次使用的内存减少。 7）增大内存

二阶椭圆偏微分方程实例求解(附matlab代码)用的是五点差分法。

针对本题，我觉得比较理想的解决方法是采用稀疏矩阵的方式定义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.1

二阶椭圆偏微分方程实例求解(附matlab代码)用的是五点差分法。

4.2当n=70,h=1/70时的结果

取n=70,h=1/70时，matlab运行的部分结果如图4.4（LU分解法）。计算时 …… 此处隐藏：2628字，全部文档内容请下载后查看。喜欢就下载吧 ……