数值计算方法上机实习题考证.docx
《数值计算方法上机实习题考证.docx》由会员分享，可在线阅读，更多相关《数值计算方法上机实习题考证.docx（38页珍藏版）》请在冰点文库上搜索。
数值计算方法上机实习题考证
---------------------------------------------------
此文档包含我们计算方法的经典算法
包含（数值计算方法上机实习题）
1．设
，
（1）由递推公式
，从
的几个近似值出发，计算
；
（2）粗糙估计
，用
，计算
；
（3）分析结果的可靠性及产生此现象的原因（重点分析原因）。
（1）解答：
n=0，
这里可以用for循环，while循环，根据个人喜好与习惯：
for循环程序：
While循环程序：
I=0.1823;I=0.1823;
forn=1:
20i=1;
I=（-5）*I+1/n;whilei<21
EndI=（-5）*I+1/i;
Ii=i+1;
fprintf（'I20=%f',I）end
èI=-2.0558e+009>>I
I20=-2055816073.851284>>I=-2.0558e+009
（2）
粗略估计I20：
Mathcad计算结果：
for循环程序：
While循环程序：
>>I=0.007998;I=0.007998;
>>forn=1:
20n=1;
I=（-0.2）*I+1/（5*n）;whilen<21
EndI=（-0.2）*I+1/（5*n）;
>>In=n+1;
I=0.0083end
>>I
I=0.0083
（3）算法误差分析：
计算在递推过程中传递截断误差和舍入误差
第一种算法：
（从1——>20）
误差放大了5n倍，算法稳定性很不好；
第二种算法：
（从20——>1）
误差在逐步缩小，算法趋近稳定，收敛。
2．求方程
的近似根，要求
，并比较计算量。
（1）在[0，1]上用二分法；
function[ti]=erfenfa（a,b）
f=@（x）（exp（x）+10*x-2）
t=（a+b）./2;
i=0;
whileabs（f（t））>0.001
iff（a）*f（t）<0
b=t;t=（a+b）/2;
elseiff（b）*f（t）<0
a=t;t=（b+a）/2;
end
i=i+1;
end
（2）取初值
，并用迭代
；
functionx=diedai（x0）%x0初值
x=x0;
fori=1:
10000
y=（2-exp（x））./10;x=y;y=（2-exp（x））./10;
ifabs（x-y）<5*10^（-4）
disp（'迭代次数'）;
2*i
break;
end
end
（3）加速迭代的结果（艾特肯Aitken加速方法）；
function[ym]=aitken（a）
func=@（x）（（2-exp（x））./10）
x
（1）=a;
wucha=1;m=1;
whilewucha>5*10^（-4）
p（m+1）=func（x（m））;
q（m+1）=func（p（m+1））;
x（m+1）=q（m+1）-（（q（m+1）-p（m+1））^2）./（q（m+1）-2*p（m+1）+x（m））;
wucha=abs（x（m+1）-x（m））;
m=m+1;
ifm>1000
break;
end
end
formatlong
y=x（m-1）;
m=m-1;
运行结果y=0.090483741803596
m=2
（4）取初值
，并用牛顿迭代法；
functionx=newtondiedai（x0）
x=x0;
fori=1:
10000
y=x-（exp（x）+10*x-2）./（exp（x）+10）;x=y;
y=x-（exp（x）+10*x-2）./（exp（x）+10）;
ifabs（x-y）<0.0001
disp（'迭代次数'）;
i
break;
end
end
（5）分析绝对误差。
|x-x*|=
|0.090525101307255-0.0905|=
0.000025101307255<1/（2*（0+1））*
10^（-（11-1）=0.5*10^（-10）
3．钢水包使用次数多以后，钢包的容积增大，数据如下：
x
2
3
4
5
6
7
8
9
y
6.42
8.2
9.58
9.5
9.7
10
9.93
9.99
10
11
12
13
14
15
16
10.49
10.59
10.60
10.8
10.6
10.9
10.76
方法一：
一个拉格朗日函数解析式和插值结果的计算程序
functionf=Language00（x,y,x0）
%求已知数据点的拉格朗日插值多项式
%已知数据点的x坐标向量:
x
%已知数据点的y坐标向量:
y
%插值点的x坐标:
x0
%求得的拉格朗日插值多项式或在x0处的插值:
f
symstl;
if（length（x）==length（y））
n=length（x）;
else
disp（'x和y的维数不相等！
'）;
return;%检错
end
h=sym（0）;
for（i=1:
n）
l=sym（y（i））;
for（j=1:
n）
ifi~=j
l=l*（t-x（j））/（x（i）-x（j））;
end
end
h=h+l;
end
simplify（h）;
if（nargin==3）
f=subs（h,'t',x0）;%计算插值点的函数值
else
f=collect（h）;
f=vpa（f,6）;%将插值多项式的系数化成6位精度的小数
end
多种拟合方法的曲线对比
function[y1x1]=lagrange（x,y,x0）%x0为拉格朗日差值函数点
nx=length（x）;ny=length（y）;
ifnx~=ny
warning（'向量x与y的长度应该相等'）
return
end
m=length（x0）;
x1=x0;
fori=1:
m;
t=0.0;
forj=1:
nx;
u=1.0;
fork=1:
nx;
ifk~=j
u=u*（x0（i）-x（k））/（x（j）-x（k））;
end
end
t=t+u*y（j）;
end
y1（i）=t;
end
xp=2:
0.1:
16;%这里的点的间距细一些.
yp=spline（x,y,xp）;
t=polyfit（x,y,4）;
disp（'最小二乘法拟合函数解析式'）;
L=poly2sym（t,'x'）
disp（'显示拉格朗日插值x0结果'）
yf=polyval（t,xp）;
subplot（321）
plot（xp,yp,'g-','LineWidth',2）;title（'样条插值曲线'）;
subplot（322）
plot（x,y,'r--','LineWidth',2）;title（'原始数据点图'）;
subplot（323）
plot（xp,yf,'b+','LineWidth',2）;title（'最小二乘法拟合'）;
subplot（324）
plot（x1,y1,'y-','LineWidth',2）;
title（'拉格朗日插值点的函数值图像'）;
subplot（325）
plot（xp,yf,'b+',xp,yp,'g-',x,y,'r--',x1,y1,'y-','LineWidth',2）;
title（'样条、拉格朗日、最小二乘法与原始数据对比图'）;
在commandwindow输入一下函数值调用
x=[2：
16]
y=[6.428.29.589.59.7109.939.9910.4910.5910.6010.810.610.910.76];
x0=[2.53.54.55.56.57.58.59.510.511.512.513.514.515.515.7]
方法二：
因为不知道其函数形式，可先plot而后确定使用哪种模型比较合适。
函数图形程序：
x=[2345678910111213141516];
y=[6.428.29.589.59.7109.939.9910.4910.5910.610.810.610.910.76];
plot（x,y）
与指数函数趋势类似（但是趋势不同，一个趋于递减，一个趋于递增，使用倒数），故拟合模型为：
两边同时取对数：
用此方程进行数据拟合：
x=[2345678910111213141516];
y=[6.428.29.589.59.7109.939.9910.4910.5910.610.810.610.910.76];
t=1./x;
s=polyfit（t,log（y）,1）
s=
-1.11072.4578
即是：
最终拟合曲线为：
将此拟合曲线和源数据进行对比：
主程序：
x=[2345678910111213141516];
y=[6.428.29.589.59.7109.939.9910.4910.5910.610.810.610.910.76];
plot（x,y,'r'）
holdon
lims1=[2,16];
fplot（'11.679*exp（-1.1107/x）',lims1）
holdoff
均方差的求解
x=[2:
16];
y=[6.428.29.589.59.7109.939.9910.4910.5910.6010.810.610.910.76];
f（x）=11.679*exp（-1.1107./x）;
c=0;
fori=1:
15
a=y（i）;
b=x（i）;
c=c+（a-f（b））^2;
end
averge=c/15
averge=
0.0594
4．设
，
，
分析下列迭代法的收敛性，并求
的近似解及相应的迭代次数。
（1）JACOBI迭代；
function[x,j]=JACOBI（A,b,x0）
n=length（b）;
l=zeros（n,n）;
u=zeros（n,n）;
x=x0;
fori=1:
n;
forj=1:
n;
ifi>j
l（i,j）=-A（i,j）;
elseifi==j
d（i,j）=A（i,j）;
else
u（i,j）=-A（i,j）;
end
end
end
disp（'增广矩阵B'）
B=[A,b]
B=inv（d）*（l+u）;
f=inv（d）*b;
forj=1:
1000
x1=B*x+f;
x=x1;x1=B*x+f;
ifabs（x1-x）<0.0001
break;
elseifj>1000
disp（'已达到迭代设置1000次或者迭代不收敛'）;
end
end
（2）GAUSS-SEIDEL迭代；
function[x,j]=Gauss（A,b,x0）
n=length（b）;
l=zeros（n,n）;
u=zeros（n,n）;
x=x0;
fori=1:
n;
forj=1:
n;
ifi>j
l（i,j）=-A（i,j）;
elseifi==j
d（i,j）=A（i,j）;
else
u（i,j）=-A（i,j）;
end
end
end
disp（'增光矩阵'）
B=[A,b]
G=inv（d-l）*u;
d1=inv（d-l）*b;
forj=1:
1000
x1=G*x+d1;
x=x1;x1=G*x+d1;
ifabs（x1-x）<0.0001
break;
elseifj>1000
disp（'迭代不收敛或迭代次数超出1000'）;
end
end
（3）SOR迭代（
）
function[x,j]=SOR（A,b,x0,w）
n=length（b）;
l=zeros（n,n）;
u=zeros（n,n）;
x=x0;
fori=1:
n;
forj=1:
n;
ifi>j
l（i,j）=-A（i,j）;
elseifi==j
d（i,j）=A（i,j）;
else
u（i,j）=-A（i,j）;
end
end
end
disp（'曾广矩阵'）
B=[A,b]
G=inv（d-w*l）*（（1-w）*d+w*u）;
d1=w*inv（d-w*l）*b;
forj=1:
1000
x1=G*x+d1;
x=x1;x1=G*x+d1;
ifabs（x1-x）<0.0001
break;
elseifj>1000
disp（'迭代不收敛或者迭代次数太少'）;
end
end
5．用逆幂迭代法求
最接近于11的特征值和特征向量，准确到
。
解答：
以下所得结果是最小特征值对应的参数结果
function[m,u,index,k]=pow_inv（A,ep,it_max）
%求矩阵最小特征值的反幂法，其中?
%?
A为矩阵；
%ep为精度要求，缺省为1e-5；
%it_max为最大迭代次数，缺省为100；m为绝对值最大的特征值；
%index，当index=1时，迭代成功，当index=0时，迭代失败
ifnargin<3
it_max=1000;end
ifnargin<2
ep=1e-5;
end
n=length（A）;index=0;k=0;m1=0;m0=0;
%?
修改移位参数，原点移位法加速收敛，为0时，即为反幂法?
I=eye（n）;
T=A-m0*I;
invT=inv（T）;
u=ones（n,1）
whilek<=it_max
v=invT*u;
[vmax,i]=max（abs（v））;m=v（i）;u=v/m;
ifabs（m-m1）index=1;
break;
end
m1=m;k=k+1;
end
m=1/m;
m=m+m0;
以下所得结果是最大特征值对应的参数结果
function[m,u,index,k]=pow（A,ep,it_max）
%求矩阵最大特征值的幂法，其中A为矩阵；
%ep为精度要求，缺省为1e-5；
%it_max为最大迭代次数，缺省为100；
%m为绝对值最大的特征值；
%index，当index=1时，迭代成功，当index=0时，迭代失败
ifnargin<3
it_max=100;
end
ifnargin<2
ep=1e-5;
end
n=length（A）;index=0;k=0;m1=0;m0=0.01;%修改移位参数，原点移位法加速收敛，为0时，即为幂法
I=eye（n）;
T=A-m0*I;
u=ones（n,1）;
whilek<=it_max
v=T*u;
[vmax,i]=max（abs（v））;
m=v（i）;u=v/m;
ifabs（m-m1）index=1;
break;
end
m=m+m0;
m1=m;
k=k+1;
end
6．用经典R-K方法求解初值问题
（1）
，
，
；
和精确解
比较，分析结论。
解答：
程序函数
*******************************************************************************
a=0;
b=10;
h=0.01;
n=（b-a）/h;
t=a;
u1=2;
u2=3;
fori=1:
n
%每个方程的k1
k（1,1）=h*（-2*u1+u2-2*sin（t））;
k（1,2）=h*（u1-2*u2+2*cos（t）-2*sin（t））;
%每个方程的k2
k（2,1）=h*（3*（u1+0.5*k（1,1））+2*（u2+0.5*k（1,2））-（2*（t+0.5*h）^2+1）*exp（2*（t+0.5*h）））;k（2,2）=h*（4*（u1+0.5*k（1,1））+（u2+0.5*k（1,2））+（（t+0.5*h）^2+2*（t+0.5*h）-4）*exp（2*（t+0.5*h）））;
%每个方程的k3
k（3,1）=h*（3*（u1+0.5*k（2,1））+2*（u2+0.5*k（2,2））-（2*（t+0.5*h）^2+1）*exp（2*（t+0.5*h）））;
k（3,2）=h*（4*（u1+0.5*k（2,1））+（u2+0.5*k（2,2））+（（t+0.5*h）^2+2*（t+0.5*h）-4）*exp（2*（t+0.5*h）））;
%每个方程的k4
k（4,1）=h*（3*（u1+k（3,1））+2*（u2+k（3,2））-（2*（t+h）^2+1）*exp（2*（t+h）））;
k（4,2）=h*（4*（u1+k（3,1））+（u2+k（3,2））+（（t+h）^2+2*（t+h）-4）*exp（2*（t+h）））;
%求每个方程的w
u1=u1+（k（1,1）+2*k（2,1）+2*k（3,1）+k（4,1））/6;
u2=u2+（k（1,2）+2*k（2,2）+2*k（3,2）+k（4,2））/6;
t=a+i*h;
end
disp（'w1='）;
disp（u1）;
disp（'误差为'）;
disp（abs（u1-（exp（5）/3-exp（-1）/3+exp
（2））））;
disp（'w2='）;
disp（u2）;
disp（'误差为'）;
disp（abs（u2-exp（5）/3-exp（-1）/3*2-exp
（2）））;
（2）
，
，
。
未做
以下是参照类似试题
题目：
用四阶R-K方法求下列初值问题的解。
1、
2
1、%用四阶R-K方法求P322的1a
disp（'P3221（a）'）;
a=0;
b=1;
h=0.02;
n=（b-a）/h;
t=a;
u1=1;
u2=1;
fori=1:
n%每个方程的k1
k（1,1）=h*（3*u1+2*u2-（2*t^2+1）*exp（2*t））;
k（1,2）=h*（4*u1+u2+（t^2+2*t-4）*exp（2*t））;%每个方程的k2
k（2,1）=h*（3*（u1+0.5*k（1,1））+2*（u2+0.5*k（1,2））-（2*（t+0.5*h）^2+1）*exp（2*（t+0.5*h）））;k（2,2）=h*（4*（u1+0.5*k（1,1））+（u2+0.5*k（1,2））+（（t+0.5*h）^2+2*（t+0.5*h）-4）*exp（2*（t+0.5*h）））;%每个方程的k3
k（3,1）=h*（3*（u1+0.5*k（2,1））+2*（u2+0.5*k（2,2））-（2*（t+0.5*h）^2+1）*exp（2*（t+0.5*h）））;
k（3,2）=h*（4*（u1+0.5*k（2,1））+（u2+0.5*k（2,2））+（（t+0.5*h）^2+2*（t+0.5*h）-4）*exp（2*（t+0.5*h）））;%每个方程的k4
k（4,1）=h*（3*（u1+k（3,1））+2*（u2+k（3,2））-（2*（t+h）^2+1）*exp（2*（t+h）））;
k（4,2）=h*（4*（u1+k（3,1））+（u2+k（3,2））+（（t+h）^2+2*（t+h）-4）*exp（2*（t+h）））;%求每个方程的w
u1=u1+（k（1,1）+2*k（2,1）+2*k（3,1）+k（4,1））/6;
u2=u2+（k（1,2）+2*k（2,2）+2*k（3,2）+k（4,2））/6;
t=a+i*h;
end
disp（'w1='）;
disp（u1）;
disp（'误差为'）;
disp（abs（u1-（exp（5）/3-exp（-1）/3+exp
（2））））;
disp（'w2='）;
disp（u2）;
disp（'误差为'）;
disp（abs（u2-exp（5）/3-exp（-1）/3*2-exp
（2）））;
2、%用四阶R-K方法求P322的2a
第二个方程组的解析参照（非题程序）自己可以修改也与最后一个题类似
disp（'P3222（a）'）;
a=0;
b=1;
h=0.1;
n=（b-a）/h;
t=a;
u1=0;
u2=0;
fori=1:
n%每个方程的k1
k（1,1）=h*（u2）;
k（1,2）=h*（-u1+2*u2+t*exp（t）-t）;%每个方程的k2
k（2,1）=h*（u2+0.5*k（1,2））;
k（2,2）=h*（-（u1+0.5*k（1,1））+2*（u2+0.5*k（1,2））+（t+0.5*h）*exp（t+0.5*h）-（t+0.5*h））;%每个方程的k3
k（3,1）=h*（u2+0.5*k（2,2））;
k（3,2）=h*（-（u1+0.5*k（2,1））+2*（u2+0.5*k（2,2））+（t+0.5*h）*exp（t+0.5*h）-（t+0.5*h））;%每个方程的k4
k（4,1）=h*（u2+k（3,2））;
k（4,2）=h*（-（u1+k（3,1））+2*（u2+k（3,2））+（t+h）*exp（t+h）-（t+h））;
%求每个方程的w
u1=u1+（k（1,1）+2*k（2,1）+2*k（3,1）+k（4,1））/6;
u2=u2+（k（1,2）+2*k（2,2）+2*k（3,2）+k（4,2））/6;
t=a+i*h;
end
disp（'y1='）;
disp（u1）;
disp（'误差为'）;
disp（abs（u1-exp
（1）/6-exp
（1）+3））;
7．用有限差分法求解边值问题（h=0.1）：
.
functionys=dbf（f,a,b,alfa,beta,h,eps）
ff=@（x,y）[y
（2）,f（y
（1）,y
（2）,x）];
xvalue=a:
h:
b;%x取值范围
n=length（xvalue）;
s0=a-0.01;%选取适当的s的初值
x0=[alfa,s0];%迭代初值
flag=0;%用于判断精度
y0=rk4（ff,a,x0,h,a,b）;
ifabs（y0（1,n）-b
升级会员 