搜档网
当前位置:搜档网 › 数值分析实验报告

数值分析实验报告

数值分析实验报告
数值分析实验报告

实验五 解线性方程组的直接方法

实验5.1 (主元的选取与算法的稳定性) 问题提出:Gauss 消去法是我们在线性代数中已经熟悉的。但由于计算机的数值运算是在一个有限的浮点数集合上进行的,如何才能确保Gauss 消去法作为数值算法的稳定性呢?Gauss 消去法从理论算法到数值算法,其关键是主元的选择。主元的选择从数学理论上看起来平凡,它却是数值分析中十分典型的问题。

实验内容:考虑线性方程组

编制一个能自动选取主元,又能手动选取主元的求解线性方程组的Gauss 消去过程。 实验要求:

(1)取矩阵??

?

??

??

?????????=????????????????=1415157,6816816816 b A ,则方程有解T x )1,,1,1(* =。取n=10计算矩阵的

条件数。让程序自动选取主元,结果如何?

(2)现选择程序中手动选取主元的功能。每步消去过程总选取按模最小或按模尽可能小的元素作为主元,观察并记录计算结果。若每步消去过程总选取按模最大的元素作为主元,结果又如何?分析实验的结果。

(3)取矩阵阶数n=20或者更大,重复上述实验过程,观察记录并分析不同的问题及消去过程中选择不同的主元时计算结果的差异,说明主元素的选取在消去过程中的作用。

(4)选取其他你感兴趣的问题或者随机生成矩阵,计算其条件数。重复上述实验,观察记录并分析实验结果。 思考题一:(Vadermonde 矩阵)设

??

??????????????????????=?

?

?

?????????????=∑∑∑∑====n i i n n i i n

i i n i i n n n n n

n n

x x x x b x x x x x x x x x x x x A 0020

10022222121102001111 ,,

其中,n k k x k ,,1,0,1.01 =+=,

(1)对n=2,5,8,计算A 的条件数;随n 增大,矩阵性态如何变化?

(2)对n=5,解方程组Ax=b ;设A 的最后一个元素有扰动10-4,再求解Ax=b (3)计算(2)扰动相对误差与解的相对偏差,分析它们与条件数的关系。

(4)你能由此解释为什么不用插值函数存在定理直接求插值函数而要用拉格朗日或牛顿插值法的原因吗?

相关MATLAB 函数提示:

zeros(m,n) 生成m 行,n 列的零矩阵

ones(m,n) 生成m 行,n 列的元素全为1的矩阵 eye(n) 生成n 阶单位矩阵

rand(m,n) 生成m 行,n 列(0,1)上均匀分布的随机矩阵 diag(x) 返回由向量x 的元素构成的对角矩阵

tril(A) 提取矩阵A 的下三角部分生成下三角矩阵

triu(A) 提取矩阵A的上三角部分生成上三角矩阵rank(A) 返回矩阵A的秩

det(A) 返回方阵A的行列式

inv(A) 返回可逆方阵A的逆矩阵

[V,D]=eig(A) 返回方阵A的特征值和特征向量norm(A,p) 矩阵或向量的p范数

cond(A,p) 矩阵的条件数

[L,U,P]=lu(A) 选列主元LU分解

R=chol(X) 平方根分解

Hi=hilb(n) 生成n阶Hilbert矩阵

5.1实验过程:

5.1.1程序:

function x=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)

for i=1:4

p=input('条件数对应的范数是p-范数:p=')

pp=cond(A,p)

end

pause

[m,n]=size(A);

nb=n+1;Ab=[A b]

r=input('请输入是否为手动,手动输入1,自动输入0:r=')

for i=1:n-1

if r==0

[pivot,p]=max(abs(Ab(i:n,i)));

ip=p+i-1;

if ip~=i

Ab([i ip],:)=Ab([ip i],:);disp(Ab); pause

end

end

if r==1

i=i

ip=input('输入i列所选元素所处的行数:ip=');

Ab([i ip],:)=Ab([ip i],:);disp(Ab); pause

end

pivot=Ab(i,i);

for k=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);

for i=n-1:-1:1

x(i)=(Ab(i,nb)-Ab(i,i+1:n)*x(i+1:n))/Ab(i,i);

end

5.1.2实验结果如下:

1.按照实验要求一:取矩阵A的阶数:n=10且自动选取主元,程序结果运行如下:

(2)现选择程序中手动选取主元的功能,观察并记录计算结果。

①选取绝对值最大的元素为主元:

程序运行开始如第一问的截图也是求范数故这里不在给出。

The answer is :1 1 1 1 1 1 1 1 1 1

②选取绝对值最小的元素为主元:

The answer is:

1.0e+003*(INF 0.007 0.0057 -0.0410 -0.0303 0.3430 0.2577 -

2.7290 -2.0463 2.7308)

⑶取矩阵A的阶数:n=20,手动选取主元:

①选取绝对值最大的元素为主元:

The answer is :1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1

②选取绝对值最小的元素为主元:

The answer is:

1.0e+007*(-Inf 0.0000 0.0000 -0.0000 -0.0000 0.0000 0.0000 -0.0003 -0.0002 0.0022

0.0016 -0.0175 -0.0131 0.1398 0.1049 -1.1185 -0.8389 8.9478 6.7109 -8.9478)

⑷修改程序如下:

function x=gaussong(n,r)

n=input('请输入矩阵A的阶数:n=')

A=hilb(n)

b=A*ones(n,1)

for i=1:4

p=input('条件数对应的范数是p-范数:p=')

pp=cond(A,p)

end

pause

[m,n]=size(A);

nb=n+1;Ab=[A b]

r=input('请输入是否为手动,手动输入1,自动输入0:r=')

for i=1:n-1

if r==0

[pivot,p]=max(abs(Ab(i:n,i)));

ip=p+i-1;

if ip~=i

Ab([i ip],:)=Ab([ip i],:);disp(Ab); pause

end

end

if r==1

i=i

ip=input('输入i列所选元素所处的行数:ip=');

Ab([i ip],:)=Ab([ip i],:);disp(Ab); pause

end

pivot=Ab(i,i);

for k=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);

for i=n-1:-1:1

x(i)=(Ab(i,nb)-Ab(i,i+1:n)*x(i+1:n))/Ab(i,i);

end

①所求范数为:

自动输入结果为:

ans =

1.0000 1.00001 .0000 1.0000 1.0002 0.9996 1.0007 0.9993 1.0004 0.9999

②选取绝对值最大的元素为主元结果为:

The answer is :

NaN NaN NaN NaN NaN Inf -Inf -Inf 281.3945 -283.3708 ③选取绝对值最小的元素为主元结果为: The answer is :

NaN NaN NaN -Inf -5.8976 -1.9243 -2.0291 -4.9972 23.4548 -11.1012 5.1.3 对实验结果进行分析:

5.1.3.1 对实验要求一的结果进行分析:对于Gauss 消去法就是用行的初等变换将原线性方程组系数矩阵转化为简单形式,从而进行求解,缺点是迭代次数可能较多,效率不高,且在消去过程中不可以将主元素很小的做除数,否则将导致其他元素数量级的严重增长和舍入误差的扩散,使得计算解不可靠。

5.1.3.2 对实验要求二的结果进行分析:通过每次选取最大或最小的主元可以发现取绝对值大的元素作为主元比取绝对值小的元素作为主元时产生的结果比较准确,即选取绝对值小的主元时结果产生了较大的误差,条件数越大产生的误差就越大。所以应尽量避免很小的数作为除数。

5.1.3.3 对实验要求三的结果进行分析:此要求是对要求一和要求二的一个延续,通过实验结果可以看出若采用很小的数作为主元迭代次数越多导致的结果越不可靠,甚至出现错误。

5.1.3.4 对实验要求四的结果进行分析:对新矩阵进行实验发现依然符合上述规律,可以知道,在进行迭代时主元的选择与算法的稳定性有密切的联系选取绝对值大的元素作为主元比绝对值小的元素作为主元时对结果产生的误差较小。条件数越大对用gauss 消去法解线性方程组时,对结果产生的误差就越大。 5.1.4实验总结:

1.在用gauss 消去法解线性方程组时,主元的选取与算法的稳定性有密切的联系,选取适当的主元有利于得出稳定的算法,

2.在算法的过程中,选取绝对值较大的主元比选取绝对值较小的主元更有利于算法的稳定,选取绝对值最大的元素作为主元时,得出的结果相对较准确较稳定。

3.条件数越小,对用这种方法得出的结果更准确。

4.在算除法的过程中要尽量避免使用较小的数做为除数,以免发生结果数量级加大,使大数吃掉小数,产生舍入误差。

实验5.2(线性代数方程组的性态与条件数的估计)

问题提出:理论上,线性代数方程组b Ax =的摄动满足

矩阵的条件数确实是对矩阵病态性的刻画,但在实际应用中直接计算它显然不现实,因为计算1-A 通常要比求解方程b Ax =还困难。

实验内容:Matlab 中提供有函数“condest ”可以用来估计矩阵的条件数,它给出的是按1-范数的条件数。首先构造非奇异矩阵A 和右端,使得方程是可以精确求解的。再人为地引进系数矩阵和右端的摄动b A ??和,使得b A ??和充分小。 实验要求:

(1)假设方程Ax=b 的解为x ,求解方程b b x

