清华大学高等数值分析 第一次实验作业

高等数值分析第二次作业 第四题

姓名：---- 学号：------- 时间：11月20日

1、构造例子说明 CG 的数值性态. 当步数 = 阶数时 CG 的解如何? 当 A 的最大特征值远大于第二个最大特征值, 最小特征值远小于第二个最小特征值时方法的收敛性如何?

答：1）构造矩阵如下

首先，构造一个非常病态的矩阵。

构造矩阵的条件数K(A)为1015，可以得到收敛曲线如下图所示

010*******

400500600

7008009001000

10

-10

10-8

10-6

10-4

10-2

10

102

104

106

算法的收敛曲线 (阶数n=1002)

迭代次数

||r k ||/||b ||

第二种方法构造矩阵（良态）：

条件数不大的矩阵也能够使迭代步数大于等于矩阵阶数，构造一个良态的矩阵如下： 使其特征值为等比数列，即

1000

510

i i =λ

由此得到的矩阵A 最大特征值为105，最小为1，的条件数是105，特征值分布如下图

10

10

1

10

2

103

10

4

10

5

10

6

A 的特征值分布情况

特征值

收敛曲线如下图，发现使用CG 法仍然可以收敛，但是收敛速度非常缓慢，即使到1000步也难以达到理想的残差，说明即使问题的条件数较小，即良态的，CG 方法也可能收敛非常慢。

0200400

60080010001200

10

-3

10

-2

10

-1

10

10

1

算法的收敛曲线 (阶数n=1001)

迭代次数

||r k ||/||b ||

从上面的结果中可以看出，当矩阵的迭代步数等于阶数时

a. 往往收敛曲线不平滑，震荡非常严重，残差的2范数没有最优性；

b. 但是方法的相对残差大体上仍然具有下降趋势，最终仍然能够收敛。矩阵的病态性延迟了方

法的收敛。

c. 条件数较小的矩阵也会出现收敛缓慢的现象，这与特征值的分布有关。

2）构造矩阵1002阶的A ，其中，A 的最大的特征值为109，最小特征值为1，其具有如下的特征值分布

10

-2

10

10

2

104

10

6

10

8

10

10

A 的特征值分布情况

特征值

图 1 A 的特征值分布

该矩阵的条件数k=1011，是个病态问题，结果如下

050100150

200250300350400

10

-8

10

-6

10

-4

10

-2

10

10

2

10

4

CG 算法的收敛曲线 (n=1002)

迭代次数

||r k ||

图 2 A 的收敛性态图

代入CG 法，得到的收敛曲线如图 2 A 的收敛性态图所示。从图中可以看出，一开始，方法残差下降很慢，等到后面，在150步左右，以另一种收敛速度下降。实验结果与理论结果相同。

3）结论：

a. 对于非常病态的问题，CG 法的迭代次数可能大于等于阶数，残差震荡非常严重，但是大体

趋势仍然是下降的，最终也能够收敛；对于极其病态的问题，CG 法也可能无能为力，无法收敛。

b. 对于最大特征值远远大于第二大特征值的情况，收敛曲线出现了两种下降速度，说明特征值

的分布对于CG 方法的收敛速度影响很大。 c.

2、对于同样的例子, 比较 CG 和 Lanczos 的计算结果.

在相同的A 下，取b=(1,1,1,1...1)T ,x 0=(0,0,0...0)T ，分别用CG 法和Lanczos 法求解。停机准则为相对误差小于10-8。将A 取以下几种：良态正定矩阵、近似奇异正定矩阵、病态正定矩阵、良态不定矩阵进行实验，结果分别如下：

良态的正定问题

取矩阵A 的特征值分布如下：

10

10

1

10

2

10

3

A 的特征值分布情况

特征值

lamda=[1001, 1000:-1:1, 1]

从特征值的分布上可以看出，A 的条件数为103，是个良态正定的问题。CG 结果如下图

204060

80100120140160180

10

-10

10

-8

10

-6

10-4

10

-2

10

