数值代数上机实验报告.doc

数值代数课程设计实验报告 姓名： 班级： 学号： 实验日期：一、实验名称 代数的数值解法二、实验环境 MATLAB7.0实验一、平方根法与改进平方根法一、实验要求： 用熟悉的计算机语言将不选主元和列主元Gasuss消元法编写成通用的子程序，然后用编写的程序求解下列方程组 用所编的程序分别求解40、84、120阶方程组的解。2、 算法描述及实验步骤GAuss程序如下：（1） 求A的三角分解：;（2） 求解得y；（3） 求解得x；列主元Gasuss消元法程序如下：1求A的列主元分解：;2求解得y；3求解得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; for i=2:39 b1(i)=15; end b1(40)=14; for i=2:83 b2(i)=15;end b2(40)=14; for i=2:119 b1(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.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000列主元GAuss消元法的解即为x12 = 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 16、 源程序：function A=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); function x=GAuss(A,b)n=length(b);x=b;%-分解-for i=1:n-1 for j=i+1:n mi=A(j,i)/A(i,i); b(j)=b(j)-mi*b(i); for k=i:n A(j,k)=A(j,k)-mi*A(i,k); end AB=A,b; endend%-回代-x(n)=b(n)/A(n,n);for i=n-1:-1:1 s=0; for j=i+1:n s=s+A(i,j)*x(j); end x(i)=(b(i)-s)/A(i,i);endfunction x=GAussZhu(A,b)n=length(b);x=b;%-选主元-for k=1:n-1 a_max=0; for i=k:n if abs(A(i,k)a_max a_max=abs(A(i,k);r=i; end end if rk for j=k:n z=A(k,j);A(k,j)=A(r,j);A(r,j)=z; end z=b(k);b(k)=b(r);b(r)=z; end %-消元- for i=k+1:n m=A(i,k)/A(k,k); for j=k:n A(i,j)=A(i,j)-m*A(k,j); end b(i)=b(i)-m*b(k); endendif abs(A(n,n)=0 return;endAbZhu=A,b;%-回代-x(n)=b(n)/A(n,n);for i=n-1:-1:1 for j=i+1:n b(i)=b(i)-A(i,j)*x(j); end x(i)=b(i)/A(i,i);end实验二、平方根法与改进平方根法一、实验要求： 用计算机语言将平方根法和改进的平方根法编成通用的子程序，然后用编写的程序求解对称正定方程组阶方程组AX=b， 二、算法描述及实验步骤： 平方根法函数程序如下： 1、求A的Cholesky分解：; 2、求解得y； 3、求解得x；改进平方根法函数程序如下： 1、求A的Cholesky分解：; 2、求解得y； 3、求解得x；三、调试过程及实验结果：clear;clc;%-方程系数- A=Sanduijiaozhen(1,10,1,100); b(1)=11; for i=2:99 b(i)=12; end b(100)=11; x1=Cholesky(A,b) x2=GJCholesky(A,b)运行结果：平方根法的解即为x1 = 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000改进平方根法解得的解即为x2 = 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 0.9999 1.0009 0.9908 1.0908 0.1010四、源程序：function x=Cholesky(A,b)n=size(A);n=n(1);% x=A-1*b; % disp(Matlab自带解即为x);%-Cholesky分解-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:n,k)*A(j,k); endend%-前代法求解Ly=b-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);%-回代法求解Lx=y- 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(平方根法的解即为);function b=GJCholesky(A,b)n=size(A);n=n(1);v=zeros(n,1); %-LDL分解-for j=1:n for i=1:j-1 v(i)=A(j,i)*A(i,i); end A(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); end B=diag(A);D=zeros(n); for i=1:n D(i,i)=B(i); A(i,i)=1; end %-前代法-A=tril(A); %得到L和Dfor 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=D*(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(改进平方根法解得的解即为);实验三、二次多项式拟合一、实验要求： 用计算机语言编制利用QR分解求解线性最小二乘问题的通用子程序，用编写的程序求解一个二次多项式使在残向量的范数最小的意义下拟合下面的数据-1 -0.75 -0.5 0 0.25 0.5 0.75 1.00 0.8125 0.75 1.00 1.3125 1.75 2.3125二、算法描述及实验步骤： QR分解求解程序如下： 1、求A的分解; 2、计算; 3、求解上三角方程得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)的散点图);运行后屏幕显示数据的散点图（略）. （3）编写下列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+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/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为求使达到最小，只需利用极值的必要条件，得到关于的线性方程组，这可以由下面的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/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.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-