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

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

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;

tic

[xA,rmA,flagA]=Arnoldi(A1,b,x0,e,m);

toc

tic

[xG,rmG,flagG]=GMRES(A1,b,x0,e,m);

toc

subplot(1,2,1);

semilogy(rmA)

title('Arnoldiê?á2?ú??')

xlabel('μü′ú2?êyk')

ylabel('log(||rk||)')

subplot(1,2,2);

semilogy(rmG)

title('GMRESê?á2?ú??')

xlabel('μü′ú2?êyk')

ylabel('log(||rk||)')

得到：

得到两种方法的收敛曲线如上所示，将计算结果整理在下表中：

结果讨论：

1.从图中可以看出，基本的Arnoldi方法经过546步收敛，基本的GMRES方法经过526

步收敛，基本的GMRES收敛速度会略快于基本的Arnoldi方法。

2.从图中可以看出，GMRES方法的的性态较Arnoldi方法更好。Arnoldi方法会有平台

和不光滑段，但是由于GMRES具有的残差最优性，GMRES方法平稳地收敛，收敛

曲线也更光滑。

（3）观察重新启动的Arnoldi和GMRES方法

在上述两个函数的基础上，编写了重新启动的Arnoldi函数（详情见附录）：

function [x,rm,flag,Maxi]=ArnoldiM(A,b,x0,tol,m,Maxm)

输入：m为给定步数,Maxm为人为限制的最大重启次数。

输出：x为方程的解，rm为残差向量，flag为解是否收敛的标志，Maxi为重启次数。

首先编写程序，计算重启次数Maxi与给定步数m的关系，为选取给定步数m给出依据：num_restartA=[];

num_restartG=[];

for m=10:10:150

tic

[xG,rmG,flagG,MaxiG]=GMRESM(A,b,x0,tol,m,Maxm);

[xA,rmA,flagA,MaxiA]=ArnoldiM(A,b,x0,tol,m,Maxm);

num_restartA=[num_restartA MaxiA];

num_restartG=[num_restartG MaxiG];

toc

end

m1=10:10:150;

plot(m1,num_restartA,'o-',m1,num_restartG,'*-')

从上述结果中看到在m=50之后，重启次数随着给定步长的增加减少很慢，进入饱和。故可以取m=50，最大步数可知在50步以，运算程序如下所示：

m=50;

Maxm=50;

tic

[xA,rmA,flagA,MaxiA]=ArnoldiM(A,b,x0,tol,m,Maxm);

toc

tic

[xG,rmG,flagG,MaxiG]=GMRESM(A,b,x0,tol,m,Maxm);

toc

m1=1:MaxiA;

m2=1:MaxiG;

plot(m1,log10(rmA),'o-',m2,log10(rmG),'*-')

title('á?????????·¨μ?ê?á2?ú??')

xlabel('????′?êyk')

ylabel('log(||rk||)')

得到的收敛曲线结果如下图所示：

得到两种方法的收敛曲线如上所示，将计算结果整理在下表中：

结果讨论：

1.重启次数随着m的增大而减小，当m增大到一定程度如50之后，减小很缓慢，当

m非常大的时候，就过度到了基本算法。重启的算法如何选择合适的m的因问题

而不同。

2.重启算法的总迭代次数（重启次数×m）要多于基本的算法。这是由于在重启动时，

从另一个我们认为更好的初值点（其实不一定好）x0重新开始计算，可能会丢掉一

些之前算过的信息。

3.重启算法与基本算法相比，虽然总迭代次数增加了，但是运算速度却能大大提高，

运行时间远远小于基本算法（尤其是重启GMRES算法）。

4.一般情况，对于同一个题目，重启的GMRES方法收敛速度要比Arnoldi方法快；

5.总的来说，对题中的矩阵A，四种方法均能稳定地收敛。

T2. 对于1 中的矩阵, 将特征值进行平移, 使得实部有正有负, 和1 的结果进行比较,

方法的收敛速度会如何? 基本的Arnoldi 算法有无峰点? 若有, 基本的GMRES 算法相应地会怎样?

Answer:

对A 的特征值进行平移，如对S 矩阵的实部统一减去一个数s ，则小单元矩阵S 对应的一

