matlab上机实验报告.docx
《matlab上机实验报告.docx》由会员分享，可在线阅读，更多相关《matlab上机实验报告.docx（29页珍藏版）》请在冰点文库上搜索。
matlab上机实验报告
上海电力学院
数值计算方法上机实习报告
院系：
能源与机械工程学院
专业年级：
动力机械及工程2012级
*******
学号：
ys**********
*******
2012年12月26日
数值计算方法上机实习题
1．设
，
（1）由递推公式
，从
的几个近似值出发，计算
；
解：
I0=
=0.1823
计算I20编辑matlab命令如下：
I=0.1823
forn=1:
1:
20,
I=-5*I+1/n;
fprintf（'%.1d%.4f\n',n,I）;
end
结果：
（2）粗糙估计
，用
，计算
；
解：
I20=
使用复合中点公式进行积分，相应的matlab程序如下：
I=0;
forh=0:
0.001:
1,
m=h+0.0005;
I=I+0.001*m^20/（5+m）;
fprintf（'%.1d%.4f\n',m,I）;
end
disp（I）;
fork=1:
20,
n=21-k;
I=0.2*（1/n-I）;
fprintf（'%.1d%.4f\n',n,I）;
end
disp（I）
结果：
程序结束时输出两个I值，第一个表示I20，第二个表示I0；
分别为I20=0.0082
I0=0.1823
（3）分析结果的可靠性及产生此现象的原因（重点分析原因）。
从上述计算中分析得到如果先得到I0，再从I0由递推公式得到I20，I20结果跟精确值相比误差很大；如果先估算I20，在从I20有递推公式得到I0，I0的结果跟精确值相比近似相等。
原因分析：
如果从I0推I20的近似值，需要用到递推公式In=-5In-1+1/n，I0本身结果是有误差的；经过递推公式计算20次，就等于误差被认为的放大5的20次方倍，所以得到的I20与其精确值相差甚远。
如果从I20推I0的近似值，需要用到In-1=0.2（1/n-In），尽管I20本身有误差，但是经过20次运算，其误差缩小到原来的0.2的20次方倍，所以得到的I0与其精确值比较相近。
2．求方程
的近似根，要求
，并比较计算量。
（1）在[0，1]上用二分法；
Matlab程序如下：
a=0;
b=1;
c=b-a;
n=0
whilec>0.0005,
x=（a+b）/2;
f=exp（x）+10*x-2;
iff>0,
b=x;
c=b-a;
elseiff<0,
a=x;
c=b-a;
else
x=x;
c=0;
end
n=n+1;
fprintf（'%.1d%.4f%.4f\n',n,x,c）;
end
结果如下：
解得到;x=0.0903
（2）取初值
，并用迭代
；
采用matlab进行迭代的程序如下：
x=0;
c=1;
n=0;
whilec>0.0005,
m=x;
m=（2-exp（m））/10;
c=abs（m-x）;
x=m;
n=n+1;
fprintf（'%.1d%.4f%.4f\n',n,x,c）;
end
结果：
解得x=0.0905
（3）加速迭代的结果；
采用matlab进行迭代的程序如下：
x=0;
n=0;
a=0;
b=1;
whileabs（a-b）>0.0005,
n=n+1;
a=x;
y=（2-exp（x））/10;
z=（2-exp（y））/10;
x=x-（y-x）^2/（z-2*y+x）;
b=x;
fprintf（'%.1d%.4f%.4f\n',n,x,abs（a-b））;
end
结果如下：
（4）取初值
，并用牛顿迭代法；
Matlab程序如下：
x=0;
a=1;
n=0;
whileabs（a）>0.0005,
n=n+1;
a=（exp（x）+10*x-2）/（exp（x）+10）;
x=x-a;
fprintf（'%.1d%.4f%.4f\n',n,x,abs（a））;
end
运行结果：
（5）分析绝对误差。
迭代次数
二分法
代数式迭代
加速迭代
牛顿迭代
X（k）
Erroe
X（k）
Erroe
X（k）
Erroe
X（k）
Erroe
1
0.5000
0.5000
0.1000
0.1000
0.0905
0.0905
0.0909
0.0909
2
0.2500
0.2500
0.0895
0.0105
0.0905
0．0000
0.0905
0.0004
3
0.1250
0.1250
0.0906
0.0012
4
0.0625
0.0625
0.0905
0.0001
5
0.0938
0.0313
6
0.0781
0.0156
7
0.0859
0.0078
8
0.0898
0.0039
9
0.0918
0.0020
10
0.0908
0.0010
11
0.0903
0.0005
我们可以看到，在运算要求到同一精度的情况下，采用
（1）的二分法运算了11次，采用
（2）的方法运算了4次，采用（3）的加速迭代法运算了2次，采用（4）的牛顿迭代法也需运算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
试从中找出使用次数和容积之间的关系，计算均方差。
（注：
增速减少，用何种模型）
解：
将使用次数x与体积y的关系用matlab采用如下程序绘制在二维坐标系:
x=[2345678910111213141516];
y=[6.428.29.589.59.7109.939.9910.4910.5910.6010.810.610.910.76];
plot（x,y,'b*-'）;
结果如下：
由数据点分布图可知，拟合曲线y=f（x）随着x的增加而上升，但上升速度由快到慢，当x趋于无穷大时，y趋于某个常数，故曲线有一水平渐进线。
根据上述特征很容易想到用Logistic模型来拟合该曲线。
设y=f（x）的形式为y=aeb/x（a>0,b<0），两边取对数，得lny=lna+b/x。
记A=lna，B=b，并引入新变量m=lny，n=1/x。
引入新变量后的数据表如下：
x
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
n=1/x
0.5000
0.3333
0.2500
0.2000
0.1667
0.1429
0.1250
0.1111
0.1000
0.0909
0.0833
0.0769
0.0714
0.0667
0.0625
m=lny
1.8594
2.1041
2.2597
2.2513
2.2721
2.3026
2.2956
2.3016
2.3504
2.3599
2.3609
2.3795
2.3609
2.3888
2.3758
Matlab拟合程序如下：
n=[0.50000.33330.25000.20000.16670.14290.12500.11110.10000.09090.08330.07690.07140.06670.0625];
m=[1.85942.10412.25972.25132.27212.30262.29562.30162.35042.35992.36092.37952.36092.38882.3758];
polyfit（n,m,1）
运行的结果：
由此可得A=2.4578B=-1.1107
a=eA=11.6791b=B=-1.1107
由此得到使用次数与容积的函数为
将统计表和函数用matlab绘制在同一坐标系内程序如下：
x1=[2345678910111213141516];
y1=[6.428.29.589.59.7109.939.9910.4910.5910.6010.810.610.910.76];
x=2:
1:
16;
y=11.679*exp（-1.1107*x.^（-1））;
holdon;
plot（x,y,'ro-'）;
plot（x1,y1,'b*-'）;
结果如图：
计算均方差s，matlab程序如下：
y=[6.428.29.589.59.7109.939.9910.4910.5910.6010.810.610.910.76];
s=0;
forn=2:
1:
16,
a=abs（11.679*exp（-1.1107*n.^（-1））-y（n-1））;
s=s+a.^2;
end
s=（s/15）.^（1/2）;
disp（s）;
运算结果均方差S=0.2438
小结：
根据已给的条件计算函数是十分困难的，但通过对离散点的分析及变化规律找出其中的规律，并通过计算来得到实际的函数是十分有用的方法。
本题就是这样做的一个典型，在n=1/x和m=lny的基础上找到了它们之间的关系并通过这种关系来拟合原函数，并最终验证计算结果。
4．设
，
，
分析下列迭代法的收敛性，并求
的近似解及相应的迭代次数。
（1）JACOBI迭代；
解matlab计算程序如下：
A=[4-10-100;-14-10-10;0-14-10-1;-10-14-10;0-10-14-1;00-10-14];
b=[0;5;-2;5;-2;6];
error=1;
D=diag（diag（A））;
L=D-tril（A）;
U=D-triu（A）;
X=zeros（size（b））;
whileerror>0.0001,
X=D\（b+L*X+U*X）;
error=norm（b-A*X）/norm（b）;
end
disp（x）;
disp（error）;
解得X=[0.9999;1.9999;0.9998;1.9999;0.9998;1.9999]error=7.0206e-05
（2）GAUSS-SEIDEL迭代；
A=[4-10-100;-14-10-10;0-14-10-1;-10-14-10;0-10-14-1;00-10-14];
b=[0;5;-2;5;-2;6];
error=1;
D=diag（diag（A））;
L=D-tril（A）;
U=D-triu（A）;
X=zeros（size（b））;
whileerror>0.0001,
X=（D-L）\（b+U*X）;
error=norm（b-A*X）/norm（b）;
end
disp（x）;
disp（error）;
解得X=[0.9998;1.9998;0.9998;1.9999;0.9999;1.9999]error=5.5892e-05
（3）SOR迭代（
）。
●N=1.334使用matlab求解程序如下：
A=[4-10-100;-14-10-10;0-14-10-1;-10-14-10;0-10-14-1;00-10-14];
b=[0;5;-2;5;-2;6];
error=1;
D=diag（diag（A））;
L=D-tril（A）;
U=D-triu（A）;
X=zeros（size（b））;
whileerror>0.0001,
n=1.334;
X=（D-n*L）\[（1-n）*D+n*U]*X+n*[（D-n*L）\b];
error=norm（b-A*X）/norm（b）;
disp（X）;
end
disp（error）;
此循环得到的X=[0.9999;2.0000;1.0000;1.9999;1.0000;2.0000]error=6.3632e-05
●N=1.95使用matlab求解程序如下：
A=[4-10-100;-14-10-10;0-14-10-1;-10-14-10;0-10-14-1;00-10-14];
b=[0;5;-2;5;-2;6];
error=1;
D=diag（diag（A））;
L=D-tril（A）;
U=D-triu（A）;
X=zeros（size（b））;
whileerror>0.0001,
n=1.95;
X=（D-n*L）\[（1-n）*D+n*U]*X+n*[（D-n*L）\b];
error=norm（b-A*X）/norm（b）;
disp（X）;
end
disp（error）;
此循环得到的X=[0.9999;2.0001;0.9999;1.9999;1.0001;1.9999]error=9.0363e-05
●N=0.95使用matlab求解程序如下：
A=[4-10-100;-14-10-10;0-14-10-1;-10-14-10;0-10-14-1;00-10-14];
b=[0;5;-2;5;-2;6];
error=1;
D=diag（diag（A））;
L=D-tril（A）;
U=D-triu（A）;
X=zeros（size（b））;
whileerror>0.0001,
n=0.95;
X=（D-n*L）\[（1-n）*D+n*U]*X+n*[（D-n*L）\b];
error=norm（b-A*X）/norm（b）;
disp（X）;
end
disp（error）;
此循环得到的X=[0.9997;1.9997;0.9997;1.9998;0.9998;1.9999]error=8.6235e-05
5．用逆幂迭代法求
最接近于11的特征值和特征向量，准确到
。
解：
matlab程序如下：
a=[631;321;111];
I=[100;010;001];
b0=a-11*I;
v0=[1;1;1];
m=max（abs（v0））;
flab=1;
whileflab>0.001,
u=v0/m;
v0=b0\u;
[tv,ti]=max（abs（v0））;
n=v0（ti）;
flab=abs（n-m）;
m=n;
end
m=1/m+11;
disp（m）;
运行结果如下：
即离11最近的特征值为7.8745；相应的特征向量u=[1.0000;0.5503;0.2271]。
6．用经典R-K方法求解初值问题
（1）
，
，
；
解:
采用经典R-K公式计算的MATLAB程序如下：
y1=2;
y2=3;
forh=0:
0.01:
10,
k1=0.01*（-2*y1+y2+2*sin（h））;
l1=0.01*（y1-2*y2+2*cos（h）-2*sin（h））;
k2=0.01*（-2*（y1+0.5*k1）+（y2+0.5*l1）+2*sin（h+0.005））;
l2=0.01*（（y1+0.5*k1）-2*（y2+0.5*l1）+2*cos（h+0.005）-2*sin（h+0.005））;
k3=0.01*（-2*（y1+0.5*k2）+（y2+0.5*l2）+2*sin（h+0.005））;
l3=0.01*（（y1+0.5*k2）-2*（y2+0.5*l2）+2*cos（h+0.005）-2*sin（h+0.005））;
k4=0.01*（-2*（y1+k3）+（y2+l3）+2*sin（h+0.01））;
l4=0.01*（（y1+k3）-2*（y2+l3）+2*cos（h+0.01）-2*sin（h+0.01））;
y1=y1+1/6*（k1+2*k2+2*k3+k4）;
y2=y2+1/6*（l1+2*l2+2*l3+l4）;
ifh==fix（h）;
fprintf（'%.1d%.4f%.4f\n',h,y1,y2）;
else
end
end
结果如下所示：
（2）
，
，
。
和精确解
比较，分析结论。
Matlab程序如下：
y1=2;
y2=3;
forh=0:
0.00001:
10,
k1=0.00001*（-2*y1+y2+2*sin（h））;
l1=0.00001*（998*y1-999*y2+999*cos（h）-999*sin（h））;
k2=0.00001*（-2*（y1+0.5*k1）+（y2+0.5*l1）+2*sin（h+0.000005））;
l2=0.00001*（998*（y1+0.5*k1）-999*（y2+0.5*l1）+999*cos（h+0.000005）-999*sin（h+0.000005））;
k3=0.00001*（-2*（y1+0.5*k2）+（y2+0.5*l2）+2*sin（h+0.005））;
l3=0.00001*（998*（y1+0.5*k2）-999*（y2+0.5*l2）+999*cos（h+0.000005）-999*sin（h+0.000005））;
k4=0.00001*（-2*（y1+k3）+（y2+l3）+2*sin（h+0.00001））;
l4=0.00001*（998*（y1+k3）-999*（y2+l3）+999*cos（h+0.00001）-999*sin（h+0.00001））;
y1=y1+1/6*（k1+2*k2+2*k3+k4）;
y2=y2+1/6*（l1+2*l2+2*l3+l4）;
ifh==fix（h）,
fprintf（'%.1d%.4f%.4f\n',h,y1,y2）;
else
end
end
结果如下：
精确解：
forx=0:
1:
10,
y1=2*exp（-x）+sin（x）;
y2=2*exp（-x）+cos（x）;
fprintf（'%.1d%.4f%.4f\n',x,y1,y2）;
end
结果;
结果分析：
四阶RungeKutta方法得到的结果已很接近精确解，证明这种迭代方法精确度很好，是一种有效的算法。
但是要注意龙格-库塔公式的推导基于泰勒展开方法，因而它要求所求的的解具有较好的光滑性质。
反之，如果解得光滑性差，那么，使用四阶龙格-库塔求得的数值解精度就不是太高，此种情况可以采用缩小步长来解决，比如上述计算。
7．用有限差分法求解边值问题（h=0.1）：
.
微分方程式可以变为
用有限差分法matlab程序如下：
h=0.1;
n=2/0.1-1;
g
（1）=1/（h.^2）;
g（n）=1/（h^2）;
fori=2:
1:
18,
g（i）=0;
end
g=[g
（1）;g
（2）;g（3）;g（4）;g（5）;g（6）;g（7）;g（8）;g（9）;g（10）;g（11）;g（12）;g（13）;g（14）;g（15）;g（16）;g（17）;g（18）;g（19）];
disp（g）;
fori=1:
1:
19,
forj=1:
1:
19,
ifi==1,
H（1,1）=2/（h.^2）+（1+（-1+0.1*i）.^2）;
H（1,2）=-1/（h.^2）;
elseifi==19,
H（19,18）=-1/（h.^2）;
H（19,19）=2/（h.^2）+（1+（-1+0.1*i）.^2）;
else
ifj==i,
H（i,j）=2/（h.^2）+（1+（-1+0.1*i）.^2）;
elseifj==i-1,
H（i,i-1）=-1/（h.^2）;
elseifj==i+1;
H（i,i+1）=-1/（h.^2）;
else
H（i,j）=0;
end
end
end
end
disp（H）;
y=H\g;
fori=1:
1:
19,
fprintf（'%.4f%.4f\n',-1+0.1*i,y（i,1））
end
运算结果为：
H矩阵为：
Y在各点的近似值为：
X
Y
8.用函数y=asin（bx）拟合数据.
x
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
y
0.6
1.1
1.6
1.8
2.0
1.9
1.7
1.3
Matlab上机程序为：
function[err,a,b]=nlfitb（x,y）
ifnargin<2,
x=[1:
8]'/10;
y=[0.61.11.61.82.01.91.71.3]';
end
beta0=[11]';
beta=nlinfit（x,y,@mymodel,beta0）;
fprintf（'Thenonlinearleastsquarefittingy=a*sin（b*x）fordata\n\n'）;
fprintf（'%6.1f',x）;
fprintf（'%6.1f',y）;
fprintf（'\n\nis\n\nty=%7.4f*sin（%7.4f*x）\n\n',beta）;
z=linspace（x
（1）,x（end）,100）;
plot（x,y,'ro',z,beta
（1）*sin（beta
（2）*z）,'b-.'）;
functionyb=mymodel（beta,xb）
yb=beta
（1）*sin（beta
（2）*xb）;
计算结果：
9．拟合形如f（x）≈（a+bx）/（1+cx）的函数的一种快速方法是将最小二乘法用于下列问题：
f（x）（1+cx）≈（a+bx），试用这一方
升级会员 