二阶椭圆偏微分方程实例求解附matlab代码Word格式.docx
《二阶椭圆偏微分方程实例求解附matlab代码Word格式.docx》由会员分享，可在线阅读，更多相关《二阶椭圆偏微分方程实例求解附matlab代码Word格式.docx（17页珍藏版）》请在冰豆网上搜索。
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问题产生
当按照老师要求，n=100，h=1/100时，运行编好的matlab程序时，会出现如图3.1的错误提示。
图3.1
3.3.2问题分析
在matlab的命令窗口输入“memory”，出现如图3.2的内存使用情况，可以得出：
MemoryusedbyMATLAB:
454MB（4.760e+008bytes）。
，若不用稀疏矩阵定义P，经过粗略计算，我发现矩阵P就要占800MB左右的内存，加上其他数据，内存消耗至少在1G以上。
但是我电脑上分配给matlab的内存只有：
454MB，即使在关闭杀毒软件等大部分应用程序后，分配给matlab的内存也刚够1G。
图3.2
3.3.3问题解决
经过上网查找资料后，我找到了如下几个解决方法。
1）尽量早的分配大的矩阵变量
2）尽量避免产生大的瞬时变量，当它们不用的时候应该及时clear。
3）尽量的重复使用变量（跟不用的clear掉一个意思）
4）将矩阵转化成稀疏形式
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分解法和高斯消去法的部分结果（这两个直接法结果完全一样），表4.3为迭代法的部分结果。
图4.1
i,j
数值解
真实值
误差
1，1
1.010050145335
1.010050167084
0.000000021749
1，2
1.020201264438
1.020201340027
0.000000075589
1,3
1.030454341388
1.030454533954
0.000000192565
5,7
1.419066053043
1.419067548593
0.000001495550
5,8
1.491822786765
1.491824697641
0.000001910877
5,9
1.568310733070
1.568312185490
0.000001452420
6,1
1.061837559575
1.061836546545
0.000001013029
6,2
1.127498750514
1.127496851579
0.000001898934
6,3
1.197219794676
1.197217363122
0.000002431554
6,4
1.271251608972
1.271249150321
0.000002458650
9,7
1.877615353384
1.877610579264
0.000004774120
9,8
2.054437020906
2.054433210644
0.000003810262
9,9
2.247910518362
2.247907986676
0.000002531686
表4.2
1.010049929223
0.000000237861
1.020200873723
0.000000466304
1.030453831277
0.000000702677
5,5
1.284024086148
1.284025416688
0.000001330540
5,6
1.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.247910400175
0.000002413498
表4.3
4.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分解法和高斯消去法的计算结果是完全相等的。
而gauseidel迭代法计算结果按个人设定的计算精度而定。
对于这种大型线性方程组来说，迭代法从计算速度和计算机存储方面来看具有超过直接法的决定性优点。
对于稀疏方程组来说，迭代法十分有效。
从本题的实际情况来看，当n=70时，LU分解的直接法的计算时间为17分钟左右，而gauseidel迭代法的计算时间为20秒左右，这与以上分析的情况一致。
但是当n越来越大时（从n=10起），固定迭代步数的迭代法解的精度越来越差，为了达到与直接法相同的精度，必须提高迭代步数。
然而这又会加大计算量，使计算速度变慢（见表5.1）。
所以迭代法是在计算精度和计算速度之间的权衡取舍。
n=50，h=1/50
迭代精度
迭代步数
迭代时间
1.00E-04
1555
13.3秒
1.00E-06
2681
21.6秒
1.00E-08
3808
29.8秒
表5.1
仅从本题计算结果来看（n=10时）：
当x，y靠近0时，直接法的数值解更准确，而当x，y靠近1时，迭代法的数值解更准确。
这与gauseidel迭代法的算法特点相符合。
因为gauseidel迭代法后面的解在迭代时要用到前面的解，所以x，y靠近1的后面的解更为精确。
六、附录
6.1主代码
主代码中解决了对系数矩阵的赋值，即写出了线性方程组。
在解线性方程组时可以选用LU分解法、高斯消去法和迭代法中的任意一种。
functionn=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为单位区间长度
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）;
fori=1:
2*n-1%使A.B下标i-1/2变为2i-1
forj=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）;
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,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））/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=[bbb（i,j）];
ifj<
=n-2
cc=[ccc（i,j）];
ifj>
=2
dd=[ddd（i,j）];
gg=[ggg（i,j）];
kk=[kkk（i,j）];
%给f赋值
ifi==1
ff（i,j）=F（i,j）+k（i,j）*1;
%边值为1
elseifi==n-1
ff（i,j）=F（i,j）+g（i,j）*A（i,2*j）;
%A中i取值无所谓，不影响
else
ff（i,j）=F（i,j）;
ifj==1
ff（i,j）=ff（i,j）+d（i,j）*1;
elseifj==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）;
%你应该懂的，坐标变换
bcd=diag（bb）-diag（cc,1）-diag（dd,-1）;
BCD=blkdiag（BCD,bcd）;
ifi<
G=[Ggg];
ifi>
K=[Kkk];
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-1）*（n-1）+j）=U（i,j）;
error=abs（x-u）;
result=[xuerror];
time=toc;
disp（'
计算时间为：
'
）;
disp（time）;
---------------------------------------------------------------'
formatlong;
计算结果为：
数值解真实值误差'
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的第一列
fork=2:
n
fori=k:
U（k,i）=A（k,i）-L（k,1:
（k-1））*U（1:
（k-1）,i）;
%U的第k行
forj=（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）
forj=（i+1）:
if（A（i,i）==0）
disp（'
对角元素为0！
%防止对角元素为0
return;
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;
x=SolveUpTriangle（A,b）;
%通用的求上三角系数矩阵线性方程组的函数
XA=A;
%消元后的系数矩阵
6.4上三角解线性方程组（LU分解法、高斯消去法要用到这个算法）
functionx=SolveUpTriangle（A,b）
N=size（A）;
n=N
（1）;
fori=n:
-1:
1
if（i<
n）
s=A（i,（i+1）:
n）*x（（i+1）:
n,1）;
else
s=0;
x（i,1）=（b（i）-s）/A（i,i）;
6.5下三角解线性方程组（LU分解法要用到这个算法）
functionx=SolveDownTriangle（A,b）
if（i>
1）
s=A（i,1:
（i-1））*x（1:
（i-1）,1）;
s=0;
6.6gauseidel迭代法
ifnargin==3
eps=1.0e-6;
%可以安自己的精度要求改变
M=10000;
elseifnargin==4
M=10000;
elseifnargin<
3
error
升级会员 