CG 算法的收敛曲线 (阶数n=1002)

迭代次数

||r k ||/||b ||

CG

Lanczos

从上图可以看出，对称正定良态问题，CG 和Lanczos 方法收敛性态几乎一模一样。 近似奇异的问题

上一个A 矩阵中的一个特征值减小，取为10-3，特征值分布为

10

-3

10

-2

10

-1

100

10

1

10

2

10

3

A 的特征值分布情况

特征值

50

100150

200250

10

-10

10

-8

10

-6

10

-4

10-2

10

10

2

CG 算法的收敛曲线 (阶数n=1002)

迭代次数

||r k ||/||b ||

CG

Lanczos

从收敛特性上看，对于近似奇异的正定矩阵（数值上不奇异），CG 和Lanczos 方法得到的结果一致。但是CG 的求解速度比Lanczos 小很多。 病态的正定问题

在以上相同的A （n=1002）的基础上，将最大特征值提高为1e15，则条件数为1e15，是个非常病态的问题，同时使用CG 法和Lanczos 方法得到的收敛曲线下

10

10

5

10

10

10

15

A 的特征值分布情况

特征值

100200

300400500600700

10

-10

10-8

10-6

10-4

10-2

100

102

104

106

CG 算法的收敛曲线 (阶数n=1002)

迭代次数

||r k ||/||b ||

CG

Lanczos

从结果中可以看出，对于这个问题，Lanczos 可能比CG 的迭代次数小，但是，程序的运行时间来看，CG 法为2.13s ，而Lanczos 为7.14s ，因此对于正定问题CG 法的效率要明显高于Lanczos 。 良态的不定问题

在以上相同的A 的基础上，增加两个负特征值-1和-10，在此基础上，使用CG 和Lanczos 计算，得到的结果如下

50

100150

200250

10

-10

10

-8

10

-6

10

-4

10-2

10

10

2

CG 算法的收敛曲线 (阶数n=1003)

迭代次数

||r k ||/||b ||

CG

Lanczos

从图中可以看出，虽然问题比较良性，但是也出现了尖峰，即近似中断现象。其中，CG 法也能够在没有中断的情况下求解不定问题，但是如果一旦中断，那么无法得到结果。但是Lanczos 中断之后，却能够继续进行下去。即使用Lanczos 解不定问题的风险更小。

总结：

1) 对于正定问题，不论良态还是病态，Lanczos 和CG 方法求解的结果非常相近，收敛曲线也近

似一致，收敛步数也相差不大，这说明两个方法的原理是一致的。但是CG 方法的速度明显快于Lanczos 法。

2) 对于不定问题，CG 方法也可能收敛，如果收敛，Lanczos 方法和CG 方法的到的结果也类似

（步数和收敛性态）。但是，如果CG 方法一旦中断，前面的所有运算都白费了，是恶性中断，而Lanczos 仍然能够继续。因此，再不定问题中，Lanczos 比CG 方法的风险小。

3、当A 只有m 个不同特征值时, 对于大的m 和小的m, 观察有限精度下 Lanczos 方法如何收敛.

答：对于不同的m 值进行实验，得到的结果如下：

当m=20时，取特征值分别为1,2,3,4...20，每个特征值的重数均为60重，得到结果为

10

101

10

2

A 的特征值分布情况

特征值

2468

101214161820

10

-16

10-14

10-12

10-10

10-8

10-6

10-4

10-2

100

CG 算法的收敛曲线 (阶数n=1200)

迭代次数

||r k ||/||b ||

CG

Lanczos

两个方法均在20步的时候收敛，即n it ≤m 。

当m=50时，取特征值分别为1,2,3,4...50，每个特征值的重数均为24重，得到的结果为

100101102

A的特征值分布情况

特征值

0510152025303540 10-10

10-8

10-6

10-4

10-2

100

CG算法的收敛曲线 (阶数n=1200)

迭代次数

|

|

r

k

|

|

/

|

|

b

|

|

CG

Lanczos

在38步左右收敛，CG法与Lanczos方法的收敛速度一致。同样也满足n it≤m。

