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

重 庆 大 学学 生 实 验 报 告实验课程名称 偏微分方程数值解期末课程设计 开课实验室 偏微分方程数值解 学 院 数统学院 年级 2011 专业班 学 生 姓 名 学 号 开 课 时 间 2013 至 2014 学年第 2 学期总 成 绩教师签名数学与统计学院制开课学院、实验室： 实验时间 ： 2014 年 6 月 30 日实验项目名 称偏微分方程期末课程设计实验项目类型验证演示综合设计其他指导教师成 绩一 实验目的自学，掌握有限元分析的基本理论，并运用有限元分析的方法求解第二章的两点边值问题，做出数值解，体会有限元和差分方法的不同之处。掌握平面上拉普拉斯方程的五点差分方法，体会与一般一维问题的不同，特别注意边界条件的处理。学会处理大型方程组数值解的压缩存储方法。二 实验内容实验一：考虑两点边值问题： (1) 边界条件为： 真解： 运用有限元的方法求解该方程的数值解，并和真解比较。实验二： 用五点差分格式近似Laplacian方程： for 正方形区域D和边间条件如下： 要求： 寻找一种数值算法，尽可能的让迭代步长变小即尽可能的让网格数N变大。三、实验原理、方法（算法）、步骤实验一：将方程化为标准形式：其中 为常数第一步：考虑从Galerkin法出发建立有限元方程。考虑一阶Sobolev 空间H1(a，b) 的子空间V: V=v|vH1a，b，va=0任取vV, 用它乘以方程 (1) 的两边, 然后在区间(a, b) 上积分 对左端分部积分得 利用边界条件va=0 和 ，代入每项的系数得 令 则上式即 方程(5) 即为边值问题(1)的变分方程第二步：网格剖分 对区间0，1 作等距剖分，令步长 h=1/n ， n为网格数。则有 xi=ih , i=0,1,2n有限元方程的解u在每一个分段的小区间上为分片的线性函数uhVh，设uh在网格点 0=x0，x1，xi，,xn=1 上分别取值为u0，u1，ui，,un取基函数：容易发现, 基函数有如下性质显然，基函数线性无关，对任意的一个uhV可表示为 因此，求uhV，使得 这相当于：求u0，u1，ui，,un，使得 改写成矩阵形式 其中 运用基函数的性质可以化简矩阵： 第三步： 由上面的推导可知刚度矩阵K 是三对角矩阵, 故可采用“追赶法”求解有限元方程. 求出各节点的近似解ui 后, 还可写出边值问题的解为： 实验二; 用中心差分代替Laplacian方程的二阶导数，舍去截断误差得到数值算法；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随着网格数的增加急剧的变大以至于电脑内存存不下，因此有必要对系数矩阵进行特别处理。观察矩阵网格数为三的系数阵可以发现整个矩阵只由三个不同的元素构成，所以我们不去存储完整的矩阵，而是只保存这个三个不同的元素，即可得全部的信息。令 在共轭梯度法的迭代计算中需要不停地计算A*x,因此只要解决了A*x的问题，则算法即可进行。对于A*x=y,由A的特殊性我们可以归纳计算如下（当网格大小为N时）：当 i=1 当 当i=N 当 当 当 当 当 当 因此这样就解决了矩阵A存储不下的尴尬情况，大大的减少了程序所占内存的大小。几乎几秒钟就能算出网格数N取好几百时的情况。四实验环境（所用软件、硬件等）及实验数据文件 MATLAB2014aWindow7 内存4G三 实验结果及实例分析实验一源程序：function myfun1(n)h=1/n;x=0:h:1;A=zeros(n,n);b=zeros(n,1);for i=2:n %生成刚性矩阵 A(i-1,i-1)=-(x(i-1)3+3*x(i-1)2*x(i)-3*x(i-1)*x(i)2-2*x(i)3+3*x(i)2*x(i+1)-3*x(i)*x(i+1)2+x(i+1)3-3*x(i-1)+3*x(i+1)/(3*h2); A(i-1,i-1)=A(1,1);endA(n,n)=-1/h+1/h2*(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/h2*(x(i+1)3/6-x(i)3/6+x(i)2*x(i+1)/2-x(i)*x(i+1)2/2); A(i,i-1)= A(i-1,i);endfor i=2:n %生成右端常数项 b(i-1)=1/(2*h*pi2)*(-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)*cos(pi*x(i+1)-sin(2*pi*x(i); endb(n)=1/(4*h*pi2)*(-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*pi2);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*pi2); plot(x,y,-r);hold on;legend(有限元,精确解);图一：N=10数值解和真解的比较图二：N=20数值解和真解的比较 实验二源程序：1没有对系数矩阵进行处理的function myElliptic1(N)h=1/(N+1);k=1/(N+1);l=k/h;x=0:h:1;y=0:k:1;g=sin(5*pi*x); a=zeros(N,N); %生成系数矩阵b=ones(N,1)*(2*(1+l2);c=ones(N-1,1)*(-l2);a=diag(b)+diag(c,1)+diag(c,-1);A=zeros(N2,N2);for i=1:N A(i-1)*N+1:i*N,(i-1)*N+1:i*N)=a;endA=A+diag(-ones(N*(N-1),1),N)+diag(-ones(N*(N-1),1),-N);B=zeros(N*N,1); %生成右端常向量B(N*(N-1)+1:N*N)=g(2:N+1);v=myconjgrad(A,B);u=zeros(N,N);u(:)=v;m=zeros(N+2,N+2);m(:,N+2)=g;m(2:N+1,2:N+1)=u;surf(x(N+2:-1:1),y,m); hold on;xlabel(x-axis);xlabel(y-axis);zlabel(Siolution);运行结果 当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); 由出错的语句可以看到，正是我们在生成系数矩阵，由于矩阵阶数太大，导致内存不足。2改进的压缩存贮矩阵 主程序：function myElliptic2(N)h=1/(N+1);k=1/(N+1);l=k/h;x=0:h:1;y=0:k:1;g=sin(5*pi*x); a1=(2*(1+l2);b1=(-l2);c1=-1; %只存储系数矩阵中的三个元素B=zeros(N*N,1);B(N*(N-1)+1:N*N)=g(2:N+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(2:N+1,2:N+1)=u; %作图surf(x(N+2:-1:1),y,m); hold on;xlabel(x-axis);xlabel(y-axis);zlabel(Siolution); 子程序：求解线性方程组的的共轭梯度法function x = myconjgrad1(N,a1,b1,c1,B) a2=a1;b2=b1;c2=c1;N2=N;b=B; tol=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*r; r = r - alpha*z; if( norm(r) tol ) return; end B = (r*r)/s; p = r + B*p; end 子程序：压缩存储矩阵的乘法function v=myfun(N,a,b,c,x)v=zeros(N*N,1);v(1)=a*x(1)+b*x(2)+c*x(N+1);for i=2:N-1; v(i)=b*x(i-1)+a*x(i)+b*x(i+1)+c*x(i+N);endv(N)=b*x(N-1)+a*x(N)+c*x(N+N);v(N+1)=c*x(1)+a*x(N+1)+b*x(N+2)+c*x(N+1+N);for i=N+2:N*N-N-1; v(i)=c*x(i-N)+b*x(i-1)+a*x(i)+b*x(i+1)+c*x(i+N);endv(N*N-N)=c*x(N*N-2*N)+b*x(N*N-N-1)+a*x(N*N-N)+c*x(N*N);v(N*N-N+1)=c*x(N*N-2*N+1)+a*x(