数值计算方法上机实习题Word文档格式.docx
《数值计算方法上机实习题Word文档格式.docx》由会员分享，可在线阅读，更多相关《数值计算方法上机实习题Word文档格式.docx（21页珍藏版）》请在冰点文库上搜索。
易知第一个算法每一步计算都把误差放大了5倍，n次计算后更是放大了5n倍，可靠性低。
第二个算法可得出
可以看出第二个算法每一步计算就把误差缩小5倍，n次后缩小了5n倍，可靠性高。
2．求方程
的近似根，要求
，并比较计算量。
（1）在[0，1]上用二分法；
（1）[0，1]上的二分法
二分法子程序：
function[root,n]=EFF3（f,x1,x2）
%第二题
（1）二分法
f1=subs（f,symvar（f）,x1）;
%函数在x=x1的值
f2=subs（f,symvar（f）,x2）;
%x=x2
n=0;
%步数
er=5*10^-4;
%误差
if（f1==0）
root=x1;
return;
elseif（f2==0）
root=x2;
elseif（f1*f2>
0）
disp（'
两端点函数值乘积大于0！
'
）;
else
while（abs（x1-x2）>
er）%循环
x3=（x1+x2）/2;
f3=subs（f,symvar（f）,x3）;
n=n+1;
if（f3==0）
root=x3;
break;
elseif（f1*f3>
x1=x3;
else
x2=x3;
end
root=（x1+x2）/2;
%while循环少一步需加上
计算根与步数程序：
fplot（@（x）exp（x）+10*x-2,[0,1]）;
gridon;
symsx;
f=exp（x）+10*x-2;
[root,n]=EFF3（f,0,1）;
root=%6.8f,n=%d\n'
root,n）;
root=0.09057617,n=11
（2）取初值
，并用迭代
（2）初值x0=0迭代
迭代法子程序：
function[root,n]=DDF（g,x0,err,max）（接下页）
%root根，n+1步数，g函数，x0初值，err误差，max最大迭代次数
X
（1）=x0;
forn=2:
max
X（n）=subs（g,symvar（g）,X（n-1））;
c=abs（X（n）-X（n-1））;
root=X（n）;
if（c<
err）
ifn==max
超出预设迭代次数'
n=n-1;
%步数减1
f=（2-exp（x））/10;
（接下页）
x0=0;
err=5*10^（-4）;
max=100;
[root,n]=DDF（f,x0,err,max）;
root=0.09051262,n=4
（3）加速迭代的结果;
（3）加速迭代
加速迭代计算程序：
g=x-（f（x）-x）^2/（f（f（x））-2*f（x）-x）;
[root,n]=DDF（g,x0,err,max）;
root=%6.8f,n=%d'
程序函数设置：
functiong=f（x）
g=（2-exp（x））/10;
root=0.09048375,n=2
（4）取初值
，并用牛顿迭代法；
（4）牛顿迭代法
牛顿迭代法子程序：
function[root,n]=NDDDFx（g,x0,err,max）
err）||（root==0）
end
牛顿迭代法计算程序：
g=exp（x）+10*x-2;
g=x-（g/diff（g））;
[root,n]=NDDDFx（g,x0,err,max）;
root=%6.8f,n=%d\n'
root,n）；
root=0.09052511,n=2
（5）分析绝对误差。
可以看到，在同一精度下，二分法运算了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
11
12
13
14
15
16
10.49
10.59
10.60
10.8
10.6
10.9
10.76
试从中找出使用次数和容积之间的关系，计算均方差。
（用
拟合）
拟合曲线程序：
x=2:
16;
y=[6.428.29.589.59.7109.939.9910.4910.5910.6010.810.610.910.76];
[f,fval]=fminsearch（@delta1,[0,0,0],optimset,x,y）;
%fminsearch为求解多元函数的最小值函数
%f为多元函数初值x0附近的极小值点
%fval为极小值
k=2:
100;
y1=（f
（1）.*k+f
（2））./（f（3）+k）;
figure1=figure（'
Color'
[111]）;
gxt=plot（x,y,'
xr'
k,y1,'
k-'
legend（'
原数据'
'
拟合曲线'
Location'
northwest'
%y为数据点连接曲线，y1为拟合曲线
title（'
函数y=（ax+b）/（x+c）的拟合'
FontSize'
14,'
FontWeight'
Bold'
xlabel（'
次数'
ylabel（'
容积'
set（gxt,'
LineWidth'
1.5）;
%计算均方差
fori=1:
y2（i）=（f
（1）.*x（i）+f
（2））./（f（3）+x（i））;
j=0;
j=j+（y（i）-y2（i））^2;
jfc=sqrt（j/15）;
拟合出的方程为:
（x+%4.4f）y=%4.4fx+%4.4f\n均方差为:
%4.8f\n'
f（3）,f
（1）,f
（2）,jfc）;
构造函数子程序：
functiondelta=delta1（f,x,y）
a=f
（1）;
b=f
（2）;
c=f（3）;
delta=0;
fork=1:
delta=delta+（（x（k）+c）*y（k）-（a*x（k）+b））^2;
（x+-0.7110）y=11.2657x+-15.5024
均方差为:
0.33165089
总结：
指标选择，因题设方程
为非线性的，要转化为线性方程故需提指标为：
其驻点方程为:
4．设
分析下列迭代法的收敛性，并求
的近似解及相应的迭代次数。
（1）JACOBI迭代；
（1）Jacobi迭代
Jacobi迭代子程序：
function[x,k]=Jacobifl2（A,b,x0,eps,max1）
%jacobi按矩阵的分量算法
%x0初值，eps误差，max1最大迭代次数
n=length（x0）;
max1
x
（1）=（b
（1）-A（1,2:
n）*x0（2:
n））/A（1,1）;
fori=2:
n
x（i）=（b（i）-A（i,1:
i-1）*x0（1:
i-1）-A（i,i+1:
n）*x0（i+1:
n））/A（i,i）;
ifnorm（x'
-x0,2）<
eps
return
x0=x'
;
Jacobi迭代计算程序：
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]'
eps=10^（-4）;
max1=200;
[X,k]=Jacobifl2（A,b,x0,eps,max1）;
近似解x=\n[%4.8f\n%4.8f\n%4.8f\n%4.8f\n%4.8f\n%4.8f]\n迭代次数n=%d\n'
X
（1）,X
（2）,X（3）,X（4）,X（5）,X（6）,k）;
近似解x=
[0.99998177
1.99995018
0.99997509
1.99996353]
迭代次数n=28
（2）GAUSS-SEIDEL迭代；
（2）GAUSS-SEIDEL迭代
GAUSS-SEIDEL迭代子程序：
function[x,k]=GSDD（A,b,x0,eps,max1）
n）'
）/A（1,1）;
i-1）*x（1:
i-1）'
-A（i,i+1:
）/A（i,i）;
ifnorm（x-x0,2）<
x0=x;
GAUSS-SEIDEL迭代计算程序：
x0=[000000];
[X,k]=GSDD（A,b,x0,eps,max1）;
[0.99996014
1.99995676
0.99996351
1.99996662
0.99997253
1.99998401]
迭代次数n=15
（3）SOR迭代（取
，找到迭代步数最少的
）。
（3）SOR迭代
SOR迭代子程序：
function[x,k,m]=SOR1（A,b,x0,eps,max1）
a=x0;
fors=1:
19
x0=a;
w（s）=s/10;
x
（1）=x0
（1）+w（s）*（b
（1）-A（1,1:
n）*x0（1:
x（i）=x0（i）+w（s）*（b（i）-A（i,1:
-A（i,i:
n）*x0（i:
k（s）=0;
whilenorm（x-x0,2）>
=eps
k（s）=k（s）+1;
if（k（s）>
=max1）
超出迭代次数上限！
\n'
Data{s}=x;
[min1,index]=min（k）;
x=Data{index};
k=min1;
m=w（index）;
\n解x=\n%4.8f\n%4.8f\n%4.8f\n%4.8f\n%4.8f\n%4.8f'
x）;
\n步数最小omega=%2.2f,步数n=%d。
m,k）;
SOR迭代计算程序：
max1=1000;
[X,k,m]=SOR1（A,b,x0,eps,max1）;
\n近似解x=\n%4.8f\n%4.8f\n%4.8f\n%4.8f\n%4.8f\n%4.8f'
\n步数最小omega=%2.1f,步数n=%d\n'
1.00000977
1.99999782
0.99999567
2.00001071
0.99999662
1.99999923
步数最小omega=1.20,步数n=9。
Jacobi迭代次数为28次；
GAUSS-SEIDEL迭代次数为15次；
SOR迭代次数在w=1.2时达到最小次数9次。
系数矩阵严格对角占优，则Jacobi迭代法与Gauss-Seidel迭代法均收敛；
又因系数矩阵对称正定，且0＜w＜2，则SOR迭代法也收敛。
三种迭代法的结果分析及比较：
（1）Jacobi法，其设计思想是将线性方程组的求解归结为对角方程组求解过程的重复，计算规则简单而易于编写计算程序。
但收敛速度慢，谱半径相对偏大
（2）Gauss-Seidel充分利用新信息进行计算，其迭代的效果比Jacobi迭代好，谱半径比Jacobi法小。
（3）SOR迭代法，计算公式简单、编制程序容易，是求解大型稀疏方程组的一种有效方法，且若松弛因子选择得当超松弛法收敛速度比Gauss-Seidel迭代和Jacobi迭代都要快。
5．用逆幂迭代法求
最接近于11的特征值和特征向量，准确到
。
逆幂迭代法子程序：
function[lamda,V,k]=NMSF2（A,v,eps,p）
[n,n]=size（A）;
A=A-p*eye（n）;
v0=v;
[tmax,tindex]=max（abs（v0））;
lamd0=v0（tindex）;
u0=v0/lamd0;
flag=0;
k=0;
while（flag==0）
V=A\u0;
[tmax,tindex]=max（abs（V））;
lamd1=V（tindex）;
u0=V/lamd1;
if（abs（（lamd0）^（-1）-（lamd1）^（-1）））<
=eps
flag=1;
lamd0=lamd1;
k=k+1;
end（接下页）
lamda=（lamd1）^（-1）+p;
V=u0'
逆幂迭代计算程序：
A=[631;
321;
111];
v=[111]'
eps=0.001;
p=11;
[lamda,U,k]=NMSF2（A,v,eps,p）;
最接近11的特征值为:
%4.8f\n特征向量：
\n%4.8f\n%4.8f\n%4.8f\n'
lamda,U
（1）,U
（2）,U（3））;
最接近于11的特征值为：
7.87313712
特征向量：
1.00000000
0.54922698
0.22545618
6．用经典R-K方法求解初值问题
（1）
，
（2）
和精确解
比较，进行误差分析得到结论，图形显示精确解和数值解。
经典四阶R-K方法
R-K程序：
f=@（x,y）[-2*y
（1）+y
（2）+2*sin（x）;
y
（1）-2*y
（2）+2*cos（x）-2*sin（x）];
g=@（x,y）[-2*y
（1）+y
（2）+2*sin（x）;
998*y
（1）-999*y
（2）+
999*cos（x）-999*sin（x）];
x=0:
0.1:
10;
Y1=2*exp（-x）+sin（x）;
Y2=2*exp（-x）+cos（x）;
x1=0:
10/10000:
Y3=2*exp（-x1）+sin（x1）;
Y4=2*exp（-x1）+cos（x1）;
N=100;
N1=10000;
[a1,b1]=ode23（f,[0:
10],[23]）;
[a2,b2]=ode23（g,[0:
1/1000:
[a3,b3]=ode45（f,[0:
[a4,b4]=ode45（g,[0:
n=length（b1）;
dbt=plot（x,b1（:
1）,'
r-'
x,b1（:
2）,'
x,Y1,'
g--'
x,Y2,'
y--'
y1'
y2'
精确解Y1'
精确解Y2'
四阶RK算法下的方程组的数值解&
方程的精确解'
Y'
X'
set（dbt,'
Q=b1（:
1）-Y1'
W=b1（:
2）-Y2'
Q1=b2（:
1）-Y3'
W1=b2（:
2）-Y4'
Q2=b3（:
W2=b3（:
Q3=b4（:
W3=b4（:
ER1=max（abs（Q））;
ER2=max（abs（W））;
ER3=max（abs（Q1））;
ER4=max（abs（W1））;
ER5=max（abs（Q2））;
ER6=max（abs（W2））;
ER7=max（abs（Q3））;
ER8=max（abs（W3））;
\n
（1）中y1在ode23下求出的最大误差为：
%4.8f\n
（1）中y2在ode23下求出的最大误差为：
%4.8f\n'
ER1,ER2）;
\n
（2）中y1在ode23下求出的最大误差为：
%4.8f\n
（2）中y2在ode23下求出的最大误差为：
ER3,ER4）;
\n
（1）中y1在ode45下求出的最大误差为：
%4.8f\n
（1）中y2在ode45下求出的最大误差为：
ER5,ER6）;
\n
（2）中y1在ode45下求出的最大误差为：
%4.8f\n
（2）中y2在ode45下求出的最大误差为：
ER7,ER8）;
计算与显示程序：
（1）中y1在ode23下求出的最大误差为：
0.00169850
（1）中y2在ode23下求出的最大误差为：
0.00207422
（2）中y1在ode23下求出的最大误差为：
0.00000583
（2）中y2在ode23下求出的最大误差为：
0.00581329
（1）中y1在ode45下求出的最大误差为：
0.00052155
（1）中y2在ode45下求出的最大误差为：
0.00045793
（2）中y1在ode45下求出的最大误差为：
0.00000276
（2）中y2在ode45下求出的最大误差为：
0.00275331
对比2种方程在ode23，ode45命令下，与精确值的误差可以发现，在此题中，ode45的误差相对于ode23的误差小，较适合用在此类钢性问题中。
（1）Ode23的题1
（2）Ode3的题2
（3）
Ode45的题1
（4）Ode45的题2
7．用有限差分法求解边值问题（h=0.1），并
升级会员 