两点边值微分方程的有限元分析

1、庆祝大学学生们实证报告实验过程名偏微分方程数值解最后过程设计课堂开始实验室偏微分方程的数值解法科学院首统大学2011专业班学生性是名字和学号第二学期，从2013年到2014年总绩效老师签字数学和统计大学系统开课学院，研究室：实验时间：2014年6月30日实验项目名称偏微分方程的最终过程设计实验项目类型验证做示范综合设计其他指导教师成就功绩一.实验目的理解有限元分析的基本理论，使用有限元分析方法解决第二章两点边值问题，进行数值分析，了解有限元和差分方法的区别。求解平面上拉普拉斯方程的五点差分方法与一般一维问题不同，特别注意边界条件的处理。学习处理大型方程数值解的压缩存储方法。二.实验内容实验1:

2、考虑两点边值问题：(1)边界条件包括：真正的解决方案：用有限元法将方程的数值解法与实际解法进行比较。实验2:使用五点差分格式近似Laplacian方程。For矩形区域d和边之间的条件如下：要求：找到一种数值算法，使迭代步长尽可能小，网格数尽可能大n。三、实验原理、方法(算法)、步骤实验1:使方程式成为标准形式：其中是常数第一步：考虑在Galerkin方法中建立有限元方程。考虑一阶Sobolev空间H1(a，b)的子空间V:v=v | vh1a，b，va=0V/v用于方程式(1)的两侧相乘，以及在间距(a，b)中积分左侧端点的点使用边界条件va=0和代替每个项目的系数逮捕令常识就是方程式(5)是

3、边值问题(1)的转移方程式第二步：网格间距0，1是等距分割，步骤h=1/n，n是网格数。都有Xi=ih，I=0，1，2n有限元方程的解u将每段社区之间分割的线性函数uh/VH，uh分别从栅格点0=x0，x1，Xi，xn=1设置为U0、u1、ui、un使用默认函数：很容易看出基本函数具有以下特性显然，基本函数是线性独立的，对于其中一个uhv，可以表示为所以如果你求uhv与查找U0、u1、ui和un相同以矩阵格式重写其中可以使用基本函数的特性简化矩阵。第三步：从上述推导中可以看出，刚度矩阵k是三对角矩阵，因此可以用“追赶法”求解有限元方程。求了每个节点的近似解ui后，边值问题的解法如下：实验2；用

4、中心差替换拉普拉斯方程的二阶导数，用舍去截断误差得到数值算法。For I=1，NJ=1，mWhere =ckh边界条件包括：(1) x=0: u0，j=0，for j=1，M、(2) x=1: uN 1，j=0，for j=1，M、(3) y=0: ui，0=0，for I=1，N、(4) y=1: ui，n 1=g (Xi)，for I=1，n以下是可替换的矩阵形式：Ax=b其中对称正定大型稀疏矩阵问题对于Ax=b，通常最有效的数值算法是共轭梯度法。但是，对于这些问题，通过使用共轭梯度法确定n在服用约120的情况下程序不运行，可以看出系数矩阵a需要特殊处理系数矩阵，因为随着网格数的增加，计算

5、机内存不足。观察矩阵栅格数为3的系数阵列您可以看到整个矩阵由三个不同的元素组成，因此您可以通过仅保存这三个不同的元素来获取所有信息，而不是保存整个矩阵。逮捕令共轭梯度法的迭代计算需要继续计算A*x，因此，只要解决了A*x的问题，算法就会继续。对于A*x=y，a的特殊性可以概括如下(如果栅格大小为n):当I=1时什么时候当I=N时什么时候什么时候什么时候什么时候什么时候什么时候这解决了无法使用矩阵a存储的尴尬情况，并大大减少了程序占用的内存大小。网格数n需要数百次时的情况可以在几秒钟内计算出来。四。实验环境(使用的软件、硬件等)和实验数据文件MATLAB2014aWindow7内存4G三.实验结

