数值分析Word格式文档下载.docx
《数值分析Word格式文档下载.docx》由会员分享，可在线阅读，更多相关《数值分析Word格式文档下载.docx（20页珍藏版）》请在冰豆网上搜索。
Q=1/P
Q=0.006737947000163260
结果如下：
故用
计算的结果更精确。
2.将不选列主元和列主元的Gauss消去法编成程序，并求解下面的84阶方程组。
最后将你的结果与方程精确解比较，并谈谈你对Gauss消去法的看法。
（Java语言）在不选列主元的Gauss消去法中可能会出现小主元，这就会引起其他元素数量级的剧增和舍入误差的增长，导致计算结果不可靠。
而选列主元的Gauss消去法有利于控制误差的增长，使计算结果更可靠。
（其精确解为x都为1）
不选列主元的Gauss程序如下:
publicclass第二题Guauss消去法
{
//构造一个处理增广矩阵的方法
publicstaticdouble[][]func（double[][]A）
{
finalintn=A.length;
for（inti=0;
i<
n-1;
i++）
{
for（intii=i+1;
ii<
ii++）
{
doublerate=A[ii][i]/A[i][i];
for（intjj=i;
jj<
n+1;
jj++）
{
A[ii][jj]=A[ii][jj]-rate*A[i][jj];
}
}
returnA;
}
publicstaticvoidmain（String[]args）
finalintn=84;
double[][]Ab=newdouble[n][n+1];
//增广矩阵
double[]b=newdouble[n];
i++）//初始化矩阵A
Ab[i][i]=6.0;
for（inti=1;
Ab[i][i-1]=8.0;
Ab[i-1][i]=1.0;
Ab[0][n]=7;
//初始化增广矩阵中的向量b
Ab[n-1][n]=14;
Ab[i][n]=15.0;
func（Ab）;
i++）//取出作用后的向量b
b[i]=Ab[i][n];
double[]x=newdouble[n];
//回带求解方程
x[n-1]=b[n-1]/Ab[n-1][n-1];
for（inti=n-2;
i>
=0;
i--）
doublesum=0;
for（intj=i+1;
j<
j++）
sum+=Ab[i][j]*x[j];
}
x[i]=（b[i]-sum）/Ab[i][i];
i++）//显示结果
System.out.print（"
x"
+（i+1）+"
="
）;
System.out.print（x[i]+"
\t"
if（（i+1）%3==0）System.out.println（）;
}
x1=1.0000000000000000x2=1.0000000000000002x3=0.9999999999999996
x4=1.0000000000000009x5=0.9999999999999982x6=1.0000000000000036
x7=0.9999999999999929x8=1.0000000000000142x9=0.9999999999999716
x10=1.0000000000000568x11=0.9999999999998863x12=1.0000000000002274
x13=0.9999999999995453x14=1.0000000000009095x15=0.999999999998181
x16=1.000000000003638x17=0.9999999999927245x18=1.000000000014551
x19=0.999999999970898x20=1.000000000058204x21=0.9999999998835918
x22=1.0000000002328164x23=0.9999999995343671x24=1.0000000009312657
x25=0.9999999981374685x26=1.000000003725063x27=0.9999999925498742
x28=1.0000000149002517x29=0.9999999701994966x30=1.0000000596010068
x31=0.9999998807979864x32=1.0000002384040272x33=0.9999995231919456
x34=1.0000009536161087x35=0.9999980927677825x36=1.000003814464435
x37=0.99999237107113x38=1.00001525785774x39=0.9999694842845201
x40=1.0000610314309597x41=0.9998779371380806x42=1.0002441257238388
x43=0.9995117485523225x44=1.0009765028953548x45=0.9980469942092913
x46=1.0039060115814138x47=0.9921879768371866x48=1.01562404632557
x49=0.9687519073490876x50=1.0624961853009154x51=0.8750076294018072
x52=1.2499847411818346x53=0.500030517694535x54=1.9999389643781136
x55=-0.999877927824961x56=4.99975585192486x57=-6.999511688949468
x58=16.99902331829793x59=-30.99804639819183x60=64.99609184276756
x61=-126.9921798710706x62=256.9843444842836x63=-510.96862793713626
x64=1024.9370117485487x65=-2046.873046994202x66=4096.742187976823
x67=-8190.468751907319x68=16383.875007629336x69=-32764.50003051746
x70=65531.00012207008x71=-131055.0004882807x72=262097.00195312407
x73=-524127.0078124981x74=1048001.0312499963x75=-2094975.1249999925
x76=4185857.499999985x77=-8355328.99999997x78=1.664512899999994E7
x79=-3.3028126999999E7x80=6.50077449999997E7x81=-1.2582143899999955E8
x82=2.34866688999999E8x83=-4.0262860699998E8x84=5.368381449999981E8
（大概从
处误差开始扩大,到了后面误差急剧增大。
）
选列主元的Gauss程序如下:
publicclass第二题列主元Guauss消去法
//构造一个用列主元Gauss消元法处理增广矩阵A（n*（n+1））的方法
publicstaticdouble[][]func（double[][]A）{
doubletem1=A[i][i];
inttem2=0;
for（intp=i+1;
p<
p++）
if（A[p][i]>
tem1）
{
tem1=A[p][i];
tem2=p;
if（tem2!
=0）
for（intq=i;
q<
q++）
doubletemm=A[i][q];
A[i][q]=A[tem2][q];
A[tem2][q]=temm;
i++）//显示结果
x1=1.0000000000000002x2=0.9999999999999998x3=1.0000000000000002
x4=0.9999999999999998x5=1.0000000000000002x6=0.9999999999999998
x7=1.0000000000000002x8=0.9999999999999998x9=1.0000000000000002
x10=0.9999999999999998x11=1.0000000000000002x12=0.9999999999999998
x13=1.0000000000000002x14=0.9999999999999998x15=1.0000000000000002
x16=0.9999999999999998x17=1.0000000000000002x18=0.9999999999999998
x19=1.0000000000000002x20=0.9999999999999998x21=1.0000000000000002
x22=0.9999999999999998x23=1.0000000000000002x24=0.9999999999999998
x25=1.0000000000000002x26=0.9999999999999998x27=1.0000000000000002
x28=0.9999999999999998x29=1.0000000000000002x30=0.9999999999999998
x31=1.0000000000000002x32=0.9999999999999998x33=1.0000000000000002
x34=0.9999999999999998x35=1.0000000000000002x36=0.9999999999999998
x37=1.0000000000000002x38=0.9999999999999998x39=1.0000000000000002
x40=0.9999999999999998x41=1.0000000000000002x42=0.9999999999999998
x43=1.0000000000000002x44=0.9999999999999998x45=1.0000000000000002
x46=0.9999999999999998x47=1.0000000000000002x48=0.9999999999999998
x49=1.0000000000000002x50=0.9999999999999998x51=1.0000000000000004
x52=0.9999999999999989x53=1.0000000000000024x54=0.9999999999999949
x55=1.0000000000000102x56=0.9999999999999793x57=1.0000000000000415
x58=0.9999999999999167x59=1.0000000000001665x60=0.9999999999996667
x61=1.0000000000006668x62=0.9999999999986662x63=1.0000000000026676
x64=0.9999999999946645x65=1.0000000000106712x66=0.9999999999786573
x67=1.0000000000426854x68=0.9999999999146294x69=1.0000000001707399
x70=0.9999999996585256x71=1.0000000006829282x72=0.999999998634227
x73=1.0000000027312126x74=0.9999999945389089x75=1.0000000109168468
x76=0.9999999781876492x77=1.0000000435393304x78=0.9999999132628243
x79=1.0000001721084117x80=0.9999996612469355x81=1.0000006556510925
x82=0.9999987761179607x83=1.000002098083496x84=0.9999972025553387
（我们从中可以看出，用列主元Gauss消去法能有效控制误差的增长，最终得出比较精确的结果。
3.用雅克比迭代法和高斯-赛德尔迭代法求解方程组
使误差小于
，并计较迭代次数。
分析：
（本题用matlab语言）首先分解系数矩阵为D,E,F三个矩阵，再分别用这两种迭代法的矩阵表示法计算出迭代矩阵B，然后验证B的谱半径是否小于1，最后进行迭代求解。
程序如下:
A=[2.52,0.95,1.25,-0.85;
0.39,1.69,-0.45,0.49;
0.55,-1.25,1.96,-0.98;
0.23,-1.15,-0.45,2.31];
b=[1.38;
-0.34;
0.67;
1.52];
D=diag（diag（A））;
E=-1.*tril（A,-1）;
F=-1.*triu（A,1）;
%Jacobi迭代法迭代矩阵
Bj=inv（D）*（E+F）;
gj=inv（D）*b;
%Gauss迭代法迭代矩阵
Bg=inv（D-E）*F;
gg=inv（D-E）*b;
%已经验证两种迭代法的迭代矩阵都小于1
%Jacobi迭代
x1=zeros（4,1）;
%迭代初始向量
x0=ones（4,1）;
jacobitimes=0;
%Jacobi迭代次数
while（norm（x1-x0）>
1.0e-3）
x0=x1;
x1=Bj*x0+gj;
jacobitimes=jacobitimes+1;
x0
jacobitimes
%Gauss迭代
%构造初始向量
gausstimes=0;
%Gauss迭代次数
x1=Bg*x0+gg;
gausstimes=gausstimes+1;
gausstimes
Jacobi迭代：
Gauss迭代：
从中可以看出，高斯-赛德尔迭代的收敛速度较快，仅为10次。
4．观察高次插值多项式的龙格现象。
给定函数
取等距结点，构造牛顿插值多项式
和
及三次样条插值函数
分别将三种插值多项式与
的曲线画在同一个坐标系上进行比较。
（本题采用matlab语言）构造牛顿多项式时主要是构造差商矩阵；
计算
时，由于是等距结点，故简化了计算。
x5=linspace（-1,1,6）;
%6个插值结点
y5=1./（1+25*x5.^2）;
%6个在插值结点处的函数值
F=zeros（6,7）;
%构造差商矩阵F，5*6阶矩阵
F（:
1）=rot90（x5,3）;
%初始化F第一列
2）=rot90（y5,3）;
%初始化F第二列
forj=3:
7%完整初始化F
fori=（j-1）:
6
F（i,j）=（F（i,j-1）-F（i-1,j-1））/（F（i,1）-F（i-（j-2）,1））;
end
p=[];
%用p向量存储F的关键系数
fori=1:
p（1,i）=F（i,i+1）;
p=vpa（p,8）;
symsx;
symsy;
y=p
（1）;
fori=2:
j=2;
q=1;
while（j<
=i）
q=q*（x-F（j-1,1））;
j=j+1;
y=y+p（i）*q;
symsN5
N5=vpa（simplify（y）,3）
x0=linspace（-1,1,200）;
y0=subs（y,x,x0）;
plot（x0,y0,'
b'
）%做N5图像
x=linspace（-1,1,1000）;
y=1./（1+25*x.^2）;
holdon
plot（x,y,'
r*'
'
markersize'
1）%画原函数图像
gridon
text（0.1,1,'
1/（1+25x^2）图像'
text（0.05,0.6,'
图像'
%作N10，其基本算法同上
x10=linspace（-1,1,11）;
%11个插值结点
y10=1./（1+25*x10.^2）;
%11个在插值结点处的函数值
F=zeros（11,12）;
%构造差商矩阵F，11*12阶矩阵
1）=rot90（x10,3）;
%初始化第一列
2）=rot90（y10,3）;
%初始化第二列
12%完整初始化F
11
symsn10
N10=vpa（simplify（y）,3）
plot（x0,y0）%作N10图像
text（0.6,1.5,'
N10图像'
%以下是做三次样条插值函数
sx=linspace（-1,1,11）;
%等距取11个插值结点
sy=1./（1+
升级会员 