m=i;
end
end
forj=1:
n
c=A(m,j);
A(m,j)=A(k,j);
A(k,j)=c;
end
forj=k:
n
sum=0;
forr=1:
(k-1)
sum=sum+A(k,r)*A(r,j);
end
u(k,j)=A(k,j)-sum;
A(k,j)=u(k,j);
end
fori=1:
n
l(i,i)=1;
end
fori=(k+1):
n
sum=0;
forr=1:
(k-1)
sum=sum+A(i,r)*u(r,k);
end
l(i,k)=(A(i,k)-sum)/u(k,k);
A(i,k)=l(i,k);
end
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
ans=
5-1010-51
-1030-3519-4
10-3546-276
-519-2717-4
1-46-41
2:
程序:
functiond=myfun(a,b,c,d,n)
fori=2:
n
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);
end
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);
end
求各段电流量程序:
fori=2:
8
a(i)=-2;
end
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:
n
X(i)=1+0.1*i;
end
fori=1:
n
forj=1:
n
A(i,j)=X(i)^(j-1);
end
end
fori=1:
n
b(i)=sum(A(i,:
));
end
求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)
cond2=
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
0.00000.00000.00000.00000.00000.00000.00000.00000.00000.0000
0.00000.00000.00000.00000.00000.00000.00000.00000.00000.0000
0.00000.00000.00000.00000.00000.00000.00000.00000.00000.0000
0.00000.00000.00000.00000.00000.00000.00000.00000.00000.0000
0.00000.00000.00000.00000.00000.00000.00000.00000.00000.0000
0.00000.00000.00000.00000.00000.00000.00000.00000.00000.0000
0.00000.00000.00000.00000.00000.00000.00000.00000.00000.0000
0.00000.00000.00000.00000.00000.00000.00000.00000.00000.0000
0.00000.00000.00000.00000.00000.00000.00000.00000.00000.0000
0.00000.00000.00000.00000.00000.00000.00000.00000.00000.0000
0.00000.00000.00000.00000.00000.00000.00000.00000.00000.0000
0.00000.00000.00000.00000.00000.00000.00000.00000.00000.0000
0.00000.00000.00000.00000.00000.00000.00000.00000.00000.0000
0.00000.00000.00000.00000.00000.00000.00000.00000.00000.0000
0.00000.00000.00000.00000.00000.00000.00000.00000.00000.0000
0.00000.00000.00000.00000.00000.00000.00000.00000.00000.0000
0.00000.00000.00000.00000.00000.00000.00000.00000.00000.0000
0.00000.00000.00000.00000.00000.00000.00000.00000.00000.0000
0.00000.00000.00000.00000.00000.00000.00000.00000.00000.0000
Columns11through20
0.00000.00000.00000.00000.00000.00000.00000.00000.00000.0000
0.00000.00000.00000.00000.00000.00000.00000.00000.00000.0000
0.00000.00000.00000.00000.00000.00000.00000.00000.00000.0000
0.00000.00000.00000.00000.00000.00000.00000.00000.00000.0000
0.00000.00000.00000.00000.00000.00000.00000.00000.00000.0000
0.00000.00000.00000.00000.00000.00000.00000.00000.00000.0000
0.00000.00000.00000.00000.00000.00000.00000.00000.00000.0000
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=
1.0e+009*
Columns1through10
0.00000.00000.00000.00000.00000.00000.00010.00020.00040.0010
Columns11through20
0.00250.00590.01320.02870.06060.12460.24940.48740.93161.7434
cond2=cond(A3,2)
cond2=
3.2395e+022
由上述运行结果可知:
它们是病态的,而且随着n的增大,矩阵的病态变得严重。
(2)
当n=5时:
x1=A1\b1'
x1=
1.0000
1.0000
1.0000
1.0000
1.0000
当n=10时:
x2=A2\b2'
x2=
1.0000
1.0000
1.0000
1.0001
0.9999
1.0000
1.0000
1.0000
1.0000
1.0000
当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
0.0000
0.0000
由运行结果可见:
x1与精确解吻合,x2与精确解稍有差异,x3与精确解差别很大。
可见随着n的增大,矩阵病态越来越严重。
(3)
当n=10时:
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)
ans=
1
cond(A)
ans=
2.9841e+003
eig(A)
ans=
0.0102
0.8431
3.8581
30.2887
(2)
A1=[107.28.16.9;7.085.076.025;8.25.899.969.01;6.985.048.979.98];
x=[1111]'
x1=A1\b
x1=
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)
ans=
0.8225
(cond(A)*norm(dA)/norm(A))/(1-cond(A)*norm(dA)/norm(A))
ans=
-1.0358
norm(inv(A))*norm(dA)
ans=
28.8964
显然,上式不成立。
显然,原因是因为dA较大,使norm(inv(A))*norm(dA)=28.8964>1
(3)
dA=0.5*1e-4*rand(4);
A1=A+dA
A1=
10.00007.00008.00007.0000
7.00005.00006.00005.0000
8.00006.000010.00009.0000
7.00005.00009.000010.0000
x1=A1\b
x1=
0.9996
1.0007
0.9998
1.0001
dx=x-x1
dx=
1.0e-003*
0.4360
-0.6743
0.1508
-0.0952
norm(dx)
ans=
8.2256e-004
根据式(2-39)知:
当dA充分小,使得||A-1||*||δA||<1时,则有:
norm(dx)/norm(x)
ans=
4.1128e-004
(cond(A)*norm(dA)/norm(A))/(1-cond(A)*norm(dA)/norm(A))
ans=
0.0122
norm(inv(A))*norm(dA)
ans=
0.0121
由计算结果可知dA充分小,使得||A-1||*||δA||=0.0121<1时,有:
第三章
1:
(1)用jacobi迭代法:
编写jacobi迭代法的m文件如下:
functionx1=jacobi(A,b,n,x,e,N)
fork=1:
N
fori=1:
n
sum=0;
forj=1:
n
if(j==i)
continue;
end
sum=sum+A(i,j)*x(j);
end
x1(i)=(b(i)-sum)/A(i,i);
end
if(norm(x1-x)break;
end
x=x1;
end
保存为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)
x1=
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)
fork=1:
N
sum=0;
forj=2:
n
sum=sum+A(1,j)*x(j);
end
x1
(1)=(b
(1)-sum)/A(1,1);
fori=2:
n-1
f=0;g=0;
forj=1:
i-1
f=f+A(i,j)*x1(j);
end
forj=i+1:
n
g=g+A(i,j)*x(j);
end
x1(i)=(b(i)-f-g)/A(i,i);
end
sum=0;
forj=1:
n-1
sum=sum+A(n,j)*x1(j);
end
x1(n)=(b(n)-sum)/A(n,n);
if(norm(x1-x)break;
end
x=x1;
end
保存为gausdel.m文件。
然后在matlab命令窗口中编程计算:
>>A=[101234;19-12-3;2-173-5;32312-1;4-3-5-115];
>>b=[12-2714-1712];
>>x=[00000];
>>x2=gausdel(A,b,5,x,1e-6,15)
x2=
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