数值分析上机试题对应参考答案.docx
《数值分析上机试题对应参考答案.docx》由会员分享，可在线阅读，更多相关《数值分析上机试题对应参考答案.docx（24页珍藏版）》请在冰豆网上搜索。
数值分析上机试题对应参考答案
一、问答题
1、什么是近似值x*有效数字？
若近似值x*的误差限是某一位的半个单位，该位到x*的第一位非零数字共有n位，就说x*有n位有效数字。
它可表示为
X=±10m×（a1+a2×10-1+…+an×10-（n-1）,其中ai（i=1,2,…,n）是0到9中的一个数字，a1≠0,m为整数，且︱x－x*︱≠
×10m-n+1
2、数值计算应该避免采用不稳定的算法，防止有效数字的损失.因此，在进行
数值运算算法设计过程中主要注意什么？
（1）简化计算过程，减少运算次数；
（2）避免两个相近的数相减；
（3）避免除数的绝对值远小于被除数的绝对值；
（4）防止大数“吃掉”小数的现象；
（5）使用数值稳定的算法，设法控制误差的传播。
3、写出“n阶阵A具有n个不相等的特征值”的等价条件（至少写3个）
（1）|A|不为零
（2）n阶矩阵A的列或行向量组线性无关
（3）矩阵A为满秩矩阵
（4）n阶矩阵A与n阶可逆矩阵B等价
4、迭代法的基本思想是什么？
就是用某种极限过程去逐步逼近线性方程组精确解得方法。
其基本思想为：
先任取一组近似解初值X0，然后按照某种迭代原则，由X0计算新的近似解X1，以此类推，可计算出X2,X3，…XK，。
。
。
如果｛X｝收敛，则取为原方程组的解。
5、病态线性方程组的主要判断方法有哪些？
（1）系数矩阵的某两行（列）几乎近似相关
（2）系数矩阵的行列式的值很小
（3）用主元消去法解线性方程组时出现小主元
（4）近似解x*已使残差向量r=b-Ax*的范数很小，但该近似解仍不符合问题要求。
6、Lagrange插值的前提条件是什么？
并写出二次Lagrange插值的基函数。
前提条件是：
二次Lagrange插值的基函数：
7、什么是数值积分的代数精度？
如果某一个求积公式对于次数不超过m的多项式均能准确地成立，但对于m+1次多项式就不准确成立，则称该求积公式具有m次代数精度（或代数精确度）。
8、写出解线性方程组的全主元消去法、列主元消去法和标度化列主元消去法各
自的优缺点。
一般来说列主元法能保证算法的稳定，所谓算法的稳定是指在运算过程中计算误差（对消去法这种直接法来说主要指由于计算机字长有限带来的舍入误差）能得到控制；全主元是较列主元更稳定的算法，但是它的计算量要比列主元算法大得多，列主元算法每做一次消元仅与同列的元素做比较，比较的次数与线性方程组的阶n是同阶的量，而全主元每做一次要与系数矩阵所有元素进行比较，计算量是与n2同阶的量。
9、求矩阵特征值的QR方法的算法步骤是什么？
对于给定的n阶实对称矩阵A与迭代次数M；
（1）令A1=A，对于k=1,2，…，M；
（2）迭代计算下一个矩阵：
Ak=QkRk（对Ak作QR分解），
（3）Ak=RQ（交换乘法次序），令k=k+1,Ak+1=QkTAKQK
（4）返回到2，直到k=M,输出A的主对角元素。
10、什么是高斯型求积公式与高斯点？
一般可设求积公式=In（ƒ），为具有n个节点的插值求积公式，且具有2n+1次最高代数精度，则称其求积点x0,x1,…,xn称为高斯点具有，相应的公式称为高斯型求积公式。
11、什么是插值计算的龙格现象？
增加插值节点，提高插值多项式的次数，可以使插值函数在更多的点与所逼近的函数取相同的值，但会使插值函数在两端发生激烈的振荡，这就是插值计算的龙格现象。
12、龙格—库塔方法的基本思想是什么？
基本思想是，在每一步内多预报几个点的斜率值，将其加权平均作为平均斜率，从而构造出更高的计算格式。
具体做法是，用函数f（x,y）在若干点上的函数值的线性组合来构造近似公式，构造时让近似公式在（xi,yi）处的泰勒展开式与解y（x）在xi处的泰勒展开式的前面几项重合，从而使近似公式达到所需要的阶数，这样既避免求偏导数，又提高了计算方法的精度。
二、计算题
1、随机生成4*4的矩阵，做列主元的三角分解，并验证等式PA=LU。
解：
>>A=rand（4）%随机生成的四位矩阵
A=
0.81470.63240.95750.9572
0.90580.09750.96490.4854
0.12700.27850.15760.8003
0.91340.54690.97060.1419
>>[L,U,P]=lu（A）
L=
1.0000000
0.99171.000000
0.8920-0.32501.00000
0.1390-0.45520.25671.0000
U=
0.91340.54690.97060.1419
0-0.44480.00240.3447
000.09250.9426
0000.6955
P=
0001
0100
1000
0010
>>P*A%验证PA=LU
ans=
0.91340.54690.97060.1419
0.90580.09750.96490.4854
0.81470.63240.95750.9572
0.12700.27850.15760.8003
>>L*U
ans=
0.91340.54690.97060.1419
0.90580.09750.96490.4854
0.81470.63240.95750.9572
0.12700.27850.15760.8003
2、用Lagrange插值和Newton插值拟合[0,2*pi]上的sin函数，并画图比较。
解：
建立拉格朗日函插值函数
function[yt,L]=LagInterp1（x,y,xt）
%拉格朗日差值
%x，y：
差值条件
%xt：
用拉格朗日插值函数要计算的自变量，可以是多个
%yt：
用拉格朗日插值函数计算出xt对应的函数值数组
%L：
拉格朗日插值多项式表达式
symst;
n=length（x）;
ny=length（y）;
ifn~=ny
error（'差值节点x与函数值y不一致'）;
end
L=0.0;
fork=1:
n
lk=1;
forj=1:
n
ifj~=k
lk=lk*（t-x（j））/（x（k）-x（j））;
end
end;
L=L+y（k）*lk;
end
simplify（L）;%简化拉格朗日插值多项式表达式
L=collect（L）;%将拉格朗日插值多项式展开
yt=subs（L,'t',xt）;%计算插值点处的函数值
建立牛顿函数：
function[yt,N]=NewtInterp（x,y,xt）
%已知数据点的牛顿插值
%x，y：
差值条件
%xt：
要计算的插值点，可以是多个
%yt：
用牛顿插值函数算出xt对应的函数值数组
%L：
牛顿插值多项式表达式
symst;
n=length（x）;
ny=length（y）;
ifn~=ny
error（'差值节点x与函数值y不一致'）;
end
a=zeros（1,n）;
N=y
（1）;
w=1;
fork=1:
n-1
yy=zeros（1,n）;%记录差商
forj=k+1:
n
yy（j）=（y（j）-y（k））/（x（j）-x（k））;
end
a（k）=yy（k+1）;
w=w*（t-x（k））;
N=N+a（k）*w;
y=yy;
end
yt=subs（N,'t',xt）;
simplify（N）;
N=collect（N）;%将插值多项式展开
N=vpa（N,6）;%系数转化为6位精度
命令：
>>x=0:
pi/10:
2*pi;
>>y=sin（x）;
>>z=0:
pi/20:
2*pi;
>>y1=LagInterp1（x,y,z）;%拉格朗日拟合法
>>y2=NewtInterp（x,y,z）;%牛顿拟合法
>>figure;
>>plot（z,sin（z）,'--r',z,y1,'-b',z,y2,'-.k'）%绘制函数图象
>>holdon
>>plot（x,y,'+'）
>>xlabel（'x'）
>>ylabel（'y'）
3、对不同初值用Jacobi迭代法解方程组Axb，其中
解：
建立雅可比函数
functiontx=jacobi（A,b,imax,x0,acc）
%利用Jacobi迭代法解线性方程组AX=b，迭代初值为x0，迭代次数为imax，精度为acc
%利用Jacobi迭代法解线性方程组AX=b，
%迭代初值为x0，
%迭代次数为imax，
%精度为acc
del=10^（-10）;%主对角元素不能太小
tx=[x0];
n=length（x0）;
fori=1:
n
dg=A（i,i）;
ifabs（dg）disp（'主对角元素太小'）;
return
end
end
fork=1:
imax
fori=1:
n
sm=b（i）;
forj=1:
n
ifj~=i
sm=sm-A（i,j）*x0（j）;
end
end
x（i）=sm/A（i,i）;
end
tx=[tx;x];
ifnorm（x-x0）return
else
x0=x;
end
end
命令：
>>A=[2,-1,0,0,0,0;-1,2,-1,0,0,0;0,-1,2,-1,0,0;0,0,-1,2,-1,0;0,0,0,-1,2,-1,;0,0,0,0,-1,2];
>>b=[1,0,1,0,0,1];
>>acc=1.0*10^（-6）;
>>imax=20;（10,20,30都可取）
>>x0=zeros（1,6）;
>>tx=jacobi（A,b,imax,x0,acc）（此命令一行一行的复制）
tx=
000000
0.500000.5000000.5000
0.50000.50000.50000.25000.25000.5000
0.75000.50000.87500.37500.37500.6250
0.75000.81250.93750.62500.50000.6875
0.90630.84381.21880.71880.65630.7500
0.92191.06251.28130.93750.73440.8281
1.03131.10161.50001.00780.88280.8672
1.05081.26561.55471.19140.93750.9414
1.13281.30271.72851.24611.06640.9688
1.15141.43071.77441.39751.10741.0332
1.21531.46291.91411.44091.21531.0537
1.23141.56471.95191.56471.24731.1077
1.28231.59172.06471.59961.33621.1237
1.29581.67352.09561.70041.36161.1681
1.33681.69572.18701.72861.43431.1808
1.34791.76192.21221.81061.45471.2171
1.38091.78002.28621.83351.51391.2274
1.39001.83362.30671.90011.53041.2569
1.41681.84842.36681.91861.57851.2652
1.42421.89182.38351.97271.59191.2893
4、用最小二乘法求一个多项式，使它与下列数据相拟合
x
-1.0
-0.5
0.0
0.5
1.0
1.5
2.0
y
-4.447
-0.452
0.551
0.048
-0.447
0.549
4.552
（1）求拟合曲线y=Pn（x）,（n=1,2,3）
（2）求拟合误差
（3）画出y=Pn（x）的图像
解：
（1）命令
>>x=[-1.0-0.50.00.51.01.52.0];
>>y=[-4.447-0.4520.5510.048-0.4470.5494.552];
>>p1=polyfit（x,y,1）%一次多项式拟合
p1=
2.0001-0.9495
>>p2=polyfit（x,y,2）%二次多项式拟合
p2=
0.00101.9991-0.9502
>>p3=polyfit（x,y,3）%三次多项式拟合
p3=
1.9991-2.9977-0.00000.5491
>>z=-1.0:
0.1:
2.0;
>>y1=polyval（p1,z）;
>>y2=polyval（p2,z）;
>>y3=polyval（p3,z）;
即拟合的多项式为：
y1=28001/14000*x-5317/5600
y2=1152921504607243/1152921504606846976*x^2+27987/14000*x-13303/14000
y3=2249/1125*x^3-8993/3000*x^2-45750853357891/1152921504606846976*x+23063/42000
（2）
（3）>>plot（x,y,'o',z,y1,'-r',z,y2,'-b',z,y3,'-k'）%绘制拟合函数图象
>>
5、分别用欧拉公式、梯形公式和改进欧拉方法求解：
取步长h=0.1.要求列表输出解的结果，并比较不同方法的精度。
解：
建立函数欧拉公式
function[x,y]=naeuler（dyfun,xspan,y0,h）
x=xspan
（1）:
h:
xspan
（2）;
y
（1）=y0;
forn=1:
length（x）-1
y（n+1）=y（n）+h*feval（dyfun,x（n）,y（n））;
end
x=x';y=y';
梯形公式的Matlab实现程序
function[x,y]=tixing（dyfun,xspan,y0,h）
x=xspan
（1）:
h:
xspan
（2）;
y
（1）=y0;
forn=1:
length（x）-1
y（n+1）=iter（dyfun,x（n+1）,y（n）,h）;
end
x=x';y=y';
functiony=iter（dyfun,x,y,h）
y0=y;e=1e-4;K=1e+4;
y=y+h*feval（dyfun,x,y）;
y1=y+2*e;k=1;
whileabs（y-y1）>e
y1=y;
y=y0+h*feval（dyfun,x,y）;
k=k+1;ifk>K,error（'迭代发散'）;end
end
改进欧拉公式的Matlab实现程序:
Function[x,y]=naeuler2（dyfun,xspan,y0,h）
x=xspan
（1）:
h:
xspan
（2）;
y
（1）=y0;
forn=1:
length（x）-1
k1=feval（dyfun,x（n）,y（n））;
y（n+1）=y（n）+h*k1;
k2=feval（dyfun,x（n+1）,y（n+1））;
y（n+1）=y（n）+h*（k1+k2）/2;
end
x=x';y=y';
>>dyfun=inline（'y+2*x/y^2'）;
>>[x,y]=naeuler（dyfun,[0,1],1,0.1）;%调用欧拉公式
>>[x,y]
ans=
01.0000
0.10001.1000
0.20001.2265
0.30001.3758
0.40001.5450
0.50001.7331
0.60001.9397
0.70002.1655
0.80002.4119
0.90002.6806
1.00002.9737
>>dyfun=inline（'y+2*x/y^2'）;
>>[x,y]=tixing（dyfun,[0,1],1,0.1）;%调用梯形公式
>>[x,y]
ans=
01.0000
0.10001.1286
0.20001.2810
0.30001.4549
0.40001.6492
0.50001.8644
0.60002.1017
0.70002.3631
0.80002.6510
0.90002.9682
1.00003.3182
>>dyfun=inline（'y+2*x/y^2'）;
>>[x,y]=naeuler2（dyfun,[0,1],1,0.1）;%调用改进欧拉公式
>>[x,y]
ans=
01.0000
0.10001.1133
0.20001.2520
0.30001.4128
0.40001.5936
0.50001.7939
0.60002.0143
0.70002.2560
0.80002.5207
0.90002.8107
1.00003.1287
>>x=0:
0.1:
1;
>>y=sqrt（1+2.*x）;%精确解
>>[x',y']
ans=
01.0000
0.10001.0954
0.20001.1832
0.30001.2649
0.40001.3416
0.50001.4142
0.60001.4832
0.70001.5492
0.80001.6125
0.90001.6733
1.00001.7321
6、分别用复化梯形公式、复化辛普森公式计算
（区间10等分），并与积分精确值作比较。
解：
建立函数复化梯形公式
functionI=T_quad（x,y）
n=length（x）;m=length（y）;
ifn~=m
error（'wrong'）;
return;
end
h=（x（n）-x
（1））/（n-1）;a=[12*ones（1,n-2）1];
I=h/2*sum（a.*y）;
建立复化辛普森公式
functionI=S_quad（x,y）
n=length（x）;m=length（y）;
ifn~=m
error（'hh'）;
return;
end
ifrem（n-1,2）~=0
I=T_quad（x,y）;
return;
end
N=（n-1）/2;h=（x（n）-x
（1））/N;a=zeros（1,n）;
fork=1:
N
a（2*k-1）=a（2*k-1）+1;a（2*k-1）=a（2*k）+4;
a（2*k+1）=a（2*k+1）+1;
end
I=h/6*sum（a.*y）;
命令
>>x=0:
0.1:
1;
>>y=sqrt（（1-exp（-x）））./exp（x）;
>>T=T_quad（x,y）%复化梯形公式求积分
T=
0.3284
>>S=S_quad（x,y）%复化辛普公式求积分
S=
0.2012
>>symsx;%定义符号变量
>>a=int（'sqrt（（1-exp（-x）））/exp（x）',0,1）%使用int函数求函数积分
a=
（2*（1-1/exp
（1））^（3/2））/3%有理说解
>>vpa（a）%转化为数值解
ans=
0.33504922214016864454877581388462
综合题第一题
命令：
v0=100%初速度
L=200%击球处离本垒的水平距离
g=9.8/0.3048%重力加速度，转为英尺/平方秒
symsa
a=30%速度与水平夹角
x=a*pi/180%转为弧度制
vx=v0*cos（x）%水平速度
vy=v0*sin（x）%垂直方向速度
t=L/vx
LH=vy*t-0.5*g*t^2
if（LH+3>35）
disp（'可以飞过围墙'）
else
disp（'不可以飞过围墙'）
end
symsx
vx=v0*cos（x）
vy=v0*sin（x）
t=L/vx
symsh
h=vy*t-0.5*g*t^2-32
fzero（inline（h）,0）
v0=
100
L=
200
g=
32.1522
a=
30
x=
0.5236
vx=
86.6025
vy=
50.0000
t=
2.3094
LH=
29.7308
不可以飞过围墙
vx=
100*cos（x）
vy=
100*sin（x）
t=
2/cos（x）
h=
200*sin（x）/cos（x）-24500/381/cos（x）^2-32
ans=
0.5372
验证性的题目：
2.2.2第1题
（1）>>fun=inline（'x^2-0.8*x+0.15'）;
[x_star,inedx,it]=bisect（fun,0.1,0.6）
>>x_star=
0.08000.0300
inedx=
0
it=
0
当index=0时，表明初始区间不是有根区间。
（2）
（3）
3.2.1第1题
>>A=[12-1;012;-364];B=[-101;022;351]
C=B'
D=A+B
E=A*B
F=A^2
G=inv（A）*B
H=A*B
B=
-101
022
351
C=
-103
025
121
D=
020
034
0115
E=
-4-14
6124
153213
F=
4-2-1
-61310
-152431
G=
-1.00000
升级会员 