求解大型稀疏线性方程组的选主元高斯-约当法

1、沈阳航空航天大学数值分析课程设计（设计程序）题目 求解大型稀疏线性方程组的选主元高斯约当法 班级 / 学号 学 生 姓 名 指 导 教 师 刘勇进 沈阳航空航天大学课 程 设 计 任 务 书课 程 名 称 数值分析课程设计 院（系） 理学院 专业 信息与计算科学 班级 24140101 学号 1 姓名 高莹 课程设计题目 求解大型稀疏线性方程组的选主元高斯约当法 课程设计时间: 2013 年 12月 16日至 2013 年 12月 30日课程设计的要求及内容：要求1. 学习态度要认真，要积极参与课程设计，锻炼独立思考能力；2. 严格遵守上机时间安排；3. 按照MATLAB编程训练的任务要求来编

2、写程序；4. 根据任务书来完成课程设计论文；5. 报告书写格式要求按照沈阳航空航天大学“课程设计报告撰写规范”；6. 报告上交时间：课程设计结束时上交报告；7. 严禁抄袭行为，一旦发现，课程设计成绩为不及格。一、 通识题1. 作出下列函数的图形：(1) .要求：a/ 在同一图形窗口画出图形；b/ 将图形窗口分成个小图形区域。(2) .(3) 及 要求：a/ 分别考虑；b/ .2. 用数值方法计算， 其中n分别取100, 200, 1000。3. 求由方程所确定的隐函数的导数。4. 求由参数方程确定的函数的二阶导数。5. 作出函数的曲线在区间上的两条水平切线。二、 求解大型稀疏线性方程组的选主元

3、高斯约当法 要求：1）掌握求解大型稀疏线性方程组的选主元高斯约当法的基本原理； 2）编写求解大型稀疏线性方程组的选主元高斯约当法的通用程序； 3）利用所编写的程序求解方程组，其中 ，指导教师 年 月 日负责教师 年 月 日学生签字 年 月 日沈阳航空航天大学课 程 设 计 成 绩 评 定 单课 程 名 称 数值分析课程设计 院（系） 理学院 专业 信息与计算科学 课程设计题目 求解大型稀疏线性方程组的选主元高斯约当法 学号 1 姓名 高莹 指导教师评语：课程设计成绩 指导教师签字 年 月 日目 录一、正文1二、总结13参考文献14附 录15一、 正文一·通识题1、作出下列函数的图形：

4、（1）.要求：a/ 在同一图形窗口画出图形；b/ 将图形窗口分成个小图形区域。题目分析：这道题是利用matlab绘图画出的图像，其中n的取值为2,3,4,5，采用plot函数绘制，x取值范围从0到1，所以选取适当节点，这样可以尽量减少误差。其次，指数函数要用点乘运算，即 y=x.n。程序运行结果如图1-1：图1-1源程序：x = 0:0.0001:1;n=2;i=1;for n=2:1:5 y = x.n; subplot(4,1,i); i=i+1; plot(x,y); title('y=xn'); xlabel('x');ylabel('y'

5、;);end(2) .题目分析：利用matlab绘图画出的图像，采用plot函数绘制，x取值为-1到1，所以选取适当节点，可以使画出来的图像更接近真实图像，同样，数与x之间需要用点乘。程序运行结果如图1-2：图1-2源程序：x = -1:0.000001:1; y = 1./(1+25.*x.2); plot(x,y); title('y=1/(1+25*x2),x from -1 to 1'); xlabel('x'); ylabel('y');（3）题目分析：利用matlab绘图画出y=arctan(x)的图像，采用plot函数绘制，x的范围

6、为-50到50，选取适当结点，可以使图像更加精准。根据题目要求，编写出程序。程序运行结果如图1-3-1：图1-3.1源程序：x = -50:0.0001:50;y = atan(x);plot(x,y);title('y=atan(x)');text(1,pi/4,'y=atan(x)');xlabel('x');ylabel('y');（4）题目分析：matlab绘图画出的图像，n的取值为5,10,20,50，采用plot函数绘制，选取适当节点，保证图像的连续性，减小误差。并利用matlab绘图图形窗口的分割，画出四个图比较n值

