《直接法求两点边值问题..doc》由会员分享，可在线阅读，更多相关《直接法求两点边值问题..doc（32页珍藏版）》请在三一办公上搜索。
1、 课程设计（论文）任务书数学与计算科学学院 学院 信息与计算科学 专业 班课程名称 科学仿真实验五 题 目 直接法求解两点边值问题（一） 任务起止日期： 2014 年 6 月 23 日 2014年 7月 6 日学 生 姓 名 学 号 指 导 教 师 教研室主任 年 月 日审查 课程设计（论文）任务一、课题内容 1. 查阅相关文献，弄清高斯消去法和矩阵的三角分解等问题；2. 编程实现并给出具体应用实例，撰写出设计论文；3. 通过对本课题的实践以期使学生程序编写、调试、用科学计算方法解决实际问题等能力得到较大提高。二、课题要求1. 使用有关算法语言完成本课题的程序设计；2. 程序必须得到合理结果，
2、并对所得结果做必要的分析；3. 设计论文正文篇幅不少于3000字；4. 提交的所有材料必须符合长沙理工大学课程设计管理规定（长理工大教20058号）的要求。纸制文档资料包括：题目的复述、问题的分析、算法的描述、原程序、测试数据及有关运行结果和有关说明。三、课题完成后应提交的材料1. 课程设计材料按以下排列顺序装订成册(1) 封面（统一到学校教材中心领取，并详细填写） (2) 任务书(3) 题目内容(4) 问题分析、算法描述等(5) 简短源程序及有关运行结果等(6) 参考文献(7) 附件（复杂的源程序打印件、有关运行结果）2.装订成册的材料装入资料袋 资料袋统一到学校教材中心领取，并详细填写四、
3、主要参考文献（由指导教师选定）1 严蔚敏 吴伟民数据结构(C语言版)M.北京：清华大学出版社，1997.42 李庆扬，王能超，易大义.数值分析M.武汉.华中科技大学出版社，2006.7.3 清华大学、北京大学计算方法编写组。计算方法M。北京。科学出版社，1980注：1. 此任务书由指导教师填写。如不够填写，可另加页。2. 此任务书最迟必须在课程设计（论文）开始前下达给学生。学生送交全部材料日期 学生（签名） 指导教师验收（签名） 直接法求解两点边值问题（二）摘要线性方程组的数值解法可以分为直接法和迭代法两类。所谓直接法，就是不考虑舍入误差，通过有限步骤四则运算即能求得线性方程组准确解的方法。如
4、克莱姆法则，但通过第一章的分析，我们知道用克莱姆法则来求解线性代数方程组并不实用，因而寻求线性方程组的快速而有效的解法是十分重要的。本章讨论计算机上常用而有效的直接解法高斯消去法和矩阵的三角分解等问题。为方便计，设所讨论的线性方程组的系数行列式不等于零。高斯消去法是解线性方程组最常用的方法之一，它的基本思想是通过逐步消元，把方程组化为系数矩阵为三角形矩阵的同解方程组，然后用回代法解此三角形方程组得原方程组的解。关键词：线性方程组；直接解法；高斯消去法DIRECT METHOD SOLVING TWO-POINT BOUNDARY VALUE PROBLEMS (2)ABSTRACTNumeri
5、cal algorithm of linear equations can be divided into two categories, direct method and iterative method. The so-called direct method, is not considered rounding error, through limited steps arithmetic which can obtain the accurate solution of linear equations method. Such as cramers rule, but throu
6、gh the analysis of the first chapter, we know that cramers rule is used to solve the linear algebraic equations is not practical, thus seeking quick and effective solutions of systems of linear equations solution is very important.This chapter discuss computer commonly used and effective direct solu
7、tion - gaussian elimination and triangle decomposition of matrices. For the convenience of meter, discussed the coefficient determinant of linear equations is not equal to zero.Gauss elimination method is one of the most commonly used method of solving linear equations, the basic idea is to pass a g
8、radual elimination, to coefficient matrix of the triangular matrix equations with solutions of the equations, then by back substitution method solving the triangle equations to the solution of the original equations.Key words: linear equations; Direct method; Gaussian elimination目录1 问题的提出12 理论基础1 2.
9、1 高斯消去法22.2 列主元消去法52.3 矩阵的三角分解法6 2.3.1 算法介绍62.3.2 定理结论72.3.3 计算公式92.4 解三对角方程组的追赶法103 问题的求解123.1 顺序消去法123.2 列主元消去法133.3 Doolittle分解法143.4 追赶法154 计算结果16参考文献20附录211 问题的提出 考虑两点边值问题：容易知道它的精确解为： 为了把微分方程离散，把区间等分，令，得到差分方程：简化为：从而离散后得到的线性方程组的系数矩阵为：对，分别用顺序消去法、列主元消去法、Doolittle分解法和追赶法求解线性方程组，然后比较与精确解的误差，对结果进行分析。
10、改变，讨论同样问题。 2 理论基础许多科学技术问题要归结为解含有多个未知量x1, x2, , xn的线性方程组： （2.1）这里aij (i, j = 1, 2, , n)为方程组的系数，bi（i = 1, 2, , n）为方程组自由项。方程组（2.1）的矩阵形式为：AX = b其中： 线性方程组的数值解法可以分为直接法和迭代法两类。所谓直接法，就是不考虑舍入误差，通过有限步骤四则运算即能求得线性方程组（2.1）准确解的方法。如克莱姆法则，但通过第一章的分析，我们知道用克莱姆法则来求解线性代数方程组并不实用，因而寻求线性方程组的快速而有效的解法是十分重要的。 本章讨论计算机上常用而有效的直接解
11、法高斯消去法和矩阵的三角分解等问题。为方便计，设所讨论的线性方程组的系数行列式不等于零。高斯（Gauss）消去法是解线性方程组最常用的方法之一，它的基本思想是通过逐步消元，把方程组化为系数矩阵为三角形矩阵的同解方程组，然后用回代法解此三角形方程组得原方程组的解。 2.1 高斯消去法一般形式的线性方程组的解法中，为叙述问题方便，将bi写成ai, n+1，i = 1, 2,n。 （2.1）如果a11 0，将第一个方程中x1的系数化为1，得： 其中： j = 1, , n + 1（记 i = 1, 2, , n; j = 1, 2, , n + 1） 从其它n 1个方程中消x1，使它变成如下形式 （
12、2.2）其中： 由方程（2.2）到（2.3）的过程中，元素起着重要的作用，特别地，把称为主元素。 如果（2.3）中，则以为主元素，又可以把方程组（2.3）化为： （2.3） 针对（2.4） 继续消元，重复同样的手段，第k步所要加工的方程组是： 设，第k步先使上述方程组中第k个方程中xk的系数化为1：然后再从其它（n - k）个方程中消xk，消元公式为： （2.4） 按照上述步骤进行n次后，将原方程组加工成下列形式： 回代公式为： （2.5） 综上所述，高斯消去法分为消元过程与回代过程，消元过程将所给方程组加工成上三角形方程组，再经回代过程求解。 由于计算时不涉及xi, i = 1, 2, ,
13、n，所以在存贮时可将方程组AX = b，写成增广矩阵（A, b）存贮。 下面，我们统计一下高斯消去法的工作量；在（2.5）第一个式子中，每执行一次需要次除法，在（2.5）第二个式子中，每执行一次需要次除法。因此在消元过程中，共需要：次乘作法。此外，回代过程共有：次乘法。汇总在一起，高斯消去法的计算量为：次乘除法。2.2 列主元消去法在列主元消去法中，未知数仍然是顺序地消去的，但是把各方程中要消去的那个未知数的系数按绝对值最大值作为主元素，然后用顺序消去法的公式求解。 由于解方程组取决于它的系数，因此可用这些系数（包括右端项）所构成的“增广矩阵”作为方程组的一种简化形式。 列主元消去法计算步骤：
14、1、输入矩阵阶数 ，增广矩阵 2、对 （1）按列选主元：选取 使 （2）如果 ，交换 的 第k行与底l 行元素（3）消元计算 ：（4）回代计算： （5） 输出解向量： 。2.3 矩阵的三角分解法下面我们进一步用矩阵理论来分析高斯消去法，从而建立矩阵的三角分解定理，而这个定理在解方程组的直接解法中起着重要作用，我们将利用它来推导某些具有特殊的系数矩阵方程组的数值解法。2.3.1 算法介绍考虑线性方程组：AX = b ，ARnn 设解此方程组用高期消去法能够完成（不进行变换两行的初等变换），由于对A施行初等变换相当于用初等矩阵左乘A，于是，高斯消去法的求解过程用矩阵理论来叙述如下：记：其中 记 于
15、是： 如此继续下去，到n 1步时有：也就是说： 记：则有：A = LU。其中L为单位下三角矩阵，U为上三角矩阵。2.3.2 定理结论 总结上述讨论得到重要性定理：定理1： （矩阵的三角分解）设A为n n实矩阵，如果解AX = b用高斯消去法能够完成（限制不进行行的交换，即），则矩阵A可分解为单位下三角矩阵L与上三角知阵U的乘积。A = LU且这种分解是唯一的。证明：由上述的讨论，存在性已得证，现在证唯一性。设 A = L1 U1 = LU其中L1，L1为单位下三角阵，U1，U为上三角阵，设存在，于是有上式右端为上三角矩阵，左边为单位下三角阵，故应为单位矩阵。即 L1 = L U1 = U 。
16、定理2： 约化主元素（i = 1, 2, , k）的充要条件是矩阵A的顺序主子式：由上述讨论知，解 AX = b的高斯消去法相当于实现了A的三角分解，如果我们能直接从矩阵A的元素得到计算L，U的元素的公式，实现A的三角分解，而不需要任何中间步骤，那么求解AX = b的问题就等价于求解两个三角形矩阵方程组：（1）Ly = b 求y（2）UX = y 求x下面来说明L、U的元素可以由A的元素直接计算确定。显然，由矩阵乘法。得到U的第一行元素；由得：即L的第一列元素。设已经求出U的第1行第r 1行元素，L的第1列第r 1列元素，由矩阵乘法可得：即可计算出U的第r行元素，L的第r列元素。2.3.3 计
17、算公式综上所述，可得到用直接三角分解法解AX = b的计算公式。（1） （3.1） （3.2）对于r = 2, 3, , n计算（2）计算U的第r行元素 （3.3）（3）计算L的第r列元素 (r n) （3.4）（4） （3.5）（5）（3.6）（1）、（2）、（3）是矩阵A的LU分解公式,称为Doolittle分解。同理，可推出矩阵A = LU分解的另一种计算公式，其中L为下三角，U为单位上三角，这种矩阵的分解公式称为矩阵的Crout分解。2.4 解三对角方程组的追赶法 在实际问题中，经常遇到以下形式的方程组： （4.1）这种方程组的系数矩阵称为三对角矩阵。以下针对这种方程组的特点提供一种简
18、便有效的算法追赶法。追赶法实际上是高斯消去法的一种简化形式，它同样分消元与回代两个过程。先将（4.1）第一个方程中x1的系数化为1：记： （4.2）有： 注意到剩下的方程中，实际上只有第二个方程中含有变量x1，因此消元手续可以简化。利用（4.2）可将第二个方程化为：这样一步一步地顺序加工（4.1）的每个方程，设第k 1个方程已经变成： （4.3）再利用（4.3）从第k个方程中消去xk-1，得：同除，得：记：则有： 这样做n 1步以后，便得到：将上式与（4.1）中第11个方程联立，即可解出：xn = yn这里：于是，通过消元过程，所给方程组（4.1）可归结为以下更为简单的形式： （4.4） 这种
19、方程组称作二对角型方程组，其系数矩阵中的非零元素集中分步在主对角线和一条次主对角线上： 对加工得到的方程组（3.4）自下而上逐步回代，即可依次求出xn，xn-1，x1，计算公式为： （4.5） 上述算法就是追赶法，它的消元过程与回代过程分别称作“追”过程与“赶”过程。综合追与赶的过程，得如下计算公式： （4.6）（4.7）3 问题的求解3.1 顺序消去法function x = gaussElim(A,b)N = length(b);if A(1,1)=0 n_temp=find(A(:,1)=0); c=A(n_temp(1),:); A(n_temp(1),:)=A(1,:); A(1,:
20、)=c;endfor j=2:Nm = A(j:N,j-1)/A(j-1,j-1) A(j:N,:) = A(j:N,:) - m*A(j-1,:); b(j:N) = b(j:N) - m*b(j-1); if A(j,j)=0 &j=N n_temp=find(A(j:end,j)=0); c=A(j-1+n_temp(1),:); A(j-1+n_temp(1),:)=A(j,:); A(j,:)=c; end endif A(N,N)=0 error(the coefficient matrix of the linear equations is sigular);endx = ze
21、ros(N,1);x(N) = b(N)/A(N,N); for j=N-1:-1:1, x(j) = (b(j)-A(j,j+1:N)*x(j+1:N)/A(j,j); end3.2 列主元消去法function x=f(A)%列主元高斯消去法m,n=size(A);B=zeros(size(A(1,:);for i=1:(m-1) for j=i:m if A(j,i)A(i,i) B=A(j,:);A(j,:)=A(i,:);A(i,:)=B; %选取主元 end end for j=(i+1):m A(j,:)=A(j,:)-A(i,:)*A(j,i)/A(i,i); % 消元法 en
22、d end%回代过程 x(m)=A(m,n)/A(m,m); for i=(m-1):-1:1 x(i)=(A(i,n)-A(i,i+1:m)*x(i+1:m)/A(i,i); end x3.3 Doolittle分解法functionx,l,u=malu(A,b)%用途：用LU分解法解方程组%格式：x,l,u=malu(A,b),A为系数矩阵，b为右端向量，x返回解向量，l返回下三角矩阵，u返回上三角矩阵format short %LU分解n=length(b);u=zeros(n,n);l=eye(n,n);u(1,:)=A(1,:);l(2:n,1)/u(1,1);for k=2:nu(
23、k,k:n)=A(k,k:n)-l(k,1:k-1)*u(1:k-1,k:n);l(k+1:n,k)=(A(k+1:n,k)-l(k+1:n,k-1)*u(l:k-1,k)/u(k,k);End %解下三角方程组Ly=by=zeros(n,1);X(n)=y(n)/u(n,n);For k=n-1:-1:1x(k)=(y(k)-u(k,k+1:n)*x(k+1：n)/u(k,k);end3.4追赶法function x,L,U=thomas(e,p,n)a=e.*ones(1,n-2);b=-(2*e+1/n).*ones(1,n-1);c=(e+1/n).*ones(1,n-2);f=(p*
24、(1/n)*(1/n).*ones(1,n-1);n=length(b);%对A进行分解a,b,c,fu(1)=b(1);for i=2:n if u(i-1)=0 l(i-1)=a(i-1)/u(i-1); u(i)=b(i)-l(i-1)*c(i-1); else break; endendL=eye(n)+diag(l,-1);U=diag(u)+diag(c,1);x=zeros(n,1);y=x;%求解Ly=by(1)=f(1);for i=2:n y(i)=f(i)-l(i-1)*y(i-1);end%求解Ux=yif u(n)=0 x(n)=y(n)/u(n);End4 计算结果
25、（1）当e=1，a=0.5，n=100时，以上方法求解线性方程组解出的结果完全相同，但与真正的结果相差很大。平均相对误差S1，S2更是高达1.1371.因为：与，相差很大。（2）当e=1，a=0.5，n=10时，以上方法求解线性方程组解出的结果完全相同，与真正的结果相差无几。平均相对误差S1，S2是0.0040.因为： 与，不在一个数量级 。（3）当e=1，a=0.5，n=200时，以上方法求解线性方程组解出的结果完全相同，与真正的结果相差无几。平均相对误差S1，S2是 0.0021。（4）当e=1，a=0.5，n=300时，以上方法求解线性方程组解出的结果完全相同，与真正的结果相差无几。平均
26、相对误差S1，S2是 0.0016。参考文献1 严蔚敏、吴伟民数据结构(C语言版)M.北京：清华大学出版社，1997.4.2 李庆扬、王能超、易大义.数值分析M.武汉.华中科技大学出版社，2006.7.3 清华大学、北京大学计算方法编写组.计算方法M.北京.科学出版社，1980.附录（1）命令文件e=1;a=0.5;n=100;x=(0+1/n):1/n:(1-1/n);y=(1-a)*(1-exp(-x./e)/(1-exp(-1/e)+a.*x; %正确结果A=finput(e,a,n); %输入矩阵A y1=f(A); % 列主元消去法结果y2,L,U=thomas(e,a,n);%追赶
27、法结果plot(x,y,x,y1,*,x,y2,+)s1=0;s2=0;for i=1:n-2 s1=s1+(y1(i)-y(i)/ y(i); s2=s2+(y2(i)-y(i)/ y(i);ends1= s1/(n-2)s2=s2/(n-2)（2） 输入矩阵Afunction A=finput(e,a,n)%输入矩阵A=zeros(n-1,n);h=1/n;for i=1:n-1 for j=1:nif j=i A(i,j)=-(2*e+h); end if j=(i-1) A(i,j)=e; end if j=(i+1) A(i,j)=e+h; end if j=n A(i,j)=a*h
28、*h; end endEnd（3） 顺序消去法function x = gaussElim(A,b)N = length(b);if A(1,1)=0 n_temp=find(A(:,1)=0); c=A(n_temp(1),:); A(n_temp(1),:)=A(1,:); A(1,:)=c;endfor j=2:Nm = A(j:N,j-1)/A(j-1,j-1) A(j:N,:) = A(j:N,:) - m*A(j-1,:); b(j:N) = b(j:N) - m*b(j-1); if A(j,j)=0 &j=N n_temp=find(A(j:end,j)=0); c=A(j-1
29、+n_temp(1),:); A(j-1+n_temp(1),:)=A(j,:); A(j,:)=c; end endif A(N,N)=0 error(the coefficient matrix of the linear equations is sigular);endx = zeros(N,1);x(N) = b(N)/A(N,N); for j=N-1:-1:1, x(j) = (b(j)-A(j,j+1:N)*x(j+1:N)/A(j,j); end（4） 列主元高斯消去法function x=f(A)%列主元高斯消去法m,n=size(A);B=zeros(size(A(1,:
30、);for i=1:(m-1) for j=i:m if A(j,i)A(i,i) B=A(j,:);A(j,:)=A(i,:);A(i,:)=B;%选取主元 end end for j=(i+1):m A(j,:)=A(j,:)-A(i,:)*A(j,i)/A(i,i); % 消元法 end end%回代过程x(m)=A(m,n)/A(m,m); for i=(m-1):-1:1 x(i)=(A(i,n)-A(i,i+1:m)*x(i+1:m)/A(i,i); end X（5）Doolittle分解法functionx,l,u=malu(A,b)%用途：用LU分解法解方程组%格式：x,l,u
31、=malu(A,b),A为系数矩阵，b为右端向量，x返回解向量，l返回下三角矩阵，u返回上三角矩阵format short %LU分解n=length(b);u=zeros(n,n);l=eye(n,n);u(1,:)=A(1,:);l(2:n,1)/u(1,1);for k=2:nu(k,k:n)=A(k,k:n)-l(k,1:k-1)*u(1:k-1,k:n);l(k+1:n,k)=(A(k+1:n,k)-l(k+1:n,k-1)*u(l:k-1,k)/u(k,k);End %解下三角方程组Ly=by=zeros(n,1);X(n)=y(n)/u(n,n);For k=n-1:-1:1x(
32、k)=(y(k)-u(k,k+1:n)*x(k+1：n)/u(k,k);End（6）追赶法求解function x,L,U=thomas(e,p,n)a=e.*ones(1,n-2);b=-(2*e+1/n).*ones(1,n-1);c=(e+1/n).*ones(1,n-2);f=(p*(1/n)*(1/n).*ones(1,n-1);n=length(b);%对A进行分解a,b,c,fu(1)=b(1);for i=2:n if u(i-1)=0 l(i-1)=a(i-1)/u(i-1); u(i)=b(i)-l(i-1)*c(i-1); else break; end end L=eye(n)+diag(l,-1); U=diag(u)+diag(c,1); x=zeros(n,1); y=x; %求解Ly=b y(1)=f(1); for i=2:n y(i)=f(i)-l(i-1)*y(i-1); end %求解Ux=y if u(n)=0 x(n)=y(n)/u(n);end