对特征值变为()i i s i αβ-±，当矩阵A 构造同上题，那么控制s 的大小，即控制了实部为负的特征值的个数。

i i i i i s S s αββα--??= ?-??

这里将上面的矩阵生成做如下改造：

clear all ; N=500 tol=1e-6; m=50; Maxm=50; A=zeros(2*N); s=1.5

for a=1:N

A(2*a-1,2*a-1)=a-s; A(2*a-1,2*a)=-a; A(2*a,2*a-1)=a; A(2*a,2*a)=a-s; end

U = orth(rand(2*N,2*N)); A1 = U'*A*U;

改变s 的值，进而改变实部为负的特征值个数，计算结果整理如下表所示：

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

第一章 绪论 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 都是对称正定阵。

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

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

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;

数值分析实验报告

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

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

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

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

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

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

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

线性代数方程组的数值解法 实验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 α--===+????=- ?? ?????? ∑∑∑#

数值分析报告 (李庆扬版)

《数值分析》作业 学院：机械学院 专业：机械工程 姓名：赵博 学号：2014520024 日期：2015年6月29日

第二章作业 问：用线性插值及二次插值计算ln0.54的近似值。 答：VB程序如下： Option Explicit Sub czfl(ByRef x() As Single, y() As Single, n As Integer, x1 As Double, f As Double) Dim i, j As Integer Dim p As Single Dim appexcel As Object Dim wbmybook As Object Dim wsmysheet As Object Set appexcel = CreateObject("excel.application") Set wbmybook = appexcel.workbooks.Add Set wsmysheet = appexcel.worksheets.Add f = 0 For i = 0 To n p = 1 For j = 0 To n If i <> j Then p = p * (x1 - x(j)) / (x(i) - x(j)) End If Next j wsmysheet.cells(i + 1, 1) = Str(p) wsmysheet.cells(i + 1, 2) = Str(p * y(i)) f = f + p * y(i) Next i wsmysheet.cells(n + 1, 3) = "最终结果" + Str(f) appexcel.Visible = True End Sub Private Sub Command1_Click(Index As Integer) Dim x() As Single

清华大学数值分析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?? ? ? ?-??

清华大学数值分析报告某实验报告材料

清华大学数值分析报告某实 验报告材料 本页仅作为文档页封面，使用时可以删除 This document is for reference only-rar21year.March

数值分析实验报告 一、 实验 题目： 考虑线性方程组b Ax =，n n R A ?∈，n R b ∈，编制一个能自动选取主元，又能手动选取主元的求解线性代数方程组的Gauss 消去过程。 （1）取矩阵???????????? ????=6816816816 A ，?? ? ?? ??? ????????=1415157 b ，则方程有解()T x 1,,1,1*?=。取10=n 计算矩阵的条件数。分别用顺序Gauss 消元、列主元Gauss 消元和完全选主元Gauss 消元方法求解，结果如何？ （2）现选择程序中手动选取主元的功能，每步消去过程都选取模最小或按模尽可能小的元素作为主元进行消元，观察并记录计算结果，若每步消去过程总选取按模最大的元素作为主元，结果又如何分析实验的结果。 （3）取矩阵阶数n=20或者更大，重复上述实验过程，观察记录并分析不同的问题及消去过程中选择不同的主元时计算结果的差异，说明主元素的选取在消去过程中的作用。 （4）选取其他你感兴趣的问题或者随机生成的矩阵，计算其条件数，重复上述实验，观察记录并分析实验的结果。

