数值分析实验作业gauss消去法的数值稳定性分析.docx
《数值分析实验作业gauss消去法的数值稳定性分析.docx》由会员分享，可在线阅读，更多相关《数值分析实验作业gauss消去法的数值稳定性分析.docx（11页珍藏版）》请在冰点文库上搜索。
数值分析实验作业gauss消去法的数值稳定性分析
实验3.1Gauss消去法的数值稳定性试验
令狐采学
实验目的:
观察和理解Gauss消元过程中出现小主元（即
很小）时引起的方程组解的数值不稳定性。
实验内容：
求解方程组
，其中
（1）
，
；
（2）
.
实验要求：
（1）计算矩阵的条件数，判断系数矩阵是良态的还是病态的。
（2）用Gauss列主元消去法求得L和U及解向量
.
（3）用不选主元的Gauss消去法求得
和
及解向量
.
（4）观察小主元并分析其对计算结果的影响.
程序如下:
计算矩阵条件数及Gauss列主元消去法:
formatlongeng
A1=[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-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;
ifA1（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）;
fori=1:
n
L1（i,i）=1;
end
L=L1
U=triu（A1,0）
forj=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）;
forj=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；
=[18.9882;3.3378;-34.747;-33.9865]
不选主元的Gauss消去法程序：
clear
formatlongeng
A1=[0.3e-1559.1431;5.291-6.130-12;11.2952;1211];
b1=[59.17;46.78;1;2];
n=4;
fork=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）;
fori=1:
n
L1（i,i）=1;
end
L=L1
U=triu（A1,0）
forj=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）;
forj=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
程序运行结果如下：
同理可得
对应的系数矩阵条件数及Gauss列主元消去法求解结果：
K2=8.994;
不选主元的Gauss消去法结果：
实验4.5三次样条插值函数的收敛性
问题提出：
多项式插值不一定收敛的，即插值的节点多，效果不一定就好。
对三次样条插值函数又如何呢？
理论上证明三次样条插值函数的收敛性是比较困难的，也超过了本课程的内容。
通过本实验可以验证这一理论结果.
实验内容：
请按一定的规则分别选择等距或者非等距的插值节点，并不断增加插值节点的个数。
考虑实验4.4中的函数或者选择其他感兴趣的函数，可用Matlab的函数“spline”作此函数的三次样条插值函数。
实验要求：
（1）随着节点个数的增加，比较被逼近函数和三次样条差值函数的误差变化情况。
分析所得结果并与拉格朗日插值多项式比较。
（2）三次样条插值函数的思想最早产生于工业部门。
作为工业迎合用的例子，考虑如下例子：
某汽车制造商根据三次样条差值函数设计车门曲线，其中一段的数据如下：
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的样条插值.
程序如下：
formatshort
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）;
holdon
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'）
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）
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文件为：
functiony=lagrange（x0,y0,x）;
n=length（x0）;m=length（x）;
fori=1:
m
z=x（i）;
s=0.0;
fork=1:
n
p=1.0;
forj=1:
n
ifj~=k
p=p*（z-x0（j））/（x0（k）-x0（j））;
end
end
s=p*y0（k）+s;
end
y（i）=s;
end
程序运行结果如下：
插值结果比较：
y
0.0385
0.0588
0.1
0.2
0.5
1
0.5
0.2
0.1
0.0588
0.0385
yy1
0.0385
-0.252
-0.0623
0.3687
0.8024
1
0.8024
0.3687
-0.0623
-0.252
0.0385
y1la
0.0385
-0.3793
-0.1101
0.4005
0.8342
1
0.8342
0.4005
-0.1101
-0.3793
0.0385
yy2
0.0385
0.0551
0.1085
0.1765
0.5342
1
0.5342
0.1765
0.1085
0.0551
0.0385
y2la
0.0385
-0.2587
0.3049
0.0726
0.5636
1
0.5636
0.0726
0.3049
-0.2587
0.0385
yy3
0.0385
0.0588
0.1
0.2
0.5
1
0.5
0.2
0.1
0.0588
0.0385
y3la
0.0385
0.0588
0.1
0.2
0.5
1
0.5
0.2
0.1
0.0588
0.0385
yy4
0.0385
0.0587
0.1
0.2016
0.5
1
0.5
0.1989
0.0993
0.0588
0.0385
y4la
0.0385
0.4312
0.1
0.2461
0.5
1
0.5
0.264
0.2387
0.0588
0.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.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
%cubicinterpolants
N=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）;
%Clampedsplineendpointconstraints
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））;
fork=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）;
fork=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;
fork=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
区间
三次样条差值
（3）实验4.4中的样条插值见上表中的y值.
升级会员 