数值分析实验报告1Word格式.docx
《数值分析实验报告1Word格式.docx》由会员分享，可在线阅读，更多相关《数值分析实验报告1Word格式.docx（33页珍藏版）》请在冰点文库上搜索。
的变化更敏感
思考题一：
（上述实验的改进）
在上述实验中我们会发现用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
利用符号函数：
（思考题一）
y=poly2sym（a）;
rr=solve（y）
8
ess=10^（-6-m）;
ve=zeros（1,21）;
a=poly（1:
20）+ve;
y=poly2sym（a）;
r=solve（y）;
数值实验结果及分析：
formatlong
-6-mn
-7
-8
-9
-10
2
3
0.923
4
5
6
7
10
11
12
13
14
15
16
17
18
19
20
-11
-12
-13
-14
讨论：
利用这种方法进行这类实验，可以很精确的扰动敏感性的一般规律。
即当对扰动项的系数越来越小时，对其多项式扰动的结果也就越来越小，即扰动敏感性与扰动项的系数成正比，扰动项的系数越大，对其根的扰动敏感性就越明显，当扰动的系数一定时，扰动敏感性与扰动的项的幂数成正比，扰动的项的幂数越高，对其根的扰动敏感性就越明显。
实验总结：
利用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:
subplot（2,3,m）%把窗口分割成2*3大小的窗口
largrang（6*m）%对largrang函数进行运行
ifm==1
title（'
longn=6'
）
elseifm==2
longn=12'
elseifm==3
longn=18'
elseifm==4
longn=24'
elseifm==5
longn=30'
elseifm==6
longn=36'
end%对每个窗口分别写上标题为插值点的个数
保存为：
chazhi.m
functionlargrang（longn）
mm=input（'
pleaseinputmm（运行第几个函数就输入mm为几）:
mm='
ifmm==1%d表示定义域的边界值
d=1;
elseifmm==2||mm==3
d=5;
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）;
x=sym（'
x'
）;
n=length（x0）;
s=0.0;
fork=1:
n
p=1.0;
forj=1:
ifj~=k
p=p*（x-x0（j））/（x0（k）-x0（j））;
s=p*y0（k）+s;
y=s;
ezplot（'
1/（1+25*x^2）'
x/（1+x^4）'
atan（x）'
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来进行问题的分析，得到比较准确的实验结规律。
实验五解线性方程组的直接方法
实验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
ifr==1
i=i
ip=input（'
输入i列所选元素所处的行数：
ip='
pivot=Ab（i,i）;
fork=i+1:
Ab（k,i:
nb）=Ab（k,i:
nb）-（Ab（k,i）/pivot）*Ab（i,i:
nb）;
disp（Ab）;
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）;
⑴取矩阵A的阶数:
n=10，自动选取主元：
>
formatlong
gauss
n=10
n=10
p=1
p=1
r=0
r=0
⑵取矩阵A的阶数:
n=10，手动选取主元：
①选取绝对值最大的元素为主元：
p=2
p=2
pp=
r=1
r=1
ans=1111111111
②选取绝对值最小的元素为主元：
pp=03
ans=
1.000000000000001.000000000000001.00000000000000
1.00000000000001
1.00000000000003
⑶取矩阵A的阶数:
n=20，手动选取主元：
1选取绝对值最大的元素为主元：
n=20
pp=
ans=11111111111111111111
2选取绝对值最小的元素为主元：
n=20.
n=20
1.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000011.00000000000006
999891.00000000000023
1.000000000000901.00000000000352
1.00000000001273
1.00000000002910
⑷将M文件中的第三行：
改为：
A=hilb（n）
①>
n=7
n=7
1.000000000000511.00000000031354
1.00000000268805
1.00000000084337
②>
1.00000000004337
1.000000001211431.00000000152825
该问题在主元选取与算出结果有着很大的关系，取绝对值大的元素作为主元比取绝对值小的元素作为主元时产生的结果比较准确，即选取绝对值小的主元时结果产生了较大的误差，条件数越大产生的误差就越大。
在gauss消去法解线性方程组时，主元的选择与算法的稳定性有密切的联系，选取绝对值大的元素作为主元比绝对值小的元素作为主元时对结果产生的误差较小。
条件数越大对用gauss消去法解线性方程组时，对结果产生的误差就越大。
对用gauss消去法解线性方程组时，主元的选取与算法的稳定性有密切的联系，选取适当的主元有利于得出稳定的算法，在算法的过程中，选取绝对值较大的主元比选取绝对值较小的主元更有利于算法的稳定，选取绝对值最大的元素作为主元时，得出的结果相对较准确较稳定。
条件数越小，对用这种方法得出的结果更准确。
在算除法的过程中要尽量避免使用较小的数做为除数，以免发生结果数量级加大，使大数吃掉小数，产生舍入误差。
实验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矩阵的条件数。
⑴
pleaseinputn:
）%输入矩阵的阶数
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）%用高斯消去法解方程组
AB=[AB];
pivot=AB（i,i）;
AB（k,i:
nb）=AB（k,i:
nb）-（AB（k,i）/pivot）*AB（i,i:
x（n）=AB（n,nb）/AB（n,n）;
x（i）=（AB（i,nb）-AB（i,i+1:
n））/AB（i,i）;
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）））
cond2.m
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%运行一次暂停
shiyan52.m
⑶
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求出解阵
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
⑷
forn=4:
n=n%n为矩阵的阶数
Hi=hilb（n）;
%生成Hilbert矩阵
cond1Hi=cond（Hi,1）%求Hilbert矩阵得三种条件数
cond2Hi=cond（Hi,2）
condinfHi=cond（Hi,inf）
⑴>
fanshu
n=6
n=
6
a=
142516881989
329385489260
144088501316
23521929232
401010073724
1437227701
x=
111111
b=
251410221157218187
data=
1.0e-005*
datb=
1.0
升级会员 