1. 算法介绍 首先，分析各种算法消去过程的计算公式， 顺序高斯消去法： 第k 步消去中，设增广矩阵B 中的元素() 0k kk a ≠（若等于零则可以判定系数 矩阵为奇异矩阵，停止计算），则对k 行以下各行计算 ()() ,1,2, ,k ik ik k kk a l i k k n a = =++，分别用ik l -乘以增广矩阵B 的第k 行并加到第 1,2, ,k k n ++行，则可将增广矩阵B 中第k 列中() k kk a 以下的元素消为零；重复此方法，从第1步进行到第n-1步，则可以得到最终的增广矩阵，即 ()()(),n n n B A b ??=? ?； 列主元高斯消去法： 第k 步消去中，在增广矩阵B 中的子方阵()()()()k k kk kn k k nk nn a a a a ??? ?? ????? 中，选取() k k i k a 使得()(k) max k k i k ik k i n a a ≤≤=，当k i k ≠时，对B 中第k 行与第k i 行交换，然后按照和顺序消 去法相同的步骤进行。重复此方法，从第1步进行第n-1步，就可以得到最终的增广矩阵，即( ) ()()111,n n n B A b ??=? ?； 完全主元高斯消去法： 第k 步消去中，在增广矩阵B 中对应的子方阵()()()()k k kk kn k k nk nn a a a a ??? ?? ????? 中，选取()k k k i j a 使得()(k) max k k k i j ij k i n k j n a a ≤≤≤≤=，若k i k ≠或k j k ≠，则对B 中第k 行与第k i 行、第k 列与第 k j 列交换，然后按照和顺序消去法相同的步骤进行即可。重复此方法，从第1 步进行到第n-1步，就可以得到最终的增广矩阵，即() ()()222,n n n B A b ??=? ?；

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

高等数值分析第二次作业 第四题 姓名：---- 学号：------- 时间：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，其具有如下的特征值分布

清华大学数值分析实验报告

数值分析实验报告 一、 实验3.1 题目： 考虑线性程组b Ax =，n n R A ?∈，n R b ∈，编制一个能自动选取主元，又能手动选取主元的求解线性代数程组的Gauss 消去过程。 （1）取矩阵????????????????=6816816816 A ，?? ? ?? ??? ????????=1415157 b ，则程有解()T x 1,,1,1*?=。取10 =n 计算矩阵的条件数。分别用顺序Gauss 消元、列主元Gauss 消元和完全选主元Gauss 消元法求解，结果如？ （2）现选择程序中手动选取主元的功能，每步消去过程都选取模最小或按模尽可能小的元素作为主元进行消元，观察并记录计算结果，若每步消去过程总选取按模最大的元素作为主元，结果又如？分析实验的结果。 （3）取矩阵阶数n=20或者更大，重复上述实验过程，观察记录并分析不同的问题及消去过程中选择不同的主元时计算结果的差异，说明主元素的选取在消去过程中的作用。 （4）选取其他你感兴趣的问题或者随机生成的矩阵，计算其条件数，重复上述实验，观察记录并分析实验的结果。

1. 算法介绍 首先，分析各种算法消去过程的计算公式， 顺序高斯消去法： 第k 步消去中，设增广矩阵B 中的元素() 0k kk a ≠（若等于零则可以判定系数 矩阵为奇异矩阵，停止计算），则对k 行以下各行计算() () ,1,2,,k ik ik k kk a l i k k n a ==++， 分别用ik l -乘以增广矩阵B 的第k 行并加到第1,2, ,k k n ++行， 则可将增广矩阵B 中第k 列中() k kk a 以下的元素消为零；重复此法，从第1步进行到第n-1步，则可以得到最终的增广矩阵，即()()(),n n n B A b ??=? ?； 列主元高斯消去法： 第k 步消去中，在增广矩阵B 中的子阵()()()()k k kk kn k k nk nn a a a a ??? ?? ????? 中，选取() k k i k a 使得()(k) max k k i k ik k i n a a ≤≤=，当k i k ≠时，对B 中第k 行与第k i 行交换，然后按照和顺序消去 法相同的步骤进行。重复此法，从第1步进行第n-1步，就可以得到最终的增广矩阵，即( ) ()()111,n n n B A b ??=? ?； 完全主元高斯消去法： 第k 步消去中，在增广矩阵B 中对应的子阵()()()()k k kk kn k k nk nn a a a a ??? ?? ????? 中，选取()k k k i j a 使得()(k) max k k k i j ij k i n k j n a a ≤≤≤≤=，若k i k ≠或k j k ≠，则对B 中第k 行与第k i 行、第k 列与第k j 列 交换，然后按照和顺序消去法相同的步骤进行即可。重复此法，从第1步进行到

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

