二阶椭圆偏微分方程实例求解(附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分解法）。计算时间为17多分钟（1043秒），误差至少在小数点后9位。

图4.4

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

五、结论分析

由于本人的电脑比较破旧，内存不是很大，再加上没有采取稀疏矩阵的存储方式，难以达到n=100，h=1/100的计算要求。但改为n=10，h=1/10或n=70，h=1/70后，计算精度也十分理想。尤其是后者，其误差至少在小数点后9位，

在比较使用哪种方法解线性方程组时，直接法中的LU分解法和高斯消去法的计算结果是完全相等的。而gauseidel迭代法计算结果按个人设定的计算精度而定。对于这种大型线性方程组来说，迭代法从计算速度和计算机存储方面来看具有超过直接法的决定性优点。对于稀疏方程组来说，迭代法十分有效。从本题的实际情况来看，当n=70时，LU分解的直接法的计算时间为17分钟左右，而gauseidel迭代法的计算时间为20秒左右，这与以上分析的情况一致。但是当n越来越大时（从n=10起），固定迭代步数的迭代法解的精度越来越差，为了达到与直接法相同的精度，必须提高迭代步数。然而这又会加大计算量，使计算速度变慢（见表5.1）

仅从本题计算结果来看（n=10时）：当x，y靠近0时，直接法的数值解更准确，而当x，y靠近1时，迭代法的数值解更准确。这与gauseidel迭代法的算法特点相符合。因为gauseidel迭代法后面的解在迭代时要用到前面的解，所以x，y靠近1的后面的解更为精确。

六、附录

6.1主代码

主代码中解决了对系数矩阵的赋值，即写出了线性方程组。在解线性方程组时可以选用LU分解法、高斯消去法和迭代法中的任意一种。

function n=Middle_Term_Bymyself(t)%t为区间划分数 tic;

%定义函数及初始化基本变量 FA=@(x,y)exp(y); FB=@(x,y)exp(x); FC=@(x,y)x+y; FD=@(x,y)x-y; FE=@(x,y)1;

FF=@(x,y)(y^2+x^2+1)*exp(x*y)-(y^2*exp(y)+x^2*exp(x))*exp(x*y); FU=@(x,y)exp(x*y);%真实值 n=t;%区间划分为n*n

h=1/n;%h为单位区间长度

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

A=zeros(2*n-1); B=zeros(2*n-1); C=zeros(n-1); D=zeros(n-1); E=zeros(n-1); F=zeros(n-1);

U=zeros(n-1);%真实值的矩阵表示

u=zeros((n-1)*(n-1),1);%真实值的数列表示 error=zeros((n-1)*(n-1),1);%误差

P=zeros((n-1)*(n-1));%最终要解的方程组的系数矩阵,即平时的A f=zeros((n-1)*(n-1),1); %即平时的b ff=zeros(n-1); %ff(i,j)=f((i-1）*9+j) BCD=[];%记录三对角部分 G=[];%记录上三角的那一斜条 K=[];%记录下三角的那一斜条 b=zeros(n-1);%表示a(i,j) c=zeros(n-1);%表示a(i,j+1) d=zeros(n-1);%表示a(i,j-1) g=zeros(n-1);%表示a(i+1,j) k=zeros(n-1);%表示a(i-1,j) %对函数进行赋值 x=zeros(2*n-1,1); y=zeros(2*n-1,1);

for i=1:2*n-1 %使A.B下标i-1/2变为2i-1 for j=1:2*n-1 x(i)=i*h/2; y(j)=j*h/2;

A(i,j)=FA(x(i),y(j)); B(i,j)=FB(x(i),y(j)); end end

x=zeros(n-1,1); y=zeros(n-1,1); for i=1:n-1 for j=1:n-1 x(i)=i*h; y(j)=j*h;

C(i,j)=FC(x(i),y(j)); D(i,j)=FD(x(i),y(j)); E(i,j)=1;

F(i,j)=FF(x(i),y(j)); U(i,j)=FU(x(i),y(j)); end end

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

%对最终要解的方程组的系数矩阵进行赋值 for i=1:n-1

bcd=zeros(n-1); bb=[]; cc=[]; dd=[]; gg=[]; kk=[]; for j=1:n-1

b(i,j)=(A(2*i+1,2*j)+A(2*i-1,2*j)+B(2*i,2*j+1)+B(2*i,2*j-1))/h^2+E(i,j); c(i,j)=(B(2*i,2*j+1)-h*D(i,j)/2)/h^2;

d(i,j)=(B(2*i,2*j-1)+h*D(i,j)/2)/h^2; g(i,j)=(A(2*i+1,2*j)-h*C(i,j)/2)/h^2;

k(i,j)=(A(2*i-1,2*j)+h*C(i,j)/2)/h^2; bb=[bb b(i,j)]; if j<=n-2

cc=[cc c(i,j)]; end if j>=2

dd=[dd d(i,j)]; end

gg=[gg g(i,j)]; kk=[kk k(i,j)];

%给f赋值 if i==1

ff(i,j)=F(i,j)+k(i,j)*1;%边值为1 elseif i==n-1

ff(i,j)=F(i,j)+g(i,j)*A(i,2*j);%A中i取值无所谓，不影响 else

ff(i,j)=F(i,j); end if j==1

ff(i,j)=ff(i,j)+d(i,j)*1;%边值为1 elseif j==n-1

ff(i,j)=ff(i,j)+c(i,j)*B(2*i,j); end

f((i-1)*(n-1)+j,1) = ff(i,j);%你应该懂的，坐标变换 end

bcd=diag(bb)-diag(cc,1)-diag(dd,-1); BCD=blkdiag(BCD,bcd); if i<=n-2

G=[G gg]; end

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

if i>=2

K=[K kk]; end end

P=BCD-diag(G,n-1)-diag(K,-n+1);

%BCD=BCD-diag(G,n-1)-diag(K,-n+1);x=Doolittle(BCD,f);这样是不是可以减少点内存消耗

x=Doolittle(P,f);%LU分解解方程组，这里也可以输入其他解方程组的方法 %x=GaussXQByOrder(P,f) %高斯消去法

%x0=ones((n-1)^2,1);x=gauseidel(P,f,x0) %gauseidel迭代法 for i=1:n-1 for j=1:n-1

u((i-1)*(n-1)+j)=U(i,j); end end

error=abs(x-u); result=[x u error]; time = toc;

disp('计算时间为：'); disp(time);

disp('---------------------------------------------------------------'); format long;

disp('计算结果为：');

disp(' 数值解真实值误差'); disp(result);

6.2LU分解解线性方程组

function [x,L,U]= Doolittle (A,b) N = size(A); n = N(1);

L = eye(n,n); %L的对角元素为1 U = zeros(n,n);

U(1,1:n) = A(1,1:n); %U的第一行 L(1:n,1) = A(1:n,1)/U(1,1); %L的第一列

for k=2:n for i=k:n

U(k,i) = A(k,i)-L(k,1:(k-1))*U(1:(k-1),i); %U的第k行 end

for j=(k+1):n

L(j,k) = (A(j,k)-L(j,1:(k-1))*U(1:(k-1),k))/U(k,k);

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

%L的第k列 end end

y = SolveDownTriangle(L,b);

x = SolveUpTriangle(U,y); %求解方程

6.3高斯消去法

function [x,XA]=GaussXQByOrder(A,b) N = size(A); n = N(1);

for i=1:(n-1) for j=(i+1):n if(A(i,i)==0)

disp('对角元素为0！'); %防止对角元素为0 return; end

l = A(j,i); m = A(i,i);

A(j,1:n)=A(j,1:n)-l*A(i,1:n)/m; %消元方程 b(j)=b(j)-l*b(i)/m; end end

x=SolveUpTriangle(A,b); %通用的求上三角系数矩阵线性方程组的函数 XA = A; %消元后的系数矩阵

6.4上三角解线性方程组（LU分解法、高斯消去法要用到这个算法） function x=SolveUpTriangle(A,b) N=size(A); n=N(1); for i=n:-1:1 if(i<n)

s=A(i,(i+1):n)*x((i+1):n,1); else s=0; end

x(i,1)=(b(i)-s)/A(i,i); end

6.5下三角解线性方程组（LU分解法要用到这个算法） function x=SolveDownTriangle(A,b)

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

N=size(A); n=N(1); for i=1:n if(i>1)

s=A(i,1:(i-1))*x(1:(i-1),1); else s=0; end

x(i,1)=(b(i)-s)/A(i,i); end

6.6gauseidel迭代法 ifnargin==3

eps= 1.0e-6;%可以安自己的精度要求改变 M =10000;

elseifnargin == 4 M =10000;

elseifnargin<3 error return; end

D=diag(diag(A)); %求A的对角矩阵 L=-tril(A,-1); %求A的下三角阵 U=-triu(A,1); %求A的上三角阵 G=(D-L)\U; f=(D-L)\b; x=G*x0+f;

n=1; %迭代次数 while norm(x-x0)>=eps x0=x;

x=G*x0+f; n=n+1; if(n>=M)

disp('Warning: 迭代次数太多，可能不收敛！'); return; end end

disp('迭代次数为：'); disp(n);