A A ?+=?+?)(,以1-范数,给出x

x x x

x -=

??的计算

结果。

(2)选择一系列维数递增的矩阵(可以是随机生成的),比较函数“condest ”所需机器时间的差别.考虑若干逆是已知的矩阵,借助函数“eig ”很容易给出cond 2(A)的数值。将它与函数“cond(A,2)”所得到的结果进行比较。

(3)利用“condest ”给出矩阵A 条件数的估计,针对(1)中的结果给出

x

x ?的理论估计,并将

它与(1)给出的计算结果进行比较,分析所得结果。注意,如果给出了cond(A)和A的估计,马上就可以给出1 A的估计。

(4)估计著名的Hilbert矩阵的条件数。

5.2 实验过程如下:

5.2.1.1 实验要求一的程序如下:

function n=jisuan(n)

a=fix(100*rand(n))+1

x=ones(n,1)

b=a*x

data=rand(n)*0.00001

datb=rand(n,1)*0.00001

A=a+data

B=b+datb

x0=get(A,B)

x1=norm(x0-x,1)/norm(x,1)

function x=get(A,B)

[m,n]=size(A);

nb=n+1;AB=[A B];

for i=1:n-1

pivot=AB(i,i);

for k=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);

for i=n-1:-1:1

x(i)=(AB(i,nb)-AB(i,i+1:n)*x(i+1:n))/AB(i,i);

End

5.2.1.2 实验要求一程序运行结果如下:

系数矩阵a为:

70 18 14 29 52 53

63 98 2 47 34 65

80 28 90 7 44 21

96 26 20 99 23 38

53 88 30 59 58 79

89 74 67 43 77 69

加扰动后的系数矩阵A为:

70.0000046 18.0000042 14.0000099 29.0000032 52.0000044 53.0000013

63.0000057 98.0000030 2.0000079 47.0000096 34.0000093 65.0000021

80.0000079 28.0000087 90.0000044 7.0000073 44.0000068 21.0000061

96.0000006 26.0000002 20.0000050 99.0000041 23.0000021 38.0000063

53.0000060 88.0000077 30.0000021 59.0000074 58.0000084 79.0000037

89.0000005 74.0000097 67.0000064 43.0000027 77.0000063 69.0000058

b值为:

236

309

270

302

367

419

加扰动后的b值为:

236.00000451

309.00000044

270.00000027

302.00000313

367.00000013

419.00000384 data的值为:

4.610951267E-06 4.1537486

04E-06

9.9008259

26E-06

3.2003557

75E-06

4.3992430

96E-06

1.3377274

85E-06

5.678287124E-

06 3.0499867

7E-06

7.8886169

22E-06

9.6009860

04E-06

9.3338010

82E-06

2.0713272

96E-06

7.942106514E-06 8.7436717

16E-06

4.3865853

38E-06

7.2663176

66E-06

6.8333232

43E-06

6.0719894

45E-06

5.91825935E-0 7 1.5009498

7E-07

4.9831130

34E-06

4.1195320

82E-06

2.1255986

43E-06

6.2988784

88E-06

6.028690857E-06

7.6795039

E-06

2.1396333

16E-06

7.4456578

31E-06

8.3923824

03E-06

3.7047682

61E-06

5.02688037E-0 7 9.7084493

93E-06

6.4349228

79E-06

2.6794725

07E-06

6.287846E

-06

5.7514777

9E-06

datb的值为:

4.514248268E-06

4.38953253E-07

2.7185123E-07

3.126850481E-06

1.28625747E-07

3.839672885E-06 xx的值为:

0.99999968449

1.0000003878

0.99999914348

1.0000006517

1.0000022156

0.9999975453

x0的值为:

1.146958549E-06 x1的值为:

6.8990e-007 5.2.1.3实验结果为:

x

x x x

x -=

??的计算结果为:6.8990e-007

5.2.2.1 实验要求二的程序如下:

function cond2(A) B=A'*A;

[V1,D1]=eig(B); [V2,D2]=eig(B^(-1));

cond2A=sqrt(max(max(D1)))*sqrt(max(max(D2))) end

for n=10:10:100

n=n A=fix(100*randn(n)); condestA=condest(A) cond2(A) condA2=cond(A,2) pause end

5.2.2.2 实验要求二的程序运行结果如下: N condestA cond2A condA2 10 1.46E+02 42.507830441 42.507830441 20 4.59E+02 1.23E+02 1.23E+02 30 4.05E+02 79.264954885 79.264954885 40 2.21E+03 4.26E+02 4.26E+02 50 3.02E+03 4.08E+02 4.08E+02 60 4.75E+03 7.78E+02 7.78E+02 70 4.69E+03 5.14E+02 5.14E+02 80 5.47E+03 4.89E+02 4.89E+02 90 5.66E+03 5.50E+02 5.50E+02 100 4.47E+03 5.06E+02 5.06E+02 5.2.3.1 实验要求三的程序如下: function bijiao(n)

a=fix(100*rand(n))+1; x=ones(n,1); b=a*x; data=rand(n)*0.00001; datb=rand(n,1)*0.00001; A=a+data; B=b+datb;

xx=geshow(A,B); x1=norm(xx-x,1)/norm(x,1)

x2=cond(A)/(1-norm(inv(A))*norm(xx-x))*(norm((xx-x))/(norm(A))+norm(datb)/norm(B)) datx=abs(x1-x2)

5.2.3.2 实验要求三的程序运行结果如下: 给出对

x

x x x

x -=

??的估计是:7.310559817408125e-007

x

x x x

x -=

??的理论结果是: 3.828481757617297e-007

结果相差: 3.482078059790828e-007

5.2.4.1 实验要求四的程序如下: for n=4:11

n=n

Hi=hilb(n); cond1Hi=cond(Hi,1) cond2Hi=cond(Hi,2) condinfHi=cond(Hi,inf) pause end

5.2.4.2 实验要求四的程序运行结果如下:

n cond1Hi cond2Hi condinfHi 4 2.84E+04 1.55E+04 2.84E+04 5 9.44E+05 4.77E+05 9.44E+05 6 2.91E+07 1.50E+07 2.91E+07 7 9.85E+08 4.75E+08 9.85E+08 8 3.39E+10 1.53E+10 3.39E+10 9 1.10E+12 4.93E+11 1.10E+12 10 3.54E+13 1.60E+13 3.54E+13 11 1.23E+15 5.22E+14

1.23E+15

讨论:

线性代数方程组的性态与条件数有着很重要的关系,既矩阵的条件数是刻画矩阵性质的一个重要的依据,条件数越大,矩阵“病态”性越严重,在解线性代数方程组的过程中较容易产生比较大的误差,则在实际问题的操作过程中,我们必须要减少对条件数来求解,把条件数较大的矩阵化成条件数较小的矩阵来进行求解。

实验总结:

在本次实验中,使我们知道了矩阵条件数对线性代数方程组求解的影响,条件数越大,对最后解的影响的越大,hilbert 矩阵是一个很”病态”的矩阵,他的条件数随着阶数的增加而增大,每增加一阶,条件数就增大一个数量级,在求解的过程中要尽量避免hilbert 矩阵

实验六 解线性方程组的迭代法

实验6.1(病态的线性方程组的求解)

问题提出:理论的分析表明,求解病态的线性方程组是困难的。实际情况是否如此,会出现怎样的现象呢?

实验内容:考虑方程组Hx=b 的求解,其中系数矩阵H 为Hilbert 矩阵, 这是一个著名的病态问题。通过首先给定解(例如取为各个分量均为1)再计算出右端b 的办法给出确定的问题。

实验要求:

(1)选择问题的维数为6,分别用Gauss 消去法、J 迭代法、GS 迭代法和SOR 迭代法求解方程组,其各自的结果如何?将计算结果与问题的解比较,结论如何?

(2)逐步增大问题的维数,仍然用上述的方法来解它们,计算的结果如何?计算的结果说明了什么?

(3)讨论病态问题求解的算法

由于本实验用到迭代法,故先给出迭代法相关资料:

迭代法也称辗转法,是一种不断用变量的旧值递推新值的过程,跟迭代法相对应的是直接法(或者称为一次解法),即一次性解决问题。迭代法又分为精确迭代和近似迭代。“二分法”和“牛顿迭代法”属于近似迭代法。迭代算法是用计算机解决问题的一种基本方法。它利用计算机运算速度快、适合做重复性操作的特点,让计算机对一组指令(或一定步骤)进行重复执行,在每次执行这组指令(或这些步骤)时,都从变量的原值推出它的一个新值。 迭代是数值分析中通过从一个初始估计出发寻找一系列近似解来解决问题(一般是解方程或者方程组)的过程,为实现这一过程所使用的方法统称为迭代法(Iterative Method )。一般可以做如下定义:对于给定的线性方程组f Bx x +=(这里的f B x ,,同为矩阵,任意线性方程组都可以变换成此形式),用公式f k Bx k x +=+)()1((括号中为上标,代表迭代k 次得到的x ,初始时k=0)逐步带入求近似解的方法称为迭代法(或称一阶定常迭代法)。如果k 趋向无穷大时)(lim k x 存在,记为*x ,称此迭代法收敛。显然*x 就是此方程组的解,否则称为迭代法发散。跟迭代法相对应的是直接法(或者称为一次解法),即一次性的快速解决问题,例如通过开方解决方程43=+x 。一般如果可能,直接解法总是优先考虑的。但当遇到复杂问题时,特别是在未知量很多,方程为非线性时,我们无法找到直接解法(例如五次以及更高次的代数方程没有解析解,参见阿贝尔定理),这时候或许可以通过迭代法寻求方程(组)的近似解。最常见的迭代法是牛顿法。其他还包括最速下降法、共轭迭代法、变尺度迭代法、最小二乘法、线性规划、非线性规划、单纯型法、惩罚函数法、斜率投影法、遗传算法、模拟退火法等等。利用迭代算法解决问题,需要做好以下三个方面的工作: (1)确定迭代变量

在可以用迭代算法解决的问题中,至少存在一个直接或间接地不断由旧值递推出新值的变量,这个变量就是迭代变量。 (2)建立迭代关系式

所谓迭代关系式,指如何从变量的前一个值推出其下一个值的公式(或关系)。迭代关系式的建立是解决迭代问题的关键,通常可以顺推或倒推的方法来完成。 (3)对迭代过程进行控制

在什么时候结束迭代过程?这是编写迭代程序必须考虑的问题。不能让迭代过程无休止地重复执行下去。迭代过程的控制通常可分为两种情况:一种是所需的迭代次数是个确定的值,可以计算出来;另一种是所需的迭代次数无法确定。对于前一种情况,可以构建一个固

定次数的循环来实现对迭代过程的控制;对于后一种情况,需要进一步分析出用来结束迭代过程的条件。

实验过程:

6.1.1.1 实验要求一的gauss实现程序如下:

Gauss消去法程序:

function x=gaussong(n,r)

A=hilb(n)

b=A*ones(n,1)

for i=1:4

p=input('条件数对应的范数是p-范数:p=')

pp=cond(A,p)

end

pause

[m,n]=size(A);

nb=n+1;Ab=[A b]

r=input('请输入是否为手动,手动输入1,自动输入0:r=')

for i=1:n-1

if r==0

[pivot,p]=max(abs(Ab(i:n,i)));

ip=p+i-1;

if ip~=i

Ab([i ip],:)=Ab([ip i],:);disp(Ab); pause

end

end

if r==1

i=i

ip=input('输入i列所选元素所处的行数:ip=');

Ab([i ip],:)=Ab([ip i],:);disp(Ab); pause

end

pivot=Ab(i,i);

for k=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);