当m=500时，取特征值分别为1,2,3,4...500，每个特征值的重数均为2重，得到的结果为

100101102103

A的特征值分布情况

特征值

020406080100120140 10-10

10-8

10-6

10-4

10-2

100

CG算法的收敛曲线 (阶数n=1000)

迭代次数

|

|

r

k

|

|

/

|

|

b

|

|

CG

Lanczos

从结果中可以看出，迭代120多步就收敛了，而且CG法与Lanczos法的结果一致。同样满足n it≤m。当m=1000时，取特征值分别为1,2,3,4...100，每个特征值的重数均为1重，得到的结果为

100101102103

A的特征值分布情况

特征值

020406080100120140160180 10-10

10-8

10-6

10-4

10-2

100

CG算法的收敛曲线 (阶数n=1000)

迭代次数

|

|

r

k

|

|

/

|

|

b

|

|

CG

Lanczos

经过多次实验，将得到的迭代次数与不同特征值个数的关系做图如下，从图中可以看出，当m很小

时，迭代次数等于不同特征值个数，随着m 增大，迭代次数也不断增大，但是增大的趋势变缓，最终，几乎不随着m 增大而增大。同时，发现，迭代次数始终小于不同特征值的个数，即与书上讲的“若A 有m 个不同的特征值，则至多m 步收敛”一致！

010*******

4005006007008009001000

20406080100120140

160180200不同特征值个数(m)

收敛时的迭代步长

不同特征值个数与收敛迭代步长的关系（停机准则=10-8)

注：图中蓝色曲线为y=x 曲线。

总结：

1) 特征值的分布对于Lanczos 方法的收敛性态影响非常大； 2) 如果A 有m 个不同的特征值，那么就会至多m 步收敛。

4、取初始近似解为零向量, 右端项 b 仅由 A 的 m 个不同个特征向量的线性组合表示时, Lanczos 方法的收敛性如何? 数值计算中方法的收敛性和 m 的大小关系如何? 答：首先，确定一个比较良态的矩阵A ，取其特征值为

lamda=[3100:10:2500, 1000:-1:1]

分布如下图

10

10

1

10

2

10

3

A 的特征值分布情况

特征值

分别取b 为

b m =Q m (1,1,1...1)m T

其中Q m 的列为A 的特征向量。

变化分别取m 的值为5,10,15,30,40,50,80,90,100,200,300,400,500,600,800,900,1000分别得到收敛曲线，下面仅列出m=5、15、40、50、90、100、400、600、900时的收敛曲线。

0246

81012141618

10

-10

10

-8

10

-6

10-4

10

-2

10

Lanczos 收敛曲线(b 由m 个特征向量组成）m=5)

迭代次数

||r k ||/||b ||

01020

3040506070

10

-10

10

-8

10

-6

10

-4

10

-2

10

10

2

Lanczos 收敛曲线(b 由m 个特征向量组成）m=15)

迭代次数

||r k ||/||b ||

01020

3040506070

10

-10

10

-8

10

-6

10

-4

10-2

10

10

2

Lanczos 收敛曲线(b 由m 个特征向量组成）m=40)

迭代次数

||r k ||/||b ||

01020

3040506070

10

-10

10

-8

10

-6

10

-4

10-2

10

10

2

Lanczos 收敛曲线(b 由m 个特征向量组成）m=50)

迭代次数

||r k ||/||b ||

02040

6080100120

10

-10

10

-8

10

-6

10

-4

10-2

10

10

2

Lanczos 收敛曲线(b 由m 个特征向量组成）m=90)

迭代次数

||r k ||/||b ||

02040

6080100120140

10

-10

10

-8

10

-6

10

-4

10-2

10

10

2

Lanczos 收敛曲线(b 由m 个特征向量组成）m=100)

迭代次数

||r k ||/||b ||

0204060

80100120140160180

10

-10

10

-8

10

-6

10-4

10

-2

10

Lanczos 收敛曲线(b 由m 个特征向量组成）m=400)

