二阶椭圆偏微分方程实例求解附matlab代码Word格式.docx

1、Pu=f其中P是阶方阵，具体如下：而u是维的列向量，具体如下：而f是三、求解过程3.1对系数矩阵的分析对上述模型的求解就是对线性方程组的求解。通过观察，我发现P是一个对角占优的矩阵，这不仅确定了解的唯一性，还保证了迭代法的收敛性。此外，还可以确定进行LU分解，若使用高斯消去法还可以省去选主元的工作。3.2matlab编程因为矩阵阶数过大，所以此题的编程难点为构造系数矩阵，即对线性方程组的赋值。我采用的方法是分块赋值。对于P的赋值，过程如下：第一步：第二步：第三步： P=BCD-diag（G,99）-diag（K,99）.P和 f的具体构造见附录6.1主代码3.3编程求解过程中的问题 3.3.1

2、问题产生当按照老师要求，n=100，h=1/100时，运行编好的matlab程序时，会出现如图3.1的错误提示。图3.13.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问题解决经过

3、上网查找资料后，我找到了如下几个解决方法。1）尽量早的分配大的矩阵变量2）尽量避免产生大的瞬时变量，当它们不用的时候应该及时clear。3）尽量的重复使用变量（跟不用的clear掉一个意思）4）将矩阵转化成稀疏形式5）使用pack命令6）如果可行的话，将一个大的矩阵划分为几个小的矩阵，这样每一次使用的内存减少。7）增大内存针对本题，我觉得比较理想的解决方法是采用稀疏矩阵的方式定义P。这样可以有效的减小P的内存消耗。但是考虑到老师的这次期中作业主要是考察我们对二阶椭圆偏微分方程的理解与实例操作，而不是旨在考察我们的matlab编程能力。因此我在此，略作偷懒，把n改成10或70（75以上内存就不够

4、用了），适当的降低精度来得到结果。四、计算结果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.010050145335 1.010050167084 0.000000021749 1，21.020201264438 1.020201340027 0.000000075589 1,31.030454341388 1.030454533954 0.000000192565 5,71.4190660530

5、43 1.419067548593 0.000001495550 5,81.491822786765 1.491824697641 0.000001910877 5,91.568310733070 1.568312185490 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.2712491

6、50321 0.000002458650 9,71.877615353384 1.877610579264 0.000004774120 9,82.054437020906 2.054433210644 0.000003810262 9,92.247910518362 2.247907986676 0.000002531686 表4.21.010049929223 0.000000237861 1.020200873723 0.000000466304 1.030453831277 0.000000702677 5,51.284024086148 1.284025416688 0.000001

7、330540 5,61.349856719969 1.349858807576 0.000002087607 1.419064868769 0.000002679825 1.491821975367 0.000002722274 1.568310332118 0.000001853372 1.061837000239 0.000000453693 1.127497727757 0.000000876177 1.197218445256 0.000001082134 1.877615007879 0.000004428615 2.054436783189 0.000003572545 2.247

8、910400175 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=1/100的计算要求。但改为n=10，h=1/10或n=70，h=1/70后，计算精度也十分理想。尤其是后者，其误差至少在小数点后9位， 在比较使用哪种方法解线性方程组时，直接法中的LU分解法和高斯消去法的计算结果是完全相等的。而gauseide

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

10、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的后面的解更为精确。六、附录6.1主代码主代码中解决了对系数矩阵的赋值，即写出了线性方程组。在解线性方程组时可以选用LU分解法、高斯消去法和迭代法中的任意一种。function n=Middle_Term_Bymyself（t）%t为区间划分数tic;%定义函数

11、及初始化基本变量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为单位区间长度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）

12、,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=;%记录下三角的那一斜条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-

13、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,1）;y=zeros（n-1,1）;n-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,

14、j）=FU（x（i）,y（j）;%对最终要解的方程组的系数矩阵进行赋值 bcd=zeros（n-1）; bb=; cc=; dd=; gg=; kk=; 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）/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=b

15、b b（i,j）; if j=2 dd=dd d（i,j）; 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）; if j=1 ff（i,j）=ff（i,j）+d（i,j）*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）;%你应该懂的，坐

16、标变换 bcd=diag（bb）-diag（cc,1）-diag（dd,-1）; BCD=blkdiag（BCD,bcd）; if i K=K kk;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迭代法 u（i

17、-1）*（n-1）+j）=U（i,j）;error=abs（x-u）;result=x u error;time = toc;disp（计算时间为：）;disp（time）;-format long;计算结果为： 数值解 真实值 误差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=

18、2:n for i=k: U（k,i） = A（k,i）-L（k,1:（k-1）*U（1:（k-1）,i）; %U的第k行 for j=（k+1）: L（j,k） = （A（j,k）-L（j,1:（k-1）,k）/U（k,k）; %L的第k列y = SolveDownTriangle（L,b）;x = SolveUpTriangle（U,y）; %求解方程6.3高斯消去法function x,XA=GaussXQByOrder（A,b）（n-1） for j=（i+1）: if（A（i,i）=0） disp（对角元素为0！ %防止对角元素为0 return; l = A（j,i）; m = A

19、（i,i）; A（j,1:n）=A（j,1:n）-l*A（i,1:n）/m; %消元方程 b（j）=b（j）-l*b（i）/m;x=SolveUpTriangle（A,b）; %通用的求上三角系数矩阵线性方程组的函数XA = A; %消元后的系数矩阵6.4上三角解线性方程组（LU分解法、高斯消去法要用到这个算法）function x=SolveUpTriangle（A,b）N=size（A）;n=N（1）;for i=n:-1:1if（i1） s=A（i,1:（i-1）*x（1:（i-1）,1）;s=0;6.6gauseidel迭代法if nargin=3 eps= 1.0e-6;%可以安自己的精度要求改变 M = 10000;elseif nargin = 4 M = 10000;elseif nargin3 error