数值代数上机实验报告

数值代数课程设计实验报告姓名： 班级：学号： 实验日期：一、 实验名称代数的数值解法二、 实验环境MATLAB7.0实验一、平方根法与改进平方根法—、实验要求：用熟悉的计算机语言将不选主元和列主元Gasuss消元法编写成通用的子程序，然后用编写的程序求解下列方程组nxn n用所编的程序分别求解40、84、120阶方程组的解。二、算法描述及实验步骤GAuss程序如下：求A的三角分解：A=LU■求解Ly=b得y；(3)求解Ux=y得x；列主元Gasuss消元法程序如下:1求A的列主元分解：PA=LU■2求解Ly=Pb得y；3求解Ux=y得x；三、调试过程及实验结果：% 方程系数 >>A1=Sanduijiaozhen(8,6,1,40);>>A2=Sanduijiaozhen(8,6,1,84);>>A3=Sanduijiaozhen(8,6,1,120);>>b1(1)=7;b2(1)=7;b3(1)=7;>>fori=2:39b1(i)=15;end>>b1(40)=14;>>fori=2:83b2(i)=15;end>>b2(40)=14;>>fori=2:119b1(i)=15;end>>b3(120)=14;% 方程解 >>x11=GAuss(A1,b1')>>x12=GAussZhu(A1,b1')>>x21=GAuss(A2,b2')>>x22=GAussZhu(A3,b3')>>x31=GAuss(A3,b3')>>x32=GAussZhu(A3,b3')运行结果：(n=40)GAuss消元法的解即为x11=1.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.0000列主元GAuss消元法的解即为x12=

