数值线性代数实验

1、数值线性代数实验题目： 数值线性代数 专业：信息与计算科学 班级： 班 姓名： 山 东 科 技 大 学2013年 1 月 16日实验报告说明学院：信息学院 专业：信息 班级10-2 姓名： 一、主要参考资料：（1） Matlab数值计算-案例分析北京航空出版 （2） Matlab数值分析 机械工业出版 二、课程设计应解决的主要问题：（1）平方根 （2）QR方法 （3）最小二乘法 三、应用软件：（1）Matlab7.0（2）数学公式编辑器四、发出日期： 课程设计完成日期： 指导教师签字： 系主任签字： 指导教师对课程设计的评语指导教师签字： 年 月 日一、问题描述先用你所熟悉的计算机语言将平方根

2、和改进的平方根法编成写通用的子程序，然后用你编写的程序求解对称正定方程组，其中（1）b随机的选取，系数矩阵位100阶矩阵（2）系数矩阵为40阶Hilbert矩阵，即系数矩阵A的第i行第j列元素为，向量b的第i个分量为。二、分析与程序1. 平方根法函数程序如下：function x,b=pingfanggenfa(A,b)n=size(A);n=n(1);x=A-1*b; disp(Matlab自带解即为x);for k=1:n A(k,k)=sqrt(A(k,k); A(k+1:n,k)=A(k+1:n,k)/A(k,k); for j=k+1:n; A(j:n,j)=A(j:n,j)-A(j

3、:n,k)*A(j,k); endend for j=1:n-1 b(j)=b(j)/A(j,j); b(j+1:n)=b(j+1:n)-b(j)*A(j+1:n,j);endb(n)=b(n)/A(n,n); A=A;for j=n:-1:2 b(j)=b(j)/A(j,j); b(1:j-1)=b(1:j-1)-b(j)*A(1:j-1,j);endb(1)=b(1)/A(1,1); disp(平方根法的解即为b);endfunction x=ave(A,b,n) 求解Ax=bL=zeros(n,n); D=diag(n,0); S=L*D; for i=1:n %L的主对角元素均为1L(

4、i,i)=1;end for i=1:n for j=1:n if (eig(A) In pingfanggenfa at 4 In qiujie at 10Matlab自带解即为x平方根法的解即为bx = 1.6035 8.9685 0.8562 1.0195 0.9375 -50.2500 -3.0000 -16.0000 24.0000 -49.5000 -30.0000 39.0000 22.0000 -64.0000 -12.0000 2.0000 10.2500 -10.5000 -1.0000 -10.8750 83.0000 46.0000 -98.0000 12.0000 -

5、69.0000 68.0000 21.0000 17.0000 -50.7188 -8.7500 -8.0000 112.0000 6.0000 -68.7500 22.0000 44.0000 -28.0000 8.0000 -44.0000 12.0000b = 1.0e+007 * 0.0000 -0.0000 0.0001 -0.0004 -0.0014 0.0424 -0.2980 1.1419 -2.7335 4.2539 -4.3018 2.7733 -1.1989 0.5406 -0.3688 0.3285 -0.4438 0.4621 -0.2513 0.0565 0.000

6、0 -0.0051 0.0071 -0.0027 -0.0031 0.0036 -0.0019 0.0009 0.0002 -0.0002 -0.0006 0.0004 0.0001 -0.0002 0.0001 0.0000 -0.0000 0.0000 -0.0000 -0.0000改进平方根法解得的解即为bb = 1.0e+024 * 0.0000 -0.0000 0.0001 -0.0012 0.0139 -0.0954 0.4208 -1.2101 2.0624 -1.0394 -3.3343 6.2567 -0.2463 -7.4594 2.8030 3.6990 0.7277 -

7、1.7484 -0.4854 -3.6010 0.2532 5.1862 -2.1299 1.4410 0.8738 -4.5654 1.0422 4.0920 -2.7764 -2.2148 -0.8953 0.3665 4.8967 1.0416 0.1281 -4.3387 -1.1902 -2.8334 8.4610 -3.6008一、问题描述先用你所熟悉的计算机语言将算法2.5.1编成写通用的子程序，然后用你编写的程序完成下面两个计算任务：（1） 估计5到20阶Hilbert矩阵的范数条件数；（2）设先随机的选取，并计算出；然后再用列主元Gauss消去法求解该方程组，假定计算解为。试

8、对n从5到30估计计算解的精度，并且与真实相对误差做比较。二、分析与程序Hilbert矩阵：Hilbert矩阵的分量满足H(i,j)=1/(i+j-1)比如3阶Hilbert矩阵是1/1 1/2 1/31/2 1/3 1/41/3 1/4 1/5程序1. 用MathCAD计算方法如下：其中表示n阶Hilbert矩阵，即是n阶Hilbert矩阵的范数条件数2. 比较估计精度与真实相对误差用MathCAD计算方法如下： 其中表示n阶矩阵，是任选的，即是，和分别是对进行LU分解得到的下三角矩阵和上三角矩阵，是用列主元Gauss消去法求得的解，是n阶矩阵计算解的精度，是n阶矩阵计算解的真实相对误差。3

9、. 估计范数条件数用MathCAD处理计算的结果如下：当矩阵阶数 时，即是m阶Hilbert矩阵的范数条件数。以矩阵阶数m为横坐标，条件范数为纵坐标，其图像如下：条件数在一定程度上刻画了扰动对方程组解的影响程度。通常，若线性方程组的系数矩阵A的条件数(A)很大，则A是病态的；反之，A是良态的。可以看出，Hilbert矩阵的条件数很大。4. 比较估计精度与真实相对误差用MathCAD处理计算的结果如下：表示i阶矩阵计算解的精度，在图中以红色的点表示，其中 ；表示j阶矩阵计算解的真实相对误差，在图中以蓝色的点表示， 。 以i，j为横坐标，为纵坐标，在同一图像上表示如下：这一方法给出了计算解相对误差

10、的相当好的估计，可以看出，真实相对误差不大于计算解的精度。5. 结论条件数在一定程度上刻画了扰动对方程组解的影响程度。通常，若线性方程组的系数矩阵A的条件数(A)很大，则A是病态的；反之，A是良态的。可以看出，Hilbert矩阵的条件数很大，故Hilbert矩阵是十分病态的。一、 问题描述用你所熟悉的计算机语言编制利用QR分解求解线性方程组和线性最小二乘问题的通用子程序，并用你编写的程序完成下面两个计算任务：（1）求解第一章上机练习题中的三个线性方程组，并将所得的计算结果与前面的计算结果相比较说明各方法的优劣；（2）求一个二次多项式使在残向量的2范数最小意义下拟合下面的数据-1-0.75-0.

11、500.250.50.751.000.81250.751.001.31251.752.3125（3）在房产估价的线性模型中，分别表示税、浴室数目、占地面积、居住面积、车库数目、房屋数目、居室数目、房龄、建筑类型、户型及壁炉数目，y代表房屋价格。先根据如下的28组数据求出模型中的参数的最小二乘结果。y25.929.527.925.929.929.930.928.984.982.935.931.531.030.930.028.936.941.940.543.937.537.944.537.938.936.945.841.0二、 分析及程序利用QR分解求解程序如下： 1、求A的分解; 2、计算; 3

12、、求解上三角方程得x；调试过程及实验结果： t=-1 -0.75 -0.5 0 0.25 0.5 0.75; y= 1.00 0.8125 0.75 1.00 1.3125 1.75 2.3125; plot(t,y,r*); legend(实验数据(ti,yi); xlabel(t), ylabel(y); title(二次多项式拟合的数据点(ti,yi)的散点图);编写下列MATLAB程序计算在处的函数值，即输入程序 syms a b c t=-1 -0.75 -0.5 0 0.25 0.5 0.75; fi=a.*t.2+ b.*t+c%运行后屏幕显示关于 的线性方程组fi = a-b+

13、c,9/16*a-3/4*b+c,1/4*a-1/2*b+c,c,1/16*a+1/4*b+c,1/4*a+1/2*b+c,9/16*a+3/4*b+c编写构造残向量2范数的MATLAB程序 y= 1.00 0.8125 0.75 1.00 1.3125 1.75 2.3125; y= 1.00 0.8125 0.75 1.00 1.3125 1.75 2.3125; fy=fi-y; fy2=fy.2; J=sum(fy.2);运行后屏幕显示误差平方和如下J=(a-b+c-1)2+(9/16*a-3/4*b+c-13/16)2+(1/4*a-1/2*b+c-3/4)2+(c-1)2+(1/1

14、6*a+1/4*b+c-21/16)2+(1/4*a+1/2*b+c-7/4)2+(9/16*a+3/4*b+c-37/16)2为求使达到最小，只需利用极值的必要条件，得到关于的线性方程组，这可以由下面的MATLAB程序完成，即输入程序 Ja1=diff(J,a); Ja2=diff(J,b); Ja3=diff(J,c); Ja11=simple(Ja1), Ja21=simple(Ja2), Ja31=simple(Ja3) 运行后屏幕显示J分别对的偏导数如下 Ja11 = 451/128*a-63/32*b+43/8*c-887/128 Ja21 = -63/32*a+43/8*b-3/

15、2*c-61/32 Ja31 = 43/8*a-3/2*b+14*c-143/8解线性方程组，输入下列程序 A=451/128, -63/32, -3/2 ;-63/32,43/8,-3/2;43/8,-3/2,14; B=887/128,61/32,143/8; C=B/A, f=poly2sym(C)运行后屏幕显示拟合函数f及其系数C如下C = 0.3081 0.8587 1.4018f = 924/2999*x2+10301/11996*x+4204/2999 故所求的拟合曲线为 源程序： t=-1 -0.75 -0.5 0 0.25 0.5 0.75; y= 1.00 0.8125 0

16、.75 1.00 1.3125 1.75 2.3125; plot(t,y,r*); legend(实验数据(ti,yi); xlabel(t), ylabel(y); title(二次多项式拟合的数据点(ti,yi)的散点图); syms a b c t=-1 -0.75 -0.5 0 0.25 0.5 0.75; fi=a.*t.2+ b.*t+c fi = a-b+c, 9/16*a-3/4*b+c, 1/4*a-1/2*b+c, c, 1/16*a+1/4*b+c, 1/4*a+1/2*b+c, 9/16*a+3/4*b+c y= 1.00 0.8125 0.75 1.00 1.312

17、5 1.75 2.3125; y= 1.00 0.8125 0.75 1.00 1.3125 1.75 2.3125; fy=fi-y; fy2=fy.2; J=sum(fy.2) J = (a-b+c-1)2+(9/16*a-3/4*b+c-13/16)2+(1/4*a-1/2*b+c-3/4)2+(c-1)2+(1/16*a+1/4*b+c-21/16)2+(1/4*a+1/2*b+c-7/4)2+(9/16*a+3/4*b+c-37/16)2 Ja1=diff(J,a); Ja2=diff(J,b); Ja3=diff(J,c); Ja11=simple(Ja1), Ja21=simple(Ja2), Ja31=simple(Ja3) Ja11