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

1、微分方程数值解法期中作业实验报告二阶椭圆偏微分方程第一边值问题姓名： 学号：班级：2013年11月19日二阶椭圆偏微分方程第一边值问题摘要对于解二阶椭圆偏微分方程第一边值问题，课本上已经给出了相应的差分方程。而留给我的难题就是把差分方程组表示成系数矩阵的形式，以及对系数进行赋值。解决完这个问题之后，我在利用matlab解线性方程组时，又出现“out of memory”的问题。因为99*99阶的矩阵太大，超出了分配给matlab的使用内存。退而求其次，当n=10，h=1/10或n=70，h=1/70时，我都得出了很好的计算结果。然而在解线性方程组时，无论是LU分解法或高斯消去法，还是gause

2、idel迭代法，都能达到很高的精度。 关键字：二阶椭圆偏微分方程 差分方程 out of memory LU分解 高斯消去法 gauseidel迭代法一、题目重述解微分方程：已知边界:求数值解, 把区域分成，n=100注：老师你给的题F好像写错了，应该把改成。二、问题分析与模型建立2.1微分方程上的符号说明Ax,y=ey Bx,y=ex Cx,y=x+y Dx,y=x-y Ex,y=1 Fx,y=2.2课本上差分方程的缺陷课本上的差分方程为：aijuij-ai-1,jui-1,j+ai,j-1ui,j-1+ai+1,jui+1,j+ai,j+1ui,j+1=Fijai-1,j=h-2Ai-1/

3、2,j+hCij2ai,j-1=h-2Bi,j-1/2+hDij2ai+1,j=h-2Ai+1/2,j-hCij2ai,j+1=h-2Bi,j+1/2-hDij2aij=h-2Ai+1/2,j+Ai-1/2,j+Bi,j-1/2+Bi,j+1/2+Eij举一个例子：当i=2,j=3时，aij=a23；当i=3,j=3时，ai-1，j=a23。但是，显然这两个a23不是同一个数，其大小也不相等。2.3差分方程的重新定义因此，为了避免2.2中赋值重复而产生的错误，我在利用matlab编程时，对这些系数变量进行了重新定义：bij=aij，cij=ai,j+1，dij=ai,j-1，gij=ai+1,

4、j，kij=ai-1,j.2.4模型建立 这里的模型建立就是把差分方程组改写成系数矩阵的形式。经过研究，我觉得写成如下的系数矩阵不仅看起来简单明了，而且在matlab编程时比较方便。系数矩阵为：Pu=f其中P是（n-1）2阶方阵，具体如下：b11c110d12b12000c1,n-2d1,n-1b1,n-1g11g12g13g1,n-1k21k22k2,3k2,n-1b21c210d22b22000c2,n-2d1,n-1b2,n-100gn-2,1gn-2,2gn-23gn-2,n-1kn-1,1kn-1,2kn-1,3kn-1,n-1bn-1,1cn-1,10dn-1,2bn-1,2000

5、cn-1,n-2dn-1,n-1bn-1,n-1而u是（n-1）2维的列向量，具体如下：u=u11u12u1,n-1u21un-1,n-1而f是（n-1）2维的列向量，具体如下：f=f11f12f1,n-1f21fn-1,n-1三、求解过程3.1对系数矩阵的分析对上述模型的求解就是对线性方程组的求解。通过观察，我发现P是一个对角占优的矩阵，这不仅确定了解的唯一性，还保证了迭代法的收敛性。此外，还可以确定进行LU分解，若使用高斯消去法还可以省去选主元的工作。3.2matlab编程因为矩阵阶数过大，所以此题的编程难点为构造系数矩阵，即对线性方程组的赋值。我采用的方法是分块赋值。对于P的赋值，过程如

6、下：第一步：bcdi=bi1-ci10-di2bi2000-ci,n-2-di,n-1bi,n-1,gi=gi1gi2gi,n-1,ki=ki1ki2ki,n-1第二步：BCD=bcd1bcd200bcdi ,G=g1g2gn-2,K=k2k3kn-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的内存使用

7、情况，可以得出：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）将矩阵转化成稀

8、疏形式5）使用pack命令6）如果可行的话，将一个大的矩阵划分为几个小的矩阵，这样每一次使用的内存减少。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分解法和高斯消去法的部分结果（这两个直接法

9、结果完全一样），表4.3为迭代法的部分结果。图4.1i,j数值解真实值误差1，11.010050145335 1.010050167084 0.000000021749 1，21.020201264438 1.020201340027 0.000000075589 1,31.030454341388 1.030454533954 0.000000192565 5,71.419066053043 1.419067548593 0.000001495550 5,81.491822786765 1.491824697641 0.000001910877 5,91.568310733070 1.568

