数值代数上机实验报告Word格式文档下载.docx
《数值代数上机实验报告Word格式文档下载.docx》由会员分享，可在线阅读，更多相关《数值代数上机实验报告Word格式文档下载.docx（13页珍藏版）》请在冰豆网上搜索。
n-1
b（j）=b（j）/A（j,j）;
b（j+1:
n）=b（j+1:
n）-b（j）*A（j+1:
n,j）;
end
b（n）=b（n）/A（n,n）;
%前代法
A=A'
;
forj=n:
-1:
2
b（1:
j-1）=b（1:
j-1）-b（j）*A（1:
j-1,j）;
b
（1）=b
（1）/A（1,1）;
%回代法
disp（'
平方根法的解即为b'
改进平方根法函数程序如下：
functionb=gaijinpinfanggenfa（A,b）
v=zeros（n,1）;
fori=1:
j-1
v（i）=A（j,i）*A（i,i）;
A（j,j）=A（j,j）-A（j,1:
j-1）*v（1:
j-1）;
A（j+1:
n,j）=（A（j+1:
n,j）-A（j+1:
n,1:
j-1））/A（j,j）;
end%LDL'
分解
B=diag（A）;
D=zeros（n）;
D（i,i）=B（i）;
A（i,i）=1;
End
A=tril（A）;
%得到L和D
A=D*（A'
改进平方根法解得的解即为b'
调用函数解题：
clear;
clc;
n=input（'
请输入矩阵维数：
'
b=zeros（n,1）;
A=zeros（n）;
fori=1:
forj=1:
A（i,j）=1/（i+j-1）;
b（i）=b（i）+1/（i+j-1）;
end%生成hilbert矩阵
[x,b]=pingfanggenfa（A,b）
b=gaijinpinfanggenfa（A,b）
运行结果：
40
Warning:
Matrixisclosetosingularorbadlyscaled.
Resultsmaybeinaccurate.RCOND=6.570692e-020.
>
Inpingfanggenfaat4
Inqiujieat10
Matlab自带解即为x
平方根法的解即为b
x=
1.6035
8.9685
0.8562
1.0195
0.9375
-50.2500
-3.0000
-16.0000
24.0000
-49.5000
-30.0000
39.0000
22.0000
-64.0000
-12.0000
2.0000
10.2500
-10.5000
-1.0000
-10.8750
83.0000
46.0000
-98.0000
12.0000
-69.0000
68.0000
21.0000
17.0000
-50.7188
-8.7500
-8.0000
112.0000
6.0000
-68.7500
44.0000
-28.0000
8.0000
-44.0000
b=
1.0e+007*
0.0000
-0.0000
0.0001
-0.0004
-0.0014
0.0424
-0.2980
1.1419
-2.7335
4.2539
-4.3018
2.7733
-1.1989
0.5406
-0.3688
0.3285
-0.4438
0.4621
-0.2513
0.0565
-0.0051
0.0071
-0.0027
-0.0031
0.0036
-0.0019
0.0009
0.0002
-0.0002
-0.0006
0.0004
改进平方根法解得的解即为b
1.0e+024*
-0.0012
0.0139
-0.0954
0.4208
-1.2101
2.0624
-1.0394
-3.3343
6.2567
-0.2463
-7.4594
2.8030
3.6990
0.7277
-1.7484
-0.4854
-3.6010
0.2532
5.1862
-2.1299
1.4410
0.8738
-4.5654
1.0422
4.0920
-2.7764
-2.2148
-0.8953
0.3665
4.8967
1.0416
0.1281
-4.3387
-1.1902
-2.8334
8.4610
-3.6008
实验二、利用QR分解解线性方程组：
利用QR分解解线性方程组Ax=b，其中
A=[16484;
41084;
881210;
441012];
b=[32263830];
求解程序如下：
定义house函数：
function[v,B]=house（x）
n=length（x）;
y=norm（x,inf）;
x=x/y;
Q=x（2:
n）'
*x（2:
n）;
v
（1）=1;
v（2:
n）=x（2:
ifn==1
Q=0;
B=0;
else
a=sqrt（x
（1）^2+Q）;
ifx
（1）<
=0
v
（1）=x
（1）-a;
else
v
（1）=-Q/（x
（1）+a）;
B=2*v
（1）^2/（Q+v
（1）^2）;
v=v/v
（1）;
进行QR分解：
b=b'
x=size（A）;
m=x
（1）;
n=x
（2）;
d=zeros（n,1）;
[v,B]=house（A（j:
m,j））;
m,j:
n）=（eye（m-j+1）-B*（v'
）*v）*A（j:
d（j）=B;
ifj<
m
m,j）=v（2:
m-j+1）;
end%QR分解
R=triu（A）;
%得到R
D=A;
I=eye（m,n）;
Q=I;
D（i,i）=1;
H=tril（D）;
M=H'
N=I-d（i）*H（1:
m,i）*M（i,1:
m）;
Q=Q*N;
end%得到Q
b=（Q'
）*b;
%Q是正交阵
b（j）=b（j）/R（j,j）;
j-1）-b（j）*R（1:
b
（1）=b
（1）/R（1,1）;
%回带法
运行结果如下：
R=
18.76179.807215.776911.0864
09.99099.33587.5341
005.99459.8013
000-0.5126
Q=
0.8528-0.4368-0.2297-0.1709
0.21320.7916-0.4594-0.3417
0.42640.38220.28440.7689
0.21320.19110.8095-0.5126
b=
1.00000000000000
1.00000000000001
0.999999999999988
实验三、Newton下山法解非线性方程组：
3x-cos（yz）-
=0，
-81
+sinz+1.06=0，
exp（-xy）+20z+
=0；
要求满足数值解
满足
或
定义所求方程组的函数：
Newtonfun.m
functionF=Newtonfun（X）
F（1,1）=3*X
（1）-cos（X
（2）*X（3））-1/2;
F（2,1）=X
（1）^2-81*（X
（2）+0.1）^2+sin（X（3））+1.06;
F（3,1）=exp（-X
（1）*X
（2））+20*X（3）+（10*pi-3）/3;
End
向量求导：
Xiangliangqiudao.m
functionJ=xiangliangqiudao（）
symsxyz
X=[x,y,z];
F=[3*X
（1）-cos（X
（2）*X（3））-1/2;
X
（1）^2-81*（X
（2）+0.1）^2+sin（X（3））+1.06;
exp（-X
（1）*X
（2））+20*X（3）+（10*pi-3）/3];
J=jacobian（F,[xyz]）;
代值函数：
Jacobi.m
functionF=Jacobi（x）
F=[3,x（3）*sin（x
（2）*x（3））,x
（2）*sin（x
（2）*x（3））;
2*x
（1）,-162*x
（2）-81/5,cos（x（3））;
-x
（2）/exp（x
（1）*x
（2））,-x
（1）/exp（x
（1）*x
（2））,20];
方程组求解：
formatlong;
%数据表示为双精度型
X1=[0,0,0]'
eps=10^（-8）;
k=1;
i=1;
X2=X1-Jacobi（X1）^（-1）*Newtonfun（X1）;
while（norm（subs（X2-X1,pi,3.1415926）,2）>
=eps）&
&
（norm（Newtonfun（X1）,2）>
=eps）
ifnorm（Newtonfun（X2）,2）<
norm（Newtonfun（X1）,2）%判断先后两次迭代的大小
X1=X2;
B=inv（Jacobi（X2））;
C=Newtonfun（X2）;
X2=X2-B*C;
i=i+1;
else
v=1/（2^k）;
%引入下山因子
X2=X2-v*B*C;
k=k+1;
j=i+k-1%迭代次数
X=X2%输出结果
j=
5
X=
0.500000000000000
-0.000000000000000
-0.523598775598299
升级会员 