第1章作业完整版本

第一题问题叙述：构造算法并编程序精确计算二次方程的根●设a≠0,b2●包括b2≈b问题分析：对于一般情况，可以直接使用求根公式x=-bclearcleara=input('a=');b=input('b=');c=input('c=');x1=-b/2*a+1/2*a*sqrt(b*b-4*a*c)x2=-b/2*a+1/2*a*sqrt(b*b-4*a*c)但是当b2≫4ac时，b2≈解决方法：（1）方法1：使用双精度计算（2）方法2：变换求根公式当b>0时，x当b<0时，x这样可以很好地避免出现减性抵销现象实例分析：以a=c=1,b=±1000000.000001为例（1）方法1：使用双精度计算（与单精度作对比）%calculateusingdoubleprecision%calculateusingdoubleprecision%forthepurposesofcomparison,alsousesingletocalculatecleara=input('a=');b=input('b=');c=input('c=');%a,b,cisthecoefficientofequationax^2+bx+c=0delta=b*b-4*a*c;%first,usingdoubled1=(-b+sqrt(delta))/2/a;d2=(-b-sqrt(delta))/2/a;d1=vpa(d1,12),d2=vpa(d2,12)%than,usingsingledelta=single(delta);a=single(a);b=single(b);c=single(c);s1=single((-b+sqrt(delta))/2/a);s2=single((-b-sqrt(delta))/2/a);s1=vpa(s1,12),s2=vpa(s2,12)对于不同的系数，计算结果如下：方程系数单精度双精度abcxxxx175-0.807417631-6.192582130-0.807417596-6.19258240311000010.0-10000.0-0.00010000-9999.999911000000.00000110.0-1000000.01.0000076e-6-1000000.01-1000000.00000110.01000000.01.0000076e-61000000.0结论：显然，在一般情况下，单精度还能保持一定准确度；但是当a=1，b=10000，c=1时，减性抵销现象很明显，单精度下的x1计算结果为0，在b=1000000.000001（2）方法2：变换求根公式当b>0时，x当b<0时，x%calculateusingdouble%calculateusingdouble%forthepurposesofcomparison,alsousesingletocalculatecleara=input('a=');b=input('b=');c=input('c=');%a,b,cisthecoefficientofequationax^2+bx+c=0delta=b*b-4*a*c;%first,usingdoubleifb>0d1=-2*c/(b+sqrt(delta));%usingnewformulad2=(-b-sqrt(delta))/2/a;elsed1=-2*c/(b-sqrt(delta));%usingnewformulad2=(-b+sqrt(delta))/2/a;endd1=vpa(d1,12),d2=vpa(d2,12)%than,usingsingledelta=single(delta);a=single(a);b=single(b);c=single(c);ifb>0s1=single(-2*c/(b+sqrt(delta)));%usingnewformulas2=single((-b-sqrt(delta))/2/a);elses1=single(-2*c/(b-sqrt(delta)));%usingnewformulas2=single((-b+sqrt(delta))/2/a);ends1=vpa(s1,12),s2=vpa(s2,12)计算结果如下：方程系数单精度双精度abcxxxx175-0.80741763-6.19258213-0.80741759-6.1925824011000019.99999975e-5-10000.0-0.00010000-9999.999911000000.0000011-9.99999975e-7-1000000.0-0.000001-1000000.01-1000000.00000119.99999975e-71000000.00.0000011000000.0结论： 在普通情况下，变换公式的提升效果不明显。但是当b2≫4ac时，变换公式能够有效地避免减性抵销，在单精度的计算结果可以看出，即使当b=1000000.000001时，第二题问题叙述：正弦函数的无穷级数展开为sin（1）分别以单精度和双精度数据类型，计算x=0.5时近似值，要求计算结果具有4位有效数字；（2）如果采用单精度数据类型要求计算结果达到机器精度，此时结果如何？（测试机器精度，满足1+问题分析：本题应使用迭代法求解。当前迭代结果对真值的误差为ϵ当前迭代结果的误差为ϵ为了使计算结果达到一定准确度，应先设定一个容限ϵs，当ϵa<ϵϵ