for i=n-1:-1:1

x(i)=(Ab(i,nb)-Ab(i,i+1:n)*x(i+1:n))/Ab(i,i);

end

6.1.1.2实验要求一的gauss实现程序运行结果如下:

6.1.2.1实验要求一的J程序如下:

function n=JD(n)

A=hilb(n);

for i=1:n;

a0(i)=1;

x(i)=0;

end

b=A*a0';

for i=1:100;

y=x;

for j=1:n;

x(j)=b(j)/A(j,j);

for k=1:j-1;

x(j)=x(j)-A(j,k)*y(k)/A(j,j);

end

for k=j+1:n;

x(j)=x(j)-A(j,k)*y(k)/A(j,j);

end

end

end

x

6.1.2.2实验要求一的J程序运行结果如下:6.1.3.1实验要求一的GS程序如下:function n=GS(n)

A2=hilb(n);

for i=1:n;

a02(i)=1;

x2(i)=0;

end

b2=A2*a02';

for i=1:100000;

for j=1:n;

x2(j)=b2(j)/A2(j,j);

for k=1:j-1;

x2(j)=x2(j)-A2(j,k)*x2(k)/A2(j,j);

end

for k=j+1:n;

x2(j)=x2(j)-A2(j,k)*x2(k)/A2(j,j);

end

end

end

x2

6.1.3.2实验要求一的GS程序运行结果如下:6.1.4.1实验要求一的SOR程序如下:function [n ss]=SOR(n,ss)

A3=hilb(n);

for i=1:n;

a03(i)=1; x3(i)=0; end

b3=A3*a03'; for i=1:100000; for j=1:n; rc=x3(j);

x3(j)=b3(j)/A3(j,j); for k=1:j-1;

x3(j)=x3(j)-A3(j,k)*x3(k)/A3(j,j); end for k=j+1:n;

x3(j)=x3(j)-A3(j,k)*x3(k)/A3(j,j); end

x3(j)=(1-ss)*rc+ss*x3(j); end end x3

6.1.4.2实验要求一的SOR 程序运行结果如下:

(注意由于运行结果迭代过程很长故不体现在实验报告中)

6.2.1实验要求二的gauss 程序运行结果如下:

x = 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.9975 6.2.2实验要求二的J 程序运行结果如下: 6.2.3实验要求二的GS 程序运行结果如下: 6.2.4实验要求二的SOR 程序运行结果如下:

实验要求二的结果分析: 选择问题的维数为20时:

1. 用Gauss 消去法求得的解与精确解相差很大,说明能否得到优秀的解取决于算法的稳定性,如果算法不够稳定产生的结果将变的无法理喻。

2. 取初始向量为0,用J 迭代方法迭代发散,无法求解;

3. 取初始向量为0,用GS 迭代方法迭代不发散,能求得解,但收敛非常缓慢,而且迭代次数越多,与准确解的偏差就越大,说明GS 迭代适合迭代次数少的,但是通常我们无法得知需要迭代的次数。

4. 取初始向量为0,用SOR 迭代方法迭代不发散,能求得解,但同样收敛非常缓慢。 总之,从上面的结果可以看出当病态问题的阶数升高时作为直接法的Gauss 消去法又能求解变成不能求解。而GS 和SOR 迭代法在阶数升高时仍能求解。但在阶数较低时直接法能求得精确解而迭代发却总存在一定的误差。可见直接法与迭代法各有各的优势与不足。 关于病态问题的求解,主要的方法是对原方程作某些预处理,以降低系数矩阵的条件数。可以采取将系数矩阵A 的每一行本别乘上适当常数的方法。即找到可逆的对角阵1D 和2D 使

方程组Ax b =化为 1

1212,D AD y D b x D y -== 理论上最好选择对角阵12,D D 满足:

1212()min ()cond D AD cond D AD =。

补充:

本实验用到了Gauss 消去法、J 迭代法、GS 迭代法、SOR 迭代,故对其重新理解学习了一下:

Gauss 消去法:

1 将其增广矩阵化为行阶梯形

2 若最后有形如[0 0 ... 0 a] (a<>0)的行则无解

3 若含有自由变量则有无穷组解

4 原方程有唯一解。采用回代求解。

至于有无穷组解的方程组的求解,需将其化为行最简形矩阵,其方法称为高斯-若尔当消元法。

J 迭代法: 一、算法理论

Jacobi 迭代格式的引出是依据迭代法的基本思想:

构造一个向量系列}{)(n x ,使其收敛至某个极限*X ,则*X 就是要求的方程组的准确解。

Jacobi 迭代

将方程组:

在假设0≠ii a ,改写成 如果引用系数矩阵

??????????=nn n n a a a a A 1111,??????????=0011 n n b b B 及向量??

??

?

?????=??????????=??????????=n n n g g g b b b x x X 111,, 方程组(1)和(2)分别可写为:b AX =及g BX X +=,这样就得到了Jacobi 迭代格式01g BX X k k +=+用Jacobi 迭代解方程组b AX =时,就可任意取初值0X 带入迭代可知式g BX X k k +=+1,然后求k k X ∞

→lim 。但是,n 比较大的时候,写方程

组(1)和(2)是很麻烦的,如果直接由A ,b 能直接得到B ,g 就是矩阵与向量的运算了,那么如何得到B ,g 呢

实际上,如果引进非奇异对角矩阵0)(≠ii a

将A 分解成:D D A A +-=要求b AX =的解,实质上就有DX X D A AX --=)(,

而D 是非奇异的,所以1-D 存在,X A D AX DX )(-+=,从而有b

D AX D X 11--+=

我们在这里不妨令A D I B 1--=,b D g 1-=就得到Jacobi 迭代格式:

g BX X k k +=+1现在考虑Jacobi 迭代法的计算程序

float a[3][3]={{10,-2,-1},{-2,10,-1},{-1,-2,-5}};float b[3]={3,15,10};分别代表X 的系数和等号右边的常数项, 即

先输入方程,运行main 函数,如果first 不为null ,则执行if 括号里的,否则执行else 里面的,最后会调用方法sum ()。在sum ()中sum=a[m][n]*x[n]+sum;y=(b[m]-sum)/a[m][m],之后运行for 循环,最后输出结果,算法结束。 二、算法框图

GS 迭代法:

1.高斯 - 塞德尔迭代法 公式的矩阵形式

首先将高斯 - 塞德尔迭代法的公式表示为矩阵形式,为此设

这里 D 是系数矩阵A 的对角部分,L -是严格下三角部分,U -是严格上三角部分,则高斯 - 塞德尔迭代法的公式可表示为 用矩阵 D 乘等式两边得

再用矩阵 1

)(--L D 乘等式两边得

其中矩阵 U L D B 1

1)(--=称为高斯—塞德尔迭代矩阵。

由此可见,高斯 - 塞德尔迭代法是一般迭代法中迭代矩阵为 U L D B 1

1)(--=的特殊情形。需要指出的是,由于矩阵 U L D 1

)(--难于计算,所以式(2)多用在理论分析中。 2.高斯—塞德尔迭代法计算框图(见图)

开始

读入[][][],b a

输出 结束

调用()

sum 方法

高斯—塞德尔迭代法计算框图

松弛法:

ii

k i

k i

n

i j k j ij i j i ij i ii k i

a r x x a k x a

b a x )

1()

(1

)(11)

