第1章作业

第一题问题叙述：构造算法并编程序精确计算二次方程的根●设a≠0,𝑏2−4𝑎𝑐0，且有方程a𝑥2+𝑏𝑥+𝑐=0●包括𝑏2≈𝑏2−4𝑎𝑐的情况（a=c=1,b=±1000000.000001）问题分析：对于一般情况，可以直接使用求根公式𝑥=−𝑏±√𝑏2−4𝑎𝑐2𝑎计算，MATLAB语言代码如下但是当b2≫4𝑎𝑐时，𝑏2≈𝑏2−4𝑎𝑐，分子是两个几乎相等的浮点数相减，在计算机运算中会出现减性抵消的情况，用该求根公式会造成较大的误差。解决方法：（1）方法1：使用双精度计算（2）方法2：变换求根公式当b0时，𝑥1=−2𝑐𝑏+√𝑏2−4𝑎𝑐,𝑥2=−𝑏−√𝑏2−4𝑎𝑐2𝑎当b0时，𝑥1=−2𝑐𝑏−√𝑏2−4𝑎𝑐,𝑥2=−𝑏+√𝑏2−4𝑎𝑐2𝑎这样可以很好地避免出现减性抵销现象cleara=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)实例分析：以𝐚=𝐜=𝟏,𝐛=±𝟏𝟎𝟎𝟎𝟎𝟎𝟎.𝟎𝟎𝟎𝟎𝟎𝟏为例（1）方法1：使用双精度计算（与单精度作对比）对于不同的系数，计算结果如下：方程系数单精度双精度abcx1x2x1x2175-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：变换求根公式当b0时，%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)𝑥1=−2𝑐𝑏+√𝑏2−4𝑎𝑐,𝑥2=−𝑏−√𝑏2−4𝑎𝑐2𝑎当b0时，𝑥1=−2𝑐𝑏−√𝑏2−4𝑎𝑐,𝑥2=−𝑏+√𝑏2−4𝑎𝑐2𝑎计算结果如下：方程系数单精度双精度abcx1x2x1x2175-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≫4𝑎𝑐时，变换公式能够有%calculateusingdouble%forthepurposesofcomparison,alsousesingletocalculatecleara=input('a=');b=input('b=');c=input('c=');%a,b,cisthecoefficientofequationax^2+bx+c=0delta=b*b-4*a*c;%first,usingdoubleifb0d1=-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);ifb0s1=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)效地避免减性抵销，在单精度的计算结果可以看出，即使当b=1000000.000001时，x1仍然能够计算出，而不是像之前一样计算结果为0.0第二题问题叙述：正弦函数的无穷级数展开为sin𝑥=∑(−1)𝑖𝑥2𝑖+1(2𝑖+1)!∞𝑖=0（1）分别以单精度和双精度数据类型，计算x=0.5时近似值，要求计算结果具有4位有效数字；（2）如果采用单精度数据类型要求计算结果达到机器精度，此时结果如何？（测试机器精度，满足1+ϵ1时的最小浮点数）问题分析：本题应使用迭代法求解。当前迭代结果对真值的误差为ϵr=(真值−近似值)真值∗100%当前迭代结果的误差为ϵa=当前近似值−前一近似值当前近似值∗100%为了使计算结果达到一定准确度，应先设定一个容限ϵs，当|ϵa|𝜖𝑠时，迭代结束。在本题中，计算结果需要具有4位有效数字，则取ϵs=(0.5∗102−𝑛)%=5∗10−4问题解答：（1）先计算双精度数据类型计算结果如下表所示：（注：真值为0.479425538604203）迭代计算值当前迭代结果误差ϵa对真值相对误差ϵr0.5000000000/-0.04291480.4791666666-0.043478265.3996276e-40.47942708335.4318305e-4-3.2220418e-60.4794255332-3.233243e-61.120106e-8%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）迭代计算值当前迭代结果误差ϵa对真值相对误差ϵr0.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变换到合适的区间，从而避免出现级数项过大的情况。%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)%Relativeerroroftruevalueend（2）首先计算单精度数据类型下的机器精度测试结果是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的差值。对于不同的数据类型和不同的值，机器精度也不一样。题目中要求采用单精度数据类型，要求计算结果达到机器精度。可以理解为，当前计算值与前一计算值的差值要小于等于当前计算值的机器精度。%amachineprecisiontestclearepsilon=single(1);while(single(epsilon+1)single(1))epsilon=single(epsilon/2);endepsilon=single(epsilon*2)%theresultofmachineprecisioninthistesteps(single(1))%themachineprecisioncalculatedbybuilt-infunction计算结果如下表所示：（注：真值为0.479425538604203）迭代计算值当前计算值机器精度eps_y对真值相对误差ϵr0.50000005.9604645e-08-0.04291480.47916672.8560558e-085.3998345e-040.