迭代次数

||r k ||/||b ||

0204060

80100120140160180

10

-10

10

-8

10

-6

10-4

10

-2

10

Lanczos 收敛曲线(b 由m 个特征向量组成）m=900)

迭代次数

||r k ||/||b ||

从图中可以看出，各个收敛性态都很好，曲线比较平滑。 通过实验，可以得到m 与收敛迭代步长的曲线如下图

0100

2003004005006007008009001000

20406080100120

140160180m （b 由m 个特征向量线性组合而成）

收敛迭代次数

收敛迭代次数与m 的关系 （b 由m 个特征向量线性组合而成）

注：图中蓝色曲线为y=x 曲线。

从图中可以看出，收敛性态与特征向量有关！随着m 的增大，迭代次数有上升趋势，即若b 仅有少数几个特征向量的线性组合而成，那么可以得到更好的收敛效果。注：图中蓝色曲线为y=x 曲线。

总结：

1) 收敛特性与特征向量有关；

2) 若b 由m 个特征向量线性组合而成，m 越小，收敛速度越快；

3) 理论上，至多m 步收敛，但是实验中并不是严格按照这个关系，说明特征向量问题比特征值

问题复杂一些。

5、构造对称不定的矩阵, 验证 Lanczos 方法的近似中断, 观察收敛曲线中的峰点个数和特征值的分布关系; 观察当出现峰点时, MINRES 方法的收敛性态怎样.

答：取b=(1,1,1,...,1,1)T ，x 0=(0,0,0,...,0,0)T ，停机准则为10-8。

首先，取特征值为

lamda=[3100:10:2500, 1000:-1:1,-1 ];

其中有一个负特征值-1，分别使用MINRES 和Lanczos 方法得到的曲线如下图

20406080

100120140160180200

10

-10

10

-8

10

-6

10

-4

10

-2

10

10

2

算法的收敛曲线 (阶数n=1001)

迭代次数

||r k ||/||b ||

Lanczos MINRES

从图中可以看出，在40步左右Lanczos 发生近似中断现象，然而MINRES 方法的相对误差单调不增，没有出现尖峰。

我们再多实验几组数据。

在前面的特征值分布基础上增加一个负特征值-10得到的收敛曲线如下图

50

100150

200250

10

-10

10

-8

10

-6

10

-4

10

-2

10

10

2

10

4

算法的收敛曲线 (阶数n=1002)

迭代次数

||r k ||/||b ||

Lanczos MINRES

可以看出，MIRES 方法的相对误差仍然具有单调不增的特性。

再将特征值取得更多，再增加10个负特征值，-10：-10：-100，得到的曲线如下

50100

150200250300350

10

-10

10

-8

10

-6

10

-4

10

-2

10

10

2

算法的收敛曲线 (阶数n=1011)

迭代次数

||r k ||/||b ||

Lanczos MINRES

通过对比可以发现，负特征值越多，Lanczos 方法的纹波越大，而MINRES 方法却仍然维持单调不减的特性。

总结：

1) 对于Lanczos 方法，随着负特征值的增多，振荡越来越严重，发生近似中断的次数越来越多； 2) 然而，对于相同的A ，Minres 方法的相对残差没有出现峰值，随着迭代数增加而单调不减。

收敛性态比Lanczos 好。

程序代码

CG 法 CGllb11.m

% CG method for solving Ax=b % By Liu libin @ 2011-11-06

% Email:llb11@https://www.sodocs.net/doc/8912702826.html,

function [x,Error,i,flag]=CGllb11(A,b,x,ErrSet,uplimit)

%deal with input data [m,n]=size(b); if m

[m,n]=size(x); if m

%CG method begins r=b-A*x;

p=r;

%Loop begin

i=1;

temp_rkrkplus=r'*r;

Error=[sqrt(temp_rkrkplus)/norm(b,2)];

while 1

temp_AP=A*p;

temp_rkrk=temp_rkrkplus;

temp_pAP=p'*temp_AP;

