二阶椭圆偏微分方程实例求解(附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问题解决

经过上网查找资料后，我找到了如 …… 此处隐藏：5424字，全部文档内容请下载后查看。喜欢就下载吧 ……