计算方法上机实习题.docx
《计算方法上机实习题.docx》由会员分享，可在线阅读，更多相关《计算方法上机实习题.docx（26页珍藏版）》请在冰点文库上搜索。
计算方法上机实习题
数值计算方法上机实习题
1．设
，
（1）由递推公式
，从
出发，计算
；
（2）
，
用
，计算
；
（3）分析结果的可靠性及产生此现象的原因（重点分析原因）。
解：
（1）程序如下:
clearall
clc
I=0.1822;%题中的已知数据
forn=1:
20;
I=（-5）*I+1/n;%由递推公式所得
end
fprintf（'I20=%f\n',I）
M=0.1823;%与I的计算结果形成对比
fori=1:
20;
M=（-5）*M+1/i;%由递推公式所得
end
fprintf（'M20=%f\n',M）
输出结果为：
I20=-11592559237.912731
M20=-2055816073.851284
（2）程序如下:
clearall
clc
I=0;%赋予I20的初始值
forn=0:
19;
I=（-1/5）*I+1/（5*（20-n））;%有递推公式得
end
fprintf（'I0=%f\n',I）
M=10000;
fori=0:
19;
M=（-1/5）*M+1/（5*（20-i））;%有递推公式得
end
fprintf（'M0=%f\n',M）
输出结果为：
I0=0.182322
M0=0.182322
（3）由输出结果可看出第一种算法为不稳定算法，第二中算法为稳定算法。
由于误差
第一种算法为正向迭代算法，每计算一步误差增长5倍，虽然所给的初始值很接近，随着n的增大，误差也越来越大。
第二种算法为倒向迭代算法，每计算一步误差缩小5倍，虽然所给的初始值之间差很多，随着n的增大，误差也越来越小。
2．求方程
的近似根，要求
，并比较计算量。
（1）在[0，1]上用二分法；
（2）取初值
，并用迭代
；
（3）加速迭代的结果；
（4）取初值
，并用牛顿迭代法；
（5）分析绝对误差。
解：
（1）程序如下（二分法）：
clearall
clc
a=0;b=1;
f=@（x）（exp（x）+10*x-2）;
%@是定义函数句柄的运算符
c=（a+b）/2;%取区间中点
i=0;%分割次数
whileabs（f（c））>5*10^（-4）
%判断f（x）的精度是否满足要求
iff（a）*f（c）<0
b=c;c=（a+b）/2;
elseiff（b）*f（c）<0
a=c;c=（b+a）/2;
end
i=i+1;
end
fprintf（'二分法运算次数为%i\n',i）
fprintf（'二分法计算结果为%f\n',c）
运行结果如下：
二分法运算次数为13
二分法计算结果为0.090515
（2）程序如下（不动点迭代）
clearall
clc
x0=0;
x=x0;
fork=1:
10000%规定迭代次数上限
y=（2-exp（x））/10;%迭代结果存到y中
ifabs（x-y）<5*10^（-4）
fprintf（'初始值x0为%i\n迭代次数为%i\n',x0,k）;
break
end
x=y;
ifk==10000;
fprintf（'迭代次数超出上限%i\n',k）;
end
end
fprintf（'迭代法计算结果为%f\n',y）;
运行结果为：
初始值x0为0
迭代次数为4
迭代法计算结果为0.090513
（3）程序如下（艾肯特迭代加速）
clearall
clc
x0=0;
x=x0;%给定迭代初值
p（1000）=0;%先定义p矩阵的长度
p
（1）=x;
fork=2:
1000;%规定迭代次数上限
y=（2-exp（x））/10;
z=（2-exp（y））/10;
x=x-（（y-x）^2）/（z-2*y+x）;
（4）程序如下（牛顿法迭代）
clearall
clc
x0=0;
x=x0;%给定迭代初值
fork=1:
1000;%规定迭代次数上限
y=x-（exp（x）+10*x-2）/（exp（x）+10）;
%牛顿计算公式
ifabs（y-x）<5*10^（-4）
fprintf（'初始值x0为%f\n迭代次数为%i\n',x0,k）;
%艾特肯加速公式
p（k）=x;%p是用来存储每一步迭代结果的矩阵
ifabs（p（k-1）-p（k））<5*10^（-4）
fprintf（'初始值x0为%f\n迭代次数为%i\n',x0,k-1）;
break
end
ifk==1000;
fprintf（'迭代次数超出上限%i\n',k）;
end
end
fprintf（'艾特肯加速迭代法计算结果为%f\n',x）;
运行结果为：
初始值x0为0.000000
迭代次数为2
艾特肯加速迭代法计算结果为0.090525
break
end
x=y;
ifk==1000;
fprintf（'迭代次数超出上限%i\n',k）;
end
end
fprintf（'牛顿迭代法计算结果为%f\n',y）;
运行结果为：
初始值x0为0.000000
迭代次数为2
牛顿迭代法计算结果为0.090525
（5）分析绝对误差
通过指令求得方程精确解的近似解：
>>solve（'exp（x）+10*x-2=0'）
ans=
1/5-lambertw（0,exp（1/5）/10）
>>1/5-lambertw（0,exp（1/5）/10）
ans=
0.090525101307255
各种计算方法的绝对误差为:
二分法的绝对误差：
1.0e-04*0.508986927450078
不动点迭代方法的绝对误差：
1.0e-04*0.121013072549997
艾特肯加速迭代的绝对误差：
1.0e-04*0.001013072550016
牛顿法的绝对误差:
1.0e-04*0.001013072550016
由上述各种算法计算出方程的值，二分法迭代了11次，收敛速度最慢，不动点迭代法迭代了4次，艾特肯迭代法迭代了2次，牛顿法迭代了2次，后两种方法的收敛速度都比较快，但计算量大。
由绝对误差可以看出二分法的误差最大，而艾特肯加速迭代和牛顿法误差最小。
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
试从中找出使用次数和容积之间的关系，计算均方差。
（用
拟合）
解：
用已知函数模型来拟合，程序如下：
先编辑函数
function[err]=f31（f,x,y）
a=f
（1）;
b=f
（2）;
c=f（3）;
err=0;
fork=1:
15
err=err+（（c+x（k））*y（k）-（a*x（k）+b））^2;
end
主函数如下：
x=[2345678910111213141516];
y=[6.428.29.589.59.7109.939.9910.4910.5910.610.810.610.910.76];
[f,fval,exitflag]=fminsearch（@f31,[0;0;0],optimset,x,y）
%fminsearch函数求解多元函数的最小值
%f返回多元函数y=f（x）在初始值x0附近的局部极小值点
%fval返回局部极小值
%exitflag返回函数输出条件值，exitflag=1说明函数收敛到一个解x；exitflag=0说明函数达到最大迭代次数；exitflag=-1说明输出函数终止算法。
fprintf（'\n拟合出来的方程式为\t（%7.4f+x）y=%7.4f+%7.4fx\n\n',f（3）,f
（2）,f
（1））;
z=linspace（x
（1）,x（end）,15）;
Y=（f
（1）*x+f
（2））./（f（3）+x）;
plot（x,y,'r*-',z,Y,'b-'）
legend（'y','Y'）;
title（'非线性的最小二乘拟合'）;
%均方差
fori=1:
15
yt（i）=（f
（1）*x（i）+f
（2））./（f（3）+x（i））;
end
c=0;
fori=1:
15
c=c+（y（i）-yt（i））^2;
end
jfc=sqrt（c/15）;
fprintf（'均方差为%f\n',jfc）
运行结果如下：
拟合出来的方程式为（-0.7110+x）y=-15.5024+11.2657x
均方差为0.331651
4．设
，
，
分析下列迭代法的收敛性，并求
的近似解及相应的迭代次数。
（1）JACOBI迭代；
（2）GAUSS-SEIDEL迭代；
（3）SOR迭代（取
，找到迭代步数最少的
）。
解：
（1）JACOBI迭代（程序如下）：
先编辑JACOBI函数：
functiontx=jacobi（A,b,imax,x0,tol）
%初始值x0,次数imax,精度tol
del=10^-10;
tx=[x0];n=length（x0）;
fori=1:
n
dg=A（i,i）;
ifabs（dg）disp（'对角元太小'）;%防止现溢出现象
return
end
end
fork=1:
imax%jacobi循环
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）;
%x
（1）到x（n）即完成一次迭代
end
tx=[tx;x];%矩阵中又加一行。
ifnorm（x-x0）return
else
x0=x;
end
end
主函数程序如下：
clearall
clc
A=[4-10-100;-14-10-10;0-14-10-1;-10-14-10;0-10-14-1;00-10-14];
b=[05-25-26];
x0=[000000];
imax=100;%迭代次数
tol=10^-4;%精度
tx=jacobi（A,b,imax,x0,tol）;
forj=1:
length（tx）
fprintf（'%d%f%f%f%f%f%f\n',j,tx（j,1）,tx（j,2）,tx（j,3）,tx（j,4）,tx（j,5）,tx（j,6））;
end
运行结果如下:
10.0000000.0000000.0000000.0000000.0000000.000000
20.0000001.250000-0.5000001.250000-0.5000001.500000
…
280.9999471.9999640.9999271.9999640.9999271.999973
290.9999821.9999500.9999751.9999500.9999751.999964
（2）GAUSS-SEIDEL迭代（程序如下）：
先编辑GAUSS-SEIDEL函数：
functiontx=gseidel（A,b,imax,x0,tol）
%初始值x0,次数imax,精度tol
del=10^-10;
tx=[x0];n=length（x0）;
%tx是个二维矩阵，存储着每一步迭代的结果
fori=1:
n
dg=A（i,i）;
ifabs（dg）disp（'对角元太小'）;
return
end
end
fork=1:
imax%G-S循环
x=x0;
fori=1:
n
sm=b（i）;
forj=1:
n
ifj~=i
sm=sm-A（i,j）*x（j）;
end
end
x（i）=sm/A（i,i）;
end
tx=[tx;x];
ifnorm（x-x0）return
else
x0=x;
end
end
主程序如下:
clearall
clc
A=[4-10-100;-14-10-10;0-14-10-1;-10-14-10;0-10-14-1;00-10-14];
b=[05-25-26];
x0=[000000];
imax=100;
tol=10^-4;
tx=gseidel（A,b,imax,x0,tol）;
forj=1:
length（tx）
fprintf（'%d%f%f%f%f%f%f\n',j,tx（j,1）,tx（j,2）,tx（j,3）,tx（j,4）,tx（j,5）,tx（j,6））
end
运行结果如下：
10.0000000.0000000.0000000.0000000.0000000.000000
20.0000001.250000-0.1875001.2031250.1132811.481445
…
150.9999171.9999100.9999241.9999310.9999431.999967
160.9999601.9999570.9999641.9999670.9999731.999984
（3）SOR迭代（程序如下）：
先编辑SOR函数：
functiontx=sor（A,b,imax,x0,tol,w）
%初始值x0,次数imax,精度tol
del=10^-10;
tx=[x0];n=length（x0）;
fori=1:
n
dg=A（i,i）;
ifabs（dg）disp（'对角元太小'）;
return
end
end
fork=1:
imax%SOR循环
x=x0;
fori=1:
n
sm=b（i）;
forj=1:
n
ifj~=i
sm=sm-A（i,j）*x（j）;
%先高斯迭代
end
end
x（i）=sm/A（i,i）;
x（i）=w*x（i）+（1-w）*x0（i）;
%比较高斯与超松弛迭代公式，补上省缺的项
end
tx=[tx;x];
ifnorm（x-x0）return
else
x0=x;
end
end
主程序如下：
clearall
clc
A=[4-10-100;-14-10-10;0-14-10-1;-10-14-10;0-10-14-1;00-10-14];
b=[05-25-26];
x0=[000000];
imax=100;
tol=10^-4;
w=1.2;%给定不同的w，w=0.1:
0.1:
1.9，找出使SOR迭代法收敛速度最快
tx=sor（A,b,imax,x0,tol,w）;
forj=1:
size（tx,1）
fprintf（'%d%f%f%f%f%f%f\n',j,tx（j,1）,tx（j,2）,tx（j,3）,tx（j,4）,tx（j,5）,tx（j,6））
end
运行结果如下:
10.0000000.0000000.0000000.0000000.0000000.000000
20.0000001.500000-0.1500001.4550000.2865001.840950
…
100.9999572.0000220.9999791.9999821.0000181.999992
111.0000101.9999980.9999962.0000110.9999971.999999
w值得范围在0~2之间，SOR迭代才会收敛，令w=0.1:
0.1:
1.9找出w*，使得SOR迭代的收敛速度最快。
w=0.1时需要迭代101次；w=1.9时需要迭代101次；w=1时需要迭代16次；w=0.9时需要迭代20次；w=1.1时需要迭代13次；w=1.3时需要迭代13次；w=1.2时需要迭代11次。
有这几次试代可得出w=1.2时迭代次数最少，SOR迭代法收敛速度最快。
总结:
由于A为严格对角占优矩阵，根据定理可知雅可比迭代和GS迭代都收敛，SOR迭代收敛的必要条件是0三种迭代算法的结果分析:
（1）JACOBI迭代是根据A分解成上三角、下三角和对角阵，将线性方程组的求解归结为对角方程组求解过程的重复，思路简单易于编程。
迭代次数比较多，收敛速度慢。
（2）GAUSS-SEIDEL迭代，是在JACOBI收敛的基础上将计算出来的未知量马上投入下一个迭代方程中使用，使得迭代的数度加快，效果更好。
但GS迭代与JACOBI迭代的收敛域并不能相互包含，所以不能相互代替。
但如果两者都收敛时，一般说GS迭代比JACOBI迭代的收敛速度快。
（3）SOR迭代为了加快收敛速度，在GS的基础上加入一修正参数即松弛因子w，使得收敛速度更快。
但选着不同的w，收敛速度不一样，为了达到最好的效果，要选着最合适的松弛因子。
5．用逆幂迭代法求
最接近于11的特征值和特征向量，准确到
。
解：
程序如下
先编辑函数文件：
function[l,v]=fanpower（A,v0,p,ep,max1）
%p为A最接近的特征值，ep为精度，max1为最大迭代次数
n=length（A）;
B=A-p*eye（n）;
k=0;err=1;
v0=v0';%求v0的转秩
while（（k<=max1）&（err>ep））
[mj]=max（abs（v0））;%找出v0中模最大的元素，j为此元素的序号
y=v0/（v0（j））;
x0=[010];
%解B*v1=y';求逆矩阵占用空间较大,用高斯迭代。
tx=gseidel（B,y',100,x0,10^-7）;
[rq]=size（tx）;%得出的tx是一个二维矩阵
v1=[tx（r,1）tx（r,2）tx（r,3）]';%提取出最后一行的三个元素
err=abs（1/norm（v1,2）-1/norm（v0,2））;
%norm（v,2）求矩阵v的2范数，即：
最大特征值
v0=v1;
k=k+1;
end
[mj]=max（abs（v0））;
l=p+（1/v0（j））;
v=y;
fprintf（'最接近给定数的特征值是l\n'）;
fprintf（'对应的特征向量为v'）;
主程序如下：
clearall
clc
A=[631;321;111];
v0=[111];
p=11;ep=0.001;max1=1000;
[l,v]=fanpower（A,v0,p,ep,max1）
运行结果如下:
最接近给定数的特征值是l
对应的特征向量为v
l=
7.8731
v=
1.0000
0.5493
0.2256
总结:
此方法为结合原点平移的反幂法，可以根据矩阵A的已知的特征值的粗略近似值p，计算出与数p最接近的特征值及特征向量。
6．用经典R-K方法求解初值问题
（1）
，
，
；
（2）
，
，
。
4解：
（1）程序如下：
clc;clear;
f=@（x,y1,y2）（-2*y1+y2+2*sin（x））;
g=@（x,y1,y2）（y1-2*y2+2*cos（x）-2*sin（x））;
h=0.1;
y1
（1）=2;y2
（1）=3;x
（1）=0;
fori=1:
100;
K1=f（x（i）,y1（i）,y2（i））;
L1=g（x（i）,y1（i）,y2（i））;
K2=f（x（i）+0.5*h,y1（i）+0.5*h*K1,y2（i）+0.5*h*L1）;
L2=g（x（i）+0.5*h,y1（i）+0.5*h*K1,y2（i）+0.5*h*L1）;
K3=f（x（i）+0.5*h,y1（i）+0.5*h*K2,y2（i）+0.5*h*L2）;
L3=g（x（i）+0.5*h,y1（i）+0.5*h*K2,y2（i）+0.5*h*L2）;
K4=f（x（i）+h,y1（i）+h*K3,y2（i）+h*L3）;
L4=g（x（i）+h,y1（i）+h*K3,y2（i）+h*L3）;
x（i+1）=x（i）+h;
y1（i+1）=y1（i）+h*（1/6）*（K1+2*K2+2*K3+K4）;
y2（i+1）=y2（i）+h*（1/6）*（L1+2*L2+2*L3+L4）;
end
x=0:
0.1:
10;
fori=1:
101;
Y1（i）=2*exp（-x（i））+sin（x（i））;
Y2（i）=2*exp（-x（i））+cos（x（i））;
end
c=y1-Y1;
d=y2-Y2;
subplot（2,1,1）
plot（x,y1,'r-',x,y2,'b--','LineWidth',4）
legend（'y1','y2'）;
title（'R--K四阶龙格库塔算法下方程组的解'）;
ylabel（'y1曲线y2曲线'）
subplot（2,1,2）
plot（x,c,'r-',x,d,'b--'）
legend（'c','d'）;
title（'误差值'）;
ylabel（'c曲线d曲线'）
[x;y1;c;y2;d]'
（2）程序如下：
clc;clear;
f=@（x,y1,y2）（-2*y1+y2+2*sin（x））;
g=@（x,y1,y2）（998*y1-999*y2+999*cos（x）-999*sin（x））;
h=0.001;
y1
（1）=2;y2
（1）=3;x
（1）=0;
fori=1:
10000;
K1=f（x（i）,y1（i）,y2（i））;
L1=g（x（i）,y1（i）,y2（i））;
K2=f（x（i）+0.5*h,y1（i）+0.5*h*K1,y2（i）+0.5*h*L1）;
L2=g（x（i）+0.5*h,y1（i）+0.5*h*K1,y2（i）+0.5*h*L1）;
K3=f（x（i）+0.5*h,y1（i）+0.5*h*K2,y2（i）+0.5*h*L2）;
L3=g（x（i）+0.5*h,y1（i）+0.5*h*K2,y2（i）+0.5*h*L2）;
K4=f（x（i）+h,y1（
升级会员 