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

重 庆 大 学学 生 实 验 报 告实验课程名称 偏微分方程数值解期末课程设计 开课实验室 偏微分方程数值解 学 院 数统学院 年级 2011 专业班 学 生 姓 名 学 号 开 课 时 间 2013 至 2014 学年第 2 学期总 成 绩教师签名数学与统计学院制开课学院、实验室： 实验时间 ： 2014 年 6 月 30 日实验项目名 称偏微分方程期末课程设计实验项目类型验证演示综合设计其他指导教师成 绩实验目的自学，掌握有限元分析的基本理论，并运用有限元分析的方法求解第二章的两点边值问题，做出数值解，体会有限元和差分方法的不同之处。掌握平面上拉普拉斯方程的五点差分方法，体会与一般一维问题的不同，特别注意边界条件的处理。学会处理大型方程组数值解的压缩存储方法。实验内容实验一：考虑两点边值问题： (1) 边界条件为： 真解： 运用有限元的方法求解该方程的数值解，并和真解比较。实验二： 用五点差分格式近似Laplacian方程： for 正方形区域D和边间条件如下： 要求： 寻找一种数值算法，尽可能的让迭代步长变小即尽可能的让网格数N变大。三、实验原理、方法（算法）、步骤实验一：将方程化为标准形式：其中 为常数第一步：考虑从Galerkin法出发建立有限元方程。考虑一阶Sobolev 空间的子空间V: 任取, 用它乘以方程 (1) 的两边, 然后在区间(a, b) 上积分 对左端分部积分得 利用边界条件和 ，代入每项的系数得 令 则上式即 方程(5) 即为边值问题(1)的变分方程第二步：网格剖分 对区间[0，1] 作等距剖分，令步长 h=1/n ， n为网格数。则有 有限元方程的解u在每一个分段的小区间上为分片的线性函数，设在网格点 上分别取值为取基函数：容易发现, 基函数有如下性质显然，基函数线性无关，对任意的一个可表示为 因此，求，使得 这相当于：求，使得 改写成矩阵形式 其中 运用基函数的性质可以化简矩阵： 第三步： 由上面的推导可知刚度矩阵K 是三对角矩阵, 故可采用“追赶法”求解有限元方程. 求出各节点的近似解ui 后, 还可写出边值问题的解为： 实验二; 用中心差分代替Laplacian方程的二阶导数，舍去截断误差得到数值算法；For Where 边界条件为：(1) x = 0: ,for j = 1,..., M,(2) x = 1: ,for j = 1,..., M,(3) y = 0: ,for i = 1,..., N,(4) y = 1: ,for i = 1,..., N.可以改写为矩阵形式为：其中对于对称正定大型稀疏矩阵问题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*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;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); A(i,i-1)= A(i-1,i);endfor i=2:n%生成右