数值分析实验报告_清华大学__线性代数方程组的数值解法

线性代数方程组的数值解法

实验1.主元的选取与算法的稳定性

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

n

n

n R

b R

A b Ax ∈∈=?,,

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

（1）取矩阵??

????

?

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

?

??

?=141515

7,68

168

16816

b A ，则方程有解T x )1,,1,1(* =。取n=10

计算矩阵的条件数。让程序自动选取主元，结果如何？

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

（3）取矩阵阶数n=20或者更大，重复上述实验过程，观察记录并分析不同的问题及消去过程中选择不同的主元时计算结果的差异，说明主元素的选取在消去过程中的作用。 （4）选取其他你感兴趣的问题或者随机生成矩阵，计算其条件数。重复上述实验，观察记录并分析实验结果。

1.1程序清单

n=input('矩阵A 的阶数:n=');

A=6*diag(ones(1,n))+diag(ones(1,n-1),1)+8*diag(ones(1,n-1),-1); b=A*ones(n,1);

p=input('计算条件数使用p-范数,p='); cond_A=cond(A,p) [m,n]=size(A); Ab=[A b];

r=input('选主元方式(0:自动;1:手动)，r=');

Ab

for i=1:n-1

switch r

case(0)

[aii,ip]=max(abs(Ab(i:n,i)));

ip=ip+i-1;

case (1)

ip=input(['第',num2str(i),'步消元，请输入第',num2str(i),'列所选元素所处的行数：']);

end;

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

aii=Ab(i,i);

for k=i+1:n

Ab(k,i:n+1)=Ab(k,i:n+1)-(Ab(k,i)/aii)*Ab(i,i:n+1);

end;

if r==1

Ab

end

end;

x=zeros(n,1);

x(n)=Ab(n,n+1)/Ab(n,n);

for i=n-1:-1:1

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

end

x

1.2运行结果

(1)n=10，矩阵的条件数及自动选主元

Cond(A,1) =2.5575×103

Cond(A,2) =1.7276×103

Cond(A,inf) =2.5575×103

程序自动选择主元（列主元）

a.输入数据

矩阵A的阶数:n=10

计算条件数使用p-范数,p=1

选主元方式(0:自动;1:手动)，r=0

b.计算结果

x=[1，1，1，1，1，1，1，1，1，1]T (2)n=10，手动选主元

a.每步消去过程总选取按模最小或按模尽可能小的元素作为主元 矩阵A 的阶数:n=10 计算条件数使用p-范数,p=1 选主元方式(0:自动;1:手动)，r=1

(1)

(1)

61786115[]8

61158

6

14A

b ???????

?=???????

?

第1步消元，请输入第1列所选元素所处的行数：1

(2)

(2)

6.0000

1.00007.00004.6667 1.0000 5.6667

8.0000

6.000015.0000[]8.0000

1.000015.00006.0000 1.00008.0000

6.0000 1.000015.00008.0000

6.0000

14.0000A

b ??

??

??

??

??

=?

?

????

?????

?

第2步消元，请输入第2列所选元素所处的行数：2 …（实际选择时，第k 步选择主元处于第k 行） 最终计算得

x=[1.000000000000000, 1.000000000000000, 1.000000000000000, 1.000000000000001, 0.999999999999998, 1.000000000000004, 0.999999999999993, 1.000000000000012, 0.999999999999979, 1.000000000000028]T

b.每步消去过程总选取按模最大的元素作为主元 矩阵A 的阶数:n=10 计算条件数使用p-范数,p=1

选主元方式(0:自动;1:手动)，r=1

(1)

(1)

61786115[]8

61158

6

14A

b ???????

?=???????

?

第1步消元，请输入第1列所选元素所处的行数：2

(2)

(2)

8.0000

6.0000 1.000015.0000-3.50000.7500-4.25000

8.0000

6.0000 1.000015.0000[]8.0000

6.000015.00008.0000

1.00006.0000 1.000015.00008.0000

6.0000

14.0000A

b ??

??

-??

??

??

=?

?

????

?????

?

第2步消元，请输入第2列所选元素所处的行数：3 …（实际选择时，第k 步选择主元处于第k+1行） 最终计算得

x=[1，1，1，1，1，1，1，1，1，1]T (3)n=20，手动选主元

a.每步消去过程总选取按模最小或按模尽可能小的元素作为主元 矩阵A 的阶数:n=20 计算条件数使用p-范数,p=1 选主元方式(0:自动;1:手动)，r=1

(1)

(1)

6

1786115[]8

61158

6

14A

b ???????

?=???????

?

第1步消元，请输入第1列所选元素所处的行数：1

(2)

(2)

6.0000 1.0000

7.00004.6667 1.0000 5.6667

8.0000

6.000015.0000[]8.0000

1.000015.00006.0000 1.00008.0000

6.0000 1.000015.00008.0000

6.0000

14.0000A

b ??

??

??

??

??

=?

?

????

?????

?

第2步消元，请输入第2列所选元素所处的行数：2 …（实际选择时，第k 步选择主元处于第k 行） 最终计算得

x=[1.000000000000000，1.000000000000000，1.000000000000000，1.000000000000001，0.999999999999998，1.000000000000004，0.999999999999993，1.000000000000014，0.999999999999972，1.000000000000057，0.999999999999886，1.000000000000227，0.999999999999547，1.000000000000902，0.999999999998209，1.000000000003524，0.999999999993179，1.000000000012732，0.999999999978173，1.000000000029102]T b.每步消去过程总选取按模最大的元素作为主元 矩阵A 的阶数:n=20 计算条件数使用p-范数,p=1 选主元方式(0:自动;1:手动)，r=1

(1)

(1)

6

1786115[]8

61158

6

14A

b ???????

?=???????

?

第1步消元，请输入第1列所选元素所处的行数：2

(2)

(2)

8.0000

6.0000 1.000015.0000-3.50000.7500-4.25000

8.0000

6.0000 1.000015.0000[]8.0000

6.000015.00008.0000

1.00006.0000 1.000015.00008.0000

6.0000

14.0000A

b ??

??

-??

??

??

=?

?

????

?????

?

第2步消元，请输入第2列所选元素所处的行数：3

…（实际选择时，第k步选择主元处于第k+1行）

最终计算得

x=[1，1，1，1，1，1，1，1，1，1，1，1，1，1，1，1，1，1，1，1]T

(4)A分别为幻方矩阵，Hilbert矩阵，pascal矩阵和随机矩阵

将计算结果列于下表：

1.3简要分析

计算(1)表明：对于同一矩阵，不同范数定义的条件数是不同的；Gauss消去法在消去过程中选择模最大的主元能够得到比较精确的解。

计算(2)表明：通过比较每次选取模最大或模最小的元素的选主元方式，可以发现，在本题给定的问题中，选取模最大的元素作为主元比取模最小的元素作为主元时产生的结果更精确。因为这样做可以避免使用较小的数做为除数，以免发生结果数量级加大，使大数吃掉小数，产生舍入误差。

计算(3)表明：首先，n=20得到与计算(2),即n=10时一样的结论，即选取模最大的元素作为主元比取模最小的元素作为主元时产生的结果更精确；其次，与计算(2) （Cond(A10×10,1) =2.5575×103）比较，Cond(A20×20,1) =2.6214×106显著增大，且计算(3)的误差也远大于计算(2)，即，矩阵的条件数越大，产生的误差也越大。

计算(4)表明： Gauss消去法在消去过程中，主元的选择与算法的稳定性有密切的联系。一般来说，选取绝对值大的元素作为主元比绝对值小的元素作为主元时的计算结果更加精确。但这并不是绝对的，一些特殊的方阵，如Pascal方矩阵，则恰恰是选择模最小的元素作为主元时计算结果最精确（选模最小的元素只是一个表象，这种选主元方法优于其他选主元方法的本质是这种选择方法能使消去过程不产生浮点数，而全是整数运算，只有在回代过程中才有可能会产生浮点数）。在系数矩阵性质未知，或者说对于绝大多数的系数矩阵来说，选择模最大的元素作为主元是一种比较稳定和精确的方法。

实验2.病态的线性方程组的求解

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

实验内容：考虑方程组Hx=b 的求解，其中系数矩阵H 为Hilbert 矩阵，

n

j i j i h h H j i n n j i ,,2,1,,1

1,)(,, =-+=

=?

这是一个著名的病态问题。通过首先给定解（例如取为各个分量均为1）再计算出右端b 的办法给出确定的问题。 实验要求：

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

（2）逐步增大问题的维数，仍然用上述的方法来解它们，计算的结果如何？计算的结果说明了什么？ （3）讨论病态问题求解的算法

2.1程序清单

clear

n=6; % Hilbert 矩阵的阶数 H=hilb(n);

cond_H=cond(H); b=H*ones(n,1); D=diag(diag(H)); U=-triu(H,1); L=-tril(H,-1); w=1.5; tol=1.0e-6; maxk=10000; Ab=[H b];

% Gauss 消去法 for i=1:n-1

[aii,ip]=max(abs(Ab(i:n,i))); ip=ip+i-1;

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

aii=Ab(i,i);

if abs(aii)<1.0e-20

disp('Singular matrix, stop!');

pause; % 如果主元绝对值小于限值，则认为A奇异，终止end

for k=i+1:n

Ab(k,i:n+1)=Ab(k,i:n+1)-(Ab(k,i)/aii)*Ab(i,i:n+1);

end;

end;

x(n)=Ab(n,n+1)/Ab(n,n);

for i=n-1:-1:1

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

end

x_Gauss=x

% Jacobi迭代法

BJ=D\(L+U);

fJ=D\b;

x(:,1)=1.2*ones(n,1);

x(:,2)=BJ*x(:,1)+fJ;

k=2;

while (norm(x(:,k)-x(:,k-1))>=tol)&&(k

k=k+1;

x(:,k)=BJ*x(:,k-1)+fJ;

end

k_J=k

x_J=x(:,k)

% GS迭代法

BG=(D-L)\U;

fG=(D-L)\b;

x(:,1)=1.2*ones(n,1);

x(:,2)=BG*x(:,1)+fG;

k=2;

while (norm(x(:,k)-x(:,k-1))>=tol)&&(k

k=k+1;

x(:,k)=BG*x(:,k-1)+fG;

end;

k_G=k

x_G=x(:,k)

% SOR迭代法

lw=(D-w*L)\((1-w)*D+w*U);

fS=(D-w*L)\b*w;

x(:,1)=1.2*ones(n,1); x(:,2)=lw*x(:,1)+fS; k=2;

while (norm(x(:,k)-x(:,k-1))>=tol)&&(k

x(:,k)=lw*x(:,k-1)+fS; end

k_SOR=k x_SOR=x(:,k)

2.2运行结果及简要分析

(1)6阶Hilbert 矩阵

以

()

(-1)

1.06

k k x

x

tol e -<=-作为收敛的标准：J 法迭代矩阵的谱半径大于1，

迭代不收敛；

GS 法和SOR 法迭代矩阵的谱半径都略小于1，迭代是收敛的，但收敛速度非常慢，在

很多次迭代之后仍与精确解有一定误差，GS 收敛速度略快于SOR 法。 (2)n 阶Hilbert 矩阵

计算从6阶到25阶的Hilbert 矩阵，为观察收敛速度，以

()

(-1)

1.05

k k x

x

tol e -<=-作

为收敛的标准。计算结果如下（仅列出有代表性的6、7、10、11、12、13、21、25阶计算结果）

得到的主要结论如下：

(1)Gauss 消去法：Gauss 消去法求得的解与精确解的误差随Hilbert 矩阵阶数的增加而增加，Hilbert 矩阵阶数不大于11时，误差较小（小于1%），对于要求不要高的工程问题，这样的误差可以接受。当阶数大于11时，误差迅速增加，当阶数为13时，误差已经超过100%，一般来说这样的近似解不可接受。当阶数为25时，误差达到38900%。即低阶Hilbert 矩阵可用Gauss 消元法求解。

(2)Jacobi 迭代方法：无论Hilbert 矩阵为多少阶，Jacobi 迭代矩阵的谱半径都大于1，迭代不稳定、不收敛。

(3)GS 迭代方法：GS 迭代矩阵的谱半径略小于1，迭代收敛，但收敛速度非常缓慢。 (4)SOR 迭代方法：取w=1.5，SOR 迭代矩阵的谱半径略小于1，迭代收敛，但收敛速度亦非常缓慢。在本实验中，多数情况SOR 法较GS 迭代更慢一些（收敛步数更大）。

从上面的结果可以看出：Hilbert 矩阵阶数较小时，可用Gauss 消元法直接求解，解的精度比迭代法高；随着阶数增加，Gauss 消元法的精度迅速下降，解变得不可靠，这时，有一些迭代法，如GS 法和SOR 法仍可继续求得比较精确的解。另外，在三种迭代法中，GS 和SOR 法相对Jacobi 法更有优势，但这两种方法的迭代矩阵谱半径已经非常接近1（病态问题），收敛速度都很慢。

(3)病态问题的求解

求解病态问题，主要的方法是对原方程进行预处理，以降低系数矩阵的条件数。例如选择非奇异矩阵n n R Q P ?∈,，使方程组b Ax =化为等价方程组Pb y PAQ =)(，原方

6

81012

14161820222426

Hilbert 矩阵的阶数

迭代步数

程的解Qy x =。原则上应使矩阵PAQ 的条件数比A 有所改善。一般P 和Q 可选择为三角形矩阵或对角矩阵。理论上最好选择对角阵1D 和2D ，满足：

1212()m in ()

c o n

d D A D c o n d D A D =。

对于1—100阶的Hilbert 矩阵，假设1

2D D D

==，从设H D 为由Hilbert 矩阵对角元

素组成的对角阵，设D 是0.5

0.4

0.3

0.2

1;23;4;;H H

H H

D D D D D D D D ----??===?=????中，使得()

co n d D A D 最小的一种。

()/()

co n d D H D co n d H 的结果如下图所示：

虽然这里取的并不是最优的D 矩阵，但在上图中，可观察到经过预处理后的DHD 的条件数相比原Hilbert 矩阵减小了，在一定程度上改善了原Hilbert 矩阵的性质。

010203040

50

60708090100

0.2

0.4

0.6

0.8

1

n

c o n

d (D H D )/c o n d (H )

李庆扬数值分析第五版习题复习资料清华大学出版社

第一章 绪论 1．设0x >,x 的相对误差为δ，求ln x 的误差。 解：近似值* x 的相对误差为* **** r e x x e x x δ-= = = 而ln x 的误差为()1 ln *ln *ln ** e x x x e x =-≈ 进而有(ln *)x εδ≈ 2．设x 的相对误差为2%，求n x 的相对误差。 解：设()n f x x =，则函数的条件数为'() | |() p xf x C f x = 又1 '()n f x nx -=Q , 1 ||n p x nx C n n -?∴== 又((*))(*)r p r x n C x εε≈?Q 且(*)r e x 为2 ((*))0.02n r x n ε∴≈ 3．下列各数都是经过四舍五入得到的近似数，即误差限不超过最后一位的半个单位，试指 出它们是几位有效数字：*1 1.1021x =,*20.031x =, *3385.6x =, *456.430x =,* 57 1.0.x =? 解：* 1 1.1021x =是五位有效数字； *20.031x =是二位有效数字； *3385.6x =是四位有效数字； *456.430x =是五位有效数字； *57 1.0.x =?是二位有效数字。 4．利用公式(2.3)求下列各近似值的误差限：(1) ***124x x x ++,(2) ***123x x x ,(3) **24/x x . 其中**** 1234,,,x x x x 均为第3题所给的数。 解：

*4 1* 3 2* 13* 3 4* 1 51()1021()1021()1021()1021()102 x x x x x εεεεε-----=?=?=?=?=? *** 124***1244333 (1)()()()() 1111010102221.0510x x x x x x εεεε----++=++=?+?+?=? *** 123*********123231132143 (2)() ()()() 111 1.10210.031100.031385.610 1.1021385.610222 0.215 x x x x x x x x x x x x εεεε---=++=???+???+???≈ ** 24**** 24422 *4 33 5 (3)(/) ()() 11 0.0311056.430102256.43056.430 10x x x x x x x εεε---+≈ ??+??= ?= 5计算球体积要使相对误差限为1，问度量半径R 时允许的相对误差限是多少？ 解：球体体积为343 V R π= 则何种函数的条件数为 2 3'4343 p R V R R C V R ππ===g g (*)(*)3(*)r p r r V C R R εεε∴≈=g 又(*)1r V ε=Q

清华大学数值分析A第一次作业

7、设y0=28，按递推公式 y n=y n?1? 1 100 783，n=1,2,… 计算y100，若取≈27.982，试问计算y100将有多大误差？ 答：y100=y99?1 100783=y98?2 100 783=?=y0?100 100 783=28?783 若取783≈27.982，则y100≈28?27.982=0.018，只有2位有效数字，y100的最大误差位0.001 10、设f x=ln?(x? x2?1)，它等价于f x=?ln?(x+ x2?1)。分别计算f30，开方和对数取6位有效数字。试问哪一个公式计算结果可靠？为什么？ 答： x2?1≈29.9833 则对于f x=ln x?2?1，f30≈?4.09235 对于f x=?ln x+2?1，f30≈?4.09407 而f30= ln?(30?2?1) ，约为?4.09407，则f x=?ln?(x+ x2?1)计算结果更可靠。这是因为在公式f x=ln?(x? x2?1)中，存在两相近数相减（x? x2?1）的情况，导致算法数值不稳定。 11、求方程x2+62x+1=0的两个根，使它们具有四位有效数字。 答：x12=?62±622?4 2 =?31±312?1 则 x1=?31?312?1≈?31?30.98=?61.98 x2=?31+312?1= 1 31+312?1 ≈? 1 ≈?0.01613

12.（1）、计算101.1?101，要求具有4位有效数字 答：101.1?101= 101.1+101≈0.1 10.05+10.05 ≈0.004975 14、试导出计算积分I n=x n 4x+1dx 1 的一个递推公式，并讨论所得公式是否计算稳定。 答：I n=x n 4x+1dx 1 0= 1 4 4x+1x n?1?1 4 x n?1 4x+1 dx= 1 1 4 x n?1 1 dx?1 4 x n?1 4x+1 dx 1 = 1 4n ? 1 4 I n?1，n=1，2… I0= 1 dx= ln5 1 记εn为I n的误差，则由递推公式可得 εn=?1 εn?1=?=(? 1 )nε0 当n增大时，εn是减小的，故递推公式是计算稳定的。

清华大学高等数值计算(李津)实践题目一(共轭梯度CG法,Lanczos算法与MINRES算法)

高等数值计算实践题目一 1. 实践目的 本次计算实践主要是在掌握共轭梯度法，Lanczos 算法与MINRES 算法的基础上，进一步探讨这3种算法的数值性质，主要研究特征值特征向量对算法收敛性的影响。 2. 实践过程 （一）生成矩阵 （1）作5个100阶对角阵i D 如下： 1D 对角元：1,1,...,20,1+0.1(-20),21,...,100j j d j d j j ==== 2D 对角元：1,1,...,20,1+(-20),21,...,100j j d j d j j ==== 3D 对角元：,1,...,80,81,81,...,100j j d j j d j ==== 4D 对角元：,1,...,40,41,41,...,60,41+(60),61,...,100j j j d j j d j d j j =====-= 5D 对角元：,1,...,100j d j j == 记i D 的最大模特征值和最小模特征值分别为1i λ和i n λ，则i D 特征值分布有如下特点： 1D 的特征值有较多接近于i n λ，并且1/i i n λλ较小， 2D 的特征值有较多接近于i n λ，并且1/i i n λλ较大， 3D 的特征值有较多接近于1i λ，并且1/i i n λλ较大， 4D 的特征值有较多接近于中间模特征值，并且1/i i n λλ较大， 5D 的特征值均匀分布，并且1/i i n λλ较大 （2）随机生成10个100阶矩阵j M ： (100(100))j M fix rand = 并作它们的QR 分解，得j Q 和j R ，这样可得50个对称的矩阵T ij j i j A Q DQ =，其中i D 的对角元就是ij A 的特征值，若它们都大于0，则ij A 正定，j Q 的列就是相应的特征向量。结合（1）可知，ij A 都是对称正定阵。

数值分析实验报告模板

数值分析实验报告模板 篇一：数值分析实验报告(一)(完整) 数值分析实验报告 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() ()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 +

清华大学贾仲孝老师高等数值分析报告第二次实验

高等数值分析第二次实验作业

T1.构造例子特征值全部在右半平面时, 观察基本的Arnoldi 方法和GMRES 方法的数值性态, 和相应重新启动算法的收敛性. Answer: （1） 构造特征值均在右半平面的矩阵A ： 根据实Schur 分解，构造对角矩阵D 由n 个块形成，每个对角块具有如下形式，对应一对特 征值i i i αβ± i i i i i S αββα-?? = ??? 这样D=diag(S 1,S 2,S 3……S n )矩阵的特征值均分布在右半平面。生成矩阵A=U T AU ，其中U 为 正交阵，则A 矩阵的特征值也均在右半平面。不妨构造A 如下所示： 2211112222 /2/2/2/2N N A n n n n ?-?? ? ? ?- ? = ? ? ? - ? ?? ? 由于选择初值与右端项：x0=zeros(2*N,1);b=ones(2*N,1); 则生成矩阵A 的过程代码如下所示： N=500 %生成A 为2N 阶 A=zeros(2*N); for a=1:N A(2*a-1,2*a-1)=a; A(2*a-1,2*a)=-a; A(2*a,2*a-1)=a; A(2*a,2*a)=a; end U = orth(rand(2*N,2*N)); A1 = U'*A*U; （2） 观察基本的Arnoldi 和GMRES 方法 编写基本的Arnoldi 函数与基本GMRES 函数，具体代码见附录。 function [x,rm,flag]=Arnoldi(A,b,x0,tol,m) function [x,rm,flag]=GMRES(A,b,x0,tol,m) 输入:A 为方程组系数矩阵，b 为右端项，x0为初值，tol 为停机准则，m 为人为限制的最大步数。 输出:x 为方程的解，rm 为残差向量，flag 为解是否收敛的标志。 外程序如下所示： e=1e-6; m=700;

数值分析实验报告之常微分方程数值解

数学与计算科学学院实验报告 实验项目名称常微分方程数值解 所属课程名称数值方法B 实验类型验证 实验日期 2013.11.11 班级 学号 姓名 成绩

【实验过程】（实验步骤、记录、数据、分析） 注：以下图形是通过Excel 表格处理数据得出，并未通过MATLAB 编程序所得！ 1、1(0)1dy y x dx y ?=-++???=? 由题可知精确解为：x y x e -=+，当x=0时，y(x)=0。 h=0.1 表1 h=0.1时三个方法与精确值的真值表 图1 h=0.1时三个方法走势图 步长 Euler 法 预估校正法 经典四阶库 精确值 0.1 1.010000 1.005000 1.004838 1.249080 0.2 1.029000 1.019025 1.018731 1.055455 0.3 1.056100 1.041218 1.040818 1.091217 0.4 1.090490 1.070802 1.070320 1.131803 0.5 1.131441 1.107076 1.106531 1.176851 0.6 1.178297 1.149404 1.148812 1.226025 0.7 1.230467 1.197211 1.196586 1.279016 0.8 1.287421 1.249975 1.249329 1.335536 0.9 1.348678 1.307228 1.306570 1.395322 1.0 1.413811 1.368541 1.367880 1.458127

h=0.05（此时将源程序中i的围进行扩大，即for(i=0;i<20;i++)） 表2 h=0.05时三个方法与精确值的真值表步长Euler法预估校正法经典四阶库精确值 0.05 1.002500 1.001250 1.001229 1.011721 0.10 1.007375 1.004877 1.004837 1.024908 0.15 1.014506 1.010764 1.010708 1.039504 0.20 1.023781 1.018802 1.018731 1.055455 0.25 1.035092 1.028885 1.028801 1.072710 0.30 1.048337 1.040915 1.040818 1.091217 0.35 1.063421 1.054795 1.054688 1.110931 0.40 1.080250 1.070436 1.070320 1.131801 0.45 1.098737 1.087752 1.087628 1.153791 0.50 1.118800 1.106662 1.106531 1.176851 0.55 1.140360 1.127087 1.126950 1.200942 0.60 1.163342 1.148954 1.148812 1.226025 0.65 1.187675 1.172193 1.172046 1.252062 0.70 1.213291 1.196736 1.196585 1.279016 0.75 1.240127 1.222520 1.222367 1.306852 0.80 1.268121 1.249485 1.249329 1.335536 0.85 1.297215 1.277572 1.277415 1.365037 0.90 1.327354 1.306728 1.306570 1.395322 0.95 1.358486 1.336900 1.336741 1.426362 1.00 1.390562 1.368039 1.367880 1.458127 图2 h=0.05时三个方法走势图

数值分析实验报告

实验一、误差分析 一、实验目的 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

数值分析实验报告

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

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

数值分析实验报告总结

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

如果 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 ，那 外，还规定其必须满足相容性: 所以

清华大学杨顶辉数值分析第6次作业

清华大学杨顶辉数值分析第6次作业

9.令*()(21),[0,1]n n T x T x x =-∈，试证*{()}n T x 是在[0,1]上带权 2 ()x x x ρ= -****0123(),(),(),()T x T x T x T x . 证明： 1 1 **2 1 1 * *20 12 2 1**20 ()()()(21)(21)211()()()()()211()22 ()()1()1()()()()()1n m n m n m n m n m n n m n m x T x T x dx x T x dx x x t x x T x T x dx t T t dt t t t T t dt t T x x x T x T x dx t T t t ρρρ---=---=-=++-= --= -???? ?令，则 由切比雪夫多项式1 01=02 m n dt m n m n ππ ≠??? =≠??==??? 所以*{()}n T x 是在[0,1]上带权2 ()x x x ρ= - *00*11* 2 2 2 2*33233()(21)1()(21)21 ()(21)2(21)188()(21)4(21)3(21)3248181 T x T x T x T x x T x T x x x x T x T x x x x x x =-==-=-=-=--=-=-=---=-+- 14.已知实验数据如下： i x 19 25 31 38 44 i y 19.0 32.3 49.0 73.3 97.8 用最小二乘法求形如2y a bx =+的经验公式，并求均方误差 解： 法方程为

清华大学杨顶辉数值分析第6次作业

9.令*()(21),[0,1]n n T x T x x =-∈，试证*{()}n T x 是在[0,1] 上带权()x ρ=的正交多项式，并求****0123(),(),(),()T x T x T x T x . 证明： 1 1 * *0 1 1 * *011**0 ()()()(21)(21)211()()()()()2()()()()()()()()n m n m n m n m n m n n m n m x T x T x dx x T x dx t x x T x T x dx t T t dt t T t dt T x x T x T x dx t T t ρρρ---=--=-== = ???? ?令，则 由切比雪夫多项式1 01=02 m n dt m n m n ππ ≠??? =≠??==??? 所以*{()}n T x 是在[0,1] 上带权()x ρ= *00*11* 22 2 2*33233()(21)1()(21)21 ()(21)2(21)188()(21)4(21)3(21)3248181 T x T x T x T x x T x T x x x x T x T x x x x x x =-==-=-=-=--=-=-=---=-+- 14.已知实验数据如下： 用最小二乘法求形如2y a bx =+的经验公式，并求均方误差 解： 法方程为

22222(1,)(1,1)(1,)(,)(,1)(,)a y x b x y x x x ?????? =???? ?????? ?? 即 5 5327271.453277277699369321.5a b ??????=???????????? 解得 0.972579 0.050035a b =?? =? 拟合公式为20.9725790.050035y x =+ 均方误差 2 4 2 2 0[]0.015023i i i y a bx σ==--=∑ 21.给出()ln f x x =的函数表如下： 用拉格朗日插值求ln 0.54的近似值并估计误差（计算取1n =及2n =） 解：1n =时，取010.5,0.6x x == 由拉格朗日插值定理有 1 100.60.5 0.693147 0.510826 0.50.(60.60.51.82321)0 1.()6047()52 j j j x x x L x f x l x ==------=-=∑ 所以1ln0.54(0.54)0.620219L ≈=- 误差为ln 0.54(0.620219)= 0.004032ε=-- 2n =时，取0120.4,0.5,0.6x x x === 由拉格朗日插值定理有

清华大学高等数值计算(李津)实践题目二(SVD计算及图像压缩)(包含matlab代码)

第1部分 方法介绍 奇异值分解（SVD ）定理： 设m n A R ?∈，则存在正交矩阵m m V R ?∈和n n U R ?∈，使得 T O A V U O O ∑??=?? ?? 其中12(,, ,)r diag σσσ∑=，而且120r σσσ≥≥≥>，(1,2, ,)i i r σ=称为A 的 奇异值，V 的第i 列称为A 的左奇异向量，U 的第i 列称为A 的右奇异向量。 注：不失一般性，可以假设m n ≥，（对于m n <的情况，可以先对A 转置，然后进行SVD 分解，最后对所得的SVD 分解式进行转置，就可以得到原来的SVD 分解式） 方法1：传统的SVD 算法 主要思想： 设()m n A R m n ?∈≥，先将A 二对角化，即构造正交矩阵1U 和1V 使得 110T B n U AV m n ?? =?? -?? 其中1200n n B δγγδ??? ???=?????? 然后，对三角矩阵T T B B =进行带Wilkinson 位移的对称QR 迭代得到:T B P BQ =。 当某个0i γ=时，B 具有形状12B O B O B ?? =? ??? ，此时可以将B 的奇异值问题分解为两个低阶二对角阵的奇异值分解问题；而当某个0i δ=时，可以适当选取'Given s 变换，使得第i 行元素全为零的二对角阵，因此，此时也可以将B 约化为两个低 阶二对角阵的奇异值分解问题。 在实际计算时，当i B δε∞≤或者() 1j j j γεδδ-≤+（这里ε是一个略大于机器精度的正数）时，就将i δ或者i γ视作零，就可以将B 分解为两个低阶二对角阵的奇异值分解问题。

数值分析实验报告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——矩阵基础

20130917题目 求证：在矩阵的LU 分解中，1 11n n T n ij i j j i j L I e e α-==+??=- ??? ∑∑ 证明： 在高斯消去过程中，假设0jj a ≠ ，若a=0，可以通过列变换使得前面的条件成立，这里不考虑这种情况。 对矩阵A 进行LU 分解，()() () ()()1 11 1111L M n M M M n ---=-=??-………… ， 其中()1n T n ij i j i j M j I e e α=+??=+ ??? ∑ ，i e 、j e 为n 维线性空间的自然基。 ()M j 是通过对单位阵进行初等变换得到， 通过逆向的变换则可以得到单位阵，由此很容易得到()M j 的逆矩阵为1n T n ij i j i j I e e α=+??- ???∑。故111n n T n ij i j n j i j L I e e I α-==+?? ??=- ? ? ????? ∏∑ 上式中的每一项均是初等变换，从右向左乘，则每乘一次相当于对右边的矩阵进行一次 向下乘法叠加的初等变换。由于最初的矩阵为单位阵，变换从右向左展开，因而每一次变换不改变已经更新的数据，既该变换是从右向左一列一列更新数据，故 11n n T n ij i j j i j L I e e α==+??=- ??? ∑∑。 数学证明：1n T ij i j i j e e α=+?? ???∑具有 ,0 00n j j A -?? ??? 和1,1000n j n j B -+-+?? ?? ? 的形式，且有 +1,-11,10000=000n j j n j n j A B --+-+???? ?????? ? 而1 1n n T ij i j j k i j e e α-==+?? ??? ∑∑具有1,1000n k n k B -+-+?? ???的形式，因此： 1 311111211121==n n n n n n T T T n ij i j n ij i j n ik i k j i j j i j k n i k n n T n i i n ik i i i k L I e e I e e I e e I e e I e ααααα---==+==+=-=+==+??????????????=---?? ? ? ? ? ? ? ? ???????????????????????=-- ? ? ?????∏∑∏∑∑∑∑∑……11211n n n T T k n ik i k k k i k e I e e α--===+????=- ?? ?????? ∑∑∑#

数值分析实验报告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次作业答案

2．定义映射22:B R R →，()B x y =，满足y Ax =，其中 0.80.40.10.4A ??=????，2,x y R ∈ 则对任意的2 ,u v R ∈ 1111119 ||()()||||||||()||||||||||||||10B u B v Au Av A u v A u v u v -=-=-≤-=- 故映射B 对一范数是压缩的 由范数定义 ||||1 ||||max |||| 1.2 x A Ax ∞∞∞===，知必然存在0 x ， 0||||1 x ∞= 使得0|||||||| 1.2 Ax A ∞∞== 设012(,)T x x x = 取 12(,0),(0,)T T u x v x ==-，则 u v x -=，有 00||()()||||||||()|||||||||| 1.21||||||||B u B v Au Av A u v Ax A x u v ∞∞∞∞∞∞∞ -=-=-===>==- 故有||()()||B u B v ∞->||||u v ∞ -，从而映射B 对无穷范数不是压缩的 4. 证明：对任意的,[,]x y a b ∈ 由拉格朗日中值定理，有 ()()'()()() 1e G x G y G x y x y e ξ ξξ-=-=-+ 其中0111b b e e e e ξξ<≤<++ 所以 |()()||()||| 11b b e e G x G y x y x y e e ξξ-=-≤-++ 故G 为[,]a b 上的压缩映射 而 ()ln(1)ln x x G x e e x =+>= 即()G x x =无根

数值计算方法实验报告

差值法实验日志 实验题目：插值法 实验目的： 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

清华大学高等数值分析实验设计及答案

高等数值分析实验一 工物研13 成彬彬2004310559 一．用CG，Lanczos和MINRES方法求解大型稀疏对称正定矩阵Ax=b 作实验中，A是利用A= sprandsym(S,[],rc,3)随机生成的一个对称正定阵，S是1043阶的一个稀疏阵 A= sprandsym(S,[],0.01,3); 检验所生成的矩阵A的特征如下： rank(A-A')=0 %即A=A’，A是对称的； rank(A)=1043 %A满秩 cond(A)= 28.5908 %A是一个“好”阵 1．CG方法 利用CG方法解上面的线性方程组 [x,flag,relres,iter,resvec] = pcg(A,b,1e-6,1043); 结果如下： Iter=35，表示在35步时已经收敛到接近真实x relres= norm(b-A*x)/norm(b)= 5.8907e-007为最终相对残差 绘出A的特征值分布图和收敛曲线： S=svd(A); %绘制特征值分布 subplot(211) plot(S); title('Distribution of A''s singular values');; xlabel('n') ylabel('singular values') subplot(212); %绘制收敛曲线 semilogy(0:iter,resvec/norm(b),'-o'); title('Convergence curve'); xlabel('iteration number'); ylabel('relative residual'); 得到如下图象：

为了观察CG方法的收敛速度和A的特征值分布的关系，需要改变A的特征值： (1).研究A的最大最小特征值的变化对收敛速度的影响 在A的构造过程中，通过改变A= sprandsym(S,[],rc,3)中的参数rc(1/rc为A的条件数)，可以达到改变A的特征值分布的目的： 通过改变rc=0.1，0.0001得到如下两幅图 以上三种情况下，由收敛定理2.2.2计算得到的至多叠代次数分别为：48，14和486，由于上实验结果可以看出实际叠代次数都比上限值要小较多。 由以上三图比较可以看出，A的条件数越大，即A的最大最小特征值的差别越大，叠代所需要的步骤就越多，收敛越慢。 (2)研究A的中间特征值的分布对于收敛特性的影响： 为了研究A的中间特征值的分布对收敛速度的影响，进行了如下实验： 固定A的条件数，即给定A的最大最小特征值，改变中间特征值得分布，再来生成A，具体的实现方法是，先将原来的生成A进行特征值分解： [U,S]=svd(A);