矩阵与数值分析课程数值实验题目Word格式文档下载.docx
《矩阵与数值分析课程数值实验题目Word格式文档下载.docx》由会员分享，可在线阅读，更多相关《矩阵与数值分析课程数值实验题目Word格式文档下载.docx（20页珍藏版）》请在冰豆网上搜索。
first（10000）
:
A=0.749900004999506
N1=13
B=0.7499000049995
N2=323
first（1000000）
A=0.749999000000522
B=0.7499990000005
分析：
在做加法运算时，按照从小到大计算的顺序得到的结果要比按从大到小计算得到的结果有效数字位数更多。
二、解线性方程组
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））<
tol
Theprocedurewassuccessful'
）
Condition||x^（k+1）-x^（k）||<
tolwasmetafterkiterations'
k=k+1;
结果输出
Theprocedurewassuccessful
tolwasmetafterkiterations
60
x=
0.79999879906731
0.59999842795870
0.39999805685009
0.199********505
（2）Gauss-Seidel迭代法程序代码：
Tg=inv（D-L）*U;
cg=inv（D-L）*b;
k+1）=Tg*x（:
k）+cg;
31
0.79999921397935
0.59999897108561
0.39999916759077
0.199********539
G—S法应用了上一步计算出来的较精确的值，所以可以用较少的迭代步数满足误差要求。
2.用Gauss列主元消去法、QR方法求解如下方程组：
（1）Gauss列主元消去法
程序代码：
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
ifc<
abs（A（i,k））
d=i;
c=abs（A（i,k））;
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
b=A（:
N+1）;
x（N）=b（N）/A（N,N）;
fork=N-1:
1
x（k）=（b（k）-A（k,k+1:
N）*x（k+1:
N））/A（k,k）;
ans=
18.00000000000000
-9.57142857142857
6.00000000000000
-0.42857142857143
（2）QR方法：
程序代码
[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
Gauss列主元消去法和QR方法结果是一致的。
三、非线性方程的迭代解法
1．用Newton迭代法求方程
的根，计算停止的条件为：
；
Newton法程序代码：
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（err<
delta）,break,end
x0
x0=0;
结果输出:
x0=-2.98650806938193
x0=1.82938360193385
Newton法收敛速度较快，但是需要一个好的初始点。
2．利用Newton迭代法求多项式
的所有实零点，注意重根的问题。
由于不知道重根的个数，采用试探法求重根。
symsx;
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'
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）;
delta）
X（j）=x1;
j=j+1;
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;
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））;
I=I2
step=n
I=
1.38654458753674
step=
53
该算法运算精度高，但是计算时间长、计算量大。
（2）Romberg公式
TOL=1e-05;
A=zeros（20,20）;
A（1,1）=（b-a）*（feval（f,a）+feval（f,b））/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）
h=h/2;
sum=0.0;
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）;
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\t\terrorestimate:
\t\t\t\t\t%.4e\n'
s,errest）;
%s\t\tnumberoffunctionevaluations:
\t%d\n'
...s,2^（i-1）+1）;
disp（s）
else
approx=A（i,i）;
approximatevalueofintegral:
1.386297441871
errorestimate:
8.7527e-006
numberoffunctionevaluations:
17
应用加速法可以使之前不收敛的式子收敛，有效的加快了收敛速度。
五、插值与逼近
1．给定
上的函数
，请做如下的插值逼近：
构造等距节点分别取
的Lagrange插值多项式；
构造分段线性取
取Chebyshev多项式
的零点：
作插值节点构造
的插值多项式
和上述的插值多项式均要求画出曲线图形（用不同的线型或颜色表示不同的曲线）。
functionLagrange
closeall;
fori=1:
3;
ifi==1
N=5;
elseifi==2
N=8;
else
N=10;
1/（1+25*x^2）'
x1=zeros（1,N+1）;
a=-1;
b=1;
N+1
x1（i）=a+（i-1）*（b-a）/N;
y1（i）=feval（f,x1（i））;
ff=0;
fori=1:
f=1;
i-1
f=f*（x-x1（j））/（x1（i）-x1（j））;
forj=i+1:
ff=f*y1（i）+ff;
f=1;
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'
LineStyle'
--'
lagrange_5=y
elseifN==8
r'
lagrange_8=y
b'
lagrange_10=y
holdon;
xlabel（'
x'
ylabel（'
y'
title（'
y=p（x）'
holdon;
Lag_Cheb（）;
x=-1:
0.01:
1;
y=1./（1+25*x.^2）;
acu=plot（x,y）;
gridon;
holdon
set（acu,'
m'
-'
legend（'
N=5'
N=8'
N=10'
Cheb,N=10'
¾
«
È
·
Ö
µ
%%%%%%%%%%%%%%%%%%%%
functionLag_Cheb（）
N=10;
x1（i）=cos（（2*i-1）*pi/（2*N））;
yy=ff;
yy=ff;
lagrange_chebshev_10=yy
cheb=ezplot（yy,[a,b]）;
gridon
YLIM（[-0.10.6]）;
set（cheb,'
g'
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．已知函数值
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）;
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;
n-1
h（i）=x（i+1）-x（i）;
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））;
A=zeros（n-2,n-2）;
A（1,1）=2;
A（1,2）=mu
（1）;
A（n-2,n-2）=lan（n-2）;
A（i,i）=2;
A（i,i-1）=mu（i）;
A（i-1,i）=lan（i）;
g
（1）=g
（1）-lan
（1）*df0;
g（n-2）=g（n-2）-mu（n-2）*dfn;
b=A\g;
m（i）=b（i-1）;
symsz
xx=x（i）:
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
升级会员 