《矩陣與數(shù)值分析》課程數(shù)值實驗報告

矩陣與數(shù)值分析實驗報告學(xué)院(系)：建設(shè)與工程學(xué)部專業(yè)：巖土工程學(xué)生姓名：陶言祺學(xué)號：21206012指導(dǎo)教師：孟兆良完成日期：2012-12-15大連理工大學(xué)Dalian

1.給定n階方程組，其中，則方程組有解。對和，分別用Gauss消去法和列主元消去法解方程組，并比較計算結(jié)果。Gauss消去法：Matlab編程（建立GS.m文件）：functionx=GS(n)A=[];b=[];fori=1:n-1A(i,i)=6;A(i,i+1)=1;A(i+1,i)=8;b(i)=15;endA(n,n)=6;b(1)=7;b(n)=14;b=b';fork=1:n-1fori=k+1:nm(i,k)=A(i,k)/A(k,k);A(i,k:n)=A(i,k:n)-m(i,k)*A(k,k:n);b(i)=b(i)-m(i,k)*b(k);endendb(n)=b(n)/A(n,n);fori=n-1:-1:1b(i)=(b(i)-sum(A(i,i+1:n).*b(i+1:n)'))/A(i,i);endclearx;x=b;disp('AX=b的解x是')end計算結(jié)果：在matlab命令框里輸出GS（10）得：>>GS(10)AX=b的解x是ans=1.00001.00001.00001.00001.00001.00001.00001.00001.00001.0000在matlab命令框里輸出GS（84）得：>>GS(84)AX=b的解x是ans=1.0e+008*0.0000………0.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00010.0002-0.00030.0007-0.00130.0026-0.00520.0105-0.02090.0419-0.08360.1665-0.33030.6501-1.25822.3487-4.02635.3684列主元消去法：Matlab編程（建立GLZX.m文件）：functionx=GLZX(n)A=[];b=[];eps=10^-2;fori=1:n-1A(i,i)=6;A(i,i+1)=1;A(i+1,i)=8;b(i)=15;endA(n,n)=6;b(1)=7;b(n)=14;b=b';fork=1:n-1[mainElement,index]=max(abs(A(k:n,k)));index=index+k-1;%indexifabs(mainElement)<epsdisp('列元素太?。?！');break;elseifindex>ktemp=A(k,:);temp1=b(k);A(k,:)=A(index,:);b(k)=b(index);A(index,:)=temp;b(index)=temp1;endfori=k+1:nm(i,k)=A(i,k)/A(k,k);%A(k,k);A(i,k:n)=A(i,k:n)-m(i,k)*A(k,k:n);b(i)=b(i)-m(i,k)*b(k);endendb(n)=b(n)/A(n,n);fori=n-1:-1:1b(i)=(b(i)-sum(A(i,i+1:n).*b(i+1:n)'))/A(i,i);endclearx;x=b;disp('AX=b的解x是')end計算結(jié)果：在matlab命令框里輸出GLZX（10）得：>>GLZX(10)AX=b的解x是ans=1111111111在matlab命令框里輸出GLZX（84）得：>>GLZX(84)AX=b的解x是ans=1.00001.00001.00001.00001.00001.0000………1.00001.00001.00001.00001.00001.00001.00001.0000分析：在n較小時，兩種方法均能得到正確解，當(dāng)n較大后，Gauss消去法計算結(jié)果嚴(yán)重偏離準(zhǔn)確值成為錯解，列主元消去法依然能得到正確解。這是因為Gauss消去法中有小主元做除數(shù)，在計算過程中的舍入誤差會對解產(chǎn)生極大影響，從而導(dǎo)致錯誤。列主元消去法則避免了這種情況。2.設(shè)有方程組，其中，選取不同的初始向量和不同的右端向量，給定迭代誤差要求，用Jacobi和Gauss-Seidel迭代法計算，觀測得出的迭代向量序列是否收斂。若收斂，記錄迭代次數(shù)，分析計算結(jié)果并得出你的結(jié)論。選定初始向量初始向量和不同的右端向量，如取。將的主對角線元素成倍增長若干次，非主對角元素不變，每次用Jacobi法計算，要求迭代誤差滿足，比較收斂速度，分析現(xiàn)象并得出你的結(jié)論。Matlab編程：Jacobi迭代法：functionx=Jacobi(b,x0)A=zeros(20);fori=1:18A(i,i)=3;A(i+1,i)=-0.5;A(i,i+1)=-0.5;A(i+2,i)=-0.25;A(i,i+2)=-0.25;endA(19,19)=3;A(20,20)=3;A(19,20)=-0.5;A(20,19)=-0.5;D=diag(diag(A));L=tril(A-D);U=triu(A-D);B=D\(L+U);f=D\b;x=x0;i=0;ep=10^-5;whilei<1000y=x;x=B*x+f;i=i+1;ifnorm(x-y,inf)<epbreak;endendifi<1000disp('迭代收斂且迭代次數(shù)為')ielsedisp('迭代不收斂')endendGauss-Seidel迭代法：functionx=G_S(b,x0)A=zeros(20);fori=1:18A(i,i)=3;A(i+1,i)=-0.5;A(i,i+1)=-0.5;A(i+2,i)=-0.25;A(i,i+2)=-0.25;endA(19,19)=3;A(20,20)=3;A(19,20)=-0.5;A(20,19)=-0.5;D=diag(diag(A));L=tril(A-D);U=triu(A-D);B=(D-L)\U;f=(D-L)\b;x=x0;i=0;ep=10^-5;whilei<1000y=x;x=B*x+f;i=i+1;ifnorm(x-y,inf)<epbreak;endendifi<1000disp('迭代收斂且迭代次數(shù)為')ielsedisp('迭代不收斂')endend計算結(jié)果：設(shè)誤差要求為1×10^-5任選了兩組不同的初始向量和不同的右端向量：b=[1,1,1,1,1,1,1,1,2,3,4,1,4,2,3,1,3,1,4,6]'x0=[2,3,1,4,4,5,5,3,1,3,1,1,1,1,1,1,1,1,1,1]'b=[1,4,5,12,42,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]'x0=[0,1,2,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1]'在matlab內(nèi)運行兩個程序得到：第一組：>>Jacobi(b,x0)迭代收斂且迭代次數(shù)為i=19>>G_S(b,x0)迭代收斂且迭代次數(shù)為i=9第二組：>>Jacobi(b,x0)迭代收斂且迭代次數(shù)為i=18>>G_S(b,x0)迭代收斂且迭代次數(shù)為i=9分析：兩組均收斂，明顯在均收斂的情況下Gauss-Seidel迭代法要比Jacobi迭代法收斂速度快。這是因為Gauss-Seidel迭代法對已經(jīng)計算出來的信息進行了充分利用。(b)Matlab編程：functionx=Jacobi2(n)A=zeros(20);x=diag(A);e=diag(eye(20));fori=1:18A(i,i)=3*n;A(i+1,i)=-0.5;A(i,i+1)=-0.5;A(i+2,i)=-0.25;A(i,i+2)=-0.25;endA(19,19)=3*n;A(20,20)=3*n;A(19,20)=-0.5;A(20,19)=-0.5;b=A*e;D=diag(diag(A));L=tril(A-D);U=triu(A-D);B=D\(L+U);f=D\b;i=0;ep=10^-6;whilei<1000y=x;x=B*x+f;i=i+1;ifnorm(x-y,inf)<epbreak;endendifi<1000disp('迭代收斂且迭代次數(shù)為')ielsedisp('迭代不收斂')endend計算結(jié)果：n代表主元增長倍數(shù)，選取不同的n得：>>Jacobi2(1)迭代收斂且迭代次數(shù)為i=20>>Jacobi2(2)迭代收斂且迭代次數(shù)為i=11>>Jacobi2(3)迭代收斂且迭代次數(shù)為i=9>>Jacobi2(4)迭代收斂且迭代次數(shù)為i=8>>Jacobi2(10)迭代收斂且迭代次數(shù)為i=6分析：隨著矩陣A主對角元增長倍數(shù)的變大，收斂速度變快，但是這種改善效果會逐漸減小。3.用迭代法求方程的全部根，要求誤差限為。首先經(jīng)過簡單計算并結(jié)合圖形得知該方程的三個單根區(qū)間為[-3,-1],[-1,0],[0,1],然后利用二分法求解。Matlab編程：functionx=f(a,b)fa=a^3+3*a^2-1;fb=b^3+3*b^2-1;iffa*fb>0disp('區(qū)間不是單根區(qū)間')return;endep=0.5*10^-8;whileabs(a-b)>epc=(a+b)/2;f1=a^3+3*a^2-1;f2=c^3+3*c^2-1;j=f1*f2;ifj==0x=c;break;elseifj<0b=c;elsea=c;endx=a;endenddisp('該區(qū)間內(nèi)根為')end計算結(jié)果：>>f(-3,-1)該區(qū)間內(nèi)根為ans=-2.879385244101286>>f(-1,0)該區(qū)間內(nèi)根為ans=-0.652703646570444>>f(0,1)該區(qū)間內(nèi)根為ans=0.5320888832211494.給定數(shù)據(jù)如下表：0123456789100.00.791.532.192.713.033.272.893.063.193.29編制程序求三次樣條插值函數(shù)在插值中點的樣條函數(shù)插值，并作點集和樣條插值函數(shù)的圖形，滿足的邊界條件為（1）（2）(1)Matlab編程：clearallx=[012345678910];f1=0.8;f2=0.2;y=[00.791.532.192.713.033.272.893.063.193.29];n=length(x);h=[];t=[];g=[];u=[];s=[];A=2*eye(n-2);yn=[];xn=[];fork=2:nhk=x(n)-x(n-1);h=[hhk];endfork=2:length(h)tk=h(k)/(h(k)+h(k-1));uk=h(k-1)/(h(k)+h(k-1));gk=3*(uk*(y(k+1)-y(k))/h(k)+tk*(y(k)-y(k-1))/h(k-1));g=[ggk];t=[t,tk];u=[uuk];endfork=2:n-2A(k-1,k)=u(k-1);A(k,k-1)=t(k);endg(1)=g(1)-t(1)*f1;g(n-2)=g(n-2)-u(n-2)*f2;m=A\g';m=m';m=[f1mf2];fork=1:n-1a1=2*y(k)/h(k)^3;a2=-2*y(k+1)/h(k)^3;a3=m(k)/h(k)^2;a4=m(k+1)/h(k)^2;b1=y(k)/h(k)^2;b2=y(k+1)/h(k)^2;c1=[1-x(k)];c2=[1-x(k+1)];c3=conv(c1,c1);c4=conv(c2,c2);sk1=b1*c4+b2*c3;sk2=(a1+a3)*conv(c1,c4)+(a2+a4)*conv(c2,c3);sk=[0sk1]+sk2;s=[s;sk];Pn=poly2str(sk,'x');Y=polyval(s(k,:),(x(k)+x(k+1))/2)endfork=1:n-1xi=[x(k):0.02:x(k+1)-0.02];yi=polyval(s(k,:),xi);yn=[ynyi];xn=[xnxi];endxn=[xnx(n)];yn=[ynpolyval(s(n-1,:),x(n))];plot(xn,yn,'x');title('3次樣條函數(shù)插值結(jié)果')計算結(jié)果：在Matlab命令框內(nèi)輸入上述程序得到插值中點的樣條函數(shù)值Y和樣條插值函數(shù)的圖形：Y=0.398564254606255Y=1.168428726968726Y=1.871470837518841Y=2.478187922956004Y=2.873277470657021Y=3.213702194414573Y=3.084413751685133Y=2.919892798844714Y=3.149765052935493Y=3.222296989411632(2)Matlab編程：clearallx=[012345678910];f1=0;f2=0;y=[00.791.532.192.713.033.272.893.063.193.29];n=length(x);h=[];t=[];g=[];u=[];s=[];A=2*eye(n);yn=[];xn=[];fork=2:nhk=x(n)-x(n-1);h=[hhk];endfork=2:length(h)tk=h(k)/(h(k)+h(k-1));uk=h(k-1)/(h(k)+h(k-1));gk=3*(uk*(y(k+1)-y(k))/h(k)+tk*(y(k)-y(k-1))/h(k-1));g=[ggk];t=[t,tk];u=[uuk];endfork=1:n-2A(k+1,k+2)=u(k);A(k+1,k)=t(k);endA(1,2)=1;A(n,n-1)=1;g0=3*(y(2)-y(1))/h(1)-h(1)*f1/2;gn=3*(y(n)-y(n-1))/h(n-2)+h(n-2)*f2/2;g=[g0,g,gn];m=A\g';m=m';fork=1:n-1a1=2*y(k)/h(k)^3;a2=-2*y(k+1)/h(k)^3;a3=m(k)/h(k)^2;a4=m(k+1)/h(k)^2;b1=y(k)/h(k)^2;b2=y(k+1)/h(k)^2;c1=[1-x(k)];c2=[1-x(k+1)];c3=conv(c1,c1);c4=conv(c2,c2);sk1=b1*c4+b2*c3;sk2=(a1+a3)*conv(c1,c4)+(a2+a4)*conv(c2,c3);sk=[0sk1]+sk2;s=[s;sk];Pn=poly2str(sk,'x');Y=polyval(s(k,:),(x(k)+x(k+1))/2)endfork=1:n-1xi=[x(k):0.02:x(k+1)-0.02];yi=polyval(s(k,:),xi);yn=[ynyi];xn=[xnxi];endxn=[xnx(n)];yn=[ynpolyval(s(n-1,:),x(n))];plot(xn,yn,'x');title('3次樣條函數(shù)插值結(jié)果')計算結(jié)果：在命令框內(nèi)輸入上述程序得到：Y=0.398428148708663Y=1.168465553874013Y=1.871459635795274Y=2.478195902944736Y=2.873256752425371Y=3.213777087353492Y=3.084134898160016Y=2.920933320005332Y=3.145881821815571Y=3.236789392730741再次輸入plot(x,y,'x');可以得到點集的圖形5.對下列數(shù)據(jù)作三次多項式擬合，取權(quán)系數(shù)，給出擬合多項式的系數(shù)、平方誤差，并作離散數(shù)據(jù)和擬合多項式的圖形。-1.0-0.50.00.51.01.52.0-4.447-0.4520.5510.048-0.4470.5494.552Matlab編程：clearallx=[-1,-0.5,0,0.5,1,1.5,2];y=[-4.447,-0.452,0.551,0.048,-0.447,0.549,4.552];a=polyfit(x,y,3)m=polyval(a,x)-y;e=norm(m,2);ep=e^2xi=[x(1):0.02:x(7)];yi=polyval(a,xi);plot(xi,yi,'x');計算結(jié)果：在Matlab命令框內(nèi)輸入上述程序得到：a=1.9991-2.9977-0.00000.5491ep=2.1762e-005數(shù)組a為擬合多項式的系數(shù)，ep為平方誤差。擬合多項式的圖形：再輸入plot(x,y,'x');得到：離散數(shù)據(jù)圖形：6.用復(fù)化梯形公式、復(fù)化Simpson公式和復(fù)化Gauss-Legendre公式計算下列積分的近似值，使絕對誤差限為，并將計算結(jié)果與精確解作比較以及比較各種算法的計算量。（1），（2）.復(fù)化梯形公式：Matlab編程：functions=tixing(f,a,b,t)ep=0.5*10^-7;n=2;s=0;whileabs(s-t)>eph=(b-a)/n;s1=0;fork=1:n-1x=a+h*k;s1=s1+feval(f,x);endn=n+1;s=h*(feval(f,a)+feval(f,b)+2*s1)/2;endn=n-1end計算結(jié)果：在matlab命令框內(nèi)輸入tixing(@(x)1/x,1,2,log(2))得到積分（1）的計算量（用分段數(shù)n刻畫）以及計算結(jié)果：n=1119ans=0.693147230473649在matlab命令框內(nèi)輸入tixing(@(x)1/(1+x^2),0,1,pi/4)得到積分（2）的計算量（用分段數(shù)n刻畫）以及計算結(jié)果：n=913ans=0.785398113411584復(fù)化Simpson公式：Matlab編程：functions=simpson(f,a,b,t)ep=0.5*10^-7;N=2;s=0;whileabs(s-t)>eph=(b-a)/(2*N);s1=0;s2=0;fork=1:Nx=a+h*(2*k-1);s1=s1+feval(f,x);endfork=1:(N-1)x=a+h*2*k;s2=s2+feval(f,x);endN=N+1;s=h*(feval(f,a)+feval(f,b)+4*s1+2*s2)/3;endN=N-1end計算結(jié)果：在matlab命令框內(nèi)輸入simpson(@(x)1/x,1,2,log(2))得到積分（1）的計算量（用分段數(shù)N刻畫）以及計算結(jié)果：（值得注意的是，simpson公式的分段數(shù)沒有刻畫出每段中點值的計算量）N=15ans=0.693147219033552在matlab命令框內(nèi)輸入simpson(@(x)1/(1+x^2),0,1,pi/4)得到積分（2）的計算量（用分段數(shù)N刻畫）以及計算結(jié)果：N=4ans=0.785398125614677復(fù)化Gauss-Legendre公式:Matlab編程：functiony=getLegendre(n,x)%返回Legendre多項式ifn==1y=x;elseifn==2y=1.5*x^2-0.5;elsey=(2-1/n)*x*getLegendre(n-1,x)-(n-1)/n*getLegendre(n-2,x);endendfunctionquad=GL(f,a,b,t)%復(fù)化Gauss-Legendre公式計算ep=0.5*10^-7;n=1;quad=0;symsx;whileabs(quad-t)>epg=getLegendre(n,x);temp=sym2poly(g);X=roots(temp);g=diff(g);g=g*g;N=length(X);A=zeros(N,1);whileN>0A(N)=2/(1-X(N)*X(N))/subs(g,'x',X(N));N=N-1;endN=length(X);T=zeros(N,1);M=zeros(1,N);T=((a+b)/2)+((b-a)/2)*X;fori=1:NM(i)=feval(f,T(i));endquad=((b-a)/2)*(M*A);n=n+1;endn=n-1end計算結(jié)果：在matlab命令框內(nèi)輸入GL(@(x)1/x,1,2,log(2))得到積分（1）的計算量（用Legendre多項式的次數(shù)n刻畫）以及計算結(jié)果：n=5ans=0.693147157853040在matlab命令框內(nèi)輸入GL(@(x)1/(1+x^2),0,1,pi/4)得到積分（2）的計算量（用Legendre多項式的次數(shù)n刻畫）以及計算結(jié)果：n=5ans=0.785398159971188分析：通過比較發(fā)現(xiàn)，精確性：復(fù)化Gauss-Legendre公式>復(fù)化Simpson公式>復(fù)化梯形公式，計算量：復(fù)化梯形公式>復(fù)化Simpson公式>復(fù)化Gauss-Legendre公式。7.用4階Runge-Kutta法求解微分方程：分別用和計算，并比較兩個近似值與精確解Matlab編程：functionR=rk4(f,a,b,ya,h)N=(b-a)/h;T=zeros(1,N+1);Y=zeros(1,N+1);y=zeros(1,N+1);T=a:h:b;Y(1)=ya;y(1)=ya;g=@(x)x^2*exp(-x)+cos(2*x);fork=1:Nk1=feval(f,T(k),Y(k));k2=feval(f,T(k)+h/2,Y(k)+h*k1/2);k3=feval(f,T(k)+h/2,Y(k)+h*k2/2);k4=feval(f,T(k)+h,Y(k)+h*k3);Y(k+1)=Y(k)+h*(k1+2*k2+2*k3+k4)/6;y(k+1)=feval(g,T(k));endR=[T'Y'y'];end計算結(jié)果：>>rk4(@(x,y)-y+cos(2*x)-2*sin(2*x)+2*x*exp(-x),0,2,1,0.1)ans=01.0000000000000001.0000000000000000.1000000000000000.9891149457346890.9891149520216010.2000000000000000.9538102071638440.9538102241260040.3000000000000000.8920092348744880.8920092547710330.4000000000000000.8039579098467340.8039579167128680.5000000000000000.6919349976895810.6919349707962980.6000000000000000.5599300266012580.5599299434705230.7000000000000000.4132941027862110.4132939417580310.8000000000000000.2583712723870450.2583710147337330.9000000000000000.1021196981119560.1021193296967981.000000000000000-0.048266907857711-0.0482673953757001.100000000000000-0.185726497573598-0.1857271059806701.200000000000000-0.303673326189410-0.3036740503876741.300000000000000-0.396309195058830-0.3963100231414661.400000000000000-0.458891377646222-0.4588922913431091.500000000000000-0.487948660815802-0.4879496362664781.600000000000000-0.481438680925772-0.4814396897284351.700000000000000-0.438841797588724-0.4388428080670581.800000000000000-0.361189039871900-0.3611900184962071.900000000000000-0.251024083624937-0.2510249965207042.000000000000000-0.112301673441313-0.112302487917161>>rk4(@(x,y)-y+cos(2*x)-2*sin(2*x)+2*x*exp(-x),0,2,1,0.05)ans=01.0000000000000001.0000000000000000.0500000000000000.9973822387440070.9973822388392770.1000000000000000.9891149517069640.9891149520216010.1500000000000000.9747024180466210.9747024185951700.2000000000000000.9538102234221290.9538102241260040.2500000000000000.9262576101293000.9262576108323360.3000000000000000.8920092542879130.8920092547710330.3500000000000000.8511664782800030.8511664782750310.4000000000000000.8039579175102620.8039579167128680.4500000000000000.7507296708922890.7507296689740730.5000000000000000.6919349741765920.6919349707962980.5500000000000000.6281234442518160.6281234390656750.6000000000000000.5599299507993260.5599299434705230.6500000000000000.4880631790988550.4880631693061170.7000000000000000.4132939543127370.4132939417580310.7500000000000000.3364434031691730.3364433875845240.8000000000000000.2583710335803170.2583710147337330.8500000000000000.1799628163369100.1799627940374300.9000000000000000.1021193555949370.1021193296967980.9500000000000000.0257442363982430.0257442068041841.0000