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

实验目的:

观察和理解Gauss消元过程中出现小主元(即组解的数值不稳定性。 实验内容：

求解方程组Axb，其中

0.310-1559.1435.291-6.130-1（1）A111.29512171032.099999999999A25110（2）

159.1746.782； ，b112210185.90000000000162b2515021. ,

(k)akk很小)时引起的方程

实验要求：

（1） 计算矩阵的条件数，判断系数矩阵是良态的还是病态的。

4x,xR12（2） 用Gauss列主元消去法求得L和U及解向量.

~4~~~x,xRU12（3） 用不选主元的Gauss消去法求得L和及解向量.

（4） 观察小主元并分析其对计算结果的影响. 程序如下:计算矩阵条件数及Gauss列主元消去法: format longeng

A1=[0.3e-15 59.14 3 1;5.291 -6.130 -1 2;11.2 9 5 2;1 2 1 1];

b1=[59.17;46.78;1;2]; n=4;

k2=cond(A1) %k2为矩阵的条件数； for k=1:n-1

a=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; if A1(k,k)~=0

A1(k+1: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); else

break end end

L1=tril(A1,0); for i=1:n

L1(i,i)=1; end L=L1

U=triu(A1,0) for j=1:n-1

b1(j)=b1(j)/L(j,j);

b1(j+1:n)=b1(j+1:n)-b1(j)*L(j+1:n,j); end

b1(n)=b1(n)/L(n,n); for j=n:-1:2

b1(j)=b1(j)/U(j,j);

b1(1:j-1)=b1(1:j-1)-b1(j)*U(1:j-1,j); end

b1(1)=b1(1)/U(1,1); x1=b1

运行结果如下： K2=68.43；

126.791018L0.47240.08930100 0.1755100.02020.492910095211.2059.1431 U002.8351.2310000.801x1=[18.9882;3.3378;-34.747;-33.9865]

不选主元的Gauss消去法程序： clear

format longeng

A1=[0.3e-15 59.14 3 1;5.291 -6.130 -1 2;11.2 9 5 2;1 2 1 1]; b1=[59.17;46.78;1;2]; n=4;

for k=1:n-1

A1(k+1: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); end

L1=tril(A1,0); for i=1:n

L1(i,i)=1; end L=L1

U=triu(A1,0) for j=1:n-1

b1(j)=b1(j)/L(j,j);

b1(j+1:n)=b1(j+1:n)-b1(j)*L(j+1:n,j); end

b1(n)=b1(n)/L(n,n); for j=n:-1:2

b1(j)=b1(j)/U(j,j);

b1(1:j-1)=b1(1:j-1)-b1(j)*U(1:j-1,j); end

b1(1)=b1(1)/U(1,1); x1=b1

程序运行结果如下：

115~17.63710L37.331015153.3331000100 2.1168100.1890100.310150~U0059.141.043101800352.91101516017.6371015 80.51~x1[23.6848;1.0005;0;0]

同理可得A2对应的系数矩阵条件数及Gauss列主元消去法求解结果： K2=8.994;

010.51L-0.3-0.410-120.40000 100.333101107002.551.5 U0062.300003.36667x2[0.4441015;1.0;1.0;;1.0]

不选主元的Gauss消去法结果：

011~-0.3L0.5-2.49981012120-0.999910000 100.40012.3 125.7495103.3667107010126~01.010U0014.99871012000~x2[1.45105;1.000;0.99994;1.000145]

实验4.5 三次样条插值函数的收敛性

问题提出：

多项式插值不一定收敛的，即插值的节点多，效果不一定就好。对三次样条插值函数又如何呢？理论上证明三次样条插值函数的收敛性是比较困难的，也超过了本课程的内容。通过本实验可以验证这一理论结果.

实验内容：

请按一定的规则分别选择等距或者非等距的插值节点，并不断增加插值节点的个数。考虑实验4.4中的函数或者选择其他感兴趣的函数，可用Matlab的函数“spline”作此函数的三次样条插值函数。

实验要求：

（1） 随着节点个数的增加，比较被逼近函数和三次样条差值函数的

误差变化情况。分析所得结果并与拉格朗日插值多项式比较。 （2） 三次样条插值函数的思想最早产生于工业部门。作为工业迎合

用的例子，考虑如下例子：某汽车制造商根据三次样条差值函数设计车门曲线，其中一段的数据如下：

xkyk'yk 0 1 2 3 4 5 6 7 8 9 10 0.0 0.79 1.53 2.19 2.71 3.03 3.27 2.89 3.06 3.19 3.29 0.8 0.2 (3)计算实验4.4的样条插值. 程序如下： format short

x1=-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);

hold on

