数值分析实验作业-gauss消去法的数值稳定性分析

1、实验3.1Gauss消去法的数值稳定性试验实验目的：Mk)观察和理解Gauss消元过程中出现小主元(即akk很小)时引起的方程组解的数值不稳定性.实验内容:求解方程组Ax=b,其中(1)Ai=0.310-1515.29111.2I159.14-6.130923-1511221bi一59.17一46.78I12a2=(2)10-35J-72.099999999999-11065012-'12b2=一815.9000000000015I1实验要求：(1)计算矩阵的条件数,判断系数矩阵是良态的还是病态的.4(2)用Gauss列主元消去法求得L和U及解向量x1,x2R.4(3)用不选主兀的Ga

2、uss消去法求得L和U及解向量x1,x2=R.(4)观察小主元并分析其对计算结果的影响.程序如下：计算矩阵条件数及Gauss列主元消去法：formatlongengA1=0.3e-1559.1431;5.291-6.130-12;11.2952;1211;b1=59.17;46.78;1;2;n=4;k2=cond(A1)%k2为矩阵的条件数；fork=1:n-1a=max(abs(A1(k:n,k);p,k=find(A1=a);B=A1(k,:);c=b1(k);A1(k,:)=A1(p,:);b1(k)=b1(p);A1(p,:)=B;b1(p)=c;ifA1(k,k)=0A1(k+1:

3、n,k)=A1(k+1:n,k)/A1(k,k);A1(k+1:n,k+1:n)=A1(k+1:n,k+1:n)-A1(k+1:n,k)*A1(k,k+1:n);elsebreakendendL1=tril(A1,0);fori=1:nL1(i,i)=1;endL=L1U=triu(A1,0)forj=1:n-1b1(j)=b1(j)/L(j,j);b1(j+1:n)=b1(j+1:n)-b1(j)*L(j+1:n,j);endb1(n)=b1(n)/L(n,n);forj=n:-1:2b1(j)=b1(j)/U(j,j);b1(1:j-1)=b1(1:j-1)-b1(j)*U(1:j-1,j

4、);endb1(1)=b1(1)/U(1,1);x1=b1运行结果如下：K2=68.43;10004826.79父10100L=0.4724-0.175510一0.0893-0.0202-0.49291_11.2952059.1431U=00-2.8351.2310000.801-Xi=18.9882;3.3378;-34.747;-33.9865不选主元的Gauss消去法程序：clearformatlongengA1=0.3e-1559.1431;5.291-6.130-12;11,2952;1211;b1=59.17;46.78;1;2;n=4;fork=1:n-1A1(k+1:n,k尸A

5、1(k+1:n,k)/A1(k,k);A1(k+1:n,k+1:n)=A1(k+1:n,k+1:n)-A1(k+1:n,k)*A1(k,k+1:n);endL1=tril(A1,0);fori=1:nL1(i,i)=1;endL=L1U=triu(A1,0)forj=1:n-1b1(j)=b1(j)/L(j,j);b1(j+1:n)=b1(j+1:n)-b1(j)*L(j+1:n,j);endb1(n)=b1(n)/L(n,n);forj=n:-1:2b1(j)=b1(j)/U(j,j);b1(1:j-1)=b1(1:j-1)-b1(j)*U(1:j-1,j);endb1(1)=b1(1)/U

6、(1,1);x1=b1程序运行结果如下:一1000一1517.63710100L=1537.33102.11681011-17.637X101580.5_15_-3.333父100.18901_O.3X10-1559.1430-1.043X1018-52.91X1015U=00-16,0001=23.6848;1.00040;0同理可得A2对应的系数矩阵条件数及Gauss列主元消去法求解结果:K2=8.994;10000.5100L=-12-0.3-0.4x1010/0.4-0.3331_10-70102.55-1.5U=0062.300003.36667一x2=0.4441015;-1.0;

7、1.0;1.0不选主元的Gauss消去法结果：10000.3100L120.5-2.4998M1010._12_0-0.9999M100.400110-701二2-_一-1.01062.312-12014.9987105.749510_0003.36675_X2=-1.4510;-1.000;0.99994;1.000145实验4.5三次样条插值函数的收敛性问题提由：多项式插值不一定收敛的,即插值的节点多,效果不一定就好.对三次样条插值函数又如何呢理论上证实三次样条插值函数的收敛性是比拟困难的,也超过了本课程的内容.通过本实验可以验证这一理论结果.实验内容：请按一定的规那么分别选择等距或者非等

8、距的插值节点,并不断增加插值节点的个数.考虑实验4.4中的函数或者选择其他感兴趣的函数,可用Matlab的函数“spline作此函数的三次样条插值函数.实验要求：(1)随着节点个数的增加,比拟被逼近函数和三次样条差值函数的误差变化情况.分析所得结果并与拉格朗日插值多项式比拟.(2)三次样条插值函数的思想最早产生于工业部门.作为工业迎合用的例子,考虑如下例子：某汽车制造商根据三次样条差值函数设计车门曲线,其中一段的数据如下：xk012345678910yk0.00.791.532.192.713.033.272.893.063.193.29'yk0.80.2(3)计算实验4.4的样条插值

9、.程序如下：formatshortx1=-1:0.5:1;y1=1./(1+25.*x1.*x1);x2=-1:0.25:1;y2=1./(1+25.*x2.*x2);x3=-1:0.1:1;y3=1./(1+25.*x3.*x3);x4=-1,-0.82,-0.6,-0.53,-0.34,-0.2,0,0.04,0.2,0.25,0.5,0.8,1;y4=1./(1+25.*x4.*x4);xx=-1:0.01:1;yy1=spline(x1,y1,xx);yy2=spline(x2,y2,xx);yy3=spline(x3,y3,xx);yy4=spline(x4,y4,xx);holdo

10、nfplot('1./(1+25.*x.*x)',-1,1,'m')plot(xx,yy1,'g')plot(xx,yy2,'b')plot(xx,yy3,'k')plot(xx,yy4,'r')直观表现不同插值节点处holdoff%比拟被逼近函数与三次样条插值函数图像,误差的变化xx=-1:0.2:1;y=1./(1+25.*xx.*xx)%函数在相应节点处的真实值;yy1=spline(x1,y1,xx)y1la=lagrange(x1,y1,xx)yy2=spline(x2,y2,xx)y21

11、a=lagrange(x2,y2,xx)yy3=spline(x3,y3,xx)y31a=lagrange(x3,y3,xx)yy4=spline(x4,y4,xx)y41a=lagrange(x4,y4,xx)其中lagrange函数对应的m文件为：functiony=lagrange(x0,y0,x);n=length(x0);m=length(x);fori=1:mz=x(i);s=0.0;fork=1:np=1.0;forj=1:nifj=kp=p*(z-x0(j)/(x0(k)-x0(j);endends=p*y0(k)+s;endy(i)=s;end程序运行结果如下：插值结果比拟:

12、y0.03850.05880.10.20.510.50.20.10.05880.0385yy10.0385-0.252-0.06230.36870.802410.80240.3687-0.0623-0.2520.03851y1la0.0385-0.3793-0.11010.40050.834210.83420.4005-0.1101-0.37930.0385yy20.03850.05510.10850.17650.534210.53420.17650.10850.05510.0385y2la0.0385-0.25870.30490.07260.563610.56360.07260.3049-0

13、.25870.03851yy30.03850.05880.10.20.510.50.20.10.05880.0385y3la0.03850.05880.10.20.510.50.20.10.05880.0385yy40.03850.05870.10.20210.510.50.19890.09930.05880.03851y41a0.03850.43120.10.24610.510.50.2640.23870.05880.0385车门曲线求解程序:x=0:10;y=0,0.79,1.53,2.19,2.71,3.03,3.27,2.89,3.06,3.19,3.29;dx0=0.8;dx10=0

14、.2;S=csfit(x,y,dx0,dx10)其中csfit函数的m文件为：functionS=csfit(X,Y,dx0,dxn)%ClampedCubicSpline%Input-Xisthe1xnabscissavector%-Yisthe1xnordinatevector%-dx0=S'(x0)firstderivativeboundarycondition%-dxn=S'(xn)firstderivativeboundarycondition%Output-S:rowsofSarethecoefficients,indescendingorder,forthe%cu

15、bicinterpolantsN=length(X)-1;H=diff(X);D=diff(Y)./H;A=H(2:N-1);B=2*(H(1:N-1)+H(2:N);C=H(2:N);U=6*diff(D);%ClampedsplineendpointconstraintsB(1)=B(1)-H(1)/2;U(1)=U(1)-3*(D(1)-dx0);B(N-1)=B(N-1)-H(N)/2;U(N-1)=U(N-1)-3*(dxn-D(N);fork=2:N-1temp=A(k-1)/B(k-1);B(k)=B(k)-temp*C(k-1);U(k)=U(k)-temp*U(k-1);en

16、dM(N)=U(N-1)/B(N-1);fork=N-2:-1:1M(k+1)=(U(k)-C(k)*M(k+2)/B(k);endM(1)=3*(D(1)-dx0)/H(1)-M(2)/2;M(N+1)=3*(dxn-D(N)/H(N)-M(N)/2;fork=0:N-1S(k+1,1)=(M(k+2)-M(k+1)/(6*H(k+1);S(k+1,2)=M(k+1)/2;S(k+1,3)=D(k+1)-H(k+1)*(2*M(k+1)+M(k+2)/6;S(k+1,4)=Y(k+1);endend程序运行结果为：S=-0.0085-0.00150.80-0.0045-0.0270.7715

17、0.79-0.0037-0.04040.70411.53-0.0409-0.05140.61232.190.1074-0.17410.38682.71-0.26850.14790.36063.030.4266-0.6575-0.14913.27-0.26790.6222-0.18442.890.0549-0.18140.25653.060.0584-0.01680.05843.19区间三次样条差值0Ex<1Si(x)=-0.0085x3+0.8x2+0.8x1<x<2S2(x)=0.0045x3-0.027x2+0.7715x+0.792<x<3S3(x)=0.0037x30.0404x2+0.7041x+1.533Mx<4S4(x)=-0.0409x3-0.0514x2+0.6123x+2.194<x<5S5(x)=0.1074x3-0.1741x2+0.3868x+