丁丽娟《数值计算方法》五章课后实验题答案源程序很详细且运行无误Word格式.docx
《丁丽娟《数值计算方法》五章课后实验题答案源程序很详细且运行无误Word格式.docx》由会员分享，可在线阅读，更多相关《丁丽娟《数值计算方法》五章课后实验题答案源程序很详细且运行无误Word格式.docx（28页珍藏版）》请在冰豆网上搜索。
u（k,j）=A（k,j）-sum;
A（k,j）=u（k,j）;
fori=1:
l（i,i）=1;
fori=（k+1）:
sum=sum+A（i,r）*u（r,k）;
l（i,k）=（A（i,k）-sum）/u（k,k）;
A（i,k）=l（i,k）;
end
求A的列主元素三角分解：
>
A=[11111;
12345;
1361015;
14102035;
15153570];
[L,U]=myfun（A）
结果：
L=
1.00000000
1.00001.0000000
1.00000.50001.000000
1.00000.75000.75001.00000
1.00000.25000.7500-1.00001.0000
U=
1.00001.00001.00001.00001.0000
04.000014.000034.000069.0000
00-2.0000-8.0000-20.5000
000-0.5000-2.3750
0000-0.2500
（2）
求矩阵的逆矩阵A-1:
inv（A）
结果为：
ans=
5-1010-51
-1030-3519-4
10-3546-276
-519-2717-4
1-46-41
（3）
检验结果：
E=diag（[11111]）
A\E
2:
程序：
functiond=myfun（a,b,c,d,n）
fori=2:
l（i）=a（i）/b（i-1）;
a（i）=l（i）;
u（i）=b（i）-c（i-1）*a（i）;
b（i）=u（i）;
y（i）=d（i）-a（i）*d（i-1）;
d（i）=y（i）;
x（n）=d（n）/b（n）;
d（n）=x（n）;
fori=（n-1）:
-1:
1
x（i）=（d（i）-c（i）*d（i+1））/b（i）;
d（i）=x（i）;
求各段电流量程序：
8
a（i）=-2;
b=[25555555];
c=[-2-2-2-2-2-2-2];
V=220;
R=27;
d=[V/R0000000];
n=8;
I=myfun（a,b,c,d,n）
运行程序得：
I=
8.14784.07372.03651.01750.50730.25060.11940.0477
3：
（1）求矩阵A和向量b的matlab程序：
function[Ab]=myfun（n）
fori=1:
X（i）=1+0.1*i;
A（i,j）=X（i）^（j-1）;
b（i）=sum（A（i,:
））;
求n=5时A1，b1及A1的2-条件数程序运行结果如下：
n=5；
[A1,b1]=myfun（n）
A1=
1.00001.10001.21001.33101.4641
1.00001.20001.44001.72802.0736
1.00001.30001.69002.19702.8561
1.00001.40001.96002.74403.8416
1.00001.50002.25003.37505.0625
b1=
6.10517.44169.043110.945613.1875
cond2=cond（A1,2）
cond2=
5.3615e+005
求n=10时A2，b2及A2的2-条件数程序运行结果如下：
n=10;
[A2,b2]=myfun（n）
A2=
1.00001.10001.21001.33101.46411.61051.77161.94872.14362.3579
1.00001.20001.44001.72802.07362.48832.98603.58324.29985.1598
1.00001.30001.69002.19702.85613.71294.82686.27498.157310.6045
1.00001.40001.96002.74403.84165.37827.529510.541414.757920.6610
1.00001.50002.25003.37505.06257.593811.390617.085925.628938.4434
1.00001.60002.56004.09606.553610.485816.777226.843542.949768.7195
1.00001.70002.89004.91308.352114.198624.137641.033969.7576118.5879
1.00001.80003.24005.832010.497618.895734.012261.2220110.1996198.3593
1.00001.90003.61006.859013.032124.761047.045989.3872169.8356322.6877
1.00002.00004.00008.000016.000032.000064.0000128.0000256.0000512.0000
b2=
1.0e+003*
0.01590.02600.04260.06980.11330.18160.28660.44510.68011.0230
cond2=cond（A2,2）
8.6823e+011
求n=20时A3，b3及A3的2-条件数程序运行结果如下：
n=20;
[A3,b3]=myfun（n）
A3=
1.0e+009*
Columns1through10
0.00000.00000.00000.00000.00000.00000.00000.00000.00000.0000
Columns11through20
0.00000.00000.00000.00000.00000.00000.00000.00000.00000.0001
0.00000.00000.00000.00000.00000.00000.00000.00010.00010.0002
0.00000.00000.00000.00000.00000.00000.00010.00010.00030.0005
0.00000.00000.00000.00000.00000.00010.00010.00030.00060.0013
0.00000.00000.00000.00000.00010.00010.00030.00070.00150.0032
0.00000.00000.00000.00010.00010.00030.00060.00140.00320.0075
0.00000.00000.00000.00010.00020.00050.00120.00290.00700.0167
0.00000.00000.00010.00010.00040.00090.00230.00580.01460.0364
0.00000.00000.00010.00020.00060.00170.00440.01130.02950.0766
0.00000.00010.00020.00040.00110.00300.00800.02150.05810.1570
0.00000.00010.00020.00070.00180.00510.01430.04000.11190.3133
0.00000.00010.00040.00100.00300.00860.02500.07260.21050.6103
0.00010.00020.00050.00160.00480.01430.04300.12910.38741.1623
b3=
0.00000.00000.00000.00000.00000.00000.00010.00020.00040.0010
0.00250.00590.01320.02870.06060.12460.24940.48740.93161.7434
cond2=cond（A3,2）
3.2395e+022
由上述运行结果可知：
它们是病态的，而且随着n的增大，矩阵的病态变得严重。
当n=5时：
x1=A1\b1'
x1=
1.0000
1.0000
当n=10时：
x2=A2\b2'
x2=
1.0001
0.9999
当n=20时：
x3=A3\b3'
x3=
1.0e+005*
0.0203
-0.1756
0.7034
-1.7228
2.8742
-3.4342
2.9927
-1.8765
0.7820
-0.1396
-0.0720
0.0745
-0.0350
0.0108
-0.0023
0.0003
-0.0000
0.0000
由运行结果可见：
x1与精确解吻合，x2与精确解稍有差异，x3与精确解差别很大。
可见随着n的增大，矩阵病态越来越严重。
A2（2,2）=A（2,2）+1e-8;
A2（10,10）=A（10,10）+1e-8;
x=A2\b2'
x=
1.0137
0.9197
1.2089
0.6844
1.3053
0.8039
1.0837
0.9771
1.0036
0.9997
比较可见，系数矩阵出现微小变动，导致解出现较大变化。
说明n=10时，系数矩阵是病态的。
4：
（1）
A=[10787;
7565;
86109;
75910];
b=[32233331]'
;
det（A）
1
cond（A）
2.9841e+003
eig（A）
0.0102
0.8431
3.8581
30.2887
A1=[107.28.16.9;
7.085.076.025;
8.25.899.969.01;
6.985.048.979.98];
x=[1111]'
x1=A1\b
0.0077
2.3117
1.0211
1.0157
dx=x1-x
dx=
-0.9923
1.3117
0.0211
0.0157
dA=A1-A
dA=
00.20000.1000-0.1000
0.08000.07000.02000
0.2000-0.1100-0.04000.0100
-0.02000.0400-0.0300-0.0200
根据式（2-39）知：
当dA充分小，使得||A-1||*||δA||<
1时,则有：
norm（dx）/norm（x）
0.8225
（cond（A）*norm（dA）/norm（A））/（1-cond（A）*norm（dA）/norm（A））
-1.0358
norm（inv（A））*norm（dA）
28.8964
显然，上式不成立。
显然，原因是因为dA较大，使norm（inv（A））*norm（dA）=28.8964>
dA=0.5*1e-4*rand（4）;
A1=A+dA
10.00007.00008.00007.0000
7.00005.00006.00005.0000
8.00006.000010.00009.0000
7.00005.00009.000010.0000
0.9996
1.0007
0.9998
dx=x-x1
1.0e-003*
0.4360
-0.6743
0.1508
-0.0952
norm（dx）
8.2256e-004
4.1128e-004
0.0122
0.0121
由计算结果可知dA充分小，使得||A-1||*||δA||=0.0121<
1时,有：
第三章
（1）用jacobi迭代法：
编写jacobi迭代法的m文件如下：
functionx1=jacobi（A,b,n,x,e,N）
N
if（j==i）
continue;
sum=sum+A（i,j）*x（j）;
x1（i）=（b（i）-sum）/A（i,i）;
if（norm（x1-x）<
e）
break;
x=x1;
保存为jacobi.m文件。
然后在matlab命令窗口中编程计算：
A=[101234;
19-12-3;
2-173-5;
32312-1;
4-3-5-115];
b=[12-2714-1712];
x=[00000];
x1=jacobi（A,b,5,x,1e-6,15）
1.0318-2.02972.9451-1.99200.9620
即用jacobi迭代法求得解为：
[1.0318-2.02972.9451-1.99200.9620]'
（2）用Gauss-Seidel迭代法解：
编写Gauss-Seidel迭代法的m文件如下：
functionx1=gausdel（A,b,n,x,e,N）
forj=2:
sum=sum+A（1,j）*x（j）;
x1
（1）=（b
（1）-sum）/A（1,1）;
fori=2:
n-1
f=0;
g=0;
i-1
f=f+A（i,j）*x1（j）;
forj=i+1:
g=g+A（i,j）*x（j）;
x1（i）=（b（i）-f-g）/A（i,i）;
sum=sum+A（n,j）*x1（j）;
x1（n）=（b（n）-sum）/A（n,n）;
保存为gausdel.m文件。
x2=gausdel（A,b,5,x,1e-6,15）
1.0055-2.00462.9921-1.99930.9950
即用Gauss-Seidel迭代法求得解为：
[1.0055-2.00462.9921-1.99930.9950]'
（3）用共轭梯度法解：
编写共轭梯度法的m文件如下：
functionx=gonger（A,b,x0,e）
r0=b-（A*x0'
）'
d0=r0;
z0=r0*d0'
/（d0*（A*d0'
x1=x0+z0*d0;
r1=b-（A*x1'
while（norm（r1
升级会员 