重 庆 大 学

运行结果 当N=100直接求解的发发还可行
当N=120时，会出现以下结果
Elliptic1(120) 错误使用 diag 内存不足。请键入 HELP MEMORY 查看选项。
出错 Elliptic1 (line 14) A=A+diag(-ones(N*(N-1),1),N)+diag(-ones(N*(N-1),1),-N);
0
for (x, y) D
正方形区域 D 和边间条件如下：
g(x) sin(5 x)
要求： 寻找一种数值算法，尽可能的让迭代步长变小即尽可能的让网格数 N 变大。
2 / 14'.
.
三、实验原理、方法（算法）、步骤 实验一： 将方程化为标准形式：
d dx
u(a)
( p(x)
,
du ) dx
q( x)u
du dx
f ( x),
a x b,
(2) (3)
p(b)u
'(b)
u(b)
,
(4)
其中
p(x) 1,q(x) 1,f(x) sin(2 x),g(x) 1
a 0,b 1, 1, 0, 2 1 4 2
第一步：考虑从Galerkin法出发建立有限元方程。
为常数
考虑一阶Sobolev 空间
二． 实验内容 实验一： 考虑两点边值问题：
y y sin(2 x) for 0 x 1 (1)
边界条件为：
y(0) 0
y(1)
0
真解：
y
sin(2 1 4
x)
2
运用有限元的方法求解该方程的数值解，并和真解比较。
实验二：
用五点差分格式近似 Laplacian 方程：
2u x 2
2u x 2
g=sin(5*pi*x);
a=zeros(N,N);
%生成系数矩阵
b=ones(N,1)*(2*(1+l^2));c=ones(N-1,1)*(-l^2);
a=diag(b)+diag(c,1)+diag(c,-1);
A=zeros(N^2,N^2);
for i=1:N
A((i-1)*N+1:i*N,(i-1)*N+1:i*N)=a;
实验项目 名称
偏微分方程期末课程设计
实验项目类型
验证 演示 综合 设计 其他
指导教师
成绩
一． 实验目的
自学，掌握有限元分析的基本理论，并运用有限元分析的方法求解第二章的两点 边值问题，做出数值解，体会有限元和差分方法的不同之处。
掌握平面上拉普拉斯方程的五点差分方法，体会与一般一维问题的不同，特别注 意边界条件的处理。学会处理大型方程组数值解的压缩存储方法。
.
重庆大学 学生实验报告
实验课程名称
偏微分方程数值解期末课程设计
开课实验室
偏微分方程数值解
学
院 数统学院 年级 2011 专业班
学生姓名
学号
开课时间
2013 至 2014 学年第 2 学期
总成绩 教师签名
数学与统计学院制
1 / 14'.
.
开课学院、实验室：
实验时间 ： 2014 年 6 月 30 日
3 / 14'.
.
第二步：网格剖分 对区间[0，1] 作等距剖分，令步长 h=1/n ， n 为网格数。则有
有限元方程的解 u 在每一个分段的小区间上为分片的线性函数
网格点
上分别取值为
，设 在
取基函数：
容易发现, 基函数有如下性质
显然，基函数线性无关，对任意的一个
可表示为
n
uh (x) uii (x) ui uh (xi ) i 1
可以发现整个矩阵只由三个不同的元素构成，所以我们不去存储完整的矩阵，而 是只保存这个三个不同的元素，即可得全部的信息。
令
a , b= 2 , c= 1,
在共轭梯度法的迭代计算中需要不停地计算 A*x,因此只要解决了 A*x 的问题， 则算法即可进行。
对于 A*x=y,由 A 的特殊性我们可以归纳计算如下（当网格大小为 N 时）： 当 i=1
当i N ^2
cxiN bxi1 axi y(i)
因此这样就解决了矩阵 A 存储不下的尴尬情况，大大的减少了程序所占内存的大
小。几乎几秒钟就能算出网格数 N 取好几百时的情况。
四．实验环境（所用软件、硬件等）及实验数据文件
MATLAB2014a
Window7 内存 4G 三． 实验结果及实例分析
h xi1
2
h2
xi1 1 (x xi1)2 dx
xi
h2
h2
i=1 n-1
=
(
xi
3 1
3xi12 xi
3xi 1 xi 2
2 xi 3
3xi2 xi1
3xi
xi
2 1
x3 i 1
3xi1
3xi 1 )
3h2
a(n ,n )
xn 1 (x xn1)2 dx
h xn1
2
h2
=
xn3
x3 n1
n
un (x) uii (x) i 1
实验二; 用中心差分代替 Laplacian 方程的二阶导数，舍去截断误差得到数值算法；
2ui 1,j 2(1 2 )ui,j 2ui 1,j ui,j 1 ui,j 1 0
For
Where
边界条件为：
(1) x = 0:
, for j = 1,..., M,
ax1 bx2 cxN1 y(1) 当 2 i N 1
当 i=N
bxi1 axi bxi1 cxN i y(i)
7 / 14'.
.
bxi1 axi cxN i y(i) 当i N 1
cxiN axi bxi1 cxiN y(i) 当 N 2 i N(N1) 1
;
A(i,i-1)= A(i-1,i);
end
for i=2:n
%生成右端常数项
b(i-1)=1/(2*h*pi^2)*(-x(i-1)*cos(2*pi*x(i))*pi+2*cos(2*pi*x(i))*pi*x(i)x(i+1)*cos(2*pi*x(i))*pi+sin(pi*x(i-1))*cos(pi*x(i-1))+sin(pi*x(i+1))*co s(pi*x(i+1))-sin(2*pi*x(i))); end b(n)=1/(4*h*pi^2)*(-2*x(n)*cos(2*pi*x(n+1))*pi+2*cos(2*pi*x(n+1))*pi*x(n +1)+2*sin(pi*x(n))*cos(pi*x(n))-sin(2*pi*x(n+1)))-2.03*pi/(1+4*pi^2); y=myconjgrad(A,b) ; %调用求解方程组的算法 u1=zeros(n+1,1); u1(2:n+1)=y; plot(x,u1,'--b');hold on; x=0:0.01:1; y=sin(2*pi*x)/(1+4*pi^2); plot(x,y,'-r');hold on; legend('有限元','精确解');
h2
=
1 h
x3 i 1
xi3
3xi2 xi1 6h2
3xi
xi
2 1
i=1 n-1
5 / 14'.
.
F (i ) -
xi sin(2
xi1
x)(x xi1)dxh
xi1 sin(2
xi
x)(xi1 x)dx h
=
xi 1cos(2
xi
)
2cos(2
xi )
xi
ccos(2
xi )
sin( h 2
A(i-1,i-1)=A(1,1); end
8 / 14'.
. A(n,n)=-1/h+1/h^2*(x(n+1)^3/3-x(n)^3/3-x(n)*x(n+1)^2+x(n)^2*x(n+1)); A(n,n)=A(1,1)/2; for i=2:n
A(i-1,i)=1/h-1/h^2*(x(i+1)^3/6-x(i)^3/6+x(i)^2*x(i+1)/2-x(i)*x(i+1)^2/2)
cxiN bxi1 axi bxi1 cxiN y(i) 当 i N (N1)
cxiN bxi1 axi cxiN y(i) 当 i N(N1) 1
cxiN axi bxi1 y(i) 当 N(N1) 2 i N ^ 2 1
cxiN bxi1 axi bxi1 y(i)
(2) x = 1:
, for j = 1,..., M,
(3) y = 0:
, for i = 1,..., N,
(4) y = 1:
,for i = 1,..., N.
可以改写为矩阵形式为：
其中
6 / 14'.
.
对于对称正定大型稀疏矩阵问题 Ax=b，一般最有效的数值算法是共轭梯度法。 但是对以上问题用共轭梯度法最终还是发现当 N 取 120 左右的时候程序就不能运 行了，由此可以看出，系数矩阵 A 随着网格数的增加急剧的变大以至于电脑内存存不 下，因此有必要对系数矩阵进行特别处理。 观察矩阵网格数为三的系数阵
的子空间V:
任取 , 用它乘以方程 (1) 的两边, 然后在区间(a, b) 上积分
b d dy