1(])1([1++=-=++

=-+-=∑∑其中∑∑≥<++--=i

j k j

ij i

j k j

ij i k i

x a x a b r )

()

1()

1(

相当于在的基础上)

(k i x 加个余项生成)

1(+k i

x ,下面令ii

k i

k i

k i

a r x x )

1()

()

1(+++=ω

,希望

通过选取ω来加速收敛,这就是松弛法。

实验七 非线性方程求根

实验7.1(迭代法、初始值与收敛性)

实验目的:初步认识非线性问题的迭代法与线性问题迭代法的差别,探讨迭代法及初始值与迭代收敛性的关系。

问题提出:迭代法是求解非线性方程的基本思想方法,与线性方程的情况一样,其构造方法可以有多种多样,但关键是怎样才能使迭代收敛且有较快的收敛速度。

实验内容:考虑一个简单的代数方程 针对上述方程,可以构造多种迭代法,如

在实轴上取初始值x 0,请分别用迭代(7.1)-(7.3)作实验,记录各算法的迭代过程。

实验要求:

(1)取定某个初始值,分别计算(7.1)-(7.3)迭代结果,它们的收敛性如何?重复选取不同的初始值,反复实验。请自选设计一种比较形象的记录方式(如利用Matlab 的图形功能),分析三种迭代法的收敛性与初值选取的关系。

(2)对三个迭代法中的某个,取不同的初始值进行迭代,结果如何?试分析迭代法对不同的初值是否有差异?

(3)线性方程组迭代法的收敛性是不依赖初始值选取的。比较线性与非线性问题迭代的差异,有何结论和问题。

实验过程:

7.1程序:

7.1.1 对于第一个迭代方程的程序: 保存为:diedai71 clear clc

a=-1.5;b=2.5; y00=0; x00=input('请输入第一个函数的初值:x00='); x=linspace(a,b,80);

y0=x; %计算直线y=x

y1=diedai7f1(x); %计算迭代函数y=f(x) clear y ; y=[y0;y1];

plot(x,y,'linewidth',1)

legend('y=x','y=f1')

title('x(n+1)=[x(n)]^2-1') %输出标题

hold on

plot([a b],[0,0],'k-',[0 0],[a b],'k-')

axis([a,b,a,b]) %画坐标轴

z=[];

for i=1:15 %画蛛网图,迭代过程为n=15次

xt(1)=x00;yt(1)=y00; %决定始点坐标

xt(2)=diedai7f1(xt(1)); %决定终点坐标

yt(2)=diedai7f1(xt(1));

diedaiplot71(xt,yt,0.6); %画蛛网图

if i<=5

pause %按任意键逐次观察前5次迭代的蛛网图end

x00=xt(2);y00=yt(2); %将本次迭代的终点作为下次的始点

z=[z,xt(1)]; %保存迭代点

end

保存为:diedai7f1

function y=diedai7f1(x)

y=(x.*x-1);

保存为:diedaiplot72

function out=diedaiplot72(x,y,p)%画一次迭代的蛛网图,改变p调节箭头的大小u(1)=0;v(1)=(y(2)-y(1)); %画出始点(x(1),y(1))终点(x(2),y(2))的有向折线段

u(2)=eps;v(2)=eps;

X=[x(1) x(1)];

Y=[y(1) y(2)];

h=quiver(X,Y,u,v,p);

set(h,'color','red');

hold on

u(1)=(x(2)-x(1));v(1)=0;

u(2)=eps;v(2)=eps;

h=quiver([x(1) x(2)],[y(2) y(2)],u,v,p);

set(h,'color','red');

plot([x(1) x(1) x(2)],[y(1) y(2) y(2)],'r.-')

7.1.1 对于第一个迭代方程的程序运行结果:

(注:由于可以估计出方程的根,故选取根附近的值开始迭代)

当x=1.5时程序运行如下图:

当x=-1时程序运行如下图:

当x=0.5时程序运行如下图:

当x=-0.5时程序运行如下图:

当x=0时程序运行如下图:

7.2.1 对于第二个迭代方程的程序:

保存为diedai72:

clear

clc

a=0.1;b=6.5; y00=0; x00=input('请输入第二个函数的初值:x00=');

x=linspace(a,b,80);

y0=x; %计算直线y=x

y1=diedai7f2(x); %计算迭代函数y=f(x)

clear y;

y=[y0;y1];

plot(x,y,'linewidth',1)

legend('y=x','y=f1')

title('x(n+1)=[x(n)]^2-1') %输出标题

hold on

plot([a b],[0,0],'k-',[0 0],[a b],'k-')

axis([a,b,a,b]) %画坐标轴

z=[];

for i=1:15 %画蛛网图,迭代过程为n=15次xt(1)=x00;yt(1)=y00; %决定始点坐标

xt(2)=diedai7f2(xt(1)); %决定终点坐标

yt(2)=diedai7f2(xt(1));

diedaiplot72(xt,yt,0.6); %画蛛网图

if i<=5

pause %按任意键逐次观察前5次迭代的蛛网图end

x00=xt(2);y00=yt(2); %将本次迭代的终点作为下次的始点

z=[z,xt(1)]; %保存迭代点

end

保存为迭代7f2:

function y=diedai7f2(x)

y=(1+1./x);

保存为:diedaiplot72

function out=diedaiplot72(x,y,p)%画一次迭代的蛛网图,改变p调节箭头的大小u(1)=0;v(1)=(y(2)-y(1)); %画出始点(x(1),y(1))终点(x(2),y(2))的有向折线段

u(2)=eps;v(2)=eps;

X=[x(1) x(1)];

Y=[y(1) y(2)];

h=quiver(X,Y,u,v,p);

set(h,'color','red');

hold on

u(1)=(x(2)-x(1));v(1)=0;

u(2)=eps;v(2)=eps;

h=quiver([x(1) x(2)],[y(2) y(2)],u,v,p);

set(h,'color','red');

plot([x(1) x(1) x(2)],[y(1) y(2) y(2)],'r.-')

7.2.2 对于第二个迭代方程的程序运行结果:

(注:由于可以估计出方程的根,故选取根附近的值开始迭代)

当x=1.5时程序运行如下图:

当x=5时程序运行如下图:

当x=2.5时程序运行如下图:

当x=0.5时程序运行如下图:

7.3.1 对于第三个迭代方程的程序:

保存为diedai73:

clear

clc

a=0;b=2; y00=0; x00=input('请输入第三个函数的初值:x00=');

x=linspace(a,b,80);

y0=x; %计算直线y=x

y1=diedai7f3(x); %计算迭代函数y=f(x)

clear y;

y=[y0;y1];

plot(x,y,'linewidth',1)

legend('y=x','y=f1')

title('x(n+1)=[x(n)]^2-1') %输出标题

hold on

plot([a b],[0,0],'k-',[0 0],[a b],'k-')

axis([a,b,a,b]) %画坐标轴

z=[];

for i=1:15 %画蛛网图,迭代过程为n=15次xt(1)=x00;yt(1)=y00; %决定始点坐标

xt(2)=diedai7f3(xt(1)); %决定终点坐标

yt(2)=diedai7f3(xt(1));

diedaiplot73(xt,yt,0.6); %画蛛网图

if i<=5

pause %按任意键逐次观察前5次迭代的蛛网图end

x00=xt(2);y00=yt(2); %将本次迭代的终点作为下次的始点

z=[z,xt(1)]; %保存迭代点

end

保存为diedai7f3:

function y=diedai7f3(x)

y=sqrt(x+1);

保存为diedaiplot73:

function out=diedaiplot73(x,y,p)%画一次迭代的蛛网图,改变p调节箭头的大小u(1)=0;v(1)=(y(2)-y(1)); %画出始点(x(1),y(1))终点(x(2),y(2))的有向折线段

u(2)=eps;v(2)=eps;

X=[x(1) x(1)];

Y=[y(1) y(2)];

h=quiver(X,Y,u,v,p);

set(h,'color','red');

hold on

u(1)=(x(2)-x(1));v(1)=0;

u(2)=eps;v(2)=eps;

h=quiver([x(1) x(2)],[y(2) y(2)],u,v,p);

set(h,'color','red');

plot([x(1) x(1) x(2)],[y(1) y(2) y(2)],'r.-')

7.3.2 对于第三个迭代方程的程序运行结果:

(注:由于可以估计出方程的根,故选取根附近的值开始迭代)

当x=1.5时程序运行如下图:

当x=0.4时程序运行如下图:

当x=0.8时程序运行如下图:

当x=1.2时程序运行如下图:

当x=1.8时程序运行如下图:

实验要求一分析:

通过上述三种迭代图像可以判断出第一个迭代方程不收敛,而第二个、第三个迭代方程收敛,同时可以看出初值的选择决定收敛速度,对于用同一个迭代方程而言,第二个迭代方程选取初始点和方程在该点的y轴差值越小,收敛速度越快,对于第三个迭代方程,大体上和第二个迭代方程相似,选取的初始点使得两条曲线间的间距越小,收敛越快

实验要求二分析:

使用迭代法求解非线性方程的根,可有很多种不同的迭代方程,对于不同的迭代方程求解可能会产生不同的结果,就像第一个迭代方程不收敛,而第二、第三个却收敛,这说明迭代方程的选取将会直接影响我们对非线性方程的求解,如果选取的迭代方程好的话,那么它至少是收敛的,然后我们自然希望收敛的速度越快越好,也就是说迭代方程的选取的好坏有两个评价标准:(1)是否收敛;(2)收敛速度的快慢:收敛速度的快慢还取决于初始值的选择,如果初始值选取的较好的话,能够很快就收敛结束。

根据上述几个原则我对第二个迭代方程进行讨论:

(1)是否收敛:由迭代图像可以很明显的看出第二个迭代方程是收敛的,因为无论初始化值选取多少总会收敛到一定的区域内,不像第一个迭代方程初始值无论选取多大都不可能收敛到一定的区域,所以它比第一个要好;

(2)收敛速度的快慢:从图像和初始值的选取可以看出,收敛速度明显取决于初始值的选择,初始值选取的越贴近解则收敛速度越快,否则就只有慢慢收敛,第三个迭代方程也是如此;

实验要求三分析:

数值分析实验报告176453

实验报告 插值法 数学实验室 数值逼近 算法设计 级 ____________________________ 号 ____________________________ 名 _____________________________ 实验项目名称 实验室 所属课程名称 实验类型 实验日期

实验概述: 【实验目的及要求】 本次实验的目的是熟练《数值分析》第二章“插值法”的相关内容,掌握三种插 多项式插值,三次样条插值,拉格朗日插值,并比较三种插值方法的 优劣。 本次试验要求编写牛顿多项式插值,三次样条插值,拉格朗日插值的程序编码,并 去实现。 【实验原理】 《数值分析》第二章“插值法”的相关内容,包括:牛顿多项式插值,三次样条插值, 拉格朗日 插值的相应算法和相关性质。 【实验环境】(使用的软硬件) 软件: MATLAB 2012a 硬件: 电脑型号:联想 Lenovo 昭阳E46A 笔记本电脑 操作系统: Win dows 8专业版 处理器:In tel ( R Core ( TM i3 CPU M 350 @2.27GHz 2.27GHz 实验内容: 【实验方案设计】 第一步,将书上关于三种插值方法的内容转化成程序语言,用 MATLA B 现; 第二步,分别用牛顿多项式插值,三次样条插值,拉格朗日插值求解不同的问题。 【实验过程】(实验步骤、记录、数据、分析) 实验的主要步骤是:首先分析问题,根据分析设计 MATLA 程序,利用程序算出 问题答案,分析所得答案结果,再得出最后结论。 实验一: 已知函数在下列各点的值为 试用4次牛顿插值多项式 P 4( x )及三次样条函数 S ( x )(自然边界条件)对数据进行插值。 用图给出{( X i , y i ), X i =0.2+0.08i , i=0 , 1, 11, 10 } , P 4 ( x )及 S ( x )。 值方法:牛顿 在MATLAB 件中

数值分析实验报告模板

数值分析实验报告模板 篇一:数值分析实验报告(一)(完整) 数值分析实验报告 1 2 3 4 5 篇二:数值分析实验报告 实验报告一 题目:非线性方程求解 摘要:非线性方程的解析解通常很难给出,因此线性方程的数值解法就尤为重要。本实验采用两种常见的求解方法二分法和Newton法及改进的Newton法。利用二分法求解给定非线性方程的根,在给定的范围内,假设f(x,y)在[a,b]上连续,f(a)xf(b) 直接影响迭代的次数甚至迭代的收敛与发散。即若x0 偏离所求根较远,Newton法可能发散的结论。并且本实验中还利用利用改进的Newton法求解同样的方程,且将结果与Newton法的结果比较分析。 前言:(目的和意义) 掌握二分法与Newton法的基本原理和应用。掌握二分法的原理,验证二分法,在选对有根区间的前提下,必是收

敛,但精度不够。熟悉Matlab语言编程,学习编程要点。体会Newton使用时的优点,和局部收敛性,而在初值选取不当时,会发散。 数学原理: 对于一个非线性方程的数值解法很多。在此介绍两种最常见的方法:二分法和Newton法。 对于二分法,其数学实质就是说对于给定的待求解的方程f(x),其在[a,b]上连续,f(a)f(b) Newton法通常预先要给出一个猜测初值x0,然后根据其迭代公式xk?1?xk?f(xk) f'(xk) 产生逼近解x*的迭代数列{xk},这就是Newton法的思想。当x0接近x*时收敛很快,但是当x0选择不好时,可能会发散,因此初值的选取很重要。另外,若将该迭代公式改进为 xk?1?xk?rf(xk) 'f(xk) 其中r为要求的方程的根的重数,这就是改进的Newton 法,当求解已知重数的方程的根时,在同种条件下其收敛速度要比Newton法快的多。 程序设计: 本实验采用Matlab的M文件编写。其中待求解的方程写成function的方式,如下 function y=f(x);

数值分析实验报告1

实验一误差分析 实验1.1(病态问题) 实验目的:算法有“优”与“劣”之分,问题也有“好”与“坏”之别。对数值方法的研究而言,所谓坏问题就是问题本身对扰动敏感者,反之属于好问题。通过本实验可获得一个初步体会。 数值分析的大部分研究课题中,如线性代数方程组、矩阵特征值问题、非线性方程及方程组等都存在病态的问题。病态问题要通过研究和构造特殊的算法来解决,当然一般要付出一些代价(如耗用更多的机器时间、占用更多的存储空间等)。 问题提出:考虑一个高次的代数多项式 显然该多项式的全部根为1,2,…,20共计20个,且每个根都是单重的。现考虑该多项式的一个扰动 其中ε(1.1)和(1.221,,,a a 的输出b ”和“poly ε。 (1(2 (3)写成展 关于α 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) n

很容易的得出对一个多次的代数多项式的其中某一项进行很小的扰动,对其多项式的根会有一定的扰动的,所以对于这类病态问题可以借助于MATLAB来进行问题的分析。 学号:06450210 姓名:万轩 实验二插值法

数值分析实验报告

实验一 误差分析 实验(病态问题) 实验目的:算法有“优”与“劣”之分,问题也有“好”与“坏”之别。对数值方法的研究而言,所谓坏问题就是问题本身对扰动敏感者,反之属于好问题。通过本实验可获得一个初步体会。 数值分析的大部分研究课题中,如线性代数方程组、矩阵特征值问题、非线性方程及方程组等都存在病态的问题。病态问题要通过研究和构造特殊的算法来解决,当然一般要付出一些代价(如耗用更多的机器时间、占用更多的存储空间等)。 问题提出:考虑一个高次的代数多项式 )1.1() ()20()2)(1()(20 1∏=-=---=k k x x x x x p Λ 显然该多项式的全部根为1,2,…,20共计20个,且每个根都是单重的。现考虑该多项式的一个扰动 )2.1(0 )(19=+x x p ε 其中ε是一个非常小的数。这相当于是对()中19x 的系数作一个小的扰动。我们希望比较()和()根的差别,从而分析方程()的解对扰动的敏感性。 实验内容:为了实现方便,我们先介绍两个Matlab 函数:“roots ”和“poly ”。 roots(a)u = 其中若变量a 存储n+1维的向量,则该函数的输出u 为一个n 维的向量。设a 的元素依次为121,,,+n a a a Λ,则输出u 的各分量是多项式方程 01121=+++++-n n n n a x a x a x a Λ 的全部根;而函数 poly(v)b = 的输出b 是一个n+1维变量,它是以n 维变量v 的各分量为根的多项式的系数。可见“roots ”和“poly ”是两个互逆的运算函数。 ;000000001.0=ess );21,1(zeros ve = ;)2(ess ve = ))20:1((ve poly roots +

数值计算实验报告

(此文档为word格式,下载后您可任意编辑修改!) 2012级6班###(学号)计算机数值方法 实验报告成绩册 姓名:宋元台 学号: 成绩:

数值计算方法与算法实验报告 学期: 2014 至 2015 第 1 学期 2014年 12月1日课程名称: 数值计算方法与算法专业:信息与计算科学班级 12级5班 实验编号: 1实验项目Neton插值多项式指导教师:孙峪怀 姓名:宋元台学号:实验成绩: 一、实验目的及要求 实验目的: 掌握Newton插值多项式的算法,理解Newton插值多项式构造过程中基函数的继承特点,掌握差商表的计算特点。 实验要求: 1. 给出Newton插值算法 2. 用C语言实现算法 二、实验内容 三、实验步骤(该部分不够填写.请填写附页)

1.算法分析: 下面用伪码描述Newton插值多项式的算法: Step1 输入插值节点数n,插值点序列{x(i),f(i)},i=1,2,……,n,要计算的插值点x. Step2 形成差商表 for i=0 to n for j=n to i f(j)=((f(j)-f(j-1)(x(j)-x(j-1-i)); Step3 置初始值temp=1,newton=f(0) Step4 for i=1 to n temp=(x-x(i-1))*temp*由temp(k)=(x-x(k-1))*temp(k-1)形成 (x-x(0).....(x-x(i-1)* Newton=newton+temp*f(i); Step5 输出f(x)的近似数值newton(x)=newton. 2.用C语言实现算法的程序代码 #includeMAX_N) { printf("the input n is larger than MAX_N,please redefine the MAX_N.\n"); return 1; } if(n<=0) { printf("please input a number between 1 and %d.\n",MAX_N); return 1; } printf("now input the (x_i,y_i)i=0,...%d\n",n); for(i=0;i<=n;i++) { printf("please input x(%d) y(%d)\n",i,i);

数值分析实验报告

实验一、误差分析 一、实验目的 1.通过上机编程,复习巩固以前所学程序设计语言及上机操作指令; 2.通过上机计算,了解误差、绝对误差、误差界、相对误差界的有关概念; 3.通过上机计算,了解舍入误差所引起的数值不稳定性。 二.实验原理 误差问题是数值分析的基础,又是数值分析中一个困难的课题。在实际计算中,如果选用了不同的算法,由于舍入误差的影响,将会得到截然不同的结果。因此,选取算法时注重分析舍入误差的影响,在实际计算中是十分重要的。同时,由于在数值求解过程中用有限的过程代替无限的过程会产生截断误差,因此算法的好坏会影响到数值结果的精度。 三.实验内容 对20,,2,1,0 =n ,计算定积分 ?+=10 5dx x x y n n . 算法1:利用递推公式 151--=n n y n y , 20,,2,1 =n , 取 ?≈-=+=1 00182322.05ln 6ln 51dx x y . 算法2:利用递推公式 n n y n y 51511-= - 1,,19,20 =n . 注意到 ???=≤+≤=10 10202010201051515611261dx x dx x x dx x , 取 008730.0)12611051(20120≈+≈y .: 四.实验程序及运行结果 程序一: t=log(6)-log(5);

n=1; y(1)=t; for k=2:1:20 y(k)=1/k-5*y(k-1); n=n+1; end y y =0.0884 y =0.0581 y =0.0431 y =0.0346 y =0.0271 y =0.0313 y =-0.0134 y =0.1920 y =-0.8487 y =4.3436 y =-21.6268 y =108.2176 y =-541.0110 y =2.7051e+003 y =-1.3526e+004 y =6.7628e+004 y =-3.3814e+005 y =1.6907e+006 y =-8.4535e+006 y =4.2267e+007 程序2: y=zeros(20,1); n=1; y1=(1/105+1/126)/2;y(20)=y1; for k=20:-1:2 y(k-1)=1/(5*k)-(1/5)*y(k); n=n+1; end 运行结果:y = 0.0884 0.0580 0.0431 0.0343 0.0285 0.0212 0.0188 0.0169

(完整版)哈工大-数值分析上机实验报告

实验报告一 题目:非线性方程求解 摘要:非线性方程的解析解通常很难给出,因此线性方程的数值解法就尤为重要。本实验采用两种常见的求解方法二分法和Newton法及改进的Newton法。 前言:(目的和意义) 掌握二分法与Newton法的基本原理和应用。 数学原理: 对于一个非线性方程的数值解法很多。在此介绍两种最常见的方法:二分法和Newton法。 对于二分法,其数学实质就是说对于给定的待求解的方程f(x),其在[a,b]上连续,f(a)f(b)<0,且f(x)在[a,b]内仅有一个实根x*,取区间中点c,若,则c恰为其根,否则根据f(a)f(c)<0是否成立判断根在区间[a,c]和[c,b]中的哪一个,从而得出新区间,仍称为[a,b]。重复运行计算,直至满足精度为止。这就是二分法的计算思想。

Newton法通常预先要给出一个猜测初值x0,然后根据其迭代公式 产生逼近解x*的迭代数列{x k},这就是Newton法的思想。当x0接近x*时收敛很快,但是当x0选择不好时,可能会发散,因此初值的选取很重要。另外,若将该迭代公式改进为 其中r为要求的方程的根的重数,这就是改进的Newton法,当求解已知重数的方程的根时,在同种条件下其收敛速度要比Newton法快的多。 程序设计: 本实验采用Matlab的M文件编写。其中待求解的方程写成function的方式,如下 function y=f(x); y=-x*x-sin(x); 写成如上形式即可,下面给出主程序。 二分法源程序: clear %%%给定求解区间 b=1.5; a=0;

%%%误差 R=1; k=0;%迭代次数初值 while (R>5e-6) ; c=(a+b)/2; if f12(a)*f12(c)>0; a=c; else b=c; end R=b-a;%求出误差 k=k+1; end x=c%给出解 Newton法及改进的Newton法源程序:clear %%%% 输入函数 f=input('请输入需要求解函数>>','s') %%%求解f(x)的导数 df=diff(f);

数值分析实验报告总结

数值分析实验报告总结 随着电子计算机的普及与发展,科学计算已成为现代科 学的重要组成部分,因而数值计算方法的内容也愈来愈广泛和丰富。通过本学期的学习,主要掌握了一些数值方法的基本原理、具体算法,并通过编程在计算机上来实现这些算法。 算法算法是指由基本算术运算及运算顺序的规定构成的完 整的解题步骤。算法可以使用框图、算法语言、数学语言、自然语言来进行描述。具有的特征:正确性、有穷性、适用范围广、运算工作量少、使用资源少、逻辑结构简单、便于实现、计算结果可靠。 误差 计算机的计算结果通常是近似的,因此算法必有误差, 并且应能估计误差。误差是指近似值与真正值之差。绝对误差是指近似值与真正值之差或差的绝对值;相对误差:是指近似值与真正值之比或比的绝对值。误差来源见表 第三章泛函分析泛函分析概要 泛函分析是研究“函数的函数”、函数空间和它们之间 变换的一门较新的数学分支,隶属分析数学。它以各种学科

如果 a 是相容范数,且任何满足 为具体背景,在集合的基础上,把客观世界中的研究对象抽 范数 范数,是具有“长度”概念的函数。在线性代数、泛函 分析及相关的数学领域,泛函是一个函数,其为矢量空间内 的所有矢量赋予非零的正长度或大小。这里以 Cn 空间为例, Rn 空间类似。最常用的范数就是 P-范数。那么 当P 取1, 2 ,s 的时候分别是以下几种最简单的情形: 其中2-范数就是通常意义下的距离。 对于这些范数有以下不等式: 1 < n1/2 另外,若p 和q 是赫德尔共轭指标,即 1/p+1/q=1 么有赫德尔不等式: II = ||xH*y| 当p=q=2时就是柯西-许瓦兹不等式 般来讲矩阵范数除了正定性,齐次性和三角不等式之 矩阵范数通常也称为相容范数。 象为元素和空间。女口:距离空间,赋范线性空间, 内积空间。 1-范数: 1= x1 + x2 +?+ xn 2-范数: x 2=1/2 8 -范数: 8 =max oo ,那 外,还规定其必须满足相容性: 所以

数值分析2016上机实验报告

序言 数值分析是计算数学的范畴,有时也称它为计算数学、计算方法、数值方法等,其研究对象是各种数学问题的数值方法的设计、分析及其有关的数学理论和具体实现的一门学科,它是一个数学分支。是科学与工程计算(科学计算)的理论支持。许多科学与工程实际问题(核武器的研制、导弹的发射、气象预报)的解决都离不开科学计算。目前,试验、理论、计算已成为人类进行科学活动的三大方法。 数值分析是计算数学的一个主要部分,计算数学是数学科学的一个分支,它研究用计算机求解各种数学问题的数值计算方法及其理论与软件实现。现在面向数值分析问题的计算机软件有:C,C++,MATLAB,Python,Fortran等。 MATLAB是matrix laboratory的英文缩写,它是由美国Mathwork公司于1967年推出的适合用于不同规格计算机和各种操纵系统的数学软件包,现已发展成为一种功能强大的计算机语言,特别适合用于科学和工程计算。目前,MATLAB应用非常广泛,主要用于算法开发、数据可视化、数值计算和数据分析等,除具备卓越的数值计算能力外,它还提供了专业水平的符号计算,文字处理,可视化建模仿真和实时控制等功能。 本实验报告使用了MATLAB软件。对不动点迭代,函数逼近(lagrange插值,三次样条插值,最小二乘拟合),追赶法求解矩阵的解,4RungeKutta方法求解,欧拉法及改进欧拉法等算法做了简单的计算模拟实践。并比较了各种算法的优劣性,得到了对数值分析这们学科良好的理解,对以后的科研数值分析能力有了极大的提高。

目录 序言 (1) 问题一非线性方程数值解法 (3) 1.1 计算题目 (3) 1.2 迭代法分析 (3) 1.3计算结果分析及结论 (4) 问题二追赶法解三对角矩阵 (5) 2.1 问题 (5) 2.2 问题分析(追赶法) (6) 2.3 计算结果 (7) 问题三函数拟合 (7) 3.1 计算题目 (7) 3.2 题目分析 (7) 3.3 结果比较 (12) 问题四欧拉法解微分方程 (14) 4.1 计算题目 (14) 4.2.1 方程的准确解 (14) 4.2.2 Euler方法求解 (14) 4.2.3改进欧拉方法 (16) 问题五四阶龙格-库塔计算常微分方程初值问题 (17) 5.1 计算题目 (17) 5.2 四阶龙格-库塔方法分析 (18) 5.3 程序流程图 (18) 5.4 标准四阶Runge-Kutta法Matlab实现 (19) 5.5 计算结果及比较 (20) 问题六舍入误差观察 (22) 6.1 计算题目 (22) 6.2 计算结果 (22) 6.3 结论 (23) 7 总结 (24) 附录

数值分析实验报告

学生实验报告实验课程名称 开课实验室 学院年级专业班 学生姓名学号 开课时间至学年学期

if(A(m,k)~=0) if(m~=k) A([k m],:)=A([m k],:); %换行 end A(k+1:n, k:c)=A(k+1:n, k:c)-(A(k+1:n,k)/ A(k,k))*A(k, k:c); %消去end end x=zeros(length(b),1); %回代求解 x(n)=A(n,c)/A(n,n); for k=n-1:-1:1 x(k)=(A(k,c)-A(k,k+1:n)*x(k+1:n))/A(k,k); end y=x; format short;%设置为默认格式显示,显示5位 (2)建立MATLAB界面 利用MA TLAB的GUI建立如下界面求解线性方程组: 详见程序。 五、计算实例、数据、结果、分析 下面我们对以上的结果进行测试,求解:

? ? ? ? ? ? ? ? ? ? ? ? - = ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? - - - - - - 7 2 5 10 13 9 14 4 4 3 2 1 13 12 4 3 3 10 2 4 3 2 1 x x x x 输入数据后点击和,得到如下结果: 更改以上数据进行测试,求解如下方程组: 1 2 3 4 43211 34321 23431 12341 x x x x ?? ???? ?? ???? ?? ???? = ?? ???? - ?? ???? - ???? ?? 得到如下结果:

数值分析实验报告资料

机电工程学院 机械工程 陈星星 6720150109 《数值分析》课程设计实验报告 实验一 函数插值方法 一、问题提出 对于给定的一元函数)(x f y =的n+1个节点值(),0,1,,j j y f x j n ==。试用Lagrange 公式求其插值多项式或分段二次Lagrange 插值多项式。 数据如下: (1 求五次Lagrange 多项式5L ()x ,计算(0.596)f ,(0.99)f 的值。(提示:结果为(0.596)0.625732f ≈, (0.99) 1.05423f ≈) 实验步骤: 第一步:先在matlab 中定义lagran 的M 文件为拉格朗日函数 代码为: function[c,l]=lagran(x,y) w=length(x); n=w-1; l=zeros(w,w); for k=1:n+1 v=1; for j=1:n+1 if(k~=j) v=conv(v,poly(x(j)))/(x(k)-x(j)); end end l(k,:)=v; end c=y*l; end

第二步:然后在matlab命令窗口输入: >>>> x=[0.4 0.55 0.65 0.80,0.95 1.05];y=[0.41075 0.57815 0.69675 0.90 1.00 1.25382]; >>p = lagran(x,y) 回车得到: P = 121.6264 -422.7503 572.5667 -377.2549 121.9718 -15.0845 由此得出所求拉格朗日多项式为 p(x)=121.6264x5-422.7503x4+572.5667x3-377.2549x2+121.9718x-15.0845 第三步:在编辑窗口输入如下命令: >> x=[0.4 0.55 0.65 0.80,0.95 1.05]; >> y=121.6264*x.^5-422.7503*x.^4+572.5667*x.^3-377.2549*x.^2+121.9718 *x-15.0845; >> plot(x,y) 命令执行后得到如下图所示图形,然后 >> x=0.596; >> y=121.6264*x.^5-422.7503*x.^4+572.5667*x.^3-377.2549*x.^2+121.9718 *x-15.084 y =0.6257 得到f(0.596)=0.6257 同理得到f(0.99)=1.0542

数值分析实验报告册

实验名称:Lagrange插值(实验一) 实验目的: 掌握Lagrange插值数值算法,能够根据给定的函数值表达求出插值多项式和函数在某一点的近似值。实验准备: 1.在开始本实验之前,请回顾教科书的相关内容; 2.需要一台准备安装Windows XP Professional操作系统和装有数学软件的计算机。 实验内容及要求 已知数据如下: 要求: 试用Lagrange插值多项式求0.5626,0.5635,0.5645 x 时的函数近似值. 实验过程: 编写Matlab函数M文件Lagrange如下: function yy=lagrange(x,y,xi) m=length(x); n=length(y); if m~=n,error('向量x与y的长度必须一致');end for k=1:length(xi) s=0; for i=1:m z=1; for j=1:n if j~=i z=z*(xi(k)-x(j))/(x(i)-x(j)); end end s=s+z*y(i); end yy=s end 在命令窗口调用函数M文件lagrange,输出结果如下: >>x=[0.56160, 0.56280, 0.56401, 0.56521]; >>y=[0.82741, 0.82659, 0.82577, 0.82495]; >>xi=[0.5626, 0.5635, 0.5645]; >>yi= lagrange (x,y,xi)

yi= 0.8628 0.8261 0.8254 实验总结(由学生填写): 教师对本次实验的评价(下面的表格由教师填写): 实验名称:曲线拟合的最小二乘方法(实验二) 实验目的: 掌握最小二乘方法,并能根据给定数据求其最小二乘一次或二次多项式,然后进行曲线拟合。实验准备: 1.在开始本实验之前,请回顾教科书的相关内容;

数值分析实验报告1

实验一 误差分析 实验1.1(病态问题) 实验目的:算法有“优”与“劣”之分,问题也有“好”与“坏”之别。对数值方法的研究而言,所谓坏问题就是问题本身对扰动敏感者,反之属于好问题。通过本实验可获得一个初步体会。 数值分析的大部分研究课题中,如线性代数方程组、矩阵特征值问题、非线性方程及方程组等都存在病态的问题。病态问题要通过研究和构造特殊的算法来解决,当然一般要付出一些代价(如耗用更多的机器时间、占用更多的存储空间等)。 问题提出:考虑一个高次的代数多项式 )1.1() ()20()2)(1()(20 1∏=-=---=k k x x x x x p 显然该多项式的全部根为1,2,…,20共计20个,且每个根都是单重的。现考虑该多项式的一个扰动 )2.1(0 )(19=+x x p ε 其中ε是一个非常小的数。这相当于是对(1.1)中19x 的系数作一个小的扰动。我们希望比较(1.1)和(1.2)根的差别,从而分析方程(1.1)的解对扰动的敏感性。 实验容:为了实现方便,我们先介绍两个Matlab 函数:“roots ”和“poly ”。 roots(a)u = 其中若变量a 存储n+1维的向量,则该函数的输出u 为一个n 维的向量。设a 的元素依次为121,,,+n a a a ,则输出u 的各分量是多项式方程 01121=+++++-n n n n a x a x a x a 的全部根;而函数 poly(v)b = 的输出b 是一个n+1维变量,它是以n 维变量v 的各分量为根的多项式的系数。可见“roots ”和“poly ”是两个互逆的运算函数。 ;000000001.0=ess );21,1(zeros ve = ;)2(ess ve =

数值分析实验报告

实验五 解线性方程组的直接方法 实验5.1 (主元的选取与算法的稳定性) 问题提出:Gauss 消去法是我们在线性代数中已经熟悉的。但由于计算机的数值运算是在一个有限的浮点数集合上进行的,如何才能确保Gauss 消去法作为数值算法的稳定性呢?Gauss 消去法从理论算法到数值算法,其关键是主元的选择。主元的选择从数学理论上看起来平凡,它却是数值分析中十分典型的问题。 实验内容:考虑线性方程组 编制一个能自动选取主元,又能手动选取主元的求解线性方程组的Gauss 消去过程。 实验要求: (1)取矩阵?? ? ?? ?? ?????????=????????????????=1415157,6816816816 b A ,则方程有解T x )1,,1,1(* =。取n=10计算矩阵的 条件数。让程序自动选取主元,结果如何? (2)现选择程序中手动选取主元的功能。每步消去过程总选取按模最小或按模尽可能小的元素作为主元,观察并记录计算结果。若每步消去过程总选取按模最大的元素作为主元,结果又如何?分析实验的结果。 (3)取矩阵阶数n=20或者更大,重复上述实验过程,观察记录并分析不同的问题及消去过程中选择不同的主元时计算结果的差异,说明主元素的选取在消去过程中的作用。 (4)选取其他你感兴趣的问题或者随机生成矩阵,计算其条件数。重复上述实验,观察记录并分析实验结果。 思考题一:(Vadermonde 矩阵)设 ?? ??????????????????????=? ? ? ?????????????=∑∑∑∑====n i i n n i i n i i n i i n n n n n n n x x x x b x x x x x x x x x x x x A 0020 10022222121102001111 ,, 其中,n k k x k ,,1,0,1.01 =+=, (1)对n=2,5,8,计算A 的条件数;随n 增大,矩阵性态如何变化? (2)对n=5,解方程组Ax=b ;设A 的最后一个元素有扰动10-4,再求解Ax=b (3)计算(2)扰动相对误差与解的相对偏差,分析它们与条件数的关系。 (4)你能由此解释为什么不用插值函数存在定理直接求插值函数而要用拉格朗日或牛顿插值法的原因吗? 相关MATLAB 函数提示: zeros(m,n) 生成m 行,n 列的零矩阵 ones(m,n) 生成m 行,n 列的元素全为1的矩阵 eye(n) 生成n 阶单位矩阵 rand(m,n) 生成m 行,n 列(0,1)上均匀分布的随机矩阵 diag(x) 返回由向量x 的元素构成的对角矩阵 tril(A) 提取矩阵A 的下三角部分生成下三角矩阵

数值分析实验报告1

实验一 误差分析 实验1.1(病态问题) 实验目的:算法有“优”与“劣”之分,问题也有“好”与“坏”之别。对数值方法的研究而言,所谓坏问题就是问题本身对扰动敏感者,反之属于好问题。通过本实验可获得一个初步体会。 数值分析的大部分研究课题中,如线性代数方程组、矩阵特征值问题、非线性方程及方程组等都存在病态的问题。病态问题要通过研究和构造特殊的算法来解决,当然一般要付出一些代价(如耗用更多的机器时间、占用更多的存储空间等)。 问题提出:考虑一个高次的代数多项式 )1.1() ()20()2)(1()(20 1∏=-=---=k k x x x x x p 显然该多项式的全部根为1,2,…,20共计20个,且每个根都是单重的。现考虑该多项式的一个扰动 )2.1(0 )(19=+x x p ε 其中ε是一个非常小的数。这相当于是对(1.1)中19x 的系数作一个小的扰动。我们希望比较(1.1)和(1.2)根的差别,从而分析方程(1.1)的解对扰动的敏感性。 实验内容:为了实现方便,我们先介绍两个Matlab 函数:“roots ”和“poly ”。 roots(a)u = 其中若变量a 存储n+1维的向量,则该函数的输出u 为一个n 维的向量。设a 的元素依次为121,,,+n a a a ,则输出u 的各分量是多项式方程 01121=+++++-n n n n a x a x a x a 的全部根;而函数 poly(v)b = 的输出b 是一个n+1维变量,它是以n 维变量v 的各分量为根的多项式的系数。可见“roots ”和“poly ”是两个互逆的运算函数。 ;000000001.0=ess );21,1(zeros ve = ;)2(ess ve =

数值分析实验报告3

实验报告 实验项目名称数值积分与数值微分实验室数学实验室 所属课程名称数值逼近 实验类型算法设计 实验日期 班级 学号 姓名 成绩

实验概述: 【实验目的及要求】 本次实验的目的是熟练《数值分析》第四章“数值积分与数值微分”的相关内容,掌握复合梯形求积公式、复合辛普森求积公式、龙贝格求积公式以及高斯-勒让德公式。 本次试验要求编写复合梯形求积公式、复合辛普森求积公式、龙贝格求积公式以及高斯-勒让德公式的程序编码,并在MATLAB软件中去实现。 【实验原理】 《数值分析》第四章“数值积分与数值微分”的相关内容,包括:复合梯形求积公式、复合辛普森求积公式、龙贝格求积公式以及高斯-勒让德公式的相应算法和相关性质。 【实验环境】(使用的软硬件) 软件: MATLAB 2012a 硬件: 电脑型号:联想 Lenovo 昭阳E46A笔记本电脑 操作系统:Windows 8 专业版 处理器:Intel(R)Core(TM)i3 CPU M 350 @2.27GHz 2.27GHz 实验内容: 【实验方案设计】 第一步,将书上关于复合梯形求积公式、复合辛普森求积公式、龙贝格求积公式以及高斯-勒让德公式的内容转化成程序语言,用MATLAB实现;第二步,分别用以上求积公式的程序编码求解不同的问题。 【实验过程】(实验步骤、记录、数据、分析) 实验的主要步骤是:首先分析问题,根据分析设计MATLAB程序,利用程序算出问题答案,分析所得答案结果,再得出最后结论。 实验:用不同数值方法计算积分 (1) 取不同的步长h.分别用复合梯形及复合辛普森求积计算积分,给出误差中关于h的函数,并与积分精确值比较两个公式的精度,是否存在一个最小的h,使得精度不能再被改善? (2) 用龙贝格求积计算完成问题(1)。 (3)用勒让德多项式确定零点,再代入计算高斯公式,使其精度达到10-4 (1)在MATLAB的Editor中建立一个M-文件,输入程序代码,实现复合梯形求积公式的程序代码如下:

数值计算方法实验报告

差值法实验日志 实验题目:插值法 实验目的: 1.掌握拉格朗日插值、牛顿插值、分段低次插值和样条插值的方法。 2.对四种插值结果进行初步分析。 实验要求: (1)写出算法设计思想; (2)程序清单; (3)运行的结果; (4)所得图形; (5)四种插值的比较; (6)对运行情况所作的分析以及本次调试程序所取的经验。如果程序未通过,应分析其原因。 实验主要步骤: 1.已知函数) f满足: (x x0.0 0.1 0.195 0.3 0.401 0.5 f(0.39894 0.39695 0.39142 0.38138 0.36812 x ) 0.35206 (1)用分段线性插值; 打开MATLAB,按以下程序输入: x0=-5:5; y0=1./(1+x0.^2); x=-5:0.1:5; y=1./(1+x.^2); y1=lagr(x0,y0,x); y2=interp1(x0,y0,x); y3=spline(x0,y0,x);

for k=1:11 xx(k)=x(46+5*k); yy(k)=y(46+5*k); yy1(k)=y1(46+5*k); yy2(k)=y2(46+5*k); yy3(k)=y3(46+5*k); end [xx;yy;yy2;yy3]' z=0*x; plot(x,z,x,y,'k--',x,y2,'r') plot(x,z,x,y,'k--',x,y1,'r') pause plot(x,z,x,y,'k--',x,y3,'r') 回车得以下图形:

(2) 拉格朗日插值。 创建M 文件,建立lagr 函数: function y=lagr1(x0,y0,x) n=length(x0);m=length(x); for i=1:m z=x(i); s=0.0; for k=1:n p=1.0; for j=1:n if j~=k p=p*(z-x0(j))/(x0(k)-x0(j)); end end s=p*y0(k)+s; end y(i)=s; end 新建一个M 文件,输入: x0=[0.0 0.1 0.195 0.3 0.401 0.5]; y0=[0.39894 0.39695 0.39142 0.38138 0.36812 0.35206]; x=0.0:0.01:0.5; y1=lagr1(x0,y0,x); 00.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

数值分析实验报告2

一、实验名称 复合梯形求积公式、复合辛普森求积公式、龙贝格求积公式及自适应辛普森积分。 二、实验目的及要求 1. 掌握复合梯形求积计算积分、复合辛普森求积计算积分、龙贝格求积计算积分和自适应辛普森积分的基本思路和步骤. 2. 培养Matlab 编程与上机调试能力. 三、实验环境 计算机,MATLAB 软件 四、实验内容 1.用不同数值方法计算积分9 4 ln 1 0-=? xdx x 。 (1)取不同的步长h 。分别用复合梯形及复合辛普森求积计算积分,给出误差中关于h 的函数,并与积分精确指比较两个公式的精度,是否存在一个最小的h ,使得精度不能再被改善。 (2)用龙贝格求积计算完成问题(1)。 (3)用自适应辛普森积分,使其精度达到10-4。 五、算法描述及实验步骤 1.复合梯形公式 将区间[a,b]划分为n 等份,分点x k =a+ah,h=(b-a)/h,k=0,1,...,n ,在每个子区间[x k ,x k +1](k=0,1,...,n-1)上采用梯形公式(),得 )]()([2 )(b f a f a b dx x f b a +-≈ ? () )]()(2)([2)]()([21 1 110b f x f b f h x f x f h T n k k k n k k n ++=+=∑∑-=+-= () ),(),(12 )(' '2b a f h a b f R n ∈-- =ηη () 其中Tn 称为复合梯形公式,Rn 为复合梯形公式的余项。 2.复合辛普森求积公式 将区间[a,b]划分为n 等份,在每个子区间[x k ,x k +1](k=0,1,...,n-1)上采用辛普森公式(),得 )]()2 (4)([6b f b a f a f a b S +++-= ()

数值分析实验报告97389

数值分析实验报告 (第二章) 实验题目: 分别用二分法、牛顿迭代法、割线法、史蒂芬森迭代法求方程 f(x)=(x2+1)(x?1)5=0 的根x=1,观察不同初始值下的收敛性,并给出结论。 问题分析: 题目有以下几点要求: 1.不同的迭代法计算根,并比较收敛性。 2.选定不同的初始值,比较收敛性。 实验原理: 各个迭代法简述 二分法:取有根区间[a,b]的重点x0,确定新的有根区间[a1,b1]的区间长度仅为[a,b]区间长度的一版。对压缩了的有根区间[a1,b1]重复以上过程,又得到新的有根区间[a2,b2],其区间长度为[a1,b1]的一半,如此反复,……,可得一系列有根区间,区间收敛到一个点即为根。 牛顿迭代法:不动点迭代法的一种特例,具有局部二次收敛的特性。迭代格式为 x n+1=x n?f(x n) f′(x n) ,n=0,1,2,… 割线法:是牛顿法的改进,具有超线性收敛的特性,收敛阶为1.618. 迭代格式为 x n+1=x n? f(x n) (n)(n?1) (x n?x n?1),n=1,2,… 史蒂芬森迭代法:采用不动点迭代进行预估校正。至少是平方收敛的。迭代格式为 y n=φ(x n) z n=φ(y n) x n+1=x n? (y n?x n)2 z n?2y n+x n 这里φ(x)可采用牛顿迭代法的迭代函数。实验内容:

1.写出该问题的f(x)函数代码如下: function py= f(x) syms k; y=(k^2+1)*(k-1)^5; yy=diff(y,k); py(1)=subs(y,k,x); py(2)=subs(yy,k,x); end 2.分别写出各个迭代法的迭代函数代码如下: 二分法: function y=dichotomie(a,b,e) i=2; m(1)=a; while abs(a-b)>e t=(a+b)/2; s1=f(a); s2=f(b); s3=f(t); if s1(1)*s3(1)<=0 b=t; else a=t; end m(i)=t; i=i+1; end y=[t,i+1,m]; end 牛顿迭代法: function y=NewtonIterative(x,e) i=2; en=2*e;m(1)=x; while abs(en)>=e s=f(x); t=x-s(1)/s(2); en=t-x; x=t; m(i)=t; i=i+1; end y=[x,i+1,m]; end 牛顿割线法: function y=Secant(x1,x2,e) i=3; m(1)=x1,m(2)=x2; while abs(x2-x1)>=e s1=f(x1); s2=f(x2); t=x2-(x2-x1)*s2(1)/(s2(1)-s1( 1)); x1=x2; x2=t; m(i)=t; i=i+1; end

数值分析实验报告77712

《数值分析》 实验报告 学院:计算机科学与软件学院姓名:XXX 班级:计算机XX班 学号:XXXXXX 实验一:舍入误差与数值稳定性

实验目的: 1、 通过上机编程,复习巩固以前所学程序设计语言; 2、 通过上机计算,了解舍入误差所引起的数值不稳定性。 3、 通过上机计算,了解运算次序对计算结果的影响,从而尽量避免大数吃小数的现象。 实验内容:用两种不同的顺序计算644834.1100001 2≈∑=-n n ,分析其误差 的变化。 实验流程图: 实验源程序:

#include #include void main() { int i; float s1=0,s2=0,d1,d2; for (i=1;i<=10000;i++) s1=s1+1.0f/(i*i); for (i=10000;i>=1;i--) s2=s2+1.0f/(i*i); d1=(float)(fabs(1.644834-s1)); d2=(float)(fabs(1.644834-s2)); printf("正向求和结果为%f\n 误差为%f\n\n",s1,d1); printf("反向求和结果为%f\n 误差为%f\n\n",s2,d2); if(d1

实验分析:第一次做数值实验,又一次使用C语言编程,没有了刚学习C语言的艰难,能够将实验步骤转换成流程图并编写出完整的实验代码,在经过多次调试、改正后得到正确的程序和结果。这个实验较简单,计算误差时如果输入数据有误差,而在计算过程中舍入误差不增长,则称此算法是稳定的,否则称此算法是数值不稳定的,减少运算次数可以减小舍入误差。在运算中,如果参加运算的数的数量级相差很大,而计算机位数有限,如不注意运算次序就可能出现大数“吃掉”小数的现象,进而影响计算结果的可靠性,所以计算过程中要注意运算次序,避免出现这种现象。 实验二:拉格朗日插值法和牛顿插值法 实验目的:分别用拉格朗日差值和牛顿插值解决数学问题,并比较各方法的优略。 1、拉格朗日插值 实验内容: x i -3.0-1.0 1.0 2.0 3.0 y i 1.0 1.5 2.0 2.0 1.0 作二次插值,并求x 1=-2,x 2 =0,x 3 =2.75时的函数近似值。

相关主题