数值分析.docx

数值分析.docx

《数值分析.docx》由会员分享，可在线阅读，更多相关《数值分析.docx（20页珍藏版）》请在冰豆网上搜索。

数值分析

数值分析实验报告

信息与计算科学02李志刚10092035

1.设

计算

，哪个接近

。

算法分析：

（matlab语言）用Taylor级数计算

。

计算的误差为δ1=

（其中ζ介于0到-5）,

计算的误差为δ2=

（其中η介于0到5），由于计算的误差较小，δ2的上界近似等于δ1的下界，故δ2误差较小。

也可以由交错级数收敛得慢得出

误差较大。

程序如下：

functionlizhigang

x=-5;

he=1;

jc=1;

n=24;

fork=1:

n;

jc=jc*k;

he=he+1/jc*x^k;

end

P=vpa（he,16）

P=0.006737963097703459

functionlzg

x=5;

he=1;

jc=1;

n=24;

fork=1:

n;

jc=jc*k;

he=he+1/jc*x^k;

end

P=vpa（he,16）;

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

{

for（intii=i+1;ii

{

doublerate=A[ii][i]/A[i][i];

for（intjj=i;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];

for（inti=0;i

{

Ab[i][i]=6.0;

}

for（inti=1;i

{

Ab[i][i-1]=8.0;

}

for（inti=1;i

{

Ab[i-1][i]=1.0;

}

Ab[0][n]=7;//初始化增广矩阵中的向量b

Ab[n-1][n]=14;

for（inti=1;i

{

Ab[i][n]=15.0;

}

func（Ab）;

for（inti=0;i

{

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

{

sum+=Ab[i][j]*x[j];

}

x[i]=（b[i]-sum）/Ab[i][i];

}

for（inti=0;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）{

finalintn=A.length;

for（inti=0;i

{

doubletem1=A[i][i];

inttem2=0;

for（intp=i+1;p

{

if（A[p][i]>tem1）

{

tem1=A[p][i];

tem2=p;

}

}

if（tem2!

=0）

{

for（intq=i;q

{

doubletemm=A[i][q];

A[i][q]=A[tem2][q];

A[tem2][q]=temm;

}

}

for（intii=i+1;ii

{

doublerate=A[ii][i]/A[i][i];

for（intjj=i;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];

for（inti=0;i

{

Ab[i][i]=6.0;

}

for（inti=1;i

{

Ab[i][i-1]=8.0;

}

for（inti=1;i

{

Ab[i-1][i]=1.0;

}

Ab[0][n]=7;//初始化增广矩阵中的向量b

Ab[n-1][n]=14;

for（inti=1;i

{

Ab[i][n]=15.0;

}

func（Ab）;

for（inti=0;i

{

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

{

sum+=Ab[i][j]*x[j];

}

x[i]=（b[i]-sum）/Ab[i][i];

}

for（inti=0;i

{

System.out.print（"x"+（i+1）+"="）;

System.out.print（x[i]+"\t"）;

if（（i+1）%3==0）System.out.println（）;

}

}

}

结果如下：

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;

end

x0

jacobitimes

%Gauss迭代

x1=zeros（4,1）;%构造初始向量

x0=ones（4,1）;

gausstimes=0;%Gauss迭代次数

while（norm（x1-x0）>1.0e-3）

x0=x1;

x1=Bg*x0+gg;

gausstimes=gausstimes+1;

end

x0

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第一列

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

end

p=[];%用p向量存储F的关键系数

fori=1:

6

p（1,i）=F（i,i+1）;

end

p=vpa（p,8）;

symsx;

symsy;

y=p

（1）;

fori=2:

6

j=2;

q=1;

while（j<=i）

q=q*（x-F（j-1,1））;

j=j+1;

end

y=y+p（i）*q;

end

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阶矩阵

F（:

1）=rot90（x10,3）;%初始化第一列

F（:

2）=rot90（y10,3）;%初始化第二列

forj=3:

12%完整初始化F

fori=（j-1）:

11

F（i,j）=（F（i,j-1）-F（i-1,j-1））/（F（i,1）-F（i-（j-2）,1））;

end

end

p=[];%用p向量存储F的关键系数

fori=1:

11

p（1,i）=F（i,i+1）;

end

p=vpa（p,8）;

symsx;

symsy;

y=p

（1）;

fori=2:

11

j=2;

q=1;

while（j<=i）

q=q*（x-F（j-1,1））;

j=j+1;

end

y=y+p（i）*q;

end

symsn10

N10=vpa（simplify（y）,3）

x0=linspace（-1,1,200）;

y0=subs（y,x,x0）;

plot（x0,y0）%作N10图像

text（0.6,1.5,'N10图像'）;

%以下是做三次样条插值函数

sx=linspace（-1,1,11）;%等距取11个插值结点

sy=1./（1+