數(shù)值代數(shù)上機(jī)實(shí)驗(yàn)報(bào)告

數(shù)值代數(shù)課程設(shè)計(jì)實(shí)驗(yàn)報(bào)告姓名： 班級：學(xué)號(hào)： 實(shí)驗(yàn)日期：一、 實(shí)驗(yàn)名稱代數(shù)的數(shù)值解法二、 實(shí)驗(yàn)環(huán)境MATLAB7.0實(shí)驗(yàn)一、平方根法與改進(jìn)平方根法—、實(shí)驗(yàn)要求：用熟悉的計(jì)算機(jī)語言將不選主元和列主元Gasuss消元法編寫成通用的子程序，然后用編寫的程序求解下列方程組nxn n用所編的程序分別求解40、84、120階方程組的解。二、算法描述及實(shí)驗(yàn)步驟GAuss程序如下：求A的三角分解：A=LU■求解Ly=b得y；(3)求解Ux=y得x；列主元Gasuss消元法程序如下:1求A的列主元分解：PA=LU■2求解Ly=Pb得y；3求解Ux=y得x；三、調(diào)試過程及實(shí)驗(yàn)結(jié)果：% 方程系數(shù) >>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')運(yùn)行結(jié)果：(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實(shí)驗(yàn)二、平方根法與改進(jìn)平方根法一、實(shí)驗(yàn)要求：用計(jì)算機(jī)語言將平方根法和改進(jìn)的平方根法編成通用的子程序，然后用編寫的程序求解對稱正定方程組100階方程組AX=b，二、 算法描述及實(shí)驗(yàn)步驟：平方根法函數(shù)程序如下：1、 求A的Cholesky分解：A-LTL；2、 求解Ly-b得y；3、 求解Lx-y得x；改進(jìn)平方根法函數(shù)程序如下：1、 求A的Cholesky分解：A-LDLT；2、 求解Ly-b得y；3、 求解DLx-y得x；三、 調(diào)試過程及實(shí)驗(yàn)結(jié)果：clear;clc;% 方程系數(shù) >>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')運(yùn)行結(jié)果：平方根法的解即為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改進(jìn)平方根法解得的解即為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('改進(jìn)平方根法解得的解即為');實(shí)驗(yàn)三、二次多項(xiàng)式擬合一、實(shí)驗(yàn)要求：用計(jì)算機(jī)語言編制利用QR分解求解線性最小二乘問題的通用子程序，用編寫的程序求解一個(gè)二次多項(xiàng)式，逾山仕使在殘向量的2范數(shù)最小的意義下擬合下面的數(shù)據(jù)t i -1-0.75-0.500.250.50.75y——i 1.000.81250.751.001.31251.752.3125二、算法描述及實(shí)驗(yàn)步驟:QR分解求解程序如下:1、求A的QR分解;2、 計(jì)算匕=Qtb；3、 求解上三角方程Rx=c1得x；三、調(diào)試過程及實(shí)驗(yàn)結(jié)果：>>t=[-1-0.75-0.500.250.50.75];>>y=[1.000.81250.751.001.31251.752.3125];>>plot(t,y,'r*');>>legend('實(shí)驗(yàn)數(shù)據(jù)(ti,yi)');>>xlabel('t'),ylabel('y');>>title('二次多項(xiàng)式擬合的數(shù)據(jù)點(diǎn)(ti,yi)的散點(diǎn)圖')；運(yùn)行后屏幕顯示數(shù)據(jù)的散點(diǎn)圖(略).(3)編寫下列MATLAB程序計(jì)算f(x)在(％,七)處的函數(shù)值，即輸入程序>>symsabc>>t=[-1-0.75-0.500.250.50.75];>>fi=a.*t.A2+b.*t+c%運(yùn)行后屏幕顯示關(guān)于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]編寫構(gòu)造殘向量2范數(shù)的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);運(yùn)行后屏幕顯示誤差平方和如下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達(dá)到最小，只需利用極值的必要條件亍=0，丁=0，—=0，得da db cc到關(guān)于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)運(yùn)行后屏幕顯示J分別對a,b,c的偏導(dǎo)數(shù)如下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)運(yùn)行后屏幕顯示擬合函數(shù)f及其系數(shù)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('實(shí)驗(yàn)數(shù)據(jù)(ti,yi)');>>xlabel('t'),ylabel('y');>>title('二次多項(xiàng)式擬合的數(shù)據(jù)點(diǎn)(ti,yi)的散點(diǎn)圖')；>>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,