7、不同时，函数的图像的走势变化。程序运行结果如图1-3.2：图1-3.2源程序：x=-50:0.0001:50;subplot(2,2,1)for n=5 p1=p1+(-1)(n-1)*(x.(2*n-1)/(2*n-1);endplot(x,p1);title('n=5');xlabel('x');ylabel('p(x)');subplot(2,2,2)for n=10 p2=p2+(-1)(n-1)*(x.(2*n-1)/(2*n-1);endplot(x,p2);title('n=10');xlabel('x

8、9;);ylabel('p(x)');subplot(2,2,3)for n=20 p3=p3+(-1)(n-1)*(x.(2*n-1)/(2*n-1);endplot(x,p3);title('n=20');xlabel('x');ylabel('p(x)');subplot(2,2,4)for n=50 p4=p4+(-1)(n-1)*(x.(2*n-1)/(2*n-1);endplot(x,p4);title('n=50');xlabel('x');ylabel('p(x)')

9、;2、用数值方法计算， 其中n分别取100, 200, 1000。题目分析：利用matlab计算级数的和，其中k的范围从1到n，编写程序，输入n值，显示结果。n分别取100,200,1000的值。程序运行如图2-1：图2-1源程序：function s=wo(n)sum=0;k=1;while(k<=n) sum=sum+1/k.2; k=k+1;ends=sum;3、求由方程所确定的隐函数的导数。题目分析：利用matlab求方程所确定的隐函数的导数，即对x,y分别求偏导数，利用diff函数，最后根据公式，求出。程序运行结果如图3-1：图3-1源程序：syms x y;f=log(x)+

10、exp(-y/x);diff(f,x);diff(f,y);ans=-diff(f,x)/diff(f,y);4、求由参数方程确定的函数的二阶导数。题目分析：利用matlab编写程序求参数方程确定的函数的二阶导数，用matlab中求导数diff函数。并运用公式，求出结果。程序运行结果如图4-1：图4-1源程序：syms x y t ;x=t-log(1+t);y=t3+t2;ans=(diff(x)*diff(y,2)-diff(x,2)*diff(y)/(diff(x)3;5、作出函数的曲线在区间上的两条水平切线。题目分析： 利用matlab绘图画出函数的曲线在区间上的两条水平切线。可以通过

11、求f(x)的一阶导数为零时x的值。通过x的值求出y值，在建立新的切线方程，即可画出曲线的水平切线。程序运行结果如图5-1：图5-1源程序：x=0:0.0001:1;y=4*x.3-6*x.2+x-2;f=inline('12*x.2-12*x+1');y1=4*fzero(f,0)3-6*fzero(f,0)2+fzero(f,0)-2;y2=4*fzero(f,1)3-6*fzero(f,1)2+fzero(f,1)-2;plot(x,y,x,y1,x,y2);title('4x3-6x2+x-2');xlabel('x');ylabel(&#

12、39;Y');二、求解大型稀疏线性方程组的选主元高斯约当法 要求：1）掌握求解大型稀疏线性方程组的选主元高斯约当法的基本原理； 2）编写求解大型稀疏线性方程组的选主元高斯约当法的通用程序； 3）利用所编写的程序求解方程组，其中 ，1）选主元高斯-约当法基本原理： 消元计算Ak的第k行对角线上方、下方的元素都进行消元计算，最后不需要回带即可求出方程组的解，这就是高斯-约当法，在进行选主元，就是选主元高斯-约当法。选主元高斯-约当消去法的步骤：设已经完成第1步到第k-1步，得到矩阵。计算过程如下：设 0；1. 选主元，即确定使；2. 如果，则方程组解不唯一，或者A接近奇异矩阵，停止计算；3

13、. 如果，则交换第行与第k行元素；4. 消元计算：当i=k时， j=k+1,k+1,.,n当时， j=k+1,k+2,.,n;i=1,2,.,n,上述过程完成后，计算解为。2）通用程序如下： function x,flag = Gau_Jor(A,B)n,m = size(A);nb=length(B);if n=m error('The rows and columns of matrix A must be equal!'); return;endif m=nb error('The columns of A must be equal the length of

14、b!' ); return;endflag='OK'% 开始计算，并赋值x=zeros(n,1);for k=1:n% 选主元 max1=0; for i=k:n if abs(A(i,k)>max1 max1 = abs(A(i,k); r=i; end end if max1<1e-10 falg='failure' return;end% 交换两行 if r>k for j=k:n z=A(k,j); A(k,j)=A(r,j); A(r,j)=z; end z=B(k); B(k)=B(r); B(r)=z;end% 消元计算

15、B(k)=B(k)/A(k,k); for j=k+1:n A(k,j)=A(k,j)/A(k,k); end for i=1:n if i=k for j=k+1:n A(i,j)=A(i,j)-A(i,k)*A(k,j); end B(i)=B(i)-A(i,k)*B(k); end end end % 输出 x for i=1:n x(i)=B(i); end 3）利用所编写的程序求解方程组，其中 ，输入矩阵A与B，调用函数Gau_Jor.m，输出结果。程序运行结果如图：程序说明： 调用函数Gau_Jor.m，并输入A与B，就会得到结果。 A为方程组的系数矩阵：B 为方程组的右端项；X为

16、方程组的解；Flag为指标向量，flag=failure表示计算失败，flag=OK表示计算成功。n,m = size(A);nb=length(B);当方程组行与列的维数不相等时，停止计算，并输出错误信息。 if n=m error('The rows and columns of matrix A must be equal!'); return;end当方程组与右端项的维数不匹配时，停止计算，并输出出错信息。二、总结随着两周对matlab的学习的深入，我对其了解也更加深入。MATLAB是美国MathWorks公司出品的商业数学软件，用于算法开发、数据可视化、数据分析以及数

17、值计算的高级技术计算语言和交互式环境，主要包括MATLAB和Simulink两大部分。MATLAB是矩阵实验室（Matrix Laboratory）的简称，MATLAB和Mathematica、Maple并称为三大数学软件它在数学类科技应用软件中在数值方面。  在这短短的两周内我从开始的一头雾水，到自己看书学习，到与同学讨论，再进行整个题目的理论分析和计算，参考课程上的代码，写出自己的代码。我也明白了学无止尽的道理，在我所查的很多参考书中，很多知识是我从没有接触过的，我对它的了解还仅限于皮毛，对它的很多功能以及函数还不是很了解，所以在这个学习的过程中我穿越在知识的海洋中，一

18、点一点吸取着它的知识。在MATLAB编程中需要很多的参考书，要尽量多的熟悉matlab自带的函数及其作用，因为matlab的自带函数特别多，基本上能够满足一般的数据和矩阵的计算，所以基本上不用你自己编函数。这一点对程序非常有帮助，可以使程序简单，运行效率高，可以节省很多时间。本次课设中我的体会是把基本的知识看过之后，然后就找一个实际的程序来动手编一下，不要等所有的知识都学好之后再去编程，你要在编程的过程中学习，程序需要什么知识再去补充，编程是一点一点积累的，所以需要做一些随手笔记什么的。在编写程序代码时，需要什么函数，需要什么模块就应该去着重看那个知识点，不要一步登天，一步一步学，如果太急于把

19、所有东西都学到，也是不好的，更是实现不了的。所以一天一天积累的，慢慢地学这个软件。 总之，通过这次课设，我了解了一些这个软件的功能，以及通过自己编写一些代码也学到了一些用法和知识。对于matlab我还有好多东西去学。人无论是做什么事最重要的是要对自己有信心，只有先对自己有信心的时候你才会发挥出你的才能，而这些能力是你以前所没发现的。并且做事还不能怕麻烦，不能做事三分钟热度，要有恒力，坚信自己的能力，相信自己，别人能做的自己也能做，因为自己并不比别人差，俗话说：世上无难事只怕有心人。所以只要我问有心去完成一件事就一定能做好的。完稿日期: 2013 年 12 月 23 日参考文献1mat

20、lab程序设计教程第二版.刘卫国.2005年.中国水利水电出版社.2数值分析第五版.李庆扬.王能超.易大义.2008年.清华大学出版社.3Matlab数值分析与应用.张德丰.2007年.国防工业出版社.4数值分析.曾繁慧.2009年.中国矿业大学出版社.附 录源程序:一、通识题1-1：x = 0:0.0001:1;n=2;i=1;for n=2:1:5 y = x.n; subplot(4,1,i); i=i+1; plot(x,y); title('y=xn'); xlabel('x'); ylabel('y');end1-2：x = -1:0

21、.000001:1; y = 1./(1+25.*x.2); plot(x,y); title('y=1/(1+25*x2),x from -1 to 1'); xlabel('x'); ylabel('y');1-3.1：x = -50:0.0001:50;y = atan(x);plot(x,y);title('y=atan(x)');text(1,pi/4,'y=atan(x)');xlabel('x');ylabel('y');1-3.2：x=-50:0.0001:50;su

22、bplot(2,2,1)for n=5 p1=p1+(-1)(n-1)*(x.(2*n-1)/(2*n-1);endplot(x,p1);title('n=5');xlabel('x');ylabel('p(x)');subplot(2,2,2)for n=10 p2=p2+(-1)(n-1)*(x.(2*n-1)/(2*n-1);endplot(x,p2);title('n=10');xlabel('x');ylabel('p(x)');subplot(2,2,3)for n=20 p3=p3+(

23、-1)(n-1)*(x.(2*n-1)/(2*n-1);endplot(x,p3);title('n=20');xlabel('x');ylabel('p(x)');subplot(2,2,4)for n=50 p4=p4+(-1)(n-1)*(x.(2*n-1)/(2*n-1);endplot(x,p4);title('n=50');xlabel('x');ylabel('p(x)');2、function s=wo(n)sum=0;k=1;while(k<=n) sum=sum+1/k.2

24、; k=k+1;ends=sum;3、syms x y;f=log(x)+exp(-y/x);diff(f,x);diff(f,y);ans=-diff(f,x)/diff(f,y);4、syms x y t ;x=t-log(1+t);y=t3+t2;ans=(diff(x)*diff(y,2)-diff(x,2)*diff(y)/(diff(x)3;5、x=0:0.0001:1;y=4*x.3-6*x.2+x-2;f=inline('12*x.2-12*x+1');y1=4*fzero(f,0)3-6*fzero(f,0)2+fzero(f,0)-2;y2=4*fzero(f,1)3-6*fze