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

Documentserialnumber【UU89WT-UU98YT-UU8CB-UUUT-UUT108】两点边值微分方程的元 分析 定性数据统计分析pdf销售业绩分析模板建筑结构震害分析销售进度分析表京东商城竞争战略分析 重庆大学学生实验 报告 软件系统测试报告下载sgs报告如何下载关于路面塌陷情况报告535n,sgs报告怎么下载竣工报告下载 实验课程名称偏微分方程数值解期末课程 设计 领导形象设计圆作业设计ao工艺污水处理厂设计附属工程施工组织设计清扫机器人结构设计 开课实验室偏微分方程数值解学院数统学院年级2011专业班学生姓名学号开课时间2013至2014学年第2学期总成绩教师签名数学与统计学院制开课学院、实验室：实验时间：2014年6月30日实验项目名称偏微分方程期末课程设计实验项目类型验证演示综合设计其他指导教师成绩实验目的自学，掌握有限元分析的基本理论，并运用有限元分析的 方法 快递客服问题件处理详细方法山木方法pdf计算方法pdf华与华方法下载八字理论方法下载 求解第二章的两点边值问题，做出数值解，体会有限元和差分方法的不同之处。掌握平面上拉普拉斯方程的五点差分方法，体会与一般一维问题的不同，特别注意边界条件的处理。学会处理大型方程组数值解的压缩存储方法。实验内容实验一：考虑两点边值问题：(1)边界条件为：真解：运用有限元的方法求解该方程的数值解，并和真解比较。实验二：用五点差分格式近似Laplacian方程：for正方形区域D和边间条件如下：要求：寻找一种数值算法，尽可能的让迭代步长变小即尽可能的让网格数N变大。三、实验原理、方法（算法）、步骤实验一：将方程化为标准形式：其中为常数第一步：考虑从Galerkin法出发建立有限元方程。考虑一阶Sobolev空间H1(a，b)的子空间V:V=v|v∈H1a，b，va=0任取v∈V,用它乘以方程(1)的两边,然后在区间(a,b)上积分对左端分部积分得利用边界条件va=0和，代入每项的系数得令则上式即方程(5)即为边值问题(1)的变分方程第二步：网格剖分对区间[0，1]作等距剖分，令步长h=1/n，n为网格数。则有xi=ih,i=0,1,2?n有限元方程的解u在每一个分段的小区间上为分片的线性函数uh∈Vh，设uh在网格点0=x0，x1，xi，,xn=1上分别取值为u0，u1，ui，,un取基函数：容易发现,基函数有如下性质显然，基函数线性无关，对任意的一个uh∈V可 表 关于同志近三年现实表现材料材料类招标技术评分表图表与交易pdf视力表打印pdf用图表说话 pdf 示为因此，求uh∈V，使得这相当于：求u0，u1，ui，,un，使得改写成矩阵形式其中运用基函数的性质可以化简矩阵：第三步：由上面的推导可知刚度矩阵K是三对角矩阵,故可采用“追赶法”求解有限元方程.求出各节点的近似解ui后,还可写出边值问题的解为：实验二;用中心差分代替Laplacian方程的二阶导数，舍去截断误差得到数值算法；Fori=1,,Nj=1,,MWhereλ=ckh边界条件为：(1)x=0:u0,j=0,forj=1,...,M,(2)x=1:uN+1,j=0,forj=1,...,M,(3)y=0:ui,0=0,fori=1,...,N,(4)y=1:ui,N+1=g(xi),fori=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实验结果及实例分析实验一源程序：functionmyfun1(n)h=1/n;x=0:h:1;A=zeros(n,n);b=zeros(n,1);fori=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*h^2);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+x(n)^2*x(n+1));A(n,n)=A(1,1)/2;fori=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)/2-x(i)*x(i+1)^2/2);A(i,i-1)=A(i-1,i);endfori=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))*cos(pi*x(i+1))-sin(2*pi*x(i)));endb(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)))*pi/(1+4*pi^2);y=myconjgrad(A,b);%调用求解方程组的算法u1=zeros(n+1,1);u1(2:n+1)=y;plot(x,u1,'--b');holdon;x=0::1;y=sin(2*pi*x)/(1+4*pi^2);plot(x,y,'-r');holdon;legend('有限元','精确解');图一：N=10数值解和真解的比较图二：N=20数值解和真解的比较实验二源程序：1．没有对系数矩阵进行处理的functionmyElliptic1(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+l^2));c=ones(N-1,1)*(-l^2);a=diag(b)+diag(c,1)+diag(c,-1);A=zeros(N^2,N^2);fori=1:NA((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);holdon;xlabel('x-axis');xlabel('y-axis');zlabel('Siolution');运行结果当N=100直接求解的发发还可行当N=120时，会出现以下结果Elliptic1(120)错误使用diag内存不足。请键入HELPMEMORY查看选项。出错Elliptic1(line14)A=A+diag(-ones(N*(N-1),1),N)+diag(-ones(N*(N-1),1),-N);由出错的语句可以看到，正是我们在生成系数矩阵，由于矩阵阶数太大，导致内存不足。2．改进的压缩存贮矩阵主程序：functionmyElliptic2(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+l^2));b1=(-l^2);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);holdon;xlabel('x-axis');xlabel('y-axis');zlabel('Siolution');子程序：求解线性方程组的的共轭梯度法functionx=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;fork=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)