10、312185490 0.000001452420 6,11.061837559575 1.061836546545 0.000001013029 6,21.127498750514 1.127496851579 0.000001898934 6,31.197219794676 1.197217363122 0.000002431554 6,41.271251608972 1.271249150321 0.000002458650 9,71.877615353384 1.877610579264 0.000004774120 9,82.054437020906 2.054433210644 0.

11、000003810262 9,92.247910518362 2.247907986676 0.000002531686 表4.2i,j数值解真实值误差1，11.010049929223 1.010050167084 0.000000237861 1，21.020200873723 1.020201340027 0.000000466304 1,31.030453831277 1.030454533954 0.000000702677 5,51.284024086148 1.284025416688 0.000001330540 5,61.349856719969 1.349858807576

12、 0.000002087607 5,71.419064868769 1.419067548593 0.000002679825 5,81.491821975367 1.491824697641 0.000002722274 5,91.568310332118 1.568312185490 0.000001853372 6,11.061837000239 1.061836546545 0.000000453693 6,21.127497727757 1.127496851579 0.000000876177 6,31.197218445256 1.197217363122 0.000001082

13、134 9,71.877615007879 1.877610579264 0.000004428615 9,82.054436783189 2.054433210644 0.000003572545 9,92.247910400175 2.247907986676 0.000002413498 表4.34.2当 n=70,h=1/70时的结果取n=70,h=1/70时，matlab运行的部分结果如图4.4（LU分解法）。计算时间为17多分钟（1043秒），误差至少在小数点后9位。图4.4五、结论分析由于本人的电脑比较破旧，内存不是很大，再加上没有采取稀疏矩阵的存储方式，难以达到n=100，h=

14、1/100的计算要求。但改为n=10，h=1/10或n=70，h=1/70后，计算精度也十分理想。尤其是后者，其误差至少在小数点后9位， 在比较使用哪种方法解线性方程组时，直接法中的LU分解法和高斯消去法的计算结果是完全相等的。而gauseidel迭代法计算结果按个人设定的计算精度而定。对于这种大型线性方程组来说，迭代法从计算速度和计算机存储方面来看具有超过直接法的决定性优点。对于稀疏方程组来说，迭代法十分有效。从本题的实际情况来看，当n=70时，LU分解的直接法的计算时间为17分钟左右，而gauseidel迭代法的计算时间为20秒左右，这与以上分析的情况一致。但是当n越来越大时（从n=10起

15、），固定迭代步数的迭代法解的精度越来越差，为了达到与直接法相同的精度，必须提高迭代步数。然而这又会加大计算量，使计算速度变慢（见表5.1）。所以迭代法是在计算精度和计算速度之间的权衡取舍。n=50，h=1/50迭代精度迭代步数迭代时间1.00E-04155513.3秒1.00E-06268121.6秒1.00E-08380829.8秒表5.1仅从本题计算结果来看（n=10时）：当x，y靠近0时，直接法的数值解更准确，而当x，y靠近1时，迭代法的数值解更准确。这与gauseidel迭代法的算法特点相符合。因为gauseidel迭代法后面的解在迭代时要用到前面的解，所以x，y靠近1的后面的解更为精

16、确。六、附录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)(y2+x2+1)*exp(x*y)-(y2*exp(y)+x2*exp(x)*exp(x*y);FU=(x,y)exp(x*y);%真实值n=t;%区间划分为n*nh=1/n;%h为单位区间长

17、度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);%最终要解的方程组的系数矩阵,即平时的Af=zeros(n-1)*(n-1),1); %即平时的bff=zeros(n-1); %ff(i,j)=f(i-1）*9+j)BCD=;%记录三对角部分G=;%记录上三角的那一斜条K=;%记

18、录下三角的那一斜条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); endendx=zeros(n-1,

19、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); endend%对最终要解的方程组的系数矩阵进行赋值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

20、*j-1)/h2+E(i,j); c(i,j)=(B(2*i,2*j+1)-h*D(i,j)/2)/h2; d(i,j)=(B(2*i,2*j-1)+h*D(i,j)/2)/h2; g(i,j)=(A(2*i+1,2*j)-h*C(i,j)/2)/h2; k(i,j)=(A(2*i-1,2*j)+h*C(i,j)/2)/h2; bb=bb b(i,j); 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,

21、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=2 K=K kk; end endP=BCD-diag(G,n-1)-diag(K,-n+1

22、);%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); endenderror=abs(x-u);result=x u error;time = toc;disp(计算时间为：

23、);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的对角元素为1U = 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) = (