if abs(temp_pAP)<1e-12

disp('恶性中断！')

break;

end

a=temp_rkrk/(temp_pAP);

x=x+a*p;

r=r-a*temp_AP;

temp_rkrkplus=r'*r;

beta=temp_rkrkplus/temp_rkrk;

p=r+beta*p;

% decide the error

Err=sqrt(temp_rkrkplus)/norm(b,2); %/(norm(b)+temp_AP);

if Err

disp('Method converge！')

disp(i)

%disp(x)

flag=1;

break;

end

Error=[Error,sqrt(temp_rkrkplus)/norm(b,2)];

if i>uplimit

flag=0;

break

end

end

End

Lanczos法 lanczosllb.m

% Lanczos method for solving Ax=b

% By Liu libin @ 2011-11-15

% Email:llb11@https://www.sodocs.net/doc/8912702826.html,

function [T,Q,x,k,Errs,tiao]=Lanczosllb(A,b,x0,Err)

%deal with input data

x=[];

tiao=[];

[m,n]=size(b);

if m

b=b';

end

[m,n]=size(x0);

if m

m=n;

x0=x0';

end

size(b)

size(A*x0)

r=b-A*x0;

r_zeros=r;

q=r/norm(r);

q0=0;

beta0=0;

T=[0];

k=1;

Q=[q];

Errs=[norm(r,2)/norm(b,2)];

%solve for tri diag matrix T

while 1

%clc

disp('k=');disp(k);

T=[T,zeros(k,1);zeros(1,k),0];

%Add a row and a collom to T

r=A*q-beta0*q0;

T(k,k)=q'*r;

r=r-T(k,k)*q;

T(k+1,k)=norm(r,2);

beta0=T(k+1,k);

T(k,k+1)=beta0;

q0=q;

q=r/beta0;

Tk=T(1:k,1:k) ;

%solve for yk

if k==1

disp('跳过此步')

tiao=[tiao,k];

else

[L,U]=LanczosLU(Tk);

bk=norm(r_zeros,2)*[1;zeros(k-1,1)];

%求Ux

Ux=zeros(k,1);

Ux(1)=bk(1);

for i=2:k

Ux(i)=bk(i)-L(i,i-1)*Ux(i-1);

end

%求ym

ym=zeros(k,1);

ym(k)=Ux(k)/U(k,k);

for i=k-1:-1:1

ym(i)=(Ux(i)-ym(i+1)*U(i,i+1))/U(i,i);

end

Errs=[Errs,abs(ym(k)*beta0)/norm(b,2)];

if abs(ym(k)*beta0)/norm(b,2)

x=x0+Q*ym;

disp('method converge, got the accurate result!') break;

else if beta0<1e-10

x=x0+Q*ym;

disp('Good Luck!')

end

end

%x=x0+Q*ym;

end

if k==1.1*m

%x=x0+Q*ym;

disp('Not Converge!')

break;

end

Q=[Q,q];

k=k+1;

end

end

MINRES法 MINRES.m

% Minres method for solving Ax=b

% By Liu libin @ 2011-11-15

% Email:llb11@https://www.sodocs.net/doc/8912702826.html,

function [T,Q,x,k,Errs]=minresllb(A,b,x0,Err)

%deal with input data

[m,n]=size(b);

if m

b=b';

end

[m,n]=size(x0);

if m

m=n;

x0=x0';

end

r=b-A*x0;

r_zeros=r;

q=r/norm(r);

q0=0;

beta0=0;

T=[0];

k=1;

Q=[q];

Errs=[norm(r_zeros,2)];

%solve for tri diag matrix T

while 1

clc

disp('k=');disp(k);

T=[T,zeros(k,1);zeros(1,k),0];

%Add a row and a collom to T

r=A*q-beta0*q0;

T(k,k)=q'*r;

r=r-T(k,k)*q;

T(k+1,k)=norm(r,2);

beta0=T(k+1,k);

T(k,k+1)=beta0;

q0=q;

q=r/beta0;

Q=[Q,q];