问题解答：（1）先计算双精度数据类型%usingiterativemethodtocalculate%usingiterativemethodtocalculatesin(0.5),doubleprecisioncleari=0;x=0.5;es=0.0005;ea=1;%initializetheerrorofcurrentiterationresultsitem=(-1)^i*x^(2*i+1)/factorial(2*i+1);y=item%yistheestimateofsin(0.5)er=(sin(0.5)-y)/sin(0.5)%Relativeerroroftruevaluewhile(abs(ea)>es)i=i+1;item=(-1)^i*x^(2*i+1)/factorial(2*i+1);y=y+itemea=item/yer=(sin(0.5)-y)/sin(0.5)%Relativeerroroftruevalueenddisp('therealvalueis:sin(0.5)=');disp(sin(0.5));计算结果如下表所示：（注：真值为0.479425538604203）迭代计算值当前迭代结果误差ϵ对真值相对误差ϵ0.5000000000/-0.04291480.4791666666-0.043478265.3996276e-40.47942708335.4318305e-4-3.2220418e-60.4794255332-3.233243e-61.120106e-8

下面计算单精度数据类型：%usingiterativemethodtocalculatesin(0.5)%usingiterativemethodtocalculatesin(0.5)，singleprecisioncleari=single(0);x=single(0.5);es=single(0.0005);ea=single(1);%initializetheerrorofcurrentiterationresultsitem=single((-1)^i*x^(2*i+1)/factorial(2*i+1));y=single(item)%yistheestimateofsin(0.5)er=(sin(0.5)-y)/sin(0.5)%Relativeerroroftruevaluewhile(single(abs(ea))>es)i=i+1;item=single((-1)^i*x^(2*i+1)/factorial(2*i+1));y=single(y+item)ea=single(item/y)er=(sin(0.5)-y)/sin(0.5)%Relativeerroroftruevalueenddisp('therealvalueis:sin(0.5)=');disp(sin(0.5));计算结果如下表所示：（注：真值为0.479425538604203）迭代计算值当前迭代结果误差ϵ对真值相对误差ϵ0.5000000/-0.04291480.4791667-0.04347835.3998345e-040.47942715.4318312e-04-3.1930326e-060.4794255-3.2332430e-063.9420897e-08结论： 可以看出，无论是单精度还是双精度，计算结果都有很高的精确度。实际上，两种计算精度都在第3次迭代就已经达到了4位有效数字的要求，可见迭代法在该问题有很高的效率。注：本题的求和级数是正负交替的级数，在求和过程中，如果某一项的值大于和值本身时，就会出现拖尾效应，会抵销大量的有效数字。在本题中，由于sin(0.5)的值较大，而且x=0.5较小，所以级数项衰减很快，不会出现拖尾效应。对于可能出现此种效应的情况，可以利用正弦函数的周期性，将x变换到合适的区间，从而避免出现级数项过大的情况。%amachineprecisiontestclearepsilon=single(1);while(single(epsilon+1)>single(1))%amachineprecisiontestclearepsilon=single(1);while(single(epsilon+1)>single(1))epsilon=single(epsilon/2);endepsilon=single(epsilon*2)%theresultofmachineprecisioninthistesteps(single(1))%themachineprecisioncalculatedbybuilt-infunction测试结果是1.1920929e-07，与MATLAB内置的eps函数的计算结果1.1920929e-07相同。在测试中发现，对于不同的类型和不同的值，机器精度是不一样的。例如double类型的1的机器精度是2.220446049250313e-16；double类型的2的机器精度是4.440892098500626e-16，而single类型的2的精度是2.3841858e-07。那么，上面这段代码求的是single类型的1的机器精度，即与单精度浮点数1最相近的浮点数与1的差值。对于不同的数据类型和不同的值，机器精度也不一样。题目中要求采用单精度数据类型，要求计算结果达到机器精度。可以理解为，当前计算值与前一计算值的差值要小于等于当前计算值的机器精度。%usingiterativemethodtocalculate%usingiterativemethodtocalculatesin(0.5)formatlongcleari=single(0);x=single(0.5);item=single((-1)^i*x^(2*i+1)/factorial(2*i+1));y=single(item)%yistheestimateofsin(0.5)er=(sin(0.5)-y)/sin(0.5)%Relativeerroroftruevalueeps_y=single(y);%themachineprecisionofywhile(single(eps_y+y)>single(y))eps_y=single(eps_y/2);endeps_y=single(eps_y*2)%theinitialvalueoftheprecisionwhile(abs(single(item))>eps_y)i=i+1;item=single((-1)^i*x^(2*i+1)/factorial(2*i+1));y=single(y+item)