数值分析计算实习题.docx
《数值分析计算实习题.docx》由会员分享，可在线阅读，更多相关《数值分析计算实习题.docx（28页珍藏版）》请在冰点文库上搜索。
数值分析计算实习题
插值法
1.下列数据点的插值
x01491625364964
y012345678
可以得到平方根函数的近似，在区间[0,64]上作图.
（1）用这9个点作8次多项式插值Ls（x）.
（2）用三次样条（第一边界条件）程序求S（x）.
从得到结果看在[0,64]上，哪个插值更精确；在区间[0,1]上，两种插值哪个更精确？
解：
（1）拉格朗日插值多项式，求解程序如下
symsxl;
x1=[01491625364964];
y1=[012345678];
n=length（x1）;
Ls=sym（0）;
fori=1:
n
l=sym（y1（i））;
fork=1:
i-1
l=l*（x-x1（k））/（x1（i）-x1（k））;
end
fork=i+1:
n
l=l*（x-x1（k））/（x1（i）-x1（k））;
end
Ls=Ls+l;
end
Ls=simplify（Ls）%为所求插值多项式Ls（x）.
输出结果为
Ls=
-24221063/63504000*x^2+95549/72072*x-1/3048192000*x^8-2168879/435456000*x^4+19/283046400*x^7+657859/10886400*x^3+33983/152409600*x^5-13003/2395008000*x^6
（2）三次样条插值，程序如下
x1=[01491625364964];
y1=[012345678];
x2=[0:
1:
64];
y3=spline（x1,y1,x2）;
p=polyfit（x2,y3,3）;%得到三次样条拟合函数
S=p
（1）+p
（2）*x+p（3）*x^2+p（4）*x^3%得到S（x）
输出结果为：
S=
2288075067923491/73786976294838206464-2399112304472833/576460752303423488*x+4552380473376713/18014398509481984*x^2+999337332656867/1125899906842624*x^3
（3）在区间[0,64]上，分别对这两种插值和标准函数作图，
plot（x2,sqrt（x2）,'b',x2,y2,'r',x2,y3,'y'）
蓝色曲线为y=
函数曲线,红色曲线为拉格朗日插值函数曲线，黄色曲线为三次样条插值曲线
可以看到蓝色曲线与黄色曲线几乎重合，因此在区间[0,64]上三次样条插值更精确。
在[0,1]区间上由上图看不出差别，不妨代入几组数据进行比较，取x4=[0:
0.2:
1]
x4=[0:
0.2:
1];
sqrt（x4）%准确值
subs（Ls,'x',x4）%拉格朗日插值
spline（x1,y1,x4）%三次样条插值
运行结果为
ans=
00.44720.63250.77460.89441.0000
ans=
00.25040.47300.67060.84551.0000
ans=
00.24290.46300.66170.84031.0000
从这几组数值上可以看出在[0,1]区间上，拉格朗日插值更精确。
数据拟合和最佳平方逼近
2.有实验给出数据表
x0.00.10.20.30.50.81.0
y1.00.410.500.610.912.022.46
试求3次、4次多项式的曲线拟合，再根据数据曲线形状，求一个另外函数的拟合曲线，用图示数据曲线及相应的三种拟合曲线。
解：
（1）三次拟合，程序如下
symx;
x0=[0.00.10.20.30.50.81.0];
y0=[1.00.410.500.610.912.022.46];
cc=polyfit（x0,y0,3）;
S3=cc
（1）+cc
（2）*x+cc（3）*x^2+cc（4）*x^3%三次拟合多项式
xx=x0
（1）:
0.1:
x0（length（x0））;
yy=polyval（cc,xx）;
plot（xx,yy,'--'）;
holdon;
plot（x0,y0,'x'）;
xlabel（'x'）;
ylabel（'y'）;
运行结果
S3=
-7455778416425083/1125899906842624+1803512222945437/140737488355328*x-655705280524945/140737488355328*x^2+4172976892199509/4503599627370496*x^3
图像如下
（2）4次多项式拟合
symx;
x0=[0.00.10.20.30.50.81.0];
y0=[1.00.410.500.610.912.022.46];
cc=polyfit（x0,y0,4）;
S3=cc
（1）+cc
（2）*x+cc（3）*x^2+cc（4）*x^3+cc（5）*x^4
xx=x0
（1）:
0.1:
x0（length（x0））;
yy=polyval（cc,xx）;
plot（xx,yy,'r'）;
holdon;
plot（x0,y0,'x'）;
xlabel（'x'）;
ylabel（'y'）;
运行结果
S3=
3248542900396215/1125899906842624-3471944732519151/281474976710656*x+4580931990070637/281474976710656*x^2-5965836931688425/1125899906842624*x^3+8491275464650307/9007199254740992*x^4
图像如下
（3）另一个拟合曲线,
新建一个M-file
程序如下：
function[C,L]=lagran（x,y）
w=length（x）;
n=w-1;
L=zeros（w,w）;
fork=1:
n+1
V=1;
forj=1:
n+1
ifk~=j
V=conv（V,poly（x（j）））/（x（k）-x（j））;
end
end
L（k,:
）=V;
end
C=y*L
在命令窗口中输入以下的命令：
x=[0.00.10.20.30.50.81.0];
y=[1.00.410.500.610.912.022.46];
cc=polyfit（x,y,4）;
xx=x
（1）:
0.1:
x（length（x））;
yy=polyval（cc,xx）;
plot（xx,yy,'r'）;
holdon;
plot（x,y,'x'）;
xlabel（'x'）;
ylabel（'y'）;
x=[0.00.10.20.30.50.81.0];
y=[0.940.580.470.521.002.002.46];%y中的值是根据上面两种拟合曲线而得到的估计数据，不是真实数据
[C,L]=lagran（x,y）;
xx=0:
0.01:
1.0;
yy=polyval（C,xx）;
holdon;
plot（xx,yy,'b',x,y,'.'）;
图像如下
解线性方程组的直接解法
3.线性方程组Ax=b的A及b为
A=
，b=
，则解x=
.用MATLAB内部函数求detA及A的所有特征值和cond（A）2.若令
A+δA=
，求解（A+δA）（x+δx）=b,输出向量x和||δx||2.从理论结果和实际计算两方面分析线性方程组Ax=b解得相对误差||δx||2/||x||2及A的相对误差||δA||2/||A||2的关系.
解：
（1）程序如下
clear;
A=[10787;7565;86109;75910];
det（A）
cond（A,2）
eig（A）
输出结果为
ans=
1
ans=
2.9841e+003
ans=
0.0102
0.8431
3.8581
30.2887
（2）程序如下
A=[1078.17.2;7.085.0465;85.989.899;6.99599.98];
b=[32233331]';
x0=[1111]';
x=A\b%扰动后方程组的解
x1=x-x0%δx的值
norm（x1,2）%δx的2-范数
运行结果为
x=
-9.5863
18.3741
-3.2258
3.5240
x1=
-10.5863
17.3741
-4.2258
2.5240
ans=
20.9322
程序如下
A=[1078.17.2;7.085.0465;85.989.899;6.99599.98];
A0=[10787;7565;86109;75910];
b=[32233331]';
x0=[1111]';
x=A\b;
x1=x-x0;
norm（x1,2）;
A1=A-A0;%δA的值
norm（x1,2）/norm（x0,2）%||δx||2/||x||2的值
norm（A1,2）/norm（A0,2）%||δA||2/||A||2的值
输出结果为
ans=
10.4661
ans=
0.0076
可见A相对误差只为0.0076，而得到的结果x的相对误差就达到了10.4661，该方程组是病态的，A的条件数为2984.1远远大于1，当A只有很小的误差就会给结果带来很大的影响。
非线性方程数值解法
4.求下列方程的实根
（1）x^2-3x+2-e^x=0;
（2）x^3+2x^2+10x-20=0.
要求：
（1）设计一种不动点迭代法，要使迭代序列收敛，然后再用斯特芬森加速迭代，计算到|x（k）-x（k-1）|<10^（-8）为止。
（2）用牛顿迭代，同样计算到|x（k）-x（k-1）|<10^（-8）。
输出迭代初值及各次迭代值和迭代次数k，比较方法的优劣。
解:
（1）先用画图的方法估计根的范围
ezplot（'x^2-3*x+2-exp（x）'）;
gridon;
可以估计到方程的根在区间（0,1）;选取迭代初值为x0=0.5；
构造不动点迭代公式x（k+1）=（x（k）^2+2-e^x（k））/3;
ψ（x）=（x^2+2-e^x）/3;
当0程序如下：
formatlong;
f=inline（'（x^2+2-exp（x））/3'）
disp（'x='）;
x=feval（f,0.5）;
disp（x）;
Eps=1E-8;
i=1;
while1
x0=x;
i=i+1;
x=feval（f,x）;
disp（x）;
ifabs（x-x0）break;
end
end
i,x
运行结果为
f=
Inlinefunction:
f（x）=（x^2+2-exp（x））/3
x=
0.20042624309996
0.27274906509837
0.25360715658413
0.25855037626494
0.25726563633509
0.25759898516219
0.25751245451483
0.25753491361525
0.25752908416796
0.25753059723833
0.25753020451046
0.25753030644564
0.25753027998767
0.25753028685501
i=
14
x=
0.25753028685501
斯特芬森加速法，程序如下:
formatlong;
f=inline（'x-（（x^2+2-exp（x））/3-x）^2/（（（（x^2+2-exp（x））/3）^2+2-exp（（x^2+2-exp（x））/3））/3-2*（x^2+2-exp（x））/3+x）'）;
disp（'x='）;
x=feval（f,0.5）;
disp（x）;
Eps=1E-8;
i=1;
while1
x0=x;
i=i+1;
x=feval（f,x）;
disp（x）;
ifabs（x-x0）break;
end
end
i,x
运行结果为x=
0.25868442756579
0.25753031771981
0.25753028543986
0.25753028543986
i=
4
x=
0.25753028543986
牛顿迭代法，程序如下：
formatlong;
x=sym（'x'）;
f=sym（'x^2-3*x+2-exp（x）'）;
df=diff（f,x）;
FX=x-f/df;
Fx=inline（FX）;
disp（'x='）;
x1=0.5;
disp（x1）;
Eps=1E-8;
i=0;
while1
x0=x1;
i=i+1;
x1=feval（Fx,x1）;
disp（x1）;
ifabs（x1-x0）break;
end
end
i,x1
运行结果如下：
x=
0.50000000000000
0.25368870241829
0.25752890079471
0.25753028543968
0.25753028543986
i=
4
x1=
0.25753028543986
（2）先用画图的方法估计根的范围
ezplot（'x^3+2*x^2+10*x-20'）;
gridon;
根大约在区间（1,2）；选取初值x0=1.5;
构造不动点迭代公式x（k+1）=（-2x（k）^2-10x（k）+20）^1/3;
ψ（x）=（-2x^2-10x+20）^1/3;
程序如下：
formatlong;
f=inline（'（-2*x^2-10*x+20）^1/3'）
disp（'x='）;
x=feval（f,1.5）;
disp（x）
Eps=1E-8;
i=1;
while1
x0=x;
i=i+1;
x=feval（f,x）;
disp（x）;
ifabs（x-x0）>1E10
break;
end
ifabs（x-x0）break;
end
end
i,x
运行结果：
f=
Inlinefunction:
f（x）=（-2*x^2-10*x+20）^1/3
x=
0.166********667
6.09259259259259
-38.38843164151806
-8.478196837919431e+002
-4.763660785374071e+005
-1.512815059604763e+011
i=
6
x=
-1.512815059604763e+011
迭代6次后x的值大得令人吃惊,表明构造的式子并不收敛.
也无法构造出收敛的不动点公式
牛顿迭代法，程序如下：
formatlong;
x=sym（'x'）;
f=sym（'x^3+2*x^2+10*x-20'）;
df=diff（f,x）;
FX=x-f/df;
Fx=inline（FX）;
disp（'x='）;
x1=0.5;
disp（x1）;
Eps=1E-8;
i=0;
while1
x0=x1;
i=i+1;
x1=feval（Fx,x1）;
disp（x1）;
ifabs（x1-x0）break;
end
end
i,x1
运行结果：
x=
1.50000000000000
1.37362637362637
1.36881481962396
1.36880810783441
1.36880810782137
i=
4
x1=
1.36880810782137
比较三种方法，牛顿法的收敛性比较好，相比不动点迭代法要构造出收敛的公式比较难，牛顿法迭代次数也较少，收敛速度快，只是对初值的要求很高，几种方法各有利弊，具体采用哪种也需因题而异。
常微分方程初值问题数值解法
5.给定初值问题
y’=-50y+50x^2+2x,
;
y（0）=1/3;
用经典的四阶R-K方法解该问题，步长分别取h=0.1,0.025,0.01,计算并打印x=0.1i（i=0,1,…,10）各点的值，与准确值y（x）=1/3e^（-50x）+x^2比较。
解：
取步长h=0.1，程序如下：
%经典的四阶R-K方法
clear;
F='-50*y+50*x^2+2*x';
a=0;b=1;
h=0.1;
n=（b-a）/0.1;
X=a:
0.1:
b;
Y=zeros（1,n+1）;
Y
（1）=1/3;
fori=1:
n
x=X（i）;
y=Y（i）;
K1=h*eval（F）;
x=x+h/2;
y=y+K1/2;
K2=h*eval（F）;
y=Y（i）+K2/2;
K3=h*eval（F）;
x=X（i）+h;
y=Y（i）+K3;
K4=h*eval（F）;
Y（i+1）=Y（i）+（K1+2*K2+2*K3+K4）/6;
end
%准确值
temp=[];
f=dsolve（'Dy=-50*y+50*x^2+2*x','y（0）=1/3','x'）;
df=zeros（1,n+1）;
fori=1:
n+1
temp=subs（f,'x',X（i））;
df（i）=double（vpa（temp））;
end
disp（'步长四阶经典R-K法准确值'）;
disp（[X',Y',df']）;
运行结果：
步长四阶经典R-K法准确值
1.0e+010*
00.000000000033330.00000000003333
0.000000000010000.000000000460550.00000000000122
0.000000000020000.000000006306250.00000000000400
0.000000000030000.000000086404940.00000000000900
0.000000000040000.000001184363000.00000000001600
0.000000000050000.000016235451100.00000000002500
0.000000000060000.000222560671340.00000000003600
0.000000000070000.003050935427780.00000000004900
0.000000000080000.041823239217400.00000000006400
0.000000000090000.573326903478090.00000000008100
0.000000000100007.859356300837710.00000000010000
%画图观察结果
figure;
plot（X,df,'k*',X,Y,'--r'）;
gridon;
title（'四阶经典R-K法解常微分方程'）;
legend（'准确值','四阶经典R-K法'）;
当x值接近1的时候，偏离准确值太大。
当步长h=0.025时，将上面程序中的h改为0.025即可，运行结果：
步长四阶经典R-K法准确值
00.333333333333330.33333333333333
0.100000000000000.103135240342880.01224598233303
0.200000000000000.044285273275990.04001513330992
0.300000000000000.051967957555070.09000010196744
0.400000000000000.093957311494390.16000000068705
0.500000000000000.160345314357570.25000000000463
0.600000000000000.248085701305570.36000000000003
0.700000000000000.356241884727580.49000000000000
0.800000000000000.484525906616270.64000000000000
0.900000000000000.632849233007510.81000000000000
1.000000000000000.801184643742061.00000000000000
图像如下：
当步长h=0.01时，将上面程序中的h改为0.01即可，运行结果：
步长四阶经典R-K法准确值
00.333333333333330.33333333333333
0.100000000000000.202357204861110.01224598233303
0.200000000000000.128817001907910.04001513330992
0.300000000000000.097991826678500.09000010196744
0.400000000000000.100949467750240.16000000068705
0.500000000000000.132********4700.25000000000463
0.600000000000000.188********1990.36000000000003
0.700000000000000.268139302784310.49000000000000
0.800000000000000.369481660283190.64000000000000
0.900000000000000.491957621994750.81000000000000
1.000000000000000.635121421679101.00000000000000
七年级英语期末考试质量分析
一、试卷分析：
本次试卷的难易程度定位在面向大多数学生。
该份试卷紧扣教材，突出重点，注重对基础知识
升级会员 