数值分析上机试题对应参考答案Word文档下载推荐.docx
《数值分析上机试题对应参考答案Word文档下载推荐.docx》由会员分享，可在线阅读，更多相关《数值分析上机试题对应参考答案Word文档下载推荐.docx（24页珍藏版）》请在冰豆网上搜索。
（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-0.44480.00240.3447
000.09250.9426
0000.6955
P=
0001
0100
1000
0010
P*A%验证PA=LU
ans=
L*U
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:
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）
%已知数据点的牛顿插值
要计算的插值点，可以是多个
用牛顿插值函数算出xt对应的函数值数组
牛顿插值多项式表达式
a=zeros（1,n）;
N=y
（1）;
w=1;
fork=1:
n-1
yy=zeros（1,n）;
%记录差商
forj=k+1:
yy（j）=（y（j）-y（k））/（x（j）-x（k））;
a（k）=yy（k+1）;
w=w*（t-x（k））;
N=N+a（k）*w;
y=yy;
yt=subs（N,'
simplify（N）;
N=collect（N）;
%将插值多项式展开
N=vpa（N,6）;
%系数转化为6位精度
命令：
x=0:
pi/10:
2*pi;
y=sin（x）;
z=0:
pi/20:
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:
dg=A（i,i）;
ifabs（dg）<
del
disp（'
主对角元素太小'
return
imax
fori=1:
sm=b（i）;
forj=1:
ifj~=i
sm=sm-A（i,j）*x0（j）;
x（i）=sm/A（i,i）;
tx=[tx;
x];
ifnorm（x-x0）<
acc
else
x0=x;
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）>
o'
-r'
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））;
x=x'
;
y=y'
梯形公式的Matlab实现程序
function[x,y]=tixing（dyfun,xspan,y0,h）
y（n+1）=iter（dyfun,x（n+1）,y（n）,h）;
x=x'
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（'
迭代发散'
改进欧拉公式的Matlab实现程序:
Function[x,y]=naeuler2（dyfun,xspan,y0,h）
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;
dyfun=inline（'
y+2*x/y^2'
[x,y]=naeuler（dyfun,[0,1],1,0.1）;
%调用欧拉公式
[x,y]
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
[x,y]=tixing（dyfun,[0,1],1,0.1）;
%调用梯形公式
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
[x,y]=naeuler2（dyfun,[0,1],1,0.1）;
%调用改进欧拉公式
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
1;
y=sqrt（1+2.*x）;
%精确解
[x'
y'
]
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）
m=length（y）;
ifn~=m
wrong'
return;
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）
hh'
ifrem（n-1,2）~=0
I=T_quad（x,y）;
N=（n-1）/2;
h=（x（n）-x
（1））/N;
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;
I=h/6*sum（a.*y）;
命令
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）%转化为数值解
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）
可以飞过围墙'
else
不可以飞过围墙'
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
200
g=
32.1522
30
x=
0.5236
vx=
86.6025
vy=
50.0000
t=
2.3094
LH=
29.7308
不可以飞过围墙
100*cos（x）
100*sin（x）
2/cos（x）
h=
200*sin（x）/cos（x）-24500/381/cos（x）^2-32
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=
当index=0时，表明初始区间不是有根区间。
（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
升级会员 