数值分析实验作业gauss消去法的数值稳定性分析说课讲解.docx
《数值分析实验作业gauss消去法的数值稳定性分析说课讲解.docx》由会员分享,可在线阅读,更多相关《数值分析实验作业gauss消去法的数值稳定性分析说课讲解.docx(13页珍藏版)》请在冰豆网上搜索。
数值分析实验作业gauss消去法的数值稳定性分析说课讲解
数值分析实验作业-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值.