6、果和案例分析实验1源程序：Function myfun1(n)h=1/n；X=0:h:1A=zeros(n，n)；B=zeros(n，1)；生成For i=2:n%刚度矩阵A (I-1，I-1)=-(x(I-1)3 * x(I-1)2 * x(I)-3 * x(I-1)A(i-1，i-1)=A(1，1)；EndA (n，n)=-1/h 1/h 2 *(x(n 1)3/3-x(n)3/3-x(n)* x(n 1)2 xA(n，n)=A(1，1)/2；For i=2:nA (I-1，I)=1/h-1/h 2 *(x(I 1)3/6-x(I)3/6 x(I)2 * x(I 1)A(i，i-1)=A(

7、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 *。Endb(n)=1/(4 * h * pi 2)*(-2 * x(n)* cos(2 * pi * x(n 1)* pi 2 * cosY=myconjgrad(A，b)；调用求解%方程的算法U1=zeros(n 1，1)；u1(23360n 1)=y；Plot (x、u1、-b)；霍尔德温；X=0:0.01:1y=sin(2 * pi * x)/(1 4 * pi 2)；Plot(x、y、-r)；霍尔德温；Legend(有

8、限元，精确解)；图1: N=10数值解法与实际解决方案的比较图2: N=20数值解法与实际解决方案的比较实验2源程序：1.未处理系数矩阵Function my ellipse1 (n)h=1/(N 1)；k=1/(N 1)；l=k/h；X=0:h:1Y=0:k:1g=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:Na(I-1)* n 13360 I * n，(I-

9、1)* n 13360 I * n)=a；EndA=a diag (-ones (n * (n-1)，1)，n) diag (-ones (n * (n-1)，1)，-nB=zeros(N*N，1)；产生%右端常数向量b(n *(n-1)1:n * n)=g(23360n 1)；V=myconjgrad(A，B)；U=zeros(N，N)；u(:)=v；M=zeros(N 2，N 2)；M(:N 2)=g；M (23360n 1，2330n 1)=u；Surf (x (n 23360-1:1)，y，m)；霍尔德温；xlabel(x-axis)；xlabel(y-axis)；z label(Si

10、olution)；运行结果N=100可以实现直接解决的头发发行时如果N=120，则结果如下Elliptic1(120)错误地使用Diag内存不足。输入HELP MEMORY视图选项。错误ellipse1 (line14)A=a diag (-ones (n * (n-1)，1)，n) diag (-ones (n * (n-1)，1)，-n如出错的语句所示，矩阵顺序太大，内存不足，因此正在生成系数矩阵。2.增强的压缩存储矩阵主程序：Function my ellipsc2 (n)h=1/(N 1)；k=1/(N 1)；l=k/h；X=0:h:1Y=0:k:1g=sin(5 * pi * x)；

11、a1=(2 *(1 l 2)；B1=(-l 2)；C1=-1；%仅存储系数矩阵的三个元素B=zeros(N*N，1)；b(n *(n-1)1:n * n)=g(23360n 1)；V=myconjgrad1(N，a1，b1，c1，B)；调用共轭梯度法计算%线性方程的解U=zeros(N，N)；u(:)=v；M=zeros(N 2，N 2)；M(:N 2)=g；M (23360n 1，2330n 1)=u；%图片Surf (x (n 23360-1:1)，y，m)；霍尔德温；xlabel(x-axis)；xlabel(y-axis)；z label(Siolution)；子程序：求解线性方程的共轭梯度法Function x=myconjgrad 1 (n、a1、B1、C1、b)A2=a1B2=b1C2=c1N2=N；B=Btol=10(-4)；%设置错误限制x=b；R=b-fun(N2、a2、b2、c2、x)；调用在%保存后压缩A*b的算法p=r；for k=1:(numel(b)2；设置%周期数Z=myfun(N2、a2、b2、c2、p)；alpha=(r * r)/(p * z)；x=x alpha* p；S=r * rr=r-alpha * z；If(norm(r) tol)ReturnEndb=(r * r)/