1111111111111111 1 11 1 11 1 11 1 11 1 11 1 11 1 11 1 11 1 1六、源程序：functionA=Sanduijiaozhen(a,b,c,n)%生成n阶以a,b,c为元素的三对角阵A=diag(b*ones(1,n),0)+diag(c*ones(1,n-1),1)+diag(a*ones(1,n-1),-1);functionx=GAuss(A,b)n=length(b);x=b;% 分解 fori=1:n-1forj=i+1:nmi=A(j,i)/A(i,i);b(j)=b(j)-mi*b(i);fork=i:nA(j,k)=A(j,k)-mi*A(i,k);endAB=[A,b];endend% 回代 x(n)=b(n)/A(n,n);fori=n-1:-1:1s=0;forj=i+1:ns=s+A(i,j)*x(j);endx(i)=(b(i)-s)/A(i,i);endfunctionx=GAussZhu(A,b)n=length(b);x=b;% 选主元 fork=1:n-1a_max=0;fori=k:nifabs(A(i,k))>a_maxa_max=abs(A(i,k));r=i;endendifr>kforj=k:nz=A(k,j);A(k,j)=A(r,j);A(r,j)=z;endz=b(k);b(k)=b(r);b(r)=z;end% 消兀 fori=k+1:nm=A(i,k)/A(k,k);forj=k:nA(i,j)=A(i,j)-m*A(k,j);endb(i)=b(i)-m*b(k);endendifabs(A(n,n))==0return;endAbZhu=[A,b];% 回代 x(n)=b(n)/A(n,n);fori=n-1:-1:1forj=i+1:nb(i)=b(i)-A(i,j)*x(j);endx(i)=b(i)/A(i,i);end实验二、平方根法与改进平方根法一、实验要求：用计算机语言将平方根法和改进的平方根法编成通用的子程序，然后用编写的程序求解对称正定方程组100阶方程组AX=b，二、 算法描述及实验步骤：平方根法函数程序如下：1、 求A的Cholesky分解：A-LTL；2、 求解Ly-b得y；3、 求解Lx-y得x；改进平方根法函数程序如下：1、 求A的Cholesky分解：A-LDLT；2、 求解Ly-b得y；3、 求解DLx-y得x；三、 调试过程及实验结果：clear;clc;% 方程系数 >>A=Sanduijiaozhen(1,10,1,100);>>b(1)=11;>>fori=2:99b(i)=12;end>>b(100)=11;>>x1=Cholesky(A,b')>>x2=GJCholesky(A,b')运行结果：平方根法的解即为x1=1.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.0000改进平方根法解得的解即为x2=1.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00000.99991.00001.00001.00001.00091.00001.00001.00000.99081.00001.00001.00001.09081.00001.00001.00000.1010四、源程序：functionx=Cholesky(A,b)n=size(A);n=n(1);% x=A"-1*b;%disp('Matlab自带解即为x');% Cholesky分解fork=1:nA(k,k)=sqrt(A(k,k));A(k+1:n,k)=A(k+1:n,k)/A(k,k);forj=k+1:n;A(j:n,j)=A(j:n,j)-A(j:n,k)*A(j,k);endend% 前代法求解Ly=b forj=1:n-1b(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);% 回代法求解L'x=y A=A';forj=n:-1:2b(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('平方根法的解即为')；functionb=GJCholesky(A,b)n=size(A);n=n(1);v=zeros(n,1);% LDL'分解 forj=1:nfori=1:j-1v(i)=A(j,i)*A(i,i);endA(j,j)=A(j,j)-A(j,1:j-1)*v(1:j-1);A(j+1:n,j)=(A(j+1:n,j)-A(j+1:n,1:j-1)*v(1:j-1))/A(j,j);endB=diag(A);D=zeros(n);fori=1:nD(i,i)=B(i);A(i,i)=1;end% 前代法 A=tril(A); %得到L和Dforj=1:n-1b(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=D*(A');forj=n:-1:2b(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('改进平方根法解得的解即为');实验三、二次多项式拟合一、实验要求：用计算机语言编制利用QR分解求解线性最小二乘问题的通用子程序，用编写的程序求解一个二次多项式，逾山仕使在残向量的2范数最小的意义下拟合下面的数据t i -1-0.75-0.500.250.50.75y——i 1.000.81250.751.001.31251.752.3125二、算法描述及实验步骤:QR分解求解程序如下:1、求A的QR分解;2、 计算匕=Qtb；3、 求解上三角方程Rx=c1得x；三、调试过程及实验结果：>>t=[-1-0.75-0.500.250.50.75];>>y=[1.000.81250.751.001.31251.752.3125];>>plot(t,y,'r*');>>legend('实验数据(ti,yi)');>>xlabel('t'),ylabel('y');>>title('二次多项式拟合的数据点(ti,yi)的散点图')；运行后屏幕显示数据的散点图(略).(3)编写下列MATLAB程序计算f(x)在(％,七)处的函数值，即输入程序>>symsabc>>t=[-1-0.75-0.500.250.50.75];>>fi=a.*t.A2+b.*t+c%运行后屏幕显示关于a,b,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]编写构造残向量2范数的MATLAB程序>>y=[1.000.81250.751.001.31251.752.3125];>>y=[1.000.81250.751.001.31251.752.3125];>>fy=fi-y;fy2=fy.A2;J=sum(fy.A2);运行后屏幕显示误差平方和如下J=(a-b+c-1)A2+(9/16*a-34*b+c-13/16)A2+(14*a-1Z2*b+c-34)A2+(c-1)A2+(1/16*a+14*b+c-21/16)A2+(1/4*a+1/2*b+c-7/4)A2+(9/16*a+34*b+c-3；/16)A2. 8J 八 6J - dJ „为求a,b,c使J达到最小，只需利用极值的必要条件亍=0，丁=0，—=0，得da db cc到关于a,b,c的线性方程组，这可以由下面的MATLAB程序完成，即输入程序>>Ja1=diff(J,a);Ja2=diff(J,b);Ja3=diff(J,c);>>Ja11=simple(Ja1),Ja21=simple(Ja2),Ja31=simple(Ja3)运行后屏幕显示J分别对a,b,c的偏导数如下Ja11=451/128*a-6332*b+438*c-88万128Ja21=-63/32*a+43/8*b-3/2*c-61/32Ja31=438*a-3/2*b+14*c-1438解线性方程组Ja11=0,Ja^1=0,Ja^1=0，输入下列程序>>A=[451/128,-6332,-3/2;-63/32,438广弛;4龄广弛,14];>>B=[88；/128,61/32,143/8];>>C=B/A,f=poly2sym(C)运行后屏幕显示拟合函数f及其系数C如下C=0.3081 0.8587 1.4018f=924/2999*xA2+10301/11996*x+4204/2999故所求的拟合曲线为f(x)=0.3081x2+0.8581x+1.4018四、源程序：>>t=[-1-0.75-0.500.250.50.75];>>y=[1.000.81250.751.001.31251.752.3125];>>plot(t,y,'r*');>>legend('实验数据(ti,yi)');>>xlabel('t'),ylabel('y');>>title('二次多项式拟合的数据点(ti,yi)的散点图')；>>symsabc>>t=[-1-0.75-0.500.250.50.75];>>fi=a.*t.”2+b.*t+cfi=[ a-b+c,9/16*a-3/4*b+c, 1/4*a-1/2*b+c,