数值计算实验报告MATLAB线性方程组的求解方法

评分签名日期湖南商学院实验报告课程名称数值分析实验名称数值计算基本概念及算法基础专业班级信科1301指导老师胡桔州组长姓名陈平学号130320015电话[phone]QQ1274839822成员姓名周红学号130320021王双学号130320022贺嘉玲学号1303200262014—2015学年度第2学期第1页共15页第1页共15页一、实验目的与要求1、通过实验，进一步了解和掌握浮点数及其运算特点，了解阶段误差与舍入误差这两种主要误差及其对结果的影响；2、通过实验，初步理解数值算法的复杂性、算法的稳定性、问题的病态性；3、通过实验，比较数值解法与纯数学解法的区别，体会算法在工程与科学计算中的地位，并形成有意识比较算法、选择算法的良好习惯；从而进一步提高设计算法的能力。二、实验环境（包括硬件、软件配置）硬件：屏幕尺寸：14英寸1366x768CPU型号：Intel酷睿i34030U...(2种)CPU主频：1.9GHz...(2种)内存容量：4GB（4GB×1）DDR3硬盘容量：500GB显卡芯片：AMDRadeonR5M230操作系统：预装Windows764bit软件：开发环境：Window7旗舰版开发工具：MATLABR2009a(MATLAB7.8)第2页共15页第2页共15页第3页共15页第3页共15页三、项目内容（包括相关准备知识）浮点数的概念，浮点数的运算，截断误差和舍入误差的区别以及它们对结果的影响。以及用MATLAB比较两个不同算法的稳定性，通过研究和构造特定算法解决病态问题。要做好实验，这些都是必须思考并掌握的技能。(1)熟悉MATLAB的主要操作命令。(2)学会简单的算法编程。(3)掌握简单的程序调适。(4)用MATLAB编程并利用函数解决问题。第4页共15页第4页共15页四、实验步骤：（要求详细对实验数据的处理、源程序代码、算法、实验结果、图表等进行描述，可以根据情况自己添加页）一、浮点数及其运算特点1、浮点数的分布functionf=fudianshu(beta,t,L,U)m=1/beta:1/beta^(t):1/beta+(beta^(t-1)-1)/beta^(t);n=length(m);s=L:U;fork=1:(U-L+1)formatratf(k,:)=beta^s(k)*m;f(k,:)plot(f,zeros(1,n),'.')axis([04-11])pauseendpausesubplot(2,1,2)hist(f(:))%其中beta为基数，t为精确度，L和U分别是指数部分的上界和下界。类似于二进制的情形。在命令窗口输入：fudianshu(2,3,-2,2)ans=1/85/323/167/32ans=1/45/163/87/16ans=1/25/83/47/8ans=15/43/27/4ans=25/237/2ans=1/85/323/167/321/45/163/87/161/25/83/47/815/43/27/4第5页共15页第5页共15页25/237/200.511.522.533.54-1-0.500.5100.511.522.533.502468分析：从图中可以看出浮点数在接近其下界时比较稠密，而接近上界时则相对比较稀疏。我们要思考的是这种分布不均匀性会给数值计算带来问题吗？假设x是精确值，p是对x的近似，则a=|x-p|;称为p近似x的绝对误差。而r=a/|x|=|x-p|/|x|,称为p近似x的相对误差。可见用浮点数近似接近上界的实数时，由于这时浮点数的间隔大，所以近似的绝对误差则会比较大；然而只要被近似的数在在下界和上界之间，相对误差在各处是一样的。所以数值计算中更多的还是相对误差的问题，因此浮点数分布的不均匀并不会影响算法噢！2、双精度浮点数下的问题--舍入误差formatlongx=4/3-1x=0.333333333333333y=3*xy=1.000000000000000z=1-yz=2.220446049250313e-016按照理论计算，得到的z值应当是0.而实际输出的非零。其实得到的z就是机器精度。如果计算中产生一个小于realmin的数则是“下溢”；如果产生了一个大于realmax的数则是“上溢”。我们再来看一个多项式的计算：x=0.988:0.001:1.012;y=x.^7-7*x.^6+21*x.^5-35*x.^4+35*x.^3-21*x.^2+7*x-1;plot(x,y)第6页共15页第6页共15页0.9850.990.99511.0051.011.015-4-3-2-1012345x10-14从画出的图来看，这完全不像一个多项式！这个现象就是舍入误差在起作用。3、浮点数不满足结合律和交换律结合律formatlonga=rand(1,20000);b=rand(1,20000);m=2^20*rand(1,20000);s1=sum(a)+sum(m.*b);s2=sum(a+m.*b);s1,s2sub=s1-s2输出结果为：s1=5.261188112279733e+009s2=5.261188112279694e+009sub=3.910064697265625e-005交换律m=100000;x=1:m;y=m:-1:1;s=m/(m+1);disp(['s=',num2str(s)])s1=sum(1./(x.*(x+1)));s2=sum(1./(y.*(y+1)));fprintf('s1=%31.30f\ns2=%31.30f\n',s1,s2)fprintf('s1-s=%31.30f\ns2-s=%31.30f\n',s1-s,s2-s)fprintf('s1-s2=%31.30f\n',s1-s2)输出结果为：s=0.99999s1=0.999990000100012160000000000000s2=0.999990000099998940000000000000s1-s=0.000000000000013100631690576847分析：实验结果表明，浮点数加法不满足交换律和结合律。这是为什么呢？首先我们要了解计算机浮点数的加法是如何进行的。假设有两个浮点数：第7页共15页第7页共15页d=+-.d1*d2*...*dt*beta1,c=+-.c1*c2*...*ct*beta2。两个浮点数相加首先要比较它们的阶码。如果阶码相同则尾数相加，相加后尾数如果大于1，阶码进位；如果阶码不等，以相对大的阶码为标准，将阶码小的浮点数移位直到阶码一致，再按阶码相等时的规则进行相加。如果一个大数加一个非常小的数，小数的个数越多，计算的误差就越大。二、算法的复杂性1、测试多项式函数计算的不同算法在运行时间上的差别①参考教材P35......tic;x=3;f=2*x^5-5*x^4-4*x^3-6*x+7;t=toct=58.987017134798052tic;x=3;x^2==x*x;x^3==x^2*x;x^4==x^3*x;x^5==x^4*x;f=2*x^5-5*x^4-4*x^3-6*x+7;t=toct=2.152971514854767e+002tic;x=3;f=((((2*x-5)*x-4)*x+3)*x-6)*x+7;t=toct=10.009796315604856x=3;tic;x=3;f=((((2*x-5)*x-4)*x+3)*x-6)*x+7;t=toct=10.840845190320568②tic;A=rand(20,20);t=toct=第8页共15页第8页共15页31.961928037591385tic;A=rand(20,20);b=det(A);t=toct=42.9568197776177172、误差传播与算法的稳定性算法一：用递推公式E1=1/e。En=1-n*En-1。n=2,3,...,20。clearep(1)=exp(-1);forn=2:9ep(n)=1-n*ep(n-1);ep(n)=round(10000*ep(n))/10000;endep'plot(ep,'b*');xlabel('算法一');%axis([0100-1010])第9页共15页第9页共15页123456789-1012345678算法一假设En=In-Ln,则有En=-n*En-1=(-1)^n!En。因此误差是逐渐放大的。算法二、递推公式E20=0，En-1=(1-En-1)/n,n=19,18,...,1。%误差传播与算法稳定性%稳定的算法clearep(500)=0;forn=500:-1:2ep(n-1)=(1-ep(n))/n;ep(n-1)=round(10000*ep(n-1))/10000;endepplot(ep,'b-');axis([0500-0.101])ep(1:9)'xlabel('算法二');%axis([0100-1010])050100150200250300350400450500-0.100.10.20.30.40.50.60.70.80.9算法二假设En=In-Pn，则有|En-1|=1/n*|En|，Ek=(-1)^(n-k)*(k!/n!)*En。因此误差是逐渐缩小的。第10页共15页第10页共15页分析：由此我们可以判断算法二比算法一计算稳定。3、病态问题的初步讨论%病态问题实验%代数多项式的根及其扰动后多项式的根p=poly(1:20);ep=zeros(1,21);ep(3)=1.0e-5;re=roots(p+ep)plot(re,'b+');holdonplot(1:20,0,'r*');holdoff0510152025-4-3-2-101234分析：从图中我们可以看出，虽然扰动项的系数非常小，从个人感觉上说对多项式的解是不会有太大的影响。但通过计算我们发现，有很多解却差异极大。即有很多解对这种扰动行为极其敏感。这个实验告诉我们，不但算法有好坏之分，被求解的问题也有好坏之分，甚至在同一问题中有好的部分和坏的部分。数值分析的大部分领域中，如线性代数方程组、矩阵特征值问题等等，都有关于病态问题的深入研究。病态问题在付出一些代价(如耗用更多机器时间，占用更多存储空间等)后，能通过研究和构造特殊的算法来解决。三、思考题1、函数sinx有幂级数展开，利用幂级数计算sinx的matlab程序为：functions=powersin(x)s=0;t=x;n=1;whiles+t~=s;s=s+t;t=-x^2/((n+1)*(n+2))*t;第11页共15页第11页共15页n=n+2;End(a)解释上述程序的终止准则(b)对于x=pi/2，11pi/2，21pi/2,计算的精度是多少？分别需要计算多少项？答：（a）终止准则为终止条件为t=0,但是t=-x^2/((n+1)*(n+2))*t始终不为0。不过在计算中，当t小于计算精度时，计算机识别为0，此时计算才会终止。（b）x=pi/2;n=23;t=-1.2539e-018;s=1.0000;x=11*pi/2;n=75;t=-2.6232e-017;s=-1.0000;x=21*pi/2;n=121;t=6.4693e-018;s=0.9999;2、已知数列{Xn}满足递推公式Xn=10X（n-1）-1，n=1,2，...若取X0=1.41（3位有效数字），问按上述递推公式，从X0计算到X10误差有多大？这个计算过程稳定吗？若不稳定请设计一种稳定的递推的算法。在M文件窗口输入functionxk=suanfayi(x);x(1)=1.41;forn=2:10x(n)=10*x(n-1)-1;disp(x(n));end输出结果为13.1000130129912989129889129888912988889129