%*********Lanzcos process completed!

%Interruption for Uplimited times

if k<=2

else

%*******start to devide Tk^ by QR method

Rk=T(1:k+1,1:k);

Qk=eye(k+1,k);

for i=1:k

Ci=Rk( i ,i)/((Rk(i,i)^2+Rk(i+1,i)^2)^0.5);

Si=Rk(i+1,i)/((Rk(i,i)^2+Rk(i+1,i)^2)^0.5);

Rk_temp=Rk;

Rk_temp( i ,:)= Ci*Rk(i,:)+Si*Rk(i+1,:);

Rk_temp(i+1,:)=-Si*Rk(i,:)+Ci*Rk(i+1,:);

Rk=Rk_temp;

%compute Qk

Qk_temp=Qk;

Qk_temp( i ,:)= Ci*Qk(i,:)+Si*Qk(i+1,:);

Qk_temp(i+1,:)=-Si*Qk(i,:)+Ci*Qk(i+1,:);

Qk=Qk_temp;

end

%****QR devision completed

Errsk=abs(Qk(k+1,1))*norm(r_zeros,2)/norm(b,2);

Errs=[Errs;Errsk];

if Errsk

gk=Qk(1:k,1)*norm(r_zeros,2);

yk=zeros(k,1);

yk( k )=gk( k )/Rk(k,k);

yk(k-1)=(gk(k-1)-yk(k)*Rk(k-1,k))/Rk(k-1,k-1);

for i=k-2:-1:1

yk(i)=(gk(i)-yk(i+1)*Rk(i,i+1)-yk(i+2)*Rk(i,i+2))/Rk(i,i);

end

if Errsk

disp('Method Converge！Congratulations！')

else

disp('Failed!')

end

x=x0+Q(:,1:k)*yk;

break;

end

end

end

end

其他附属过程函数

A和b的生成矩阵 MakeMatrix.m

clc

clear

%build Matrix A

%lamda

global A

global b

beta=0.5;

%lamda=[1e9, 1000:-1:1, 0.01];%正定的病态问题

lamda=[3100:10:2500, 1000:-1:1,-1,-10:-10:-100];%正定的良态问题

%lamda=[1001, 1000:-1:1, 0.001];%近似奇异问题

%lamda=[1e15, 1000:-1:1,1];%病态正定问题

%lamda=[1001, 1000:-1:1,-1,-10,-100];%病态正定问题

%lamda=[1001, 1000:-1:1, -0.1, -10];%不定的良态问题

%lamda=[ 100:-1:1 , -100:0.5: -1,0 ];%近似不定的良态问题

%lamda=[100*ones(1,210),-10*ones(1,210),1*ones(1,210),1e3*ones(1,210)] %lamda=[1:1100, 1e12,0.001];

% lamda=[1e9,1,2,3,4,5,6,5,8,9,10,11,4,78];

%lamda=[100,50,20,10,1];

n=length(lamda)

semilogx(lamda,ones(n,1),'bo')

Q=orth(randn(n,n));

title('A的特征值分布情况')

xlabel('特征值')

set(gca, 'ytick',[])

b=ones(n,1);

B=mdiagllb(lamda);

A=Q'*B*Q;

disp('矩阵生成完毕！可以进行后续计算！')

disp(A(1:5))

对称三对角矩阵的LU分解 LanczosLU.m

function [L,U]=LanczosLU(A)

[m,n]=size(A);

L=eye(n,n);

U(1,1)=A(1,1);

%L(2,1)=A(2,1)/A(1,1);

U(1,2)=A(1,2);

for k=2:m-1

L(k,k-1)=A(k,k-1)/U(k-1,k-1);

U(k,k)=A(k,k)-L(k,k-1)*U(k-1,k);

U(k,k+1)=A(k,k+1)-L(k,k-1)*U(k-1,k+1);

end

k=m;

L(m,m-1)=A(m,m-1)/U(m-1,m-1);

U(k,k)=A(k,k)-L(k,k-1)*U(k-1,k);

End