fplot('1./(1+25.*x.*x)',[-1,1],'m') plot(xx,yy1,'g') plot(xx,yy2,'b') plot(xx,yy3,'k') plot(xx,yy4,'r')

hold off %比较被逼近函数与三次样条插值函数图像，直观表现不同插值节点处误差的变化 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) y2la=lagrange(x2,y2,xx) yy3=spline(x3,y3,xx) y3la=lagrange(x3,y3,xx) yy4=spline(x4,y4,xx) y4la=lagrange(x4,y4,xx)

其中lagrange函数对应的m文件为： function y=lagrange(x0,y0,x); n=length(x0);m=length(x); for i=1:m z=x(i); s=0.0; for k=1:n p=1.0; for j=1:n if j~=k

p=p*(z-x0(j))/(x0(k)-x0(j)); end end

s=p*y0(k)+s; end y(i)=s; end

程序运行结果如下：

插值结果比较： 0.0588 0.1 0.2 0.5 1 0.5 0.2 0.1 0.0588 -0.252 -0.0623 0.3687 0.8024 1 0.8024 0.3687 -0.0623 -0.252 -0.3793 -0.1101 0.4005 0.8342 1 0.8342 0.4005 -0.1101 -0.3793 0.0551 0.1085 0.1765 0.5342 1 0.5342 0.1765 0.1085 0.0551 -0.2587 0.3049 0.0726 0.5636 1 0.5636 0.0726 0.3049 -0.2587 0.0588 0.1 0.2 0.5 1 0.5 0.2 0.1 0.0588 0.0588 0.1 0.2 0.5 1 0.5 0.2 0.1 0.0588 0.0587 0.1 0.2016 0.5 1 0.5 0.1989 0.0993 0.0588 0.4312 0.1 0.2461 0.5 1 0.5 0.264 0.2387 0.0588 车门曲线求解程序: 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.2;

S=csfit(x,y,dx0,dx10)

其中csfit函数的m文件为: function S=csfit(X,Y,dx0,dxn) %Clamped Cubic Spline

%Input -X is the 1xn abscissa vector % -Y is the 1xn ordinate vector

% -dx0=S'(x0) first derivative boundary condition % -dxn=S'(xn) first derivative boundary condition

%Output-S: rows of S are the coefficients, in descending order, for the % cubic interpolants N=length(X)-1; H=diff(X); D=diff(Y)./H;

y yy1 y1la yy2 y2la yy3 y3la yy4 y4la 0.0385 0.0385 0.0385 0.0385 0.0385 0.0385 0.0385 0.0385 0.0385 0.0385 0.0385 0.0385 0.0385 0.0385 0.0385 0.0385 0.0385 0.0385 A=H(2:N-1);

B=2*(H(1:N-1)+H(2:N)); C=H(2:N); U=6*diff(D);

%Clamped spline endpoint constraints B(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)); for k=2:N-1

temp=A(k-1)/B(k-1); B(k)=B(k)-temp*C(k-1); U(k)=U(k)-temp*U(k-1); end

M(N)=U(N-1)/B(N-1); for k=N-2:-1:1

M(k+1)=(U(k)-C(k)*M(k+2))/B(k); end

M(1)=3*(D(1)-dx0)/H(1)-M(2)/2; M(N+1)=3*(dxn-D(N))/H(N)-M(N)/2; for k=0:N-1

S(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); end end

程序运行结果为: S= -0.0085 -0.0015 0.8 0 -0.0045 -0.027 0.7715 0.79 -0.0037 -0.0404 0.7041 1.53 -0.0409 -0.0514 0.6123 2.19 0.1074 -0.1741 0.3868 2.71 -0.2685 0.1479 0.3606 3.03 0.4266 -0.6575 -0.1491 3.27 -0.2679 0.6222 -0.1844 2.89 0.0549 -0.1814 0.2565 3.06 0.0584 -0.0168 0.0584 3.19 区间 三次样条差值 0x1 S1(x)-0.0085x30.8x20.8x 1x2 S2(x)0.0045x30.027x20.7715x0.79 S3(x)0.0037x30.0404x20.7041x1.53 2x3 3x4 4x5 5x6 S4(x)0.0409x30.0514x20.6123x2.19 S5(x)0.1074x30.1741x20.3868x2.71 S6(x)0.2685x30.1479x20.3606x3.03 S7(x)0.4266x30.6575x20.1491x3.27 6x7 7x8 8x9 S8(x)0.2679x30.6222x20.1844x2.89 S9(x)0.0549x30.1814x20.2565x3.06 S9(x)0.0584x30.0168x20.0584x3.19 9x10

(3)实验4.4中的样条插值见上表中的y值.