我的数值计算上机报告-上海电力学院
相关主题
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
2.0000
1.0000
2.0000
n =
17
可见w取1.334时，n最小，迭代次数最少。
上机心得：从题中可以看出，系数矩阵严格对角占优，则Jacobi迭代法与Gauss-Seidel迭代法均收敛；又因系数矩阵对称正定，且0＜w＜2，则SOR迭代法也收敛。
三种迭代法的结果各有特色：（1）Jacobi法，其设计思想是将线性方程组的求解归结为对角方程组求解过程的重复，计算规则简单而易于编写汁算程序。但，我们通常要设置最大迭代次数N，以防止迭代过程不收敛或者收敛速度过于缓慢。（2）Gauss-Seidel充分利用新信息进行计算，其迭代的效果比Jacobi迭代好。（3）SOR迭代法， 计算公式简单、编制程序容易，是求解大型稀疏方程组的一种有效方法，且超松弛法收敛速度比Gauss-Seidel迭代和Jacobi迭代都要快。 从对松弛因子三种不同选择的分析发现，若松弛因子选择得当，能显著地提高收敛速度。
b=[0;5;-2;5;-2;6];
esp=0.0001;
x0=[0;0;0;0;0;0];
L=-tril(A,-1);
U=-triu(A,1);
D=diag(diag(A));
B=inv(D)*(L+U);
f=inv(D)*b;
J=B*x0+f;
n=1;
while norm(J-x0)>=esp
x0=J;
if abs(x(i)-x(i-1))<r
return;
end
end
disp(x(i))
disp(i)
其运算结果为：
0.0905
5
(3)加速迭代：
构造序列{ }，使
= ，其中 是x的一个近似值，经计算：
=0； =0.1； =0.0894829； =0.090639； =0.0905126；
=0.090484； =0.090524449； =0.090525536。
plot(x,y,'*')
holdon;
plot(x,f(x),'-')
可以得到散点与拟合曲线图如下所示：
sum=(11.6791*exp(-1.1107/x(1))-y(1))^2;
for i=1:14
sum=sum+(11.6791*exp(-1.1107/x(i+1))-y(i+1))^2 ;
r=0.00005;
x=0;
for i=1:20
x=x-(exp(x)+10*x-2)/(exp(x)+10);
if abs((exp(x)+10*x-2)/(exp(x)+10))<r;
return;
end
end
disp(x)
disp(i)
运算结果为：
0.0905
2
（5）
不同的计算方式是有收敛速度上的不同的。由本题可以看到，在要求同一运算精度的情况下，采用二分法运算了12次，采用题设的迭代方法运算了5次，采用加速迭代法运算了3次，采用牛顿迭代法则只需运算2次。
function y=sor(a,b,w,x0)
D=diag(diag(a));
U=-triu(a,1);
L=-tril(a,-1);
lw=(D-w*L)\((1-w)*D+w*U);
f=(D-w*L)\b*w;
y=lw*x0+f;n=1;
while norm(y-x0)>10^(-4)
x0=y;
y=lw*x0+f;n=n+1;
1.8594
2.1041
2.2597
2.2513
2.2721
2.3026
2.2956
2.3016
10
11
12
13
14
15
16
0.1000
0.0909
0.0833
0.0769
0.0714
0.0667
0.0625
2.3504
2.3599
2.3609
2.3795
2.3609
2.3888
2.3758
程序：
t=[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];
I1=0.18232;
I2=0.1823;
I3=0.18
for n=1:1:20
I1=-5*I1+1/n;
I2=-5*I2+1/n;
I3=-5*I3+1/n;
fprintf(' %.4f\n',I1,I2,I3)
end
I1
I2
I3
结果为：I1 = -1.4847e+008；I2 = -2.0558e+009；I3 = -2.2140e+011
k=1;
while abs(b-a)>r
k=k+1;
x=(a+b)/2;
y=exp(x)+10*x-2;
if y1*y<0
b=x
else
a=x
end
end
disp(x)
disp(k)
其运算结果为：
b =
0.5000
b =
0.2500
b =
0.1250
a =
0.0625
b =
0.0938
a =
0.0781
U1 =
0 0.2500 0 0.2500 0 0
0 0 0.2500 0 0.2500 0
0 0 0 0.2500 0 0.2500
0 0 0 0 0.2500 0
0 0 0 0 0 0.2500
0 0 0 0 0 0
1.0000
2.0000
1.0000
2.0000
1.0000
2.0000
15
（3）SOR迭代法(ω依次取1.334,1.95,0.95)，程序如下：
z=[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];
polyfit(t,z,1)
运行结果：
ans =
-1.1107 2.4578
由此可得A=2.4578，B=-1.1107， ，b=B=-1.1107
for i=1:3
w=c(i);
sor(a,b,w,x0);
end
运行结果分别为：
y =
1.0000
2.0000
1.0000
2.0000
1.0000
2.0000
n =
13
y =
1.0000
2.0000
1.0000
2.0000
1.0000
2.0000
n =
241
y =
1.0000
2.0000
1.0000
b=[0;5;-2;5;-2;6];
x0=[0;0;0;0;0;0]
I=eye(6,6)
L=-tril(A,-1);
U=-triu(A,1);
D=diag(diag(A));
L1=inv(D)*L
U1=inv(D)*U
B=inv(D)*(L+U);
f=inv(D)*b;
J=inv(I-L1)*U1*x0+inv(I-L1)*f;
方程即为
（2）校核方差
x=[2 3 4 5 6 7 8 9 10 11 12 13 14 15 16];
y=[6.42 8.20 9.58 9.50 9.70 10.00 9.93 9.99 10.49 10.59 10.60 10.80 10.60 10.90 10.76];
f(x)=11.6791*exp(-1.1107./x);
J=B*x0+f;
n=n+1;
end
disp(J);
disp(n)
运行的结果如下：
1.0000
2.0000
1.0000
2.0000
1.0000
2.0000
28
（2）用Gauss-Seidel迭代，程序如下：
esp=0.0001;
A=[4,-1,0,-1,0,0;-1,4,-1,0,-1,0;0,-1,4,-1,0,-1;-1,0,-1,4,-1,0;0,-1,0,-1,4,-1;0,0,-1,0,-1,4];
fc=sum^0.5;
end
fc
运行结果：
fc =
0.9441
所以均方差是0.9441
4．设 ， ，
分析下列迭代法的收敛性，并求 的近似解及相应的迭代次数。
（1）JACOBI迭代；
（2）GAUSS-SEIDEL迭代；
（3）SOR迭代（ ）。
解：（1）
A=[4,-1,0,-1,0,0;-1,4,-1,0,-1,0;0,-1,4,-1,0,-1;-1,0,-1,4,-1,0;0,-1,0,-1,4,-1;0,0,-1,0,-1,4];
编程如下：
x=0;A=0;B=1;n=0
whileabs(B-A)>5*1e-4
A=x;
y=exp(x)+10*x-2;
z=exp(y)+10*y-2;
x=x-(y-x)^2/(z-2*y+x);
B=x;
n=n+1
end
x
n
结果为：x =0.0995;n=3
(4)用牛顿法，取初值x=0，编写程序代码为下:
5．用逆幂迭代法求 最接近于11的特征值和特征向量，准确到 。
解：A=[6 3 1;3 2 1;1 1 1];
[x,lamda]=eig(A),
v=[1 1 1]';
q=11,
B=inv(A-q*eye(3));
for k=1:20
z=B*v;
m=max(z);
v=z/m;
fprintf('%6d %f %f %f %f\n',k,v(1),v(2),v(3),m^-1+q),
end
y
n
以文件名sor.m保存。
程序：
a=[4 -1 0 -1 0 0;-1 4 -1 0 -1 0;0 -1 4 -1 0 -1;-1 0 -1 4 -1 0;0 -1 0 -1 4 -1;0 0 -1 0 -1 4];
b=[0 5 -2 5 -2 6]';
x0=[0 0 0 0 0 0]';
c=[1.334 1.95 0.95];
数据点分布图
设y=f(x)具有指数形式 （a>0,b<0）。对此式两边取对数，得 。记A=lna，B=b，并引入新变量z=lny，t=1/x。引入新变量后的数据表如下
x
2
3
4
5
6
7
8
9
t=1/x
0.5000
0.3333
0.2500
0.2000
0.1667
0.1429
0.1250
0.1111
z=lny
0 0 1 0 0 0
0 0 0 1 0 0
0 0 0 0 1 0
0 0 0 0 0 1
L1 =
0 0 0 0 0 0
0.2500 0 0 0 0 0
0 0.2500 0 0 0 0
0.2500 0 0.2500 0 0 0
0 0.2500 0 0.2500 0 0
0 0 0.2500 0 0.2500 0
分析：精度不同时计算出来的结果是不同的
（2）
>> i=0.009167 %粗糙估计i20值
i =
0.0092
>> for n=20:-1:1
i=-i/5+1/5/n;
end
>> disp(i)
0.1823
（3）对两种递推公式的误差进行分析：设第一递推式中开始时的误差为 ，递推过程的舍入误差不计。并记 ，则有 。由于 ，因此递推式不可靠。而在第二种递推式中 ，误差在缩小，所以此递推式是可靠的。
n=1;
while norm(J-x0)>=esp
x0=J;
J=inv(I-L1)*U1*x0+inv(I-L1)*f;
n=n+1;
end
if n>400
disp('该方程不收敛')
Biblioteka Baiduelse
disp(J);
disp(n)
end
运行结果如下：
x0 =
0
0
0
0
0
0
I =
1 0 0 0 0 0
0 1 0 0 0 0
数值计算方法上机
实习报告
专业年级：
学生姓名：
学号：
2012年12月25日
1．设 ，
（1）由递推公式 ，从 的几个近似值出发，计算 ；
（2）给出 的近似值，用 ，计算 ；
（3）分析结果的可靠性及产生此现象的原因(重点分析原因)。
解：（1）
ln6-ln5=0.18232,
将 代入递推公式 中，MATLAB程序如下:
a =
0.0859
a =
0.0898
b =
0.0918
b =
0.0908
a =
0.0903
0.0903
12
(2)用题设迭代法，取初值x(1)=0，并用迭代x(i)=(2-exp(x(i-1)))/10，程序代码如下:
x(1)=0;
r=0.0005;
for i=2:20
x(i)=(2-exp(x(i-1)))/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
试从中找出使用次数和容积之间的关系，计算均方差。（注：增速减少，用何种模型）
解：（1）由数据点分布图可知，拟合曲线y=f(x)随着x的增加而上升，但上升速度由快到慢，当x趋于无穷大时，y趋于某个常数，故曲线有一水平渐进线。根据上述特征很容易想到用Logistic模型来拟合该曲线。
2．求方程 的近似根，要求 ，并比较计算量。
（1）在[0，1]上用二分法；
（2）取初值 ，并用迭代 ；
（3）加速迭代的结果；
（4）取初值 ，并用牛顿迭代法；
（5）分析绝对误差。
解：（1）用二分法，编写程序代码为下:
r=0.0005;
a=0;
b=1;
y1=exp(a)+10*a-2;
y2=exp(b)+10*b-2;
1.0000
2.0000
n =
17
可见w取1.334时，n最小，迭代次数最少。
上机心得：从题中可以看出，系数矩阵严格对角占优，则Jacobi迭代法与Gauss-Seidel迭代法均收敛；又因系数矩阵对称正定，且0＜w＜2，则SOR迭代法也收敛。
三种迭代法的结果各有特色：（1）Jacobi法，其设计思想是将线性方程组的求解归结为对角方程组求解过程的重复，计算规则简单而易于编写汁算程序。但，我们通常要设置最大迭代次数N，以防止迭代过程不收敛或者收敛速度过于缓慢。（2）Gauss-Seidel充分利用新信息进行计算，其迭代的效果比Jacobi迭代好。（3）SOR迭代法， 计算公式简单、编制程序容易，是求解大型稀疏方程组的一种有效方法，且超松弛法收敛速度比Gauss-Seidel迭代和Jacobi迭代都要快。 从对松弛因子三种不同选择的分析发现，若松弛因子选择得当，能显著地提高收敛速度。
b=[0;5;-2;5;-2;6];
esp=0.0001;
x0=[0;0;0;0;0;0];
L=-tril(A,-1);
U=-triu(A,1);
D=diag(diag(A));
B=inv(D)*(L+U);
f=inv(D)*b;
J=B*x0+f;
n=1;
while norm(J-x0)>=esp
x0=J;
if abs(x(i)-x(i-1))<r
return;
end
end
disp(x(i))
disp(i)
其运算结果为：
0.0905
5
(3)加速迭代：
构造序列{ }，使
= ，其中 是x的一个近似值，经计算：
=0； =0.1； =0.0894829； =0.090639； =0.0905126；
=0.090484； =0.090524449； =0.090525536。
plot(x,y,'*')
holdon;
plot(x,f(x),'-')
可以得到散点与拟合曲线图如下所示：
sum=(11.6791*exp(-1.1107/x(1))-y(1))^2;
for i=1:14
sum=sum+(11.6791*exp(-1.1107/x(i+1))-y(i+1))^2 ;
r=0.00005;
x=0;
for i=1:20
x=x-(exp(x)+10*x-2)/(exp(x)+10);
if abs((exp(x)+10*x-2)/(exp(x)+10))<r;
return;
end
end
disp(x)
disp(i)
运算结果为：
0.0905
2
（5）
不同的计算方式是有收敛速度上的不同的。由本题可以看到，在要求同一运算精度的情况下，采用二分法运算了12次，采用题设的迭代方法运算了5次，采用加速迭代法运算了3次，采用牛顿迭代法则只需运算2次。
function y=sor(a,b,w,x0)
D=diag(diag(a));
U=-triu(a,1);
L=-tril(a,-1);
lw=(D-w*L)\((1-w)*D+w*U);
f=(D-w*L)\b*w;
y=lw*x0+f;n=1;
while norm(y-x0)>10^(-4)
x0=y;
y=lw*x0+f;n=n+1;
1.8594
2.1041
2.2597
2.2513
2.2721
2.3026
2.2956
2.3016
10
11
12
13
14
15
16
0.1000
0.0909
0.0833
0.0769
0.0714
0.0667
0.0625
2.3504
2.3599
2.3609
2.3795
2.3609
2.3888
2.3758
程序：
t=[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];
I1=0.18232;
I2=0.1823;
I3=0.18
for n=1:1:20
I1=-5*I1+1/n;
I2=-5*I2+1/n;
I3=-5*I3+1/n;
fprintf(' %.4f\n',I1,I2,I3)
end
I1
I2
I3
结果为：I1 = -1.4847e+008；I2 = -2.0558e+009；I3 = -2.2140e+011
k=1;
while abs(b-a)>r
k=k+1;
x=(a+b)/2;
y=exp(x)+10*x-2;
if y1*y<0
b=x
else
a=x
end
end
disp(x)
disp(k)
其运算结果为：
b =
0.5000
b =
0.2500
b =
0.1250
a =
0.0625
b =
0.0938
a =
0.0781
U1 =
0 0.2500 0 0.2500 0 0
0 0 0.2500 0 0.2500 0
0 0 0 0.2500 0 0.2500
0 0 0 0 0.2500 0
0 0 0 0 0 0.2500
0 0 0 0 0 0
1.0000
2.0000
1.0000
2.0000
1.0000
2.0000
15
（3）SOR迭代法(ω依次取1.334,1.95,0.95)，程序如下：
z=[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];
polyfit(t,z,1)
运行结果：
ans =
-1.1107 2.4578
由此可得A=2.4578，B=-1.1107， ，b=B=-1.1107
for i=1:3
w=c(i);
sor(a,b,w,x0);
end
运行结果分别为：
y =
1.0000
2.0000
1.0000
2.0000
1.0000
2.0000
n =
13
y =
1.0000
2.0000
1.0000
2.0000
1.0000
2.0000
n =
241
y =
1.0000
2.0000
1.0000
b=[0;5;-2;5;-2;6];
x0=[0;0;0;0;0;0]
I=eye(6,6)
L=-tril(A,-1);
U=-triu(A,1);
D=diag(diag(A));
L1=inv(D)*L
U1=inv(D)*U
B=inv(D)*(L+U);
f=inv(D)*b;
J=inv(I-L1)*U1*x0+inv(I-L1)*f;
方程即为
（2）校核方差
x=[2 3 4 5 6 7 8 9 10 11 12 13 14 15 16];
y=[6.42 8.20 9.58 9.50 9.70 10.00 9.93 9.99 10.49 10.59 10.60 10.80 10.60 10.90 10.76];
f(x)=11.6791*exp(-1.1107./x);
J=B*x0+f;
n=n+1;
end
disp(J);
disp(n)
运行的结果如下：
1.0000
2.0000
1.0000
2.0000
1.0000
2.0000
28
（2）用Gauss-Seidel迭代，程序如下：
esp=0.0001;
A=[4,-1,0,-1,0,0;-1,4,-1,0,-1,0;0,-1,4,-1,0,-1;-1,0,-1,4,-1,0;0,-1,0,-1,4,-1;0,0,-1,0,-1,4];
fc=sum^0.5;
end
fc
运行结果：
fc =
0.9441
所以均方差是0.9441
4．设 ， ，
分析下列迭代法的收敛性，并求 的近似解及相应的迭代次数。
（1）JACOBI迭代；
（2）GAUSS-SEIDEL迭代；
（3）SOR迭代（ ）。
解：（1）
A=[4,-1,0,-1,0,0;-1,4,-1,0,-1,0;0,-1,4,-1,0,-1;-1,0,-1,4,-1,0;0,-1,0,-1,4,-1;0,0,-1,0,-1,4];
编程如下：
x=0;A=0;B=1;n=0
whileabs(B-A)>5*1e-4
A=x;
y=exp(x)+10*x-2;
z=exp(y)+10*y-2;
x=x-(y-x)^2/(z-2*y+x);
B=x;
n=n+1
end
x
n
结果为：x =0.0995;n=3
(4)用牛顿法，取初值x=0，编写程序代码为下:
5．用逆幂迭代法求 最接近于11的特征值和特征向量，准确到 。
解：A=[6 3 1;3 2 1;1 1 1];
[x,lamda]=eig(A),
v=[1 1 1]';
q=11,
B=inv(A-q*eye(3));
for k=1:20
z=B*v;
m=max(z);
v=z/m;
fprintf('%6d %f %f %f %f\n',k,v(1),v(2),v(3),m^-1+q),
end
y
n
以文件名sor.m保存。
程序：
a=[4 -1 0 -1 0 0;-1 4 -1 0 -1 0;0 -1 4 -1 0 -1;-1 0 -1 4 -1 0;0 -1 0 -1 4 -1;0 0 -1 0 -1 4];
b=[0 5 -2 5 -2 6]';
x0=[0 0 0 0 0 0]';
c=[1.334 1.95 0.95];
数据点分布图
设y=f(x)具有指数形式 （a>0,b<0）。对此式两边取对数，得 。记A=lna，B=b，并引入新变量z=lny，t=1/x。引入新变量后的数据表如下
x
2
3
4
5
6
7
8
9
t=1/x
0.5000
0.3333
0.2500
0.2000
0.1667
0.1429
0.1250
0.1111
z=lny
0 0 1 0 0 0
0 0 0 1 0 0
0 0 0 0 1 0
0 0 0 0 0 1
L1 =
0 0 0 0 0 0
0.2500 0 0 0 0 0
0 0.2500 0 0 0 0
0.2500 0 0.2500 0 0 0
0 0.2500 0 0.2500 0 0
0 0 0.2500 0 0.2500 0
分析：精度不同时计算出来的结果是不同的
（2）
>> i=0.009167 %粗糙估计i20值
i =
0.0092
>> for n=20:-1:1
i=-i/5+1/5/n;
end
>> disp(i)
0.1823
（3）对两种递推公式的误差进行分析：设第一递推式中开始时的误差为 ，递推过程的舍入误差不计。并记 ，则有 。由于 ，因此递推式不可靠。而在第二种递推式中 ，误差在缩小，所以此递推式是可靠的。
n=1;
while norm(J-x0)>=esp
x0=J;
J=inv(I-L1)*U1*x0+inv(I-L1)*f;
n=n+1;
end
if n>400
disp('该方程不收敛')
Biblioteka Baiduelse
disp(J);
disp(n)
end
运行结果如下：
x0 =
0
0
0
0
0
0
I =
1 0 0 0 0 0
0 1 0 0 0 0
数值计算方法上机
实习报告
专业年级：
学生姓名：
学号：
2012年12月25日
1．设 ，
（1）由递推公式 ，从 的几个近似值出发，计算 ；
（2）给出 的近似值，用 ，计算 ；
（3）分析结果的可靠性及产生此现象的原因(重点分析原因)。
解：（1）
ln6-ln5=0.18232,
将 代入递推公式 中，MATLAB程序如下:
a =
0.0859
a =
0.0898
b =
0.0918
b =
0.0908
a =
0.0903
0.0903
12
(2)用题设迭代法，取初值x(1)=0，并用迭代x(i)=(2-exp(x(i-1)))/10，程序代码如下:
x(1)=0;
r=0.0005;
for i=2:20
x(i)=(2-exp(x(i-1)))/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
试从中找出使用次数和容积之间的关系，计算均方差。（注：增速减少，用何种模型）
解：（1）由数据点分布图可知，拟合曲线y=f(x)随着x的增加而上升，但上升速度由快到慢，当x趋于无穷大时，y趋于某个常数，故曲线有一水平渐进线。根据上述特征很容易想到用Logistic模型来拟合该曲线。
2．求方程 的近似根，要求 ，并比较计算量。
（1）在[0，1]上用二分法；
（2）取初值 ，并用迭代 ；
（3）加速迭代的结果；
（4）取初值 ，并用牛顿迭代法；
（5）分析绝对误差。
解：（1）用二分法，编写程序代码为下:
r=0.0005;
a=0;
b=1;
y1=exp(a)+10*a-2;
y2=exp(b)+10*b-2;