数值分析几十个程序概要.docx
- 文档编号:17906414
- 上传时间:2023-08-04
- 格式:DOCX
- 页数:35
- 大小:144.13KB
数值分析几十个程序概要.docx
《数值分析几十个程序概要.docx》由会员分享,可在线阅读,更多相关《数值分析几十个程序概要.docx(35页珍藏版)》请在冰点文库上搜索。
数值分析几十个程序概要
第一题关于舍入误差累计的效果模拟。
x*是在(0,1)上服从均匀分布的随机数,对x*取5位有效数字得到x,将产生的k(比如取为10000)个x*(x)相加得到X*(X),研究X*-X的分布情况以及X*-X和k的关系。
将得到的结果用图形表示出来。
k=1;
kk=ones(1,1000);
xc=ones(1,1000);
whilek<=1000
r=rand(1,k);
v=vpa(r,5);
x0=sum(r);
x1=sum(v);
cha=x0-x1;
xc(k)=cha;
kk(k)=k;
k=k+1;
end;
plot(kk,xc,'rh')
可以看出随着处理的数据的增大误差也越来越大,但是还分布在横轴的两侧
第二题研究产生各种特定矩阵的方法(阶数在10-100),比如对称阵,三对角阵,正定矩阵,正交阵,对角占优矩阵,说明如何生成。
(1)生成对称矩阵
k=10;
a=rand(10,10);
fori=1:
k%先生成一个下三角矩阵;
forj=1:
k
ifj>i
a(i,j)=0;
end;
end;
end;
aa=a';%加上他本身的转置;
g=a+aa
g=
1.97600.03770.88520.91330.79620.09870.26190.33540.67970.1366
0.03770.21350.65380.49420.77910.71500.90370.89090.33420.6987
0.88520.65381.48810.50000.47990.90470.60990.61770.85940.8055
0.91330.49420.50001.77300.02870.48990.16790.97870.71270.5005
0.79620.77910.47990.02870.14290.52160.09670.81810.81750.7224
0.09870.71500.90470.48990.52161.60070.45380.43240.82530.0835
0.26190.90370.60990.16790.09670.45380.79850.52690.41680.6569
0.33540.89090.61770.97870.81810.43240.52690.74480.19810.4897
0.67970.33420.85940.71270.81750.82530.41680.19811.88550.4177
0.13660.69870.80550.50050.72240.08350.65690.48970.41771.9982
(2)
三对角矩阵
k=zeros(5,5);
a=rand(1,5);
b=rand(1,4);
c=rand(1,4);
a=diag(a);
b=diag(b);
c=diag(c);
k(1:
4,2:
5)=c;
aa=k;%将得到的结果赋值给其他变量因为此时的K已经不是以前的了必须再次初始化才能不会影响以后的操作(如下一步)
k=zeros(5,5);
k(2:
5,1:
4)=b;
bb=k;
aa+bb+a%将上述的三个矩阵复合得到三对角;
ans=
0.50050.8175000
0.07140.47110.722400
00.52160.05960.14990
000.09670.68200.6596
0000.81810.0424
也是三对角
k=6;
a=zeros(6,6);
a(1,1:
2)=rand(1,2);
a(6,5:
6)=rand(1,2);
i=2;%i不能忘了初始化
whilei>1&i<6
a(i,(i-1):
(i+1))=rand(1,3);
i=i+1;
end;
a
a=
0.13890.69630000
0.53030.86110.4849000
00.39350.67140.741300
000.52010.34770.15000
0000.58610.26210.0445
00000.09380.5254
(3)正定矩阵
a=rand(5,5);
b=a';
c=a*b
c=
0.35050.65540.61040.57830.7651
0.65541.62871.40200.91471.5272
0.61041.40201.47581.18231.7376
0.57830.91471.18231.43371.7436
0.76511.52721.73761.74362.3737
矩阵的转置和他本身的乘积是正定的
(4)正交矩阵
随机生成一个方阵,然后利用schmidt正交化对列向量进行正交化单位化得到一个正交矩阵;
a=rand(5,5);
a=vpa(a,7);
b=zeros(5,5);
b(:
1)=a(:
1);
fori=2:
5
sum=zeros(5,1);
forj=1:
(i-1)
sum=sum+(dot(a(:
i),b(:
j))/dot(b(:
j),b(:
j)))*b(:
j);
end;
b(:
i)=a(:
i)-sum;
end;%完成对矩阵的正交化schmidt;
fork=1:
5
b(:
k)=b(:
k)/sqrt(dot(b(:
k),b(:
k)));
end;%完成对矩阵的单位化;
b=vpa(b,7)
b*b'
b=
[0.4662487,-0.5697477,0.4464938,0.5031049,0.07435341]
[0.4330088,0.5708757,-0.2023893,0.4915742,-0.451661]
[0.6966242,-0.2311083,-0.434368,-0.5221337,0.002112758]
[0.2120276,0.3050953,0.7502452,-0.4773735,-0.266848]
[0.2547048,0.4505489,0.09021349,0.06878292,0.8480929]
ans=
[1.000000000000000658769429433138,0.00000000000000020501909475855146103181737845679,-0.00000000000000021440423350092585100221961008361,0.0000000000000005427542082813480164677324711907,-0.00000000000000011700162166438432633834107099474]
[0.00000000000000020501909475855146103181737845679,1.0000000000000003207231381623452,0.00000000000000058082238049994760442870183594865,-0.00000000000000073536380504684284367618389622692,0.00000000000000040560052001944880242963606302516]
[-0.00000000000000021440423350092585100221961008361,0.00000000000000058082238049994760442870183594865,1.0000000000000016447931562774053,0.00000000000000063268757033543168173701894986961,-0.00000000000000013226852520735450182808565768223]
[0.0000000000000005427542082813480164677324711907,-0.00000000000000073536380504684284367618389622692,0.00000000000000063268757033543168173701894986961,0.99999999999999936362009174092183,0.00000000000000025569754497442931226184064248466]
[-0.00000000000000011700162166438432633834107099474,0.00000000000000040560052001944880242963606302516,-0.00000000000000013226852520735450182808565768223,0.00000000000000025569754497442931226184064248466,0.99999999999999902311091880723728]
(5)对角占优矩阵
随机生成一个矩阵,将每行的元素求和然后复赋值给对角元素,得到对角占优矩阵;
a=4*rand(10,10);
fori=1:
6
a(i,i)=sum(abs(a(i,1:
6)));%对角占优矩阵式对角线元素的绝对值之和大于其他元素的绝对值之和,所以这里取了绝对值;
end
a
a=
12.38281.38000.27413.06663.63121.75010.90742.73880.39202.7565
1.38598.84160.87220.08733.89831.28341.78411.77240.98272.8715
2.23003.709912.00571.57240.47950.53621.06491.74272.46292.2361
1.19913.02431.65699.50472.07590.53831.83643.17211.21982.1334
0.63631.15302.64480.816911.76263.22381.73163.26223.06793.5029
2.66102.42473.13392.64912.548115.51601.03853.00851.06891.5724
2.73683.06410.99153.65903.81563.77700.53493.15700.15801.8323
3.16963.38462.21770.02763.78773.95341.67692.00511.18620.8330
1.39453.60790.91832.98573.86641.63962.02742.22072.22553.0291
1.00032.38280.02773.19870.26941.48481.29732.52303.87632.1868
下面也是三对角占优
k=6;
a=zeros(6,6);
a(1,1:
2)=rand(1,2);
a(6,5:
6)=rand(1,2);
i=2;%i不能忘了初始化哦亲
whilei>1&i<6
a(i,(i-1):
(i+1))=rand(1,3);
i=i+1;
end;
fori=1:
6
a(i,i)=sum(a(i,:
));
end;%在三对角矩阵的基础上将每一行的数相加然后复制给对角元素则得到的矩阵式对角占优的;
a
a=
1.02920.86200000
0.88431.62710.1548000
00.19991.35550.748700
000.82561.93410.31850
0000.53410.73570.1117
00000.98991.5043
第三题比较高斯消去法和带有列主元的高斯消去法的效果。
要求构造一些方程组,方程组的精确解已知,用上面两种方法分别求解,说明选主元的效果。
高斯消去法
A=[1231;2756;149-3];%增广矩阵;
x=zeros(1,3);
[m,n]=size(A);
fori=1:
(m-1)
forj=(i+1):
m
A(j,:
)=A(j,:
)-A(i,:
)*A(j,i)/A(i,i);%消去;
end
end
A;
x(m)=A(m,n)/A(m,m);%回代未知数的最后一个分量单独计算出其他的如下;
fori=(m-1):
-1:
1
x(i)=(A(i,n)-A(i,i+1:
m)*x(i+1:
m)')/A(i,i);
end
x
x=
21-1
高斯列主元法
A=[1231;2756;149-3];%增广矩阵;
[m,n]=size(A);
fori=1:
(m-1)
temp=max(abs(A(i:
m,i)));%当前处理的矩阵的第一列的绝对值最大的元素;
[a,b]=find(abs(A(i:
m,i))==temp);%找到最大元素所在的位置;
tempo=A(a
(1)+i-1,:
);
A(a
(1)+i-1,:
)=A(i,:
);
A(i,:
)=tempo;%交换两行;
forj=(i+1):
m
A(j,:
)=A(j,:
)-A(i,:
)*A(j,i)/A(i,i);%消去;
end
A;
end
x(m)=A(m,n)/A(m,m);%回代;
fori=(m-1):
-1:
1
x(i)=(A(i,n)-A(i,i+1:
m)*x(i+1:
m)')/A(i,i);
end
x
A=
2756
0-3/21/2-2
01/213/2-6
A=
2756
0-3/21/2-2
0020/3-20/3
x=
21-1
下面构造一些已知解的矩阵,然后分别用以上两种方法进行求解,比较它们运行的结果:
高斯消去
A=rand(500,500);
yizhij=rand(500,1);
b=A*yizhij;
A=[A,b];%增广矩阵;
x=zeros(1,500);
[m,n]=size(A);
fori=1:
(m-1)
forj=(i+1):
m
A(j,:
)=A(j,:
)-A(i,:
)*A(j,i)/A(i,i);%消去;
end
end
A;
x(m)=A(m,n)/A(m,m);%回代未知数的最后一个分量单独计算出其他的如下;
fori=(m-1):
-1:
1
x(i)=(A(i,n)-A(i,i+1:
m)*x(i+1:
m)')/A(i,i);
end
dot(x'-yizhij,x'-yizhij)
ans=
4.8354e-020
高斯列主元法
A=rand(500,500);
yizhij=rand(500,1);
b=A*yizhij;
A=[A,b];%增广矩阵;
[m,n]=size(A);
fori=1:
(m-1)
temp=max(abs(A(i:
m,i)));%当前处理的矩阵的第一列的绝对值最大的元素;
[a,b]=find(abs(A(i:
m,i))==temp);%找到最大元素所在的位置;
tempo=A(a
(1)+i-1,:
);
A(a
(1)+i-1,:
)=A(i,:
);
A(i,:
)=tempo;%交换两行;
forj=(i+1):
m
A(j,:
)=A(j,:
)-A(i,:
)*A(j,i)/A(i,i);%消去;
end
A;
end
x(m)=A(m,n)/A(m,m);%回代;
fori=(m-1):
-1:
1
x(i)=(A(i,n)-A(i,i+1:
m)*x(i+1:
m)')/A(i,i);
end
dot(x'-yizhij,x'-yizhij)
ans=
6.1236e-023
可以看出运用高斯列主元法求解比高斯消去法更加精确。
第四题采用紧凑格式实现矩阵的LU分解,并使用这个方法求解线性方程组。
A=[2426;49615;26918;6151840];
N=size(A);
n=N
(1);
L=zeros(4,4);
U=eye(4,4);%U的对角元素为1;
L(1:
4,1)=A(1:
4,1);%L的第一列;
U(1,1:
4)=A(1,1:
4)/L(1,1);%U的第一行;
fork=2:
4
fori=k:
4
L(i,k)=A(i,k)-L(i,1:
(k-1))*U(1:
(k-1),k);
%L的第k列;
end
forj=(k+1):
4
U(k,j)=(A(k,j)-L(k,1:
(k-1))*U(1:
(k-1),j))/L(k,k);
%U的第k行;
end
end
L
U
b=[9232247]';
y=inv(L)*b;
x=inv(U)*y%求解方程;
L=
2000
4100
2230
6361
U=
1213
0123
0012
0001
x=
1/2
2
3
-1
第五题设计对称正定线性方程组,使用平方根法计算方程组的解。
A=[2426;49615;26918;6151840];
b=[9232247]';
n=length(b);%方程个数n
L=zeros(n,n);
L(1,1)=sqrt(A(1,1));
L(2:
n,1)=A(2:
n,1)/L(1,1);
forj=2:
n-1
L(j,j)=sqrt(A(j,j)-sum(L(j,1:
j-1).^2));
fori=j+1:
n
L(i,j)=(A(i,j)-sum(L(i,1:
j-1).*L(j,1:
j-1)))/L(j,j);
end
end
L(n,n)=sqrt(A(n,n)-sum((L(n,1:
n-1)).^2));
L
y=inv(L)*b;
x=inv(L')*y%求解方程;
L=
1393/985000
3363/1189100
1393/98521351/7800
4756/112131351/3901
x=
1/2
2
3
-1
第六题追赶法求解三对角线性方程组三对角线性方程组的追赶法编程,使得算法计算次数达到最小。
a=[2.0-1000;
-12-100;
0-12-10;
00-12-1;
000-12];
b=[10000];
y=zeros(1,5);
x=zeros(1,5);
L=zeros(5,5);U=zeros(5,5);
U(1,2)=a(1,2)/a(1,1);
fori=2:
4%利用书上给出的步骤编写;
U(i,(i+1))=a(i,(i+1))/(a(i,i)-a(i,(i-1))*U((i-1),i));
end;
y
(1)=b
(1)/a(1,1);
fori=2:
5
y(i)=(b(i)-a(i,(i-1))*y(i-1))/(a(i,i)-a(i,(i-1))*U((i-1),i));
end;
x(5)=y(5);
i=4;
whilei>=1
x(i)=y(i)-U(i,(i+1))*x(i+1);
i=i-1;
end;
y
x=vpa(x,4)
y=
1/21/31/41/51/6
x=
[0.8333,0.6667,0.5,0.3333,0.1667]
第七题调用Matlab的cond()函数分析随机生成的非奇异阵的的阶数和条件数之间的关系
i=11;
kk=zeros(1,90);
con=zeros(1,90);
whilei<=100
kk(i-10)=i;
k=1;
chengji=1;
whilek<=50
a=rand(i,i);%生成i阶的方阵,对么一个方阵作50次的求解条件数的运算,最后求取平均值;
c=cond(a);
chengji=chengji*c;%把每次一运行出来的结果乘在一起,最后取几何平均数;
k=k+1;
end;
c=chengji.^(1/50);
con(i-10)=c;
i=i+1;
end;
plot(kk,con,'rh')
holdon;
p=polyfit(kk,con,2);
con=polyval(p,kk);
plot(kk,con)
可以看出随着非奇异矩阵阶数的不断增大条件数呈现多项式函数增大
第八题编写线性方程组的Jacobi方法和高斯赛德尔迭代法的程序。
使用图形反映谱半径和残差之间的关系。
Jacobi方法
高斯赛德尔迭代法
A=[8-32;411-1;6312];
b=[203336]';
e=0.000001;
n=max(size(A));
fori=1:
n
ifA(i,i)==0
'对角元素为0不能求解'
return;
end;
end;
x=ones(n,1);
k=0;%迭代次数的初始化;
mk=50;%迭代次数的最大值;
r=1;%前后项之差的无穷范数;
whilek<=mk&r>e
x0=x;
fori=1:
n
sum=0;
forj=1:
i-1
sum=sum+A(i,j)*x0(j);
end;
forj=i+1:
n
sum=sum+A(i,j)*x0(j);
end;
x(i)=(b(i)-sum)/A(i,i);
end;%以上是迭代公式的等价翻译;
r=norm(x-x0,inf);
k=k+1;
end;
ifk>mk
'迭代失败'
else
'迭代成功'
end;
x
k
ans=
迭代成功
x=
3.0000
2.0000
1.0000
k=
15
A=[8-32;411-1;6312];
b=[203336]'
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 数值 分析 几十 程序 概要
![提示](https://static.bingdoc.com/images/bang_tan.gif)