《数值分析》实验报告书要点.docx
《《数值分析》实验报告书要点.docx》由会员分享，可在线阅读，更多相关《《数值分析》实验报告书要点.docx（18页珍藏版）》请在冰点文库上搜索。
《数值分析》实验报告书要点
《数值分析》实验报告
实验一、误差分析
误差问题是数值分析的基础，又是数值分析中一个困难的课题。
在实际计算中，如果选用了不同的算法，由于舍入误差的影响，将会得到截然不同的结果。
因此，选取算法时注重分析舍入误差的影响，在实际计算中是十分重要的。
同时，由于在数值求解过程中用有限的过程代替无限的过程会产生截断误差，因此算法的好坏会影响到数值结果的精度。
一、实验目的
1、通过上机编程，复习巩固以前所学程序设计语言及上机操作指令；
2、通过上机计算，了解误差、绝对误差、误差界、相对误差界的有关概念；
3、通过上机计算，了解舍入误差所引起的数值不稳定性。
二、实验任务
对
，计算定积分
.
算法1：
利用递推公式
取
.
算法2：
利用递推公式
.
注意到
取
.
思考：
从计算结果看，哪个算法是不稳定的，哪个算法是稳定的。
算法1：
t=log（6.0）-log（5.0）;
n=0;
y=zeros（1,21）;
y
（1）=t;
fork=2:
21
y（k）=1/k-5*y（k-1）;
n=n+1;
end
y（1:
6）
y（7:
21）
运行结果：
ans=0.1823-0.41162.3914-11.706958.7343-293.5049
算法2：
y=zeros（21,1）;
n=1;
y1=（1/105+1/126）/20;
fork=21:
-1:
2
y（k-1）=1/（5*k）-y（k）/5;
n=n+1;
end
运行结果：
y=0.08840.05800.04310.03430.02850.02430.02120.01880.01690.01540.01410.01300.01200.01120.01050.00990.00930.00890.00810.00950
由数据对比可知，算法2较为稳定。
实验二、插值法
插值法是函数逼近的一种重要方法，它是数值积分、微分方程数值解等数值计算的基础与工具，其中多项式插值是最常用和最基本的方法。
拉格朗日插值多项式的优点是表达式简单明确，形式对称，便于记忆，它的缺点是如果想要增加插值节点，公式必须整个改变，这就增加了计算工作量。
而牛顿插值多项式对此做了改进，当增加一个节点时只需在原牛顿插值多项式基础上增加一项，此时原有的项无需改变，从而达到节省计算次数、节约存储单元、应用较少节点达到应有精度的目的。
一、实验目的
1、理解插值的基本概念，掌握各种插值方法，包括拉格朗日插值和牛顿插值等，注意其不同特点；
2、通过实验进一步理解并掌握各种插值的基本算法。
二、实验任务
1、已知函数表
0.561600.562800.564010.56521
0.827410.826590.825770.82495
用二次拉格朗日插值多项式求
时的函数近似值。
2、已知函数表
0.40.550.650.80.9
0.410750.578150.696750.888111.02652
用牛顿插值多项式求
和
。
1.
function[y,R]=lagranzi（X,Y,x,M）
x=0.5635;
M=2;
X=[0.56160,0.56280,0.56401,0.56521];
Y=[0.82741,0.82659,0.82577,0.82495];
n=length（X）;
m=length（x）;
fori=1:
m
z=x（i）;s=0.0;
fork=1:
n
p=1.0;
q1=1.0;
c1=1.0;
forj=1:
n
ifj~=k
p=p*（z-X（j））/（X（k）-X（j））;
end
q1=abs（q1*（z-X（j）））;
c1=c1*j;
end
s=p*Y（k）+s;
end
y（i）=s;
end
R=M.*q1./c1;
运行结果：
ans=0.8261
2.
N3（0.596）
function[y,R]=newcz（X,Y,x,M）
x=0.596;
M=3;
X=[0.4,0.65,0.9];
Y=[0.41075,0.69675,1.02652];
n=length（X）;m=length（x）;
fort=1:
m
z=x（t）;A=zeros（n,n）;A（:
1）=Y';
s=0.0;p=1.0;q1=1.0;c1=1.0;
forj=2:
n
fori=j:
n
A（i,j）=（A（i,j-1）-A（i-1,j-1））/（X（i）-X（i-j+1））;
end
q1=abs（q1*（z-X（j-1）））;c1=c1*j;
end
C=A（n,n）;q1=abs（q1*（z-X（n）））;
fork=（n-1）:
-1:
1
C=conv（C,poly（X（k）））;d=length（C）;C（d）=C（d）+A（k,k）;
end
y（k）=polyval（C,z）;
end
R=M*q1/c1;
运行结果：
ans=0.6313
N4（0.895）
function[y,R]=newcz（X,Y,x,M）
x=0.895;
M=4;
X=[0.4,0.55,0.65,0.8,0.9];
Y=[0.41075,0.57815,0.69675,0.88811,1.02652];
n=length（X）;m=length（x）;
fort=1:
m
z=x（t）;A=zeros（n,n）;A（:
1）=Y';
s=0.0;p=1.0;q1=1.0;c1=1.0;
forj=2:
n
fori=j:
n
A（i,j）=（A（i,j-1）-A（i-1,j-1））/（X（i）-X（i-j+1））;
end
q1=abs（q1*（z-X（j-1）））;c1=c1*j;
end
C=A（n,n）;q1=abs（q1*（z-X（n）））;
fork=（n-1）:
-1:
1
C=conv（C,poly（X（k）））;d=length（C）;C（d）=C（d）+A（k,k）;
end
y（k）=polyval（C,z）;
end
R=M*q1/c1;
运行结果：
ans=1.0194
实验三、解线性方程组的直接法
解线性方程组的直接法是指经过有限步运算后能求得方程组精确解的方法。
但由于实际计算中舍入误差是客观存在的，因而使用这类方法也只能得到近似解。
目前较实用的直接法是古老的高斯消去法的变形，即主元素消去法及矩阵的三角分解法。
引进选主元的技巧是为了控制计算过程中舍入误差的增长，减少舍入误差的影响。
一般说来，列主元消去法及列主元三角分解法是数值稳定的算法，它具有精确度较高、计算量不大和算法组织容易等优点，是目前计算机上解中、小型稠密矩阵方程组可靠而有效的常用方法。
一、实验目的
1、了解求线性方程组的直接法的有关理论和方法；
2、会编制列主元消去法、LU分解法的程序；
3、通过实际计算，进一步了解各种方法的优缺点，选择合适的数值方法。
二、实验任务
1、用列主元高斯消去法求解方程组
.
2、用矩阵直接三角分解法求解方程组
，其中
，
.
1.
主程序：
function[RA,RB,n,X]=liezhu（A,b）
B=[Ab];n=length（b）;RA=rank（A）;（'求矩阵的秩'）
RB=rank（B）;zhica=RB-RA;
ifzhica>0,
disp（'请注意：
因为RA~=RB，所以此方程组无解.'）
return
end
ifRA==RB
ifRA==n
disp（'请注意：
因为RA=RB=n，所以此方程组有唯一解.'）
X=zeros（n,1）;C=zeros（1,n+1）;
forp=1:
n-1
[Y,j]=max（abs（B（p:
n,p）））;C=B（p,:
）;
B（p,:
）=B（j+p-1,:
）;B（j+p-1,:
）=C;
fork=p+1:
n
m=B（k,p）/B（p,p）;
B（k,p:
n+1）=B（k,p:
n+1）-m*B（p,p:
n+1）;
end
end
b=B（1:
n,n+1）;A=B（1:
n,1:
n）;X（n）=b（n）/A（n,n）;
forq=n-1:
-1:
1X（q）=（b（q）-sum（A（q,q+1:
n）*X（q+1:
n）））/A（q,q）;
end
else
disp（'请注意：
因为RA=RBend
end
计算程序：
A=[0.1012.3043.555;-1.3473.7124.623;-2.8351.0725.643];
b=[1.183;2.137;3.035];[RA,RB,n,X]=liezhu（A,b）
运行结果：
ans=求矩阵的秩
请注意：
因为RA=RB=n，所以此方程组有唯一解.
RA=3
RB=3
n=3
X=-0.3982
0.0138
0.3351
2.
程序：
functionX=LUjfcz（A,b）
[n,n]=size（A）;
X=zeros（n,1）;
Y=zeros（n,1）;
C=zeros（1,n）;
r=1:
n;
forp=1:
n-1
[max1,j]=max（abs（A（p:
n,p）））;
C=A（p,:
）;
A（p,:
）=A（j+p-1,:
）;
A（j+p-1,:
）=C;
g=r（p）;
r（p）=r（j+p-1）;
r（j+p-1）=g;
ifA（p,p）==0
disp（'A是奇异阵,方程组无唯一解'）;
break;
end
fork=p+1:
n
H=A（k,p）/A（p,p）;
A（k,p）=H;
A（k,p+1:
n）=A（k,p+1:
n）-H*A（p,p+1:
n）;
end
end
Y
（1）=b（r
（1））;
fork=2:
n
Y（k）=b（r（k））-A（k,1:
k-1）*Y（1:
k-1）;
end
X（n）=Y（n）/A（n,n）;
fori=n-1:
-1:
1
X（i）=（Y（i）-A（i,i+1:
n）*X（i+1:
n））/A（i,i）;
end
End
计算程序：
A=[1,2,-12,8;5,4,7,-2;-3,7,9,5;6,-12,-8,3];b=[27;4;11;49];X=LUjfcz（A,b）
运行结果：
X=3.0000
-2.0000
1.0000
5.0000
实验四、解线性方程组的迭代法
解线性方程组的迭代法是用某种极限过程去逐步逼近线性方程组精确解的方法，即是从一个初始向量
出发，按照一定的迭代格式产生一个向量序列
，使其收敛到方程组
的解。
迭代法的优点是所需计算机存储单元少，程序设计简单，原始系数矩阵在计算过程中始终不变等。
但迭代法存在收敛性及收敛速度问题。
迭代法是解大型稀疏矩阵方程组的重要方法。
一、实验目的
1、熟悉迭代法的有关理论和方法；
2、会编制雅可比迭代法、高斯-塞德尔迭代法的程序；
3、注意所用方法的收敛性及其收敛速度问题。
二、实验任务
1、用雅可比迭代法解方程组
.
注意：
若用高斯-塞德尔迭代法则发散。
2、用高斯-塞德尔迭代法解方程组
.
注意：
若用雅可比迭代法则发散。
1.
主程序：
functionX=jacdd（A,b,X0,P,wucha,max1）
[nm]=size（A）;
fork=1:
max1
k
forj=1:
m
X（j）=（b（j）-A（j,[1:
j-1,j+1:
m]）*X0（[1:
j-1,j+1:
m]））/A（j,j）;
end
X
djwcX=norm（X'-X0,P）;xdwcX=djwcX/（norm（X',P）+eps）;X0=X';
if（djwcXdisp（'请注意：
雅可比迭代收敛,此方程组的精确解jX和近似解X如下：
'）
return
end
end
if（djwcX>wucha）&（xdwcX>wucha）
disp（'请注意：
雅可比迭代次数已经超过最大迭代次数max1'）
End
计算程序：
A=[12-2;111;221];b=[7;2;5];X0=[000]';X=jacdd（A,b,X0,inf,0.01,100）
运行结果：
k=1
X=725
k=2
X=13-10-13
k=3
X=12-1
k=4
X=12-1
2.
主程序：
functionX=gsdddy（A,b,X0,P,wucha,max1）
D=diag（diag（A））;U=-triu（A,1）;L=-tril（A,-1）;dD=det（D）;
ifdD==0
disp（'请注意：
因为对角矩阵D奇异，所以此方程组无解.'）
else
disp（'请注意：
因为对角矩阵D非奇异，所以此方程组有解.'）
iD=inv（D-L）;B2=iD*U;f2=iD*b;jX=A\b;X=X0;
[nm]=size（A）;
fork=1:
max1
X1=B2*X+f2;djwcX=norm（X1-X,P）;
xdwcX=djwcX/（norm（X,P）+eps）;
if（djwcXreturn
else
k,X1',k=k+1;X=X1;
end
end
if（djwcXdisp（'请注意：
高斯-塞德尔迭代收敛,此A的分解矩阵D,U,L和方程组的精确解jX和近似解X如下：
'）
else
disp（'请注意：
高斯-塞德尔迭代的结果没有达到给定的精度，并且迭代次数已经超过最大迭代次数max1,方程组的精确解jX和迭代向量X如下：
'）
X=X';jX=jX'
end
end
X=X';D,U,L,jX=jX'
计算程序：
A=[10.90.9;0.910.9;0.90.91];b=[1.9;2.0;1.7];
X0=[000]';X=gsdddy（A,b,X0,inf,0.001,100）
运行结果：
k=1
ans=1.90000.2900-0.2710
k=2
ans=1.88290.5493-0.4890
k=3
ans=1.84570.7789-0.6622
k=4
ans=1.79490.9805-0.7979
k=5
ans=1.73561.1560-0.9025
k=6
ans=1.67181.3076-0.9815
k=7
ans=1.60651.4375-1.0396
k=8
ans=1.54191.5479-1.0808
k=9
ans=1.47961.6411-1.1086
k=10
ans=1.42081.7191-1.1259
k=11
ans=1.36611.7838-1.1349
k=12
ans=1.31601.8370-1.1377
k=13
ans=1.27061.8804-1.1359
k=14
ans=1.23001.9153-1.1308
k=15
ans=1.19391.9432-1.1234
k=16
ans=1.16221.9651-1.1145
k=17
ans=1.13451.9820-1.1049
k=18
ans=1.11061.9949-1.0949
k=19
ans=1.09002.0044-1.0850
k=20
ans=1.07252.0112-1.0754
k=21
ans=1.05772.0159-1.0662
k=22
ans=1.04532.0188-1.0577
k=23
ans=1.03502.0204-1.0499
k=24
ans=1.02652.0210-1.0428
k=25
ans=1.01962.0209-1.0364
k=26
ans=1.01402.0202-1.0308
k=27
ans=1.00952.0191-1.0258
k=28
ans=1.00602.0178-1.0214
k=29
ans=1.00322.0164-1.0176
k=30
ans=1.00122.0148-1.0144
k=31
ans=0.99962.0133-1.0116
k=32
ans=0.99852.0118-1.0093
X=0.9985
2.0118
-1.0093
升级会员 