数值分析实验报告.docx
《数值分析实验报告.docx》由会员分享，可在线阅读，更多相关《数值分析实验报告.docx（33页珍藏版）》请在冰豆网上搜索。
数值分析实验报告
实验一误差分析
实验1.1（病态问题）
实验目的：
算法有“优”与“劣”之分，问题也有“好”与“坏”之别。
对数值方法的研究而言，所谓坏问题就是问题本身对扰动敏感者，反之属于好问题。
通过本实验可获得一个初步体会。
数值分析的大部分研究课题中，如线性代数方程组、矩阵特征值问题、非线性方程与方程组等都存在病态的问题。
病态问题要通过研究和构造特殊的算法来解决，当然一般要付出一些代价（如耗用更多的机器时间、占用更多的存储空间等）。
问题提出：
考虑一个高次的代数多项式
显然该多项式的全部根为1,2,…,20共计20个，且每个根都是单重的。
现考虑该多项式的一个扰动
其中
是一个非常小的数。
这相当于是对（1.1）中
的系数作一个小的扰动。
我们希望比较（1.1）和（1.2）根的差别，从而分析方程（1.1）的解对扰动的敏感性。
实验容：
为了实现方便，我们先介绍两个Matlab函数：
“roots”和“poly”。
其中若变量a存储n+1维的向量，则该函数的输出u为一个n维的向量。
设a的元素依次为
，则输出u的各分量是多项式方程
的全部根；而函数
的输出b是一个n+1维变量，它是以n维变量v的各分量为根的多项式的系数。
可见“roots”和“poly”是两个互逆的运算函数。
上述简单的Matlab程序便得到（1.2）的全部根，程序中的“ess”即是（1.2）中的
。
实验要求：
（1）选择充分小的ess，反复进行上述实验，记录结果的变化并分析它们。
如果扰动项的系数
很小，我们自然感觉（1.1）和（1.2）的解应当相差很小。
计算中你有什么出乎意料的发现？
表明有些解关于如此的扰动敏感性如何？
（2）将方程（1.2）中的扰动项改成
或其它形式，实验中又有怎样的现象出现？
（3）（选作部分）请从理论上分析产生这一问题的根源。
注意我们可以将方程（1.2）写成展开的形式，
同时将方程的解x看成是系数
的函数，考察方程的某个解关于
的扰动是否敏感，与研究它关于
的导数的大小有何关系？
为什么？
你发现了什么现象，哪些根关于
的变化更敏感？
思考题一：
（上述实验的改进）
在上述实验中我们会发现用roots函数求解多项式方程的精度不高，为此你可以考虑用符号函数solve来提高解的精确度，这需要用到将多项式转换为符号多项式的函数poly2sym,函数的具体使用方法可参考Matlab的帮助。
实验过程：
程序：
a=poly（1:
20）;
rr=roots（a）;
forn=2:
21
n
form=1:
9
ess=10^（-6-m）;
ve=zeros（1,21）;
ve（n）=ess;
r=roots（a+ve）;
-6-m
s=max（abs（r-rr））
end
end
利用符号函数：
（思考题一）
a=poly（1:
20）;
y=poly2sym（a）;
rr=solve（y）
forn=2:
21
n
form=1:
8
ess=10^（-6-m）;
ve=zeros（1,21）;
ve（n）=ess;
a=poly（1:
20）+ve;
y=poly2sym（a）;
r=solve（y）;
-6-m
s=max（abs（r-rr））
end
end
数值实验结果与分析：
formatlong
-6-mn
-7
-8
-9
-10
2
2.331
1.158
1.748
0.047
3
1.424
0.964
0.741
0.411
4
0.536
0.061
0.834
0
5
0.871
0.844
0
0
6
0
0
0
0
7
0
0
0
0
8
0
0
0
0
9
0
0
0
0
10
0
0
0
0
11
0
0
0
0
12
0
0
0
0
13
0
0
0
0
14
0
0
0
0
15
0
0
0
0
16
0
0
0
0
17
0
0
0
0
18
0
0
0
0
19
0
0
0
0
20
0
0
0
0
21
0
0
0
0
-6-mn
-11
-12
-13
-14
2
0.380
0.280
0.598
0
3
0.546
0
0
0
4
0
0
0
0
5
0
0
0
0
6
0
0
0
0
7
0
0
0
0
8
0
0
0
0
9
0
0
0
0
10
0
0
0
0
11
0
0
0
0
12
0
0
0
0
13
0
0
0
0
14
0
0
0
0
15
0
0
0
0
16
0
0
0
0
17
0
0
0
0
18
0
0
0
0
19
0
0
0
0
20
0
0
0
0
21
0
0
0
0
讨论：
利用这种方法进行这类实验，可以很精确的扰动敏感性的一般规律。
即当对扰动项的系数越来越小时，对其多项式扰动的结果也就越来越小，即扰动敏感性与扰动项的系数成正比，扰动项的系数越大，对其根的扰动敏感性就越明显，当扰动的系数一定时，扰动敏感性与扰动的项的幂数成正比，扰动的项的幂数越高，对其根的扰动敏感性就越明显。
实验总结：
利用MATLAB来进行病态问题的实验，虽然其得出的结果是有误差的，但是可以很容易的得出对一个多次的代数多项式的其中某一项进行很小的扰动，对其多项式的根会有一定的扰动的，所以对于这类病态问题可以借助于MATLAB来进行问题的分析。
学号：
06450210
：
万轩
实验二插值法
实验2.1（多项式插值的振荡现象）
问题提出：
考虑一个固定的区间上用插值逼近一个函数。
显然拉格朗日插值中使用的节点越多，插值多项式的次数就越高。
我们自然关心插值多项式的次数增加时，L（x）是否也更加靠近被逼近的函数。
龙格给出了一个极著名例子。
设区间[-1，1]上函数
f（x）=1／（1+25x^2）
实验容：
考虑区间[-1，1]的一个等距划分，分点为：
x（i）=-1+2i/n,i=0,1,2…,n
泽拉格朗日插值多项式为：
L（x）=∑l（i）（x）/（1+25x（j）^2）i=0,1,…n
其中l（i）（x）,i=0,1,…n,n是n次拉格朗日插值基函数。
实验要求：
⑴选择不断增大的分点数目n=2,3…，画出f（x）与插值多项式函数L（x）在[-1，1]上的图象，比较分析实验结果。
（2）选择其它的函数，例如定义在区间[-5，5]上的函数
h（x）=x/（1+x^4）,g（x）=arctanx
重复上述的实验看其结果如何。
（3）区间[a,b]上切比雪夫点的定义为：
xk=（b+a）/2+（（b-a）/2）cos（（2k-1）π/（2（n+1）））,k=1,2,^,n+1
以x1,x2^x（n+1）为插值节点构造上述各函数的拉格朗日插值多项式，比较其结果。
实验过程：
程序：
多项式插值的震荡现象（实验2.1）
form=1:
6
subplot（2,3,m）%把窗口分割成2*3大小的窗口
largrang（6*m）%对largrang函数进行运行
ifm==1
title（'longn=6'）
elseifm==2
title（'longn=12'）
elseifm==3
title（'longn=18'）
elseifm==4
title（'longn=24'）
elseifm==5
title（'longn=30'）
elseifm==6
title（'longn=36'）
end%对每个窗口分别写上标题为插值点的个数
end
保存为：
chazhi.m
functionlargrang（longn）
mm=input（'pleaseinputmm（运行第几个函数就输入mm为几）:
mm='）
ifmm==1%d表示定义域的边界值
d=1;
elseifmm==2||mm==3
d=5;
end
x0=linspace（-d,d,longn）;%x的节点
ifmm==1
y0=1./（1.+25.*x0.^2）;
elseifmm==2
y0=x0./（1.+x0.^4）;
elseifmm==3
y0=atan（x0）;
end
x=sym（'x'）;n=length（x0）;s=0.0;
fork=1:
n
p=1.0;
forj=1:
n
ifj~=k
p=p*（x-x0（j））/（x0（k）-x0（j））;
end
end
s=p*y0（k）+s;
end
y=s;
ifmm==1
ezplot（'1/（1+25*x^2）'）
elseifmm==2
ezplot（'x/（1+x^4）'）
elseifmm==3
ezplot（'atan（x）'）
end
holdon
ezplot（y,[-d,d]）
holdoff
保存为：
largrang.m
数值实验结果与分析：
对于第一个函数f（x）=1/（1+25x2）
对于第二个函数h（x）=x/（1+x4）
对于第三个函数g（x）=arctan（x）
讨论：
通过对三个函数得出的largrang插值多项式并在数学软件中的运行，得出函数图象，说明了对函数的支点不是越多越好，而是在函数的两端而言支点越多，而largrang插值多项式不是更加靠近被逼近的函数，反而更加远离函数，在函数两端的跳动性更加明显，argrang插值多项式对函数不收敛。
实验总结：
利用MATLAB来进行函数的largrang插值多项式问题的实验，虽然其得出的结果是有误差的，但是增加支点的个数进行多次实验，可以找出函数的largrang插值多项式的一般规律，当支点增加时，largrang插值多项式对函数两端不收敛，不是更加逼近，而是更加远离，跳动性更强。
所以对于函数的largrang插值多项式问题可以借助于MATLAB来进行问题的分析，得到比较准确的实验结规律。
学号：
06450210
：
万轩
实验五解线性方程组的直接方法
实验5.1（主元的选取与算法的稳定性）
问题提出：
Gauss消去法是我们在线性代数中已经熟悉的。
但由于计算机的数值运算是在一个有限的浮点数集合上进行的，如何才能确保Gauss消去法作为数值算法的稳定性呢？
Gauss消去法从理论算法到数值算法，其关键是主元的选择。
主元的选择从数学理论上看起来平凡，它却是数值分析中十分典型的问题。
实验容：
考虑线性方程组
编制一个能自动选取主元，又能手动选取主元的求解线性方程组的Gauss消去过程。
实验要求：
（1）取矩阵
，则方程有解
。
取n=10计算矩阵的条件数。
让程序自动选取主元，结果如何？
（2）现选择程序中手动选取主元的功能。
每步消去过程总选取按模最小或按模尽可能小的元素作为主元，观察并记录计算结果。
若每步消去过程总选取按模最大的元素作为主元，结果又如何？
分析实验的结果。
（3）取矩阵阶数n=20或者更大，重复上述实验过程，观察记录并分析不同的问题与消去过程中选择不同的主元时计算结果的差异，说明主元素的选取在消去过程中的作用。
（4）选取其他你感兴趣的问题或者随机生成矩阵，计算其条件数。
重复上述实验，观察记录并分析实验结果。
实验过程：
程序：
建立M文件：
functionx=gauss（n,r）
n=input（'请输入矩阵A的阶数:
n='）
A=diag（6*ones（1,n））+diag（ones（1,n-1）,1）+diag（8*ones（1,n-1）,-1）
b=A*ones（n,1）
p=input（'条件数对应的数是p-数：
p='）
pp=cond（A,p）
pause
[m,n]=size（A）;
nb=n+1;Ab=[Ab]
r=input（'请输入是否为手动，手动输入1，自动输入0：
r='）
fori=1:
n-1
ifr==0
[pivot,p]=max（abs（Ab（i:
n,i）））;
ip=p+i-1;
ifip~=i
Ab（[iip],:
）=Ab（[ipi],:
）;disp（Ab）;pause
end
end
ifr==1
i=i
ip=input（'输入i列所选元素所处的行数：
ip='）;
Ab（[iip],:
）=Ab（[ipi],:
）;disp（Ab）;pause
end
pivot=Ab（i,i）;
fork=i+1:
n
Ab（k,i:
nb）=Ab（k,i:
nb）-（Ab（k,i）/pivot）*Ab（i,i:
nb）;
end
disp（Ab）;pause
end
x=zeros（n,1）;x（n）=Ab（n,nb）/Ab（n,n）;
fori=n-1:
-1:
1
x（i）=（Ab（i,nb）-Ab（i,i+1:
n）*x（i+1:
n））/Ab（i,i）;
end
数值实验结果与分析：
⑴取矩阵A的阶数:
n=10，自动选取主元：
>>formatlong
>>gauss
请输入矩阵A的阶数:
n=10
n=10
条件数对应的数是p-数：
p=1
p=1
pp=2.0000e+003
请输入是否为手动，手动输入1，自动输入0：
r=0
r=0
⑵取矩阵A的阶数:
n=10，手动选取主元：
①选取绝对值最大的元素为主元：
>>gauss
请输入矩阵A的阶数:
n=10
n=10
条件数对应的数是p-数：
p=2
p=2
pp=1.3903e+003
请输入是否为手动，手动输入1，自动输入0：
r=1
r=1
ans=1111111111
②选取绝对值最小的元素为主元：
>>gauss
请输入矩阵A的阶数:
n=10
n=10
条件数对应的数是p-数：
p=2
p=2
pp=1.3903e+003
请输入是否为手动，手动输入1，自动输入0：
r=1
r=1
ans=
1.0001.0001.000
1.0001.0001.000
0.9991.0010.998
1.003
⑶取矩阵A的阶数:
n=20，手动选取主元：
1选取绝对值最大的元素为主元：
>>gauss
请输入矩阵A的阶数:
n=20
条件数对应的数是p-数：
p=1
p=1
pp=2.0000e+006
ans=11111111111111111111
2选取绝对值最小的元素为主元：
>>gauss
请输入矩阵A的阶数:
n=20.
n=20
条件数对应的数是p-数：
p=2
p=2
pp=1.1683e+006
请输入是否为手动，手动输入1，自动输入0：
r=1
r=1
ans=
1.0001.0001.0001.0001.0001.0001.0010.9971.006
0.9891.0230.955
1.0900.8211.352
0.3181.2730.817
1.910
⑷将M文件中的第三行：
A=diag（6*ones（1,n））+diag（ones（1,n-1）,1）+diag（8*ones（1,n-1）,-1）
改为：
A=hilb（n）
①>>gauss
请输入矩阵A的阶数:
n=7
n=7
条件数对应的数是p-数：
p=1
p=1
pp=9.0030e+008
请输入是否为手动，手动输入1，自动输入0：
r=1
r=1
ans=
1.0510.2511.354
0.1331.8050.181
1.337
②>>gauss
请输入矩阵A的阶数:
n=7
n=7
条件数对应的数是p-数：
p=2
p=2
pp=4.7072e+008
请输入是否为手动，手动输入1，自动输入0：
r=1
r=1
ans=
0.8691.3370.299
1.1430.0381.825
0.491
该问题在主元选取与算出结果有着很大的关系，取绝对值大的元素作为主元比取绝对值小的元素作为主元时产生的结果比较准确，即选取绝对值小的主元时结果产生了较大的误差，条件数越大产生的误差就越大。
讨论：
在gauss消去法解线性方程组时，主元的选择与算法的稳定性有密切的联系，选取绝对值大的元素作为主元比绝对值小的元素作为主元时对结果产生的误差较小。
条件数越大对用gauss消去法解线性方程组时，对结果产生的误差就越大。
实验总结：
对用gauss消去法解线性方程组时，主元的选取与算法的稳定性有密切的联系，选取适当的主元有利于得出稳定的算法，在算法的过程中，选取绝对值较大的主元比选取绝对值较小的主元更有利于算法的稳定，选取绝对值最大的元素作为主元时，得出的结果相对较准确较稳定。
条件数越小，对用这种方法得出的结果更准确。
在算除法的过程中要尽量避免使用较小的数做为除数，以免发生结果数量级加大，使大数吃掉小数，产生舍入误差。
学号：
06450210
：
万轩
实验5.2（线性代数方程组的性态与条件数的估计）
问题提出：
理论上，线性代数方程组
的摄动满足
矩阵的条件数确实是对矩阵病态性的刻画，但在实际应用中直接计算它显然不现实，因为计算
通常要比求解方程
还困难。
实验容：
Matlab中提供有函数“condest”可以用来估计矩阵的条件数，它给出的是按1-数的条件数。
首先构造非奇异矩阵A和右端，使得方程是可以精确求解的。
再人为地引进系数矩阵和右端的摄动
，使得
充分小。
实验要求：
（1）假设方程Ax=b的解为x，求解方程
，以1-数，给出
的计算结果。
（2）选择一系列维数递增的矩阵（可以是随机生成的），比较函数“condest”所需机器时间的差别.考虑若干逆是已知的矩阵，借助函数“eig”很容易给出cond2（A）的数值。
将它与函数“cond（A,2）”所得到的结果进行比较。
（3）利用“condest”给出矩阵A条件数的估计，针对
（1）中的结果给出
的理论估计，并将它与
（1）给出的计算结果进行比较，分析所得结果。
注意，如果给出了cond（A）和
的估计，马上就可以给出
的估计。
（4）估计著名的Hilbert矩阵的条件数。
实验过程：
程序：
⑴
n=input（'pleaseinputn:
n='）%输入矩阵的阶数
a=fix（100*rand（n））+1%随机生成一个矩阵a
x=ones（n,1）%假设知道方程组的解全为1
b=a*x%用矩阵a和以知解得出矩阵b
data=rand（n）*0.00001%随即生成扰动矩阵data
datb=rand（n,1）*0.00001%随即生成扰动矩阵datb
A=a+data
B=b+datb
xx=geshow（A,B）%解扰动后的解
x0=norm（xx-x,1）/norm（x,1）%得出
的理论结果
保存为：
fanshu.m
functionx=geshow（A,B）%用高斯消去法解方程组
[m,n]=size（A）;
nb=n+1;AB=[AB];
fori=1:
n-1
pivot=AB（i,i）;
fork=i+1:
n
AB（k,i:
nb）=AB（k,i:
nb）-（AB（k,i）/pivot）*AB（i,i:
nb）;
end
end
x=zeros（n,1）;
x（n）=AB（n,nb）/AB（n,n）;
fori=n-1:
-1:
1
x（i）=（AB（i,nb）-AB（i,i+1:
n）*x（i+1:
n））/AB（i,i）;
end
保存为：
geshow.m
⑵
functioncond2（A）%自定义求二阶条件数
B=A'*A;
[V1,D1]=eig（B）;
[V2,D2]=eig（B^（-1））;
cond2A=sqrt（max（max（D1）））*sqrt（max（max（D2）））
end
保存为：
cond2.m
formatlong
forn=10:
10:
100
n=n%n为矩阵的阶
A=fix（100*randn（n））;%随机生成矩阵A
condestA=condest（A）%用condest求条件数
cond2（A）%用自定义的求条件数
condA2=cond（A,2）%用cond求条件数
pause%运行一次暂停
end
保存为：
shiyan52.m
⑶
n=input（'pleaseinputn:
n='）%输入矩阵的阶数
a=fix（100*rand（n））+1;%随机生成一个矩阵a
x=ones（n,1）;%假设知道方程组的解全为1
b=a*x;%用矩阵a和以知解得出矩阵b
data=rand（n）*0.00001;%随即生成扰动矩阵data
datb=rand（n,1）*0.00001;%随即生成扰动矩阵datb
A=a+data;
B=b+datb;
xx=geshow（A,B）;%利用第一小问的geshow.m求出解阵
x0=norm（xx-x,1）/norm（x,1）%得出
的理论结果
x00=cond（A）/（1-norm（inv（A））*norm（xx-x））*（norm（（xx-x））/（norm（A））+norm（datb）/norm（B））%得出
的估计值
datx=abs（x0-x00）%求两者之间的误差
保存为：
sy5_2.m
⑷
formatlong
forn=4:
11
n=n%n为矩阵的阶数
Hi=hilb（n）;%生成Hilbert矩阵
cond1Hi=cond（Hi,1）%求Hilbert矩阵得三种条件数
cond2Hi=cond（Hi,2）
condinfHi=cond（Hi,inf）
pause
end
数值实验结果与分析：
⑴>>fanshu
pleaseinputn:
n=6
n=
6
a=
142516881989
329385489260
144088501316
23521929232
401010073724
1437227701
x=
111111
b=
251410221157218187
data=
1.0e-005*
0.9100.0500.5900.5740.9470.783
0.2500.0220.9020.9110.7810.561
0.4630.3280.6930.1900.1210.222
0.8760.5610.0110.6950.4870.009
0.1700.6140.6720.60
升级会员 