清华大学贾仲孝老师高等数值分析第二次实验
清华大学贾仲孝老师(贾哥)高等数值分析证明题汇总
又因为A可对角化， A= -1 则p(A)=p( ) -1 则 rk
n
pPk ,p (0)1
min
p( ) -1r0
取r0 i xi Xy, 其中y =(1, 2 , 3 k-1, k ,0,0 0,0)T
i 1
则 rk
pPk ,p (0)1
i 2
n
i 2 i 2 k 1 (vk , Avk ) i 1 1 1 n 1 2 2k (vk , vk ) i i
i 1
n
1212 k 1 2 2 2 2 k 1 i 2 i 2 k 1
2 2k 1 1
又因为A可对角化， A= -1 则p(A)=p( ) -1 则 rk
pPk , p (0)1
min
p( ) -1r0 ( ) r0
pPk ,p (0) 1 1 i n
pPk , p (0)1
min
p( )
即 rk ( ) r0 可取pk (t )
n
2 2 i 2 i 2 k
2 2k i 3
i 3 n
1
2 2 k ) 2 …… 1 1 12 2 2 ( 2 ) 2 k …… 1 2 2 ( 2 ) 2 k 2 2 2 ( 2 ) 2 k 1 1 1 =o ( 2 ) 2 k 2 1 1 121 2 2 (
证毕 15． 叙述幂法 解，见讲义 P80。
5 / 6
结语： 最后三道题为四个月后补录，准确性已经无法保证，仅供参考，祝君顺利。
6 / 6
12 1 C2T X 2T A X 2C 2 1 (12 1)1 C2T ( X 2T A X 2 )C 2 C2 1 X 2 C2 =o(e 2 k ) （ o C2 ）
清华大学贾仲孝老师高等数值分析报告第二次实验
高等数值分析第二次实验作业T1.构造例子特征值全部在右半平面时, 观察基本的Arnoldi 方法和GMRES 方法的数值性态, 和相应重新启动算法的收敛性. Answer:（1） 构造特征值均在右半平面的矩阵A ：根据实Schur 分解，构造对角矩阵D 由n 个块形成，每个对角块具有如下形式，对应一对特征值i i i αβ±ii i i i S αββα-⎛⎫= ⎪⎝⎭这样D=diag(S 1,S 2,S 3……S n )矩阵的特征值均分布在右半平面。
生成矩阵A=U TAU ，其中U 为正交阵，则A 矩阵的特征值也均在右半平面。
不妨构造A 如下所示：2211112222/2/2/2/2N NA 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:NA(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; endU = 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);toctic[xG,rmG,flagG]=GMRES(A1,b,x0,e,m);tocsubplot(1,2,1);semilogy(rmA)title('ArnoldiÊÕÁ²ÇúÏß')xlabel('µü´ú²½Êýk')ylabel('log(||rk||)')subplot(1,2,2);semilogy(rmG)title('GMRESÊÕÁ²ÇúÏß')xlabel('µü´ú²½Êýk')ylabel('log(||rk||)')得到：得到两种方法的收敛曲线如上所示，将计算结果整理在下表中：结果讨论：1.从图中可以看出，基本的Arnoldi方法经过546步收敛，基本的GMRES方法经过526步收敛，基本的GMRES收敛速度会略快于基本的Arnoldi方法。
清华大学 李津老师 数值分析第二次实验作业
就不再赘述了。 二、实际计算 生成十个不同的(最好属不同类型或有不同性质的)的 m n 矩阵，这里 m, n 100 ， 用你选择的算法对其做 SVD,比较不同方法的效果(比如计算小气一直和对应左右奇异向量的 误差，效率等)，计算时间和所需存储量等，根据结果提出对算法的认识。 1.误差 在实验中，我们取 m=200，n=100，利用 orth()函数生成了正交矩阵������、������，再生成 了不同奇异值分布的奇异值矩阵������，再通过������ = ������������������，计算出不同的待分解矩阵。 各矩阵奇异值分不如下表所示 序号 1 2 3 4 5 6 7 8 9 10 奇异值个数 10 20 30 40 50 60 70 80 90 100 奇异值分布 10 → 1 20 → 1 30 → 1 40 → 1 50 → 1 60 → 1 70 → 1 80 → 1 90 → 1 100 → 1
−1 −1 −1 −1 −1 −1 −1 −1 −1 −1
经过 matlab 计算，我们得到了两种算法对奇异值的估计误差表，如下所示
序号 svd 1 2 3 4 5 1.2135e-29 1.0336e-28 3.9976e-28 7.9768e-28 1.5711e-27
i 1
100
i
i
2.5232e-27 4.6720e-27 5.8535e-27 8.7958e-27 9.8885e-27
r
i 1
ui uis r
2 2
i 1
vi vis r
2 2
lansvd
svd
lansvd
2 1.6 2.1333 2 2.16 1.8 2.3429 2.15 2 2
贾忠孝高等数值分析第一次上机作业
【运行结果】 当 m=300 时，下图表示实验结果，其中 k 为停机循环次数
当 m=995 时，停机循环次果表明当 m 越小时，收敛速度越快，在有限精度下，停机循环次数越小。这 与 Lanczos 方法最多 m 步找到准确解这一结论相吻合。
4. 取初始近似解为零向量, 右端项 b 仅由 A 的 m 个不同个特征向量的线性组 合表示时, Lanczos 方法的收敛性如何? 数值计算中方法的收敛性和m 的大小关 系如何? 【程序】 clc clear clf n=1000; rs=10^(-10); %停机准则 m=10; disp('Lanczos 算法') %定义对称正定阵与 b
1 1000, 2 59.8858, 999 10.0961, 1000 0.01,
经过 k=100，基本收敛达到精度要求。
【结论】 后者的收敛速度要大于前者。说明当 A 的最大特征值远大于第二个最大特征 值, 最小特征值远小于第二个最小特征值时，收敛速度只与 2. 对于同样的例子, 比较 CG 和 Lanczos 的计算结果. 【程序】
高等数值分析实验作业一
机械系 光辉 2014310205
1. 构造例子说明 CG 的数值性态 . 当步数 = 阶数时 CG 的解如何 ? 当 A 的最 大特征值远大于第二个最大特征值, 最小特征值远小于第二个最小特征值时 方 法的收敛性如何? 问题 1 【程序】 clc clear n=1000; %定义对称正定阵与 b lans=sort(rand(n,1),'descend'); D=diag(lans); Q=orth(rand(1000,1000)); A=Q*D*Q'; b=rand(n,1); %初始情况 x0=eye(n,1); r0=b-A*x0; p(:,1)=r0; %执行第一步 af(1)=dot(r0,r0)/dot(p(:,1),A*p(:,1)); x(:,1)=x0+af(1)*p(:,1); r(:,1)=b-A*x(:,1); bat(2)=dot(r(:,1),r(:,1))/dot(r0,r0); p(:,2)=r(:,1)+bat(2)*p(:,1); %执行 k 步 K=n; for k=2:1:K k af(k)=dot(r(:,k-1),r(:,k-1))/dot(A*p(:,k),p(:,k)); x(:,k)=x(:,k-1)+af(k)*p(:,k); r(:,k)=b-A*x(:,k);
清华大学高等数值分析实验设计及答案
高等数值分析实验一工物研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步时已经收敛到接近真实xrelres= 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，由于上实验结果可以看出实际叠代次数都比上限值要小较多。
清华大学数值分析实验报告
数值分析实验报告一、 实验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 ikik k kka l i k k n a ==++,分别用ik l -乘以增广矩阵B 的第k 行并加到第1,2,,k k n ++行，则可将增广矩阵B 中第k 列中()k kka 以下的元素消为零；重复此方法，从第1步进行到第n-1步,则可以得到最终的增广矩阵，即()()(),n n n B Ab ⎡⎤=⎣⎦； 列主元高斯消去法：第k 步消去中，在增广矩阵B 中的子方阵()()()()k kkkknk k nknn a a a a ⎡⎤⎢⎥⎢⎥⎢⎥⎣⎦中,选取()k k i k a 使得()(k)max k k i k ik k i na a ≤≤=，当k i k ≠时，对B 中第k 行与第k i 行交换,然后按照和顺序消去法相同的步骤进行。
高等数值分析第二章答案
第二章习题参考答案1.解： 由于20Ax b−≥，极小化2b Ax −与极小化22Ax b −是等价的。
令22()(,)(,)2(,)x Ax b Ax Ax b b Ax b ϕ=−=+−，对于任意的n R y x ∈,和实数α，)()(),()()(,*222*2****x Ay a x Ay Ay a x ay x b Ax x ϕϕϕϕ≥+=+=+=则有满足若这表示处达到极小值。
在*)(x x ϕ反之，若必有处达到极小，则对任意在nR y x ay x ∈+*)(ϕ0),(2),(2),(20)(**0*=−=+−=+=Ay b Ax Ay Ay a Ay b Ax daay x d a 即ϕ故有 b Ax =*成立。
以上证明了求解,22b Ax b Ax −=等价于极小化即。
等价于极小化2b Ax b Ax −= 推导最速下降法过程如下：),/(),(0),(),(,0),,2)(222)()(11k T k T k T k k T k T k T k k T k k k T k k kT k T k T T x x k r AA r AA r AA r a r AA r AA a r AA r r aA x da dx a r aA x x r A Ax b A Ax A b A x grad x x k==+−=++==−=−=−++=最终得到得出（由取得极小值。
使求出取的负梯度方向，且下降最快的方向是该点在ϕϕϕ给出的算法如下：1）)(000Ax b A r A R x T T n −=∈，计算给定； 2）L ,2,1,0=k 对于）转到否则数。
为一事先给定的停机常则停止；其中若2),/(),(10,11kT k k k k T k k k k k k k k k r A p Ax b r r A a x x Ap Ap p p a k k r =−=+==+=>≤−−εε2.证明 1） 正定性由对称正定矩阵的性质，(),0x Ax ≥（当且仅当x ＝0时取等号），所以 ()12,0Axx Ax =≥（当且仅当x ＝0时取等号）2） 齐次性()()()121122,(),,AA xx A x x Ax x Ax x αααααα⎡⎤====⎣⎦3）o1方法（一）A 是对称正定矩阵，得到(,())0x y A x y λλ++≥，把它展开如下2(,)(,)(,)(,)0y Ay x Ay y Ax x Ax λλλ+++≥考虑到(,)(,)(,)x Ay Ax y y Ax ==，把上式看成关于λ的一元二次方程，则式子等价于24(,)4(,)(,)0x Ay x Ax y Ay ∆=−≤因此1/21/2(,)(,)(,)x Ay x Ax y Ay ≤所以1/21/221/21/2((,)(,))(,)(,)2(,)(,)(,)(,)2(,)(,)(,)(,)(,)((),())x Ax y Ay x Ax y Ay x Ax y Ay x Ax y Ay x Ay x Ax y Ay x Ay y Ax x y A x y +=++≥++=+++=++两边开平方即可得到AA A x yx y +≤+因此，1/2(,)A x Ax x =是一种向量范数。
清华大学高等数值分析_第三次作业答案
高等数值分析第三章作业参考答案1.考虑线性方程组Ax=b,其中A是对称正定矩阵.用Galerkin原理求解方程K=L=Span(v),这里v是一个固定的向量.e0=x∗−x0,e1=x∗−x1证明(e1,Ae1)=(e0,Ae0)−(r,v)2/(Av,v),(∗)其中r=b−Ax0.v应当取哪个向量在某种意义上是最佳的?证明.令x1=x0+αv,那么r1=r−αAv,e1=e0−αv.由Galerkin原理,有(r1,v)=0,因此α=(r,v)/(Av,v).注意到r1=Ae1,r=Ae,有(Ae1,v)=0.于是(e1,Ae1)=(e0−αv,Ae1)=(e0,Ae1)=(e0,Ae0)−α(e0,Av)=(e0,Ae0)−α(r,v)即(∗)式成立.由(∗)式知当v=e0时, e1 A=0最小,即近似解与精确解的误差在A范数意义下最小,算法一步收敛(但是实际中这个v不能精确找到);在最速下降意义下v=r时最佳.2.求证：考虑线性方程组Ax=b,其中A是对称正定矩阵.取K=L=Span(r,Ar).用Galerkin方法求解,其中r是上一步的残余向量.(a)用r和满足(r,Ap)=0的p向量构成K中的一组基.给出计算p的公式.解.设p=r+αAr,(r,Ap)=0等价于(Ar,p)=0.解得α=−(Ar,r)/(Ar,Ar).(b)写出从x0到x1的计算公式.解.设x1=x0+β1r+β2p,那么r1=r−β1Ar−β2Ap,再由Galerkin原理,有(r1,r)=(r1,p)=0,解得β1=(r,r)/(Ar,r),β2=(r,p)/(Ap,p).(c)该算法收敛吗?解.该算法可描述为：(1)选初始x0∈R n,计算初始残差r0=b−Ax0,ε>0为停机准则；(2)对k=0,1,2,...直到 r k <εαk=−(r k,Ar k) (Ar k,Ar k);p k=r k+αAr k;βk=(r k,r k) (Ar k,r k);γk=(r k,p k) (Ap k,p k);r k+1=r k−βk Ar k−γk Ap k;x k+1=x k+βk r k+γk p k.此算法本质上是由CG迭代一步就重启得到的,所以是收敛的,下面给出证法.设用此算法得到的x k+1=x k+¯p1(A)r k,那么e k+1 A=minp1∈P1e k+p1(A)r k A≤ e k+¯p1(A)r k A= e k−¯p1(A)Ae k A≤max1≤i≤n|˜p(λi)| e k A其中0<λ1≤...≤λn为A的特征值,˜p(t)=1−t¯p1(t)是过(0,1)点的二次多项式.当˜p满足˜p(λ1)=˜p(λn)=−˜p(λ1+λn2)时可使max1≤i≤n|˜p(λi)|达到最小.经计算可得min ˜p max1≤i≤n|˜p(λi)|≤(λ1−λn)2(λ1−λn)2+8λ1λn<1故若令κ=λ1/λn,则e k+1 A≤(κ−1)2κ2+6κ+1e k A,方法收敛.3.考虑方程组D1−F−E−D2x1x2=b1b2,其中D1,D2是m×m的非奇异矩阵.取L1=K1=Span{e1,e2,···,e m},L2= K2=Span{e m+1,e m+2,···,e n}.依次用(L1,K1),(L2,K2)按讲义46和47页公式Az∗=r0r0−Az m⊥LW T AV y m=W T r0x m=x0+V(W T AV)−1W T r0各进行一步计算.写出一个程序不断按这个方法计算下去,并验证算法收敛性.用L i=AK i重复上述各步骤.解.对任意给定x0=x(0)1x(0)2,令r=b−Ax0,V1=[e1,e2,...,e m],V2=[e m,e m+1,...,e n].对L i=K i情形,依次用(L1,K1),(L2,K2)各进行一步计算：(L1,K1)(L2,K2)z(1) 1=V1y1z(2)1=V2y2r0−Az(1)1⊥L1r0−Az(2)1⊥L2(V T1AV1)y1=V T1r0,D1y1=V T1r0(V T2AV2)y2=V T2r0,−D2y2=V T2r0x(1)1=x(1)+V1D−11V T1r0x(2)1=x(2)−V2D−12V T2r0得如下算法：(1)选初始x0∈R n,计算初始残差r0=b−Ax0,ε>0为停机准则；(2)对k=1,2,...直到 r k <ε求解D1y1=r k−1(1:m);求解−D2y2=r k−1(m+1:n);x k=x k−1+V1y1+V2y2;r k=r k−1−AV1y1−AV2y2.收敛性:r k=r k−1−AD−11−D−12rk−1=0−F D−12ED−11rk−1Br k−1算法收敛⇔ρ(B)<1⇔ρ(ED−11F D−12)<1.对L i=AK i情形,依次用(L1,K1),(L2,K2)各进行一步计算:(L1,K1)(L2,K2)z(1) 1=V1y1∈K1z(2)1=V2y2∈K2r0−Az(1)1⊥L1=AK1r0−Az(2)1⊥L2=AK2(V T1A T AV1)y1=V T1A T r0(V T2A T AV2)y2=V T2A T r0(D T1D1+E T E)y1=V T1A T r0(D T2D2+F T F)y2=V T2A T r0x(1) 1=x(1)+(D T1D1+E T E)−1V T1A T r0x(2)1=x(2)+(D T2D2+F T F)−1V T2A T r0得如下算法：(1)选初始x0∈R n,计算初始残差r0=b−Ax0,ε>0为停机准则；(2)对k=1,2,...直到 r k <ε求解(D T1D1+E T E)y1=(A T r k−1)(1:m);求解(D T2D2+F T F)y2=(A T r k−1)(m+1:n);x k=x k−1+V1y1+V2y2;r k=r k−1−AV1y1−AV2y2.收敛性:r k=r k−1−A(D T1D1+E T E)−1(D T2D2+F T F)−1A T rk−1(I−B)r k−1算法收敛⇔ρ(I−B)<1⇔0<λ(B)<2.4.令A=3−2−13−2...............−2−13,b=1...2用Galerkin原理求解Ax=b.取x0=0,V m=W m=(e1,e2,···,e m).对不同的m,观察 b−Ax m 和 x m−x∗ 的变化,其中x∗为方程的精确解.解.对于 b−Ax m 和 x m−x∗ ,都是前n−1步下降趋势微乎其微,到第n步突然收敛。
《数值分析简明教程》第二版（王能超编著）课后习题答案高等教育出版社
《数值分析简明教程》第⼆版（王能超编著）课后习题答案⾼等教育出版社0.1算法1、（p.11，题1）⽤⼆分法求⽅程013=--x x 在[1,2]内的近似根，要求误差不超过10-3.【解】由⼆分法的误差估计式311*10212||-++=≤=-≤-εk k k a b x x ，得到100021≥+k .两端取⾃然对数得96.812ln 10ln 3≈-≥k ，因此取9=k ，即⾄少需2、（p.11，题2）证明⽅程210)(-+=x e x f x在区间[0,1]内有唯⼀个实根；使⽤⼆分法求这⼀实根，要求误差不超过21021-?。
【解】由于210)(-+=x e x f x ，则)(x f 在区间[0,1]上连续，且012010)0(0<-=-?+=e f ，082110)1(1>+=-?+=e e f ，即0)1()0(⼜010)('>+=x e x f ，即)(x f 在区间[0,1]上是单调的，故)(x f 在区间[0,1]内有唯⼀实根.由⼆分法的误差估计式211*1021212||-++?=≤=-≤-εk k k a b x x ，得到1002≥k .两端取⾃然对数得6438.63219.322ln 10ln 2=?≈≥k ，因此取7=k ，即⾄少需⼆分0.2误差1．（p.12，题8）已知e=2.71828…，试问其近似值7.21=x ，71.22=x ，x 2=2.71，718.23=x 各有⼏位有效数字？并给出它们的相对误差限。
【解】有效数字：因为11102105.001828.0||-?=<=- x e ，所以7.21=x 有两位有效数字；因为12102105.000828.0||-?=<=- x e ，所以71.22=x 亦有两位有效数字；因为3310210005.000028.0||-?=<=- x e ，所以718.23=x 有四位有效数字；%85.17.205.0||111=<-=x x e r ε； %85.171.205.0||222=<-=x x e r ε； %0184.0718.20005.0||333=<-=x x e r ε。
哈工大数值分析实验报告
产生逼近解 x*的迭代数列{xk}，这就是割线法的思想
.1.
数值分析实验报告
宋俊霖
拟 Newton 法以 X 0 为初始 Nhomakorabea似，利用递推关系
X k 1 X k H k F ( X k ) H k 1 ( F ( X k 1 ) F ( X k )) X k 1 X k H k 1 H k H k , k 0,1,...
其中 r 为要求的方程的根的重数，这就是改进的 Newton 法，当求解已知重数的方程的根 时，在同种条件下其收敛速度要比 Newton 法快的多。 割线法通常预先给出两个猜测初值 x0 , x1 ，然后根据迭代公式
xk 1 xk f ( xk ) ( xk xk 1 ) f ( xk ) f ( xk 1 )
产生近似于方程组 F ( X ) 0 的解 X * 的迭代序列 { X k } ， 利用这个递推关系式就是拟 Newton 法， 实际计算时， 只要选择较好的初始近似 X 0 和初始矩阵 H 0 ， 一般可得到较好的近似解。
1.2
Python 程序设计
二分法源程序： import math def f(x): return math.sin(x)-0.5*x**2 def binary_method(f, x0, x1, eps, maxi): delta=100 k=0 while delta > eps: x2 = (x1+x0)/2 if f(x0)*f(x2)>0: x0=x2 else : x1=x2 delta=abs(x1-x0) k+=1 print('Root is at: ', x2) print('f(x) at root is: ', f(x2)) x0=1.0 x1=2.0 binary_method(f, x0, x1, 1e-5,100)
数学实验第二次测验题及参考答案（09级）
数学实验第二次测验题及参考答案（09级）数学实验第二次测验题及参考答案一、写出下列MATLAB指令的运算结果.1. A=[1;2;3]; transpose(A)1 2 31 2 32. A=[1, 2, 3 ; 4, 5, 6 ]; B=A([1 2], [1,3]) , d=size(A)B =1 34 6d =2 33. a=1:3; b=linspace(1,3,3); x=sum(a.*b), y=cross(a, b)x = 14y = 0 0 04. A=[1,2, 3; 4,5,6; 7,8,9]; B=ones(3); C=A-BC =0 1 23 4 56 7 85. v=[1, 2, 3]; A=diag(v); E=eig(A), D=det(A)E=123D =66. x=[1,2,3,4,5]; [mean(x), median(x), range(x),sum(x), prod(x)]3 34 15 1207. x=[2,3,4]; a=cumsum(x) ,b=sort(x)a =2 5 9b =2 3 48．format rat; v=[1, 2, 3]; A=diag(v); inv(A)ans =1 0 00 1/2 00 0 1/39. [m,v]=normstat(1,4) % 求参数为1,4的正态分布的均值与方差m =1, v =16二、写出下列MATLAB指令的实验目的.1. dsolve('x*Dy+y-exp(-x)=0', 'y(1)=2*exp(1)', 'x')求微分方程0=-+'-x e y y x 在初始条件e y x 2|1==下的特解.2. u=[1,2,3],v=[0,3,2], w=[5, 2, 1]; dot(w, cross(u, v))计算向量u, v, w 的混合积.3. A=[1 2 3; 2 2 5; 3 5 1]; b=[1;2;3]; det(A); inv(A)*b利用逆矩阵解线性方程组=++=++=++3532522132321 321321x x x x x x x x x .4. A=[0 0 1; 0 1 1; 1 1 1; 1 0 0]; rref(A)求向量组)1,0,0(1=α,)1,1,0(2=α,)1,1,1(3=α,)0,0,1(4=α的秩.或对矩阵A 做行初等变换。
高等数值分析第二次实验作业(贾忠孝)
3. 对 1 中的例子固定特征值的实部, 变化虚部, 比较收敛性. A 矩阵每个小块有如下形式：
A i i
其特征值为
i i
i ii ，对其次对角元素 i 做如下处理：
i ' k i
其特征值变为 果：
i ik i 。在本题中分别设定 k 的取值为 0.5，1，2，4。可以得到下面的结
结果分析： （1） 当 k 越大，既即特征根虚部越大时，基本 Arnoldi 方法和 GMRES 方法的收敛速度会 越来越慢，Arnoldi 方法的振荡程度会随着 k 值的增大而增大，振荡范围也增大。 当 k 大到 10 或者 10 以上时，两法都只在最后一步（第 1000 步）收敛，而且都是 急剧收敛。 （2） GMRES 法的残量一直比 Arnoldi 法的残量小，且相对残差的曲线平滑。.
GMRES算 法 计 算 m个 特 征 向 量 组 成 的 b的 收 敛 曲 线 4 2 0 -2 GMRES, GMRES, GMRES, GMRES, m m m m = = = = 10 50 100 1000
log(||rk||)
-4 -6 -8 -10 -12 -14
0
50 Iteration Times
100
150
结果分析： （1） 当 m 较小时，收敛速度很快，但其随着 m 值的增大而慢慢变小。但是一直到最后 一步之前，收敛的幅度都不大，往往是在最后一步急剧收敛； （2） 当 m 较大时，收敛的速度随 m 的增大而微弱减小，而且其在步数比较小的时候， 收敛速度较快，越往后走反而慢。 （3） Arnoldi 法比 GMRES 法收敛曲线的趋势是基本一样的，而且 Arnoldi 法收敛比 GMRES 法要快，相对残量： 1.构造例子特征值全部在右半平面时,观察基本的 Arnoldi 方法和 GMRES 方法的数值性态,和 相应重新启动算法的收敛性. 解：构造 1000 阶符合条件的矩阵 A。根据实 Schur 分解，构造如下形式的矩阵：
数值分析--清华李庆杨五版第二章_插值法
xi
(i=0,1,…,n )
处与 f ( xi ) 相等,在其它点x就用 (x) 的值作为f(x) 的近似值。这一过程称为插值，点x称为插值点。换 句话说, 插值就是根据被插函数给出的函数表“插出” 所要点的函数值。用 (x) 的值作为f(x)的近似值,不仅希
望 (x) 能较好地逼近f(x),而且还希望它计算简单 。由
l0 ( x) c( x x1 )( x x2 )
类似地可以构造出满足条件： l1 ( x1 ) 1, l1 ( x0 ) 0 , l1 ( x2 ) 0 的插值多项式
( x x0 )( x x 2 ) l1 ( x) ( x1 x0 )( x1 x 2 )
a n x0 n a n 1 x 0 n 1 a1 x0 a 0 f ( x0 ) n n 1 a n x1 a n 1 x1 a1 x1 a 0 f ( x1 ) a x n a x n 1 a x a f ( x ) n 1 n 1 n 0 n n n
j 0 j k n
于是
Ak
1
(x
j 0 j k
n
k
xj)
代入上式,得
(x x
l k ( x)
j 0 jk n
n
j
)
j 0 jk n
x xj xk x j
(x
j 0 jk
k
xj)
称
l k (x) 为关于基点
x i 的n次插值基函数(i=0,1,…,n)
l k ( xi ) ki
1 0
(i k ) (i k )
l 0 ( x) 与 l1 ( x) 称为线性插值基函数。且有
清华大学高等数值分析 第一次实验作业
高等数值分析实验作业一
附件：主要算法代码
CG 法 CG.m function [x,Error,i,flag]=CG(A,b,x,ErrSet,uplimit) [m,n]=size(b); if m<n
b=b'; end [m,n]=size(x); if m<n
x=x'; end r=b-A*x; p=r; i=1; temp_rkrkplus=r'*r; Error=sqrt(temp_rkrkplus)/norm(b,2); while 1
Lanzcos算法的收敛曲线 (阶数n=1002)
100 m=10 m=50 m=100 m=400 m=800
10-5
||rk||/||b||
-10
10
-15
10 0
20
40
60
80
100
120
140
160
迭代次数
图5 对不同的m，Lanczos法求解Ax=b的收敛曲线
高等数值分析实验作业一
结论：如果 A 只有 m 个不同的特征值，则 Lanczos 方法至多 m 步就可以找到精 确解。实验中，在 m 较大的时候，算法收敛较快，迭代次数远小于 m。当 m 较 小时，可能需要接近于 m 步才能找到准确解。
2. 对于同样的例子，比较 CG 和 Lanczos 的计算结果。 解：（1）构造 1002 阶正定矩阵 A：
（2）Lanczos法求解Ax=b，A良态： 利用matlab编程实现CG算法。b = ones(1002,1)，x0 = zeros(1002, 1)。计算 每一步迭代的残差rk相对于初始残差的2范数。相对残差2范数的对数值与 迭代步数的关系曲线如图3所示：
清华大学贾仲孝高等数值分析第二次作业
1. 解：由于求解Ax b =等价于极小化2Ax b -，相当于极小化泛函：()()1,2x Ax b Ax b ϕ=-- 从任一k x 出发，沿着()x ϕ在k x 点下降最快的方向搜索下一个近似点1k x +，使得()1k x ϕ+在该方向上达到极小值，最速下降方向为：()()T T k k x A Ax b A r ϕ-∇=--=令1,T k k k k k k p A r x x p α+==+，需要寻找合适的k α使得()()11min k k Rx x αϕϕ++∈=，则()()()()()()1,=T T T T k k k k k k k k d x x p p A A x p b p A p r d ϕϕααα+=∇=+-- 令()0d x d ϕα=，可得： ()()()(),,,,k k k k k k k k k Ap r p p Ap Ap Ap Ap α==则()22,0k k d Ap Ap d ϕα=>，因此上式的k α即为所求 于是得到极小化泛函()()1,2x Ax b Ax b ϕ=--的最速下降法： 1) 选取初值0n x R ∈，计算00r b Ax =-2) k=0,1,2,……若k r ε≤，则停止；ε为事先给定的停机场数；否则，k=k+1()()11;,;,;;T k k k k k k k k k k k k k k k p A r p p Ap Ap x x p r r Ap ααα++===+=-2. 解：()()111k k k k AAAx x x r x p A x x α***----=+-=-其中()1p A A α=-，设A 的特征根为120n λλλ≥≥≥> ，则有()11max k i k AAi nx x p x x λ**-≤≤-≤-当120αλ<<时，()1max 1i i np λ≤≤<，因此此方法收敛当()111n αλαλ-=--即12n αλλ=+时，()1max i i n p λ≤≤取最小值11n nλλλλ-+，此时收敛最快3. 解：设x *为方程组Ax b =精确解，k k e x x *=-，则()1,TT Tk k k E e e -=，则原迭代法等价于()110k k I A I E E I βαβ+⎡++-⎤=⎢⎥⎣⎦令()10I A I B I βαβ⎡++-⎤=⎢⎥⎣⎦，则迭代法收敛等价于()1B ρ<，即()1,1i B i nλ<≤≤，令0B I λ-=，即 ()10n I A ββλαλλ⎛⎫+--+-= ⎪⎝⎭ 当0λ≠时，则有10I A ββλαλ⎛⎫+--+= ⎪⎝⎭设120n μμμ≥≥≥> 是A 的特征根，则有()210101112i i i ββλαμλλβαμλβλβαμβ+--+=-+++=<⇔++<+<则有：()()()11112,1211,0,1211,0i i B i ni n ρβαμβββαμββαμ<⇔++<+<≤≤+⇔<-<<≤≤+⇔<-<< 4.5. 证明：反证法。
贾忠孝历年考题
2012高等数值分析（贾仲孝老师）
首先庆贺自己研究生的第一门也是最后一门考试顺利结束，哦也，再也没有考试啦！
贾哥延续了考题沿用往年题的优良传统，详情见拍照，全部是历年考题，基本上是按着去年考题来的，只有一道题第二个小问号改了一下而已。
1. 用householder和givens变换做QR分解，由于矩阵特殊，非常简单。
如果拿不准，不妨用GS方法做一遍验证一下，因为不同的QR分解只是符号有差异而已，GS还是比householder简单很多的。
2.1证明rayleigh商最大值等于A的最大特征值，将x拆成各个特征向量之和就容易证。
2.2幂法求一个特征值，一步收敛
3. 第三次作业第三题，也是去年的原题，基本上都不用想，直接默写就行了
4.1 去年考题。
注意到Ak*R（k-1）=R（k-1）*A（k-1），那么就类似冒泡算法把Ak移到最右边变成一个A
4.2 有点小恶心，去年这个问号问的是Qk第一列是特征向量x1，只需要两边同乘以e1即可，但今年问的是最后一行是特征向量xn，顿时就不会证了，当时打眼一看觉得貌似A也不能说明是可逆的，就没往反幂法这个方向去，但是后来想想其实最多也就是lamdaN=0，其他不为0，也许可以分情况讨论下？后来有同学说假设可逆用了反幂法也没得到什么结论...不知道...还好就是5分，丢了就丢了吧
5 考过多年的经典背诵题，默写rayleigh ritz方法和贾哥定理，以及Arnolid和精化Arnolid算法。
清华数学科学系推免指导教师名单
代数学、实变函数与泛函分析、偏微分方程与常微分 方程、数值分析（四选二)
算法方法、数值分析、实变函数与泛函分析、拓扑与 几何、数理统计（五选二）
计算数学及几何 7 数学
图像
孙剑
算法方法、数值分析、实变函数与泛函分析、拓扑与 几何、数理统计（五选二）
盛世清北—专注清华大学考研|保研|考博辅导
3 数学 基础数学
潘日新
代数学、实变函数与泛函分析、拓扑与几何、偏微分 方程与常微分方程、复分析 （五选二)
4 数学
基础数学
代数学、实变函数与泛函分析、拓扑与几何、偏微分 A.K.M.Masood-UL-Alam
方程与常微分方程、复分析 （五选二)
5 数学 应用数学 计算数学及几何
6 数学 图像
江宁 丘成桐
计算数学及几何 算法方法、数值分析、实变函数与泛函分析、拓扑与几何、数理统计（五选
3 数学
图像
二）
盛世清北—专注清华大学考研|保研|考博辅导
代数学
代数学
几何学
实变函数与泛函分析
偏微分方程
代数学
几何学
实变函数与泛函分析
偏微分方程
代数学
几何学
实变函数与泛函分析
清华大学高等数值分析 第二次实验作业
400
600
800
1000
1200
1400
1600
1800
2000
迭代次数
图4
m步的重启动GMRES法求解Ax=b的收敛曲线
结论： m步重启GMRES方法快于m步重启Arnoldi方法， 随m增加， 迭代次数减小， 但都大于未重启算法的次数。当m=40时两种方法计算时间最短，此外，m步重 启动 Arnoldi 方法的收敛曲线有峰点和波纹，收敛速度均匀性差， m 步重启动 GMRES方法的收敛曲线很平滑，单调下降，收敛速度均匀性好。（图4是五条曲 线， 只是由于m=20和m=80两条曲线比较靠近， 看起来像四条， 放大后才能看清） 2.对于 1 中的矩阵，将特征值进行平移，使得实部有正有负，和 1 的结果进行比 较，方法的收敛速度会如何？基本的 Arnoldi 算法有无峰点？若有，基本的 GMRES 算法相应地会怎样？ 解： （1）欲将特征值进行平移，使得实部有正有负，可以将矩阵 A 做如下变换：
10
0
特征值虚部按不同比例因子 k变化的 GMRES算法的收敛曲线 (阶数 n=1000）
k=0.2 k=0.5 k=2 k=5
10
-5
||rk||/||b||
10
-10
10
-15
0
100
200
300
400
500
600
700
800
900
1000
迭代次数
图8 特征值虚部按不同比例因子k变化的GMRES法求解求解Ax=b的收敛曲线
图7 特征值虚部按不同比例因子k变化的Arnoldi法求解求解Ax=b的收敛曲线
（3）GMRES法求解求解A′x=b: 利用matlab编程实现GMRES算法。b = ones(1000,1)，x0 = zeros(1000, 1)。 计算每一步迭代的残差rk相对于初始残差的2范数。相对残差2范数的对数 值与迭代步数的关系曲线如图8所示：
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
高等数值分析第二次实验作业T1. 构造例子特征值全部在右半平面时, 观察基本的Arnoldi 方法和GMRES 方法的数值性态, 和相应重新启动算法的收敛性. Answer:（1） 构造特征值均在右半平面的矩阵A ：根据实Schur 分解，构造对角矩阵D 由n 个块形成，每个对角块具有如下形式，对应一对特征值i i i αβ±&ii i i i S αββα-⎛⎫=⎪⎝⎭这样D=diag(S 1,S 2,S 3……S n )矩阵的特征值均分布在右半平面。
生成矩阵A=U TAU ，其中U 为正交阵，则A 矩阵的特征值也均在右半平面。
不妨构造A 如下所示：2211112222/2/2/2/2N NA 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:NA(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; endU = 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); toctic[xG,rmG,flagG]=GMRES(A1,b,x0,e,m); tocsubplot(1,2,1);semilogy(rmA)。
title('ArnoldiÊÕÁ²ÇúÏß') xlabel('µü´ú²½Êýk')ylabel('log(||rk||)')subplot(1,2,2);semilogy(rmG)title('GMRESÊÕÁ²ÇúÏß') xlabel('µü´ú²½Êýk')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:150tic[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、endm1=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);toctic】[xG,rmG,flagG,MaxiG]=GMRESM(A,b,x0,tol,m,Maxm);tocm1=1:MaxiA;m2=1:MaxiG;plot(m1,log10(rmA),'o-',m2,log10(rmG),'*-')title('Á½ÖÖÖØÆôËã·¨µÄÊÕÁ²ÇúÏß')xlabel('ÖØÆô´ÎÊýk')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 ii s S s αββα--⎛⎫= ⎪-⎝⎭这里将上面的矩阵生成做如下改造：clear all ; N=500 tol=1e-6; m=50;@Maxm=50;A=zeros(2*N); s=for a=1:NA(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;endU = orth(rand(2*N,2*N));A1 = U'*A*U;《改变s的值，进而改变实部为负的特征值个数，计算结果整理如下表所示：、*<结果讨论：1. 当开始有负半平面的特征值时，个数比较少时（如小于100时），随着个数的增加，基本的Arnoldi 和基本的GMRES 迭代次数明显变多，收敛速度明显变慢；当负半平面特征值个数继续增加（如大于100时），两者的迭代次数减少，收敛速度明显变快，并且将比第一题的收敛速度快很多，迭代次数少很多。
2. 当开始有负半平面的特征值时，个数比较少时（如小于100时），随着个数的增加，Arnoldi出现峰点，峰点个数与其特征值的分布有关系，但整体仍收敛。
这是由于负特征值的存在，使得海森贝格矩阵H 发生近似奇异而发生近似中断而引起的。
然而，GMRES 的残差始终平稳下降，当Aronldi 出现尖峰时，GMRES 的残差不变具有最优性。
T3. 对1 中的例子固定特征值的实部, 变化虚部, 比较收敛性。
Answer:首先对T1的代码进行简单修改，实部不变，虚部进行一定的放大（引入系数k ）。
[每个对角块具有如下形式，对应一对特征向量'i i i αβ+''i i ii S αββα⎛⎫-= ⎪⎝⎭上式中的'i i k ββ=相应代码如下所示（以Arnoldi 为例），取k=,,1,2,5这五种情况。
N=500A1=zeros(2*N); k1=;for a=1:NA1(2*a-1,2*a-1)=a;A1(2*a-1,2*a)=-k1*a;$A1(2*a,2*a-1)=k1*a;A1(2*a,2*a)=a;endU = orth(rand(2*N,2*N));A1 = U'*A1*U;%篇幅所限，此处省去了A2~A5，同理可得。
k5=5;for a=1:NA5(2*a-1,2*a-1)=a;A5(2*a-1,2*a)=-k5*a;A5(2*a,2*a-1)=k5*a;(A5(2*a,2*a)=a;endU = orth(rand(2*N,2*N));A5 = U'*A5*U;x0=zeros(2*N,1);b=ones(2*N,1);e=1e-6;m=1000;[xA1,rmA1,flagA1]=Arnoldi(A1,b,x0,e,m);[xA2,rmA2,flagA2]=Arnoldi(A2,b,x0,e,m);》[xA3,rmA3,flagA3]=Arnoldi(A3,b,x0,e,m);[xA4,rmA4,flagA4]=Arnoldi(A4,b,x0,e,m);[xA5,rmA5,flagA5]=Arnoldi(A5,b,x0,e,m);m1=1:size(rmA1,2);m2=1:size(rmA2,2);m3=1:size(rmA3,2);m4=1:size(rmA4,2);m5=1:size(rmA5,2);plot(m1,log10(rmA1),m2,log10(rmA2),m3,log10(rmA3),m4,log10(rmA4),m5,log10(rmA5)) title('»ù±¾Arnoldi·½·¨ÊÕÁ²ÇúÏß')xlabel('µü´ú´ÎÊýk')$ylabel('log(||rk||)')hg1 = legend('k=', 'k=','k=1','k=2','k=5');计算结果如下所示：结果讨论：1.随着比例因子k的增大，虚部被放大，Arnoldi方法和GMRES方法的收敛速度均减慢，迭代次数增加；2.同时，随着比例因子k的增大，Arnoldi过程跳动越来越严重，峰点个数增加，发生近似中断的次数增加；而GMRES方法的残差始终保持单调平稳下降。