矩阵与数值分析课程数值实验题目.docx
《矩阵与数值分析课程数值实验题目.docx》由会员分享，可在线阅读，更多相关《矩阵与数值分析课程数值实验题目.docx（20页珍藏版）》请在冰豆网上搜索。
矩阵与数值分析课程数值实验题目
《矩阵与数值分析》课程数值实验题目
一、设
，分别编制从小到大和从大到小的顺序计算
，
，
，并指出有效位数。
解：
程序代码如下：
function[A,B]=first（n）
formatlongg
T=1/2*（1+1/2-1/n-1/（n+1））;
k=0;
S=0;
fori=2:
n
S=S+1/（i^2-1）;
end
A=S;
x=abs（A-T）;
i=0;
while1
ifx<1/2*10^（-i）
i=i+1;
else
break
end
end
N1=k+i-1;
display（'从大到小:
'）,A,N1
S=0;
fori=n:
-1:
2
S=S+1/（i^2-1）;
end
B=S;
x=abs（B-T）;
i=0;
while1
ifx<1/2*10^（-i）
i=i+1;
else
break
end
end
N2=k+i-1;
display（'从小到大:
'）,B,N2
运行结果：
>>first（100）
从大到小:
A=0.740049504950495
N1=15
从小到大:
B=0.740049504950495
N2=15
>>first（10000）
从大到小:
:
A=0.749900004999506
N1=13
从小到大:
B=0.7499000049995
N2=323
>>first（1000000）
从大到小:
A=0.749999000000522
N1=13
从小到大:
B=0.7499990000005
N2=15
分析：
在做加法运算时，按照从小到大计算的顺序得到的结果要比按从大到小计算得到的结果有效数字位数更多。
二、解线性方程组
1．分别Jacobi迭代法和Gauss-Seidel迭代法求解线性方程组。
迭代法计算停止的条件为：
。
解：
（1）Jacobi迭代法程序代码：
clc;clear;
A=[-2100;1-210;01-21;001-2];
b=[-1000]';
N=100;
n=size（A,1）;
D=diag（diag（A））;
L=tril（-A,-1）;
U=triu（-A,1）;
Tj=inv（D）*（L+U）;
cj=inv（D）*b;
tol=1e-06;
k=1;
formatlong
x=zeros（n,1）;
whilek<=N
x（:
k+1）=Tj*x（:
k）+cj;
disp（k）;disp（'x='）;disp（x（:
k+1））;
ifnorm（x（:
k+1）-x（:
k））disp（'Theprocedurewassuccessful'）
disp（'Condition||x^（k+1）-x^（k）||disp（k）;disp（'x='）;disp（x（:
k+1））;
break
end
k=k+1;
end
结果输出
Theprocedurewassuccessful
Condition||x^（k+1）-x^（k）||60
x=
0.79999879906731
0.59999842795870
0.39999805685009
0.199********505
（2）Gauss-Seidel迭代法程序代码：
clc;clear;
A=[-2100;1-210;01-21;001-2];
b=[-1000]';
N=100;
n=size（A,1）;
D=diag（diag（A））;
L=tril（-A,-1）;
U=triu（-A,1）;
Tg=inv（D-L）*U;
cg=inv（D-L）*b;
tol=1e-06;
k=1;
x=zeros（n,1）;
whilek<=N
x（:
k+1）=Tg*x（:
k）+cg;
disp（k）;disp（'x='）;disp（x（:
k+1））;
ifnorm（x（:
k+1）-x（:
k））disp（'Theprocedurewassuccessful'）
disp（'Condition||x^（k+1）-x^（k）||disp（k）;disp（'x='）;disp（x（:
k+1））;
break
end
k=k+1;
end
结果输出
Theprocedurewassuccessful
Condition||x^（k+1）-x^（k）||31
x=
0.79999921397935
0.59999897108561
0.39999916759077
0.199********539
分析：
G—S法应用了上一步计算出来的较精确的值，所以可以用较少的迭代步数满足误差要求。
2.用Gauss列主元消去法、QR方法求解如下方程组：
解：
（1）Gauss列主元消去法
程序代码：
clc;clear;
formatlong
A=[1212;253-2;-2-235;1323];
b=[47-10]';
N=length（A）;
x=zeros（N,1）;
y=zeros（N,1）;
c=0;
d=0;
A（:
N+1）=b;
fork=1:
N-1
fori=k:
4
ifcd=i;c=abs（A（i,k））;
end
end
y=A（k,:
）;
A（k,:
）=A（d,:
）;
A（d,:
）=y;
fori=k+1:
N
c=A（i,k）;
forj=1:
N+1
A（i,j）=A（i,j）-A（k,j）*c/A（k,k）;
end
end
end
b=A（:
N+1）;
x（N）=b（N）/A（N,N）;
fork=N-1:
-1:
1
x（k）=（b（k）-A（k,k+1:
N）*x（k+1:
N））/A（k,k）;
end
结果输出
ans=
18.00000000000000
-9.57142857142857
6.00000000000000
-0.42857142857143
（2）QR方法：
程序代码
clc;clear;
formatlong
A=[1212;253-2;-2-235;1323];
b=[47-10]';
[Q,R]=qr（A）
y=Q'*b;
x=R\y
结果输出
Q=
-0.31622776601684-0.04116934847963-0.751646028002830.57735026918962
-0.63245553203368-0.49403218175557-0.150********056-0.57735026918963
0.63245553203368-0.74104827263336-0.22549380840085-0.00000000000000
-0.31622776601684-0.452862833275940.601316822402260.57735026918963
R=
-3.16227766016838-6.00832755431992-0.948683298050512.84604989415154
0-2.42899156029822-4.65213637819829-4.158********272
00-0.67648142520255-0.52615221960200
.0414********
x=
17.99999999999989
-9.57142857142851
5.99999999999997
-0.42857142857143
分析：
Gauss列主元消去法和QR方法结果是一致的。
三、非线性方程的迭代解法
1．用Newton迭代法求方程
的根，计算停止的条件为：
；
解：
Newton法程序代码：
clc;clear;
symsx
f=inline（'exp（x）+2^（-x）+2*cos（x）-6'）;
df=inline（'exp（x）-2^（-x）*log
（2）-2*sin（x）'）;
x0=-3;
delta=1e-06;
while1
x1=x0-feval（f,x0）/feval（df,x0）;
err=abs（x1-x0）;
x0=x1;
y=feval（f,x0）;
if（errend
x0
x0=0;
while1
x1=x0-feval（f,x0）/feval（df,x0）;
err=abs（x1-x0）;
x0=x1;
y=feval（f,x0）;
if（errend
x0
结果输出:
x0=-2.98650806938193
x0=1.82938360193385
分析：
Newton法收敛速度较快，但是需要一个好的初始点。
2．利用Newton迭代法求多项式
的所有实零点，注意重根的问题。
解：
由于不知道重根的个数，采用试探法求重根。
clc;clear;
symsx;
delta=1e-06;
f=inline（'x^4-5.4*x^3+10.56*x^2-8.954*x+2.7951'）;
df=inline（'4*x^3-16.2*x^2+21.12*x-8.954'）;
u=inline（'（x^4-5.4*x^3+10.56*x^2-8.954*x+2.7951）/（4*x^3-16.2*x^2+21.12*x-8.954）'）;
du=inline（'1-（x^4-27/5*x^3+264/25*x^2-4477/500*x+27951/10000）/（4*x^3-81/5*x^2+528/25*x-4477/500）^2*（12*x^2-162/5*x+528/25）'）;
j=1;
fori=0:
1:
3
x0=i;
while1
x1=x0-feval（u,x0）/feval（du,x0）;
err=abs（x1-x0）;
x0=x1;
y=feval（f,x0）;
if（errX（j）=x1;
j=j+1;
break
end
end
end
X
运行结果：
X=.0999********99823023722.100000000000011.10000116729456
分析：
该多项式的解有重根，Newton法的收敛速度较慢，此时为了提高收敛阶，采取试探法求解，效果良好。
若直接采用matlab内置函数求解：
solve（'x^4-5.4*x^3+10.56*x^2-8.954*x+2.7951=0'）
得到四个跟分别为：
2.1、1.1、1.1、1.1
四、数值积分
分别用复化梯形公式和Romberg公式计算积分
的近似值，要求误差不超过
，并给出节点个数。
解：
（1）复化梯形公式
程序代码：
Clc;clear;
f='1/x';a=2;b=8;eps=1.0e-5;
n=1;
h=（b-a）/2;
I1=2;
I2=（subs（sym（f）,findsym（sym（f））,a）+subs（sym（f）,findsym（sym（f））,b））/h;
whileabs（I2-I1）>eps
n=n+1;
h=（b-a）/n;
I1=I2;
I2=0;
fori=0:
n-1
x=a+h*i;
x1=x+h;
I2=I2+（h/2）*（subs（sym（f）,findsym（sym（f））,x）+...
subs（sym（f）,findsym（sym（f））,x1））;
end
end
I=I2
step=n
结果输出
I=
1.38654458753674
step=
53
分析：
该算法运算精度高，但是计算时间长、计算量大。
（2）Romberg公式
程序代码：
clc;clear;
f=inline（'1/x'）;
a=2;b=8;
TOL=1e-05;
A=zeros（20,20）;
A（1,1）=（b-a）*（feval（f,a）+feval（f,b））/2;
h=（b-a）/2;
A（2,1）=A（1,1）/2+h*feval（f,a+h）;
A（2,2）=（4*A（2,1）-A（1,1））/3;
errest=abs（A（2,2）-A（1,1）/2）;
i=2;
while（errest>TOL）
i=i+1;
h=h/2;
sum=0.0;
forj=1:
2:
2^（i-1）-1
sum=sum+feval（f,a+j*h）;
end;
A（i,1）=A（i-1,1）/2+h*sum;
forj=2:
i
power=4^（j-1）;
A（i,j）=（power*A（i,j-1）-A（i-1,j-1））/（power-1）;
end;
errest=abs（A（i,i）-A（i-1,i-1））/2^（i-1）;
end;
if（nargout==0）
s=sprintf（'\t\tapproximatevalueofintegral:
\t%.12f\n',A（i,i））;
s=sprintf（'%s\t\terrorestimate:
\t\t\t\t\t%.4e\n',s,errest）;
s=sprintf（'%s\t\tnumberoffunctionevaluations:
\t%d\n',...s,2^（i-1）+1）;
disp（s）
else
approx=A（i,i）;
end
运行结果：
approximatevalueofintegral:
1.386297441871
errorestimate:
8.7527e-006
numberoffunctionevaluations:
17
分析：
应用加速法可以使之前不收敛的式子收敛，有效的加快了收敛速度。
五、插值与逼近
1．给定
上的函数
，请做如下的插值逼近：
构造等距节点分别取
，
，
的Lagrange插值多项式；
构造分段线性取
的Lagrange插值多项式；
取Chebyshev多项式
的零点：
，
作插值节点构造
的插值多项式
和上述的插值多项式均要求画出曲线图形（用不同的线型或颜色表示不同的曲线）。
解：
程序代码
functionLagrange
clc;clear;closeall;
fori=1:
3;
ifi==1
N=5;
elseifi==2
N=8;
else
N=10;
end
f=inline（'1/（1+25*x^2）'）;
x1=zeros（1,N+1）;
a=-1;
b=1;
fori=1:
N+1
x1（i）=a+（i-1）*（b-a）/N;
y1（i）=feval（f,x1（i））;
end
symsx
ff=0;
fori=1:
N+1
f=1;
forj=1:
i-1
f=f*（x-x1（j））/（x1（i）-x1（j））;
end
forj=i+1:
N+1
f=f*（x-x1（j））/（x1（i）-x1（j））;
end
ff=f*y1（i）+ff;
f=1;
end
ff=collect（ff,x）;
ff=vpa（ff,4）;
y=ff;
p=ezplot（y,[a,b]）;grid
YLIM（[-0.10.6]）;
ifN==5
set（p,'Color','black'）;
set（p,'LineStyle','--'）;
lagrange_5=y
elseifN==8
set（p,'Color','r'）;
set（p,'LineStyle','--'）;
lagrange_8=y
else
set（p,'Color','b'）;
set（p,'LineStyle','--'）
lagrange_10=y
end
holdon;
xlabel（'x'）;ylabel（'y'）;
title（'y=p（x）'）;holdon;
end
Lag_Cheb（）;
x=-1:
0.01:
1;
y=1./（1+25*x.^2）;
acu=plot（x,y）;gridon;holdon
set（acu,'Color','m'）;
set（acu,'LineStyle','-'）;
legend（'N=5','N=8','N=10','Cheb,N=10','¾«È·Öµ'）;
%%%%%%%%%%%%%%%%%%%%
functionLag_Cheb（）
f=inline（'1/（1+25*x^2）'）;
N=10;
x1=zeros（1,N+1）;
a=-1;
b=1;
fori=1:
N+1
x1（i）=cos（（2*i-1）*pi/（2*N））;
y1（i）=feval（f,x1（i））;
end
symsx
ff=0;
fori=1:
N+1
f=1;
forj=1:
i-1
f=f*（x-x1（j））/（x1（i）-x1（j））;
end
forj=i+1:
N+1
f=f*（x-x1（j））/（x1（i）-x1（j））;
end
ff=f*y1（i）+ff;
f=1;
end
ff=collect（ff,x）;
ff=vpa（ff,4）;
yy=ff;
ff=collect（ff,x）;
yy=ff;
lagrange_chebshev_10=yy
cheb=ezplot（yy,[a,b]）;gridon
YLIM（[-0.10.6]）;
set（cheb,'Color','g'）;
set（cheb,'LineStyle',':
'）;
结果输出
lagrange_5=.5673+1.202*x^4-1.731*x^2
lagrange_8=1.+53.69*x^8-102.8*x^6+61.37*x^4-13.20*x^2
lagrange_10=1.-220.9*x^10+494.9*x^8-381.4*x^6+123.4*x^4-16.86*x^2
lagrange_chebshev_10=.7413-5.359*x^10-.4e-2*x^9+18.96*x^8-.1321*x^7-
25.78*x^6+.10*x^5+16.81*x^4+.5e-2*x^3-5.336*x^2+.5288e-3*x
分析：
（1）拉格朗日插值法随着插值节点增多而更加接近真实曲线，但是还是会在区间
［-1，-0.4］及［0.4，1］上，由于高次振荡产生容格现象。
（2）分段插值法则优于前者，因为分段插值在每个小区间上都可以最大限度的满足边界条件。
（3）利用契比雪夫迭代法可以更接近真实曲线。
2．已知函数值
0
1
2
3
4
5
6
7
8
9
10
2.51
3.30
4.04
4.70
5.22
5.54
5.78
5.40
5.57
5.70
5.80
和边界条件：
． 求三次样条插值函数
并画出其图形．
解：
程序实现代码：
functionSpline3_1（x,y,df0,dfn）
formatshort
x=[012345678910];
y=[2.513.304.044.705.225.545.785.405.575.705.80];
plot（x,y,'g*','MarkerSize',15）;holdon;
df0=1;
dfn=0;
n=length（x）;
h=zeros（n-1,1）;
lan=zeros（n-2,1）;
mu=zeros（n-2,1）;
g=zeros（n-2,1）;
m=zeros（n,1）;
m
（1）=df0;
m（n）=dfn;
fori=1:
n-1
h（i）=x（i+1）-x（i）;
end
fori=1:
n-2
mu（i）=h（i）/（h（i+1）+h（i））;
lan（i）=h（i+1）/（h（i+1）+h（i））;
g（i）=3*（mu（i）*（y（i+2）-y（i+1））/h（i+1）+lan（i）*（y（i+1）-y（i））/h（i））;
end
A=zeros（n-2,n-2）;
A（1,1）=2;
A（1,2）=mu
（1）;
A（n-2,n-2）=lan（n-2）;
fori=2:
n-2
A（i,i）=2;
A（i,i-1）=mu（i）;
A（i-1,i）=lan（i）;
end
g
（1）=g
（1）-lan
（1）*df0;
g（n-2）=g（n-2）-mu（n-2）*dfn;
b=A\g;
fori=2:
n-1
m（i）=b（i-1）;
end
symsz
fori=1:
n-1
xx=x（i）:
0.01:
x（i+1）;
sx1=y（i）*（xx-x（i+1））.^2.*（h（i）+2*（xx-x（i）））/h（i）^3;
sx2=y（i+1）*（xx-x（i））.^2.*（h（i）-2*（xx-x（i+1）））/h（i）^3;
sx3=m（i）*（xx-x（i+1））.^2.*（xx-x（i））/h（i）^2;
sx4=m（i+1）*（xx-x（i））.^2.*（xx-x（i+1））/h（i）^2;
sx=sx1+sx2+sx3+sx4;
z1=y（i）*（z-x（i+1））.^2.*（h（i）+2*（z-x（i）））/h（i）^3;
z2=y（i+1）*（z-x（i））.^2.*（h（i）-2
升级会员 