实验控制文件（调用这一系列函数） Fortest.m

clc

global A

global b

[m,n]=size(A);

n_Matrix=m;

%调用CG

uplimit=2000;

tic

[y,error,n]=CGllb11(A,b,zeros(n_Matrix,1),1e-8,uplimit); toc

figure(2)

semilogy([1:length(error)],error,'b.-')

norm(A*y-b)

stringsCG=' 算法的收敛曲线 (阶数n=';

stringsCG=[stringsCG,num2str(m),')'];

title(stringsCG)

xlabel('迭代次数')

ylabel('||r_k||/||b||')

grid on

%调用Lanczos

tic

[~,~,x,k,Errs]=Lanczosllb(A,b,zeros(n_Matrix,1),1e-8) ; toc

disp('最大的误差为')

disp(max(A*x-b))

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

第一章 绪论 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;

数值分析实验报告

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

清华大学杨顶辉数值分析第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 === 由拉格朗日插值定理有

数值分析实验报告总结

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

如果 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 =+的经验公式，并求均方误差 解： 法方程为

清华大学杨顶辉数值分析第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 =无根

数值分析实验报告

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

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 ?? ???? ?? ???? ?? ???? = ?? ???? - ?? ???? - ???? ?? 得到如下结果：

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

高等数值分析实验一 工物研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);

数值分析实验报告资料

机电工程学院 机械工程 陈星星 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

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

线性代数方程组的数值解法 实验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=');

清华大学高等数值计算(李津)实践题目二(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 分解为两个低阶二对角阵的奇异值分解问题。

清华大学高等数值分析作业李津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 =

数值计算方法实验报告

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

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

11. 解：计算中保留5位有效数字，第一步，选取作为主元，则 消去，得 第二步，选择 作为主元，则 消去，得 回代计算得到方程的解为 12. （1）证明： 先证明 的对称性，易得 再证的正定性，只要证明的顺序主子式 21311 0.250253.9960.0020.0005005 3.996l l = =-==-(2)(2) 3.9960 5.562547.4178(|)00.61077 1.00100.474700 2.0028 2.00200.40371A b ?? ? =--- ? ??? 320.61077 0.30496 2.0028 l -= =-(3)(3) 3.9960 5.562547.4178(|)0 2.0028 2.00200.40371000.390470.35158A b ?? ? = ? ?--?? 1231.92729,0.69847,0.90040 x x x ==-=(2)(2)11111111 11111111(2)(2)22 ,,2,3,...,,2,3,...,A ===A ij ij i j ji ji j i i j i j ij ji i j j i T ij ji a a l a a a l a i n j n a a a a a a l a l a a a a a A =-=-=====由于对称正定，则 ，则 ，即(2) 22()2()0i ii n nn a A a a ?? ? ? ?→ ? ? ?? ? O O

易得将作Gauss 消去，最终得到 由于这种变换不改变矩阵的行列式，则 由于A 对称正定，则，因此，即的顺序主子式大于零 综上， 对称正定。 （2）证明： 只需证明 由于 则 由于A 严格对角占优，则 则严格对角占优。 13. 解：显然A 对称，，则A 为对称正定矩阵，用平方 根法求得下三角矩阵L 为 由得 ，再由 得 (2)(3)() 22313...,2,3,...,i i ii a a a i n -=?=(2)(2) 2 ||||0,2,3,...,n ii ij j j i a a i n =≠->=∑(2)11ij ij i j a a l a =-(2)(2)11112211112222 11|||||||||||||||||||||| n n ii ij ii i i ij i j j j j i j i n n n n i ii ij i j ii ij j j j j j j i j i a a a l a a l a a a a l a a a a a ==≠≠====≠≠-=---≥ --=--∑∑∑∑∑∑12 1111 2 11(2)(2)21||||||,||||,1 || ||||||||0 n j n n j ii ij j j j j i n n ii ij ii ij j j j i j i a a a a a a a a a a ===≠==≠≠>><->->∑∑∑∑∑则 4L=12233?? ? ? ?-??