第二次实验作业 清华大学 高等数值分析 贾仲孝老师作业 第一题：构造例子特征值全部在右半平面时，观察基本的Arnoldi 方法和GMRES 方法的数值性态，和相应重新启动算法的收敛性。 答： 1、计算初始条件 1) 矩阵A 的生成 根据实Schur 分解，构造矩阵如下形式 11112222 /2/2/2/2n n A n n n n ?-?? ? ? ?- ?= ? ? ? - ? ?? ? 其中，A 由n/2个块形成，每个对角块具有如下形式，对应一对特征向量i i i αβ+ i i i i A α ββα-?? = ??? 、 这里，取n=1000，得到矩阵A 。经过验证，A 的特征值分布均在右半平面，如下图所示 50 100 150 200250300 350 400 450 500 -500-400-300-200-1000100200 300400500 复平面中A 的特征值分布情况 实部 Im(x) 虚部 R e (x ) 特征值 2) b 的初值为 b=(1,1,1,1..1)T 3) 迭代初值为 x 0=0 4) 停机准则为 ε=10-6

2、基本的Arnoldi 和GMRES 方法 代入前面提到的初始值A 、b 、x0，得到的收敛结果如下 100200 300400500600 10-7 10 -6 10 -5 10 -4 10 -3 10 -2 10 -1 10 两种基本算法的||r k ||收敛曲线 (阶数n=1000) 迭代次数 ||r k ||/||b || 基本的Arnoldi 算法基本的GMRES 算法 结果讨论：从图中可以看出，基本的Arnoldi 方法经过554步收敛，基本的GMRES 方法经过535步收敛。这是由于GMRES 具有残差最优性，会略快于Arnoldi 方法，但是，由于两种方法的基本原理近似，GMRES 方法不会实质性的提速。 此外，从收敛曲线上看，由于特征值均处在右半平面，收敛曲线平滑，收敛速度（收敛因子）比较均匀。 3、重启动的GMRES 和Arnoldi 算法 对上述A 、b 、x0使用重启动的Arnoldi 和GMRES 算法。变化m 经过多次实验，得到如下结果： 1. 总迭代步长与m 之间的关系 50 100 150 9001000110012001300140015001600 170018001900两种重启算法总 迭代步长 与m 的关系曲线比较 总迭代次数 m 重启GMRES 重启Aronldi 从表中数据可以看出，结果相差并不大 2. 重启次数与m 之间的关系

清华大学高等数值分析作业李津3——特征值问题

20131112题目1 问：证明若H 为有0特征值的上Hess 阵则一步QR 迭代即可算得H 的0特征值并可收缩（降阶）。 证明：H 为上Hess 阵，进行一次QR 迭代有11H Q R = ，1Q 为正交阵，1R 是上三角矩阵且 1R 对角元中0元均在下方（可以通过given‘s 变换实现）。由于H 有0特征值，故1R 对 角元中至少有一个0元（由矩阵行列式为0即可说明），因而可以通过一步QR 迭代找到H 的0特征值。 设1Q ，1R 有如下形式： () 11=,q n n Q Q - ，110n R R -?? '= ? ??? () 1 1111111q =,q 000n n n n n n n R R Q R H R Q Q -----???? '''== ? ? ? ??? ?? 1H 也是一个上Hess 矩阵，且有如下的形式： 1* ******* **00H ?? ????=? ???????? ? 因此可以收缩降阶，且不需额外操作。

20131114题目1 求：对 210 110.1 00.11 H ?? ?? = ?? ?? ?? ，取 33 1 h=为位移，完成一步隐式位移QR迭代。 解： S1：位移 110 100.1 00.10 H I ?? ??-= ?? ?? ?? S2：利用given’s 变换作QR分解 220 22 110 0100.10 22220 00.10 00100.10 ?? ?? ? ???? ? ???? -=- ? ???? ?? ???? ???? ?? ?? ?? ???? 100 000 220 00.10 00 ? ?? ? ?? ??? ? ??? ? -= ?? ? ??? ? ???? ?? ???? ?? ?????? 00 100100 2222 0000 2222 001001 00 2 T T Q ??? ???? - ??? ??????? ?? ??? ?? =-= ??? ?? ???? ?? ???? ???????? ??? ??? ???? ? ? =? ?? ?? ?? ?? ??