清华大学高等数值分析_第二次实验作业
数值分析第二次大作业——编写ln函数实验报告
为一个整体，将 ln(a)进行 Taylor 级数分解，进而按式计算可以求得 ln(x)的值。由于我们选 取级数的前 n 项和近似 ln(a)，则有：
Rn (a
1)
(1)n
(a 1)n1 n 1
…；
|
Rn
(a
1)
||
(a 1)n1 n 1
|
。
2
自 92 乔晖 2009011414
精度为 20 位，即 h4 1020 。
int * GetNumber() {return number;}
//获得大数中的整数数组
bool GetSgn() const {return sgn;}
//获得大数的符号
bool IsZero();
//判断大数是否为 0
CBigNumber Reverse() {sgn = !sgn;return *this;} //取相反数
数值分析——编写 ln 函数
实验报告
自 92 乔晖 2009011414
自 92 乔晖 2009011414
数值分析——编写 ln 函数实验报告
目录
一、 实验任务 ................................................... 2
二、 编译环境 ................................................... 2
1u
x 1
35
5
自 92 乔晖 2009011414
数值分析——编写 ln 函数实验报告
其前 n 项和，近似可得 ln(x)的值。
基于以上原理，可以设计迭代算法，累加 2u2i1 得到前 n 项和。每次计算后，判断计 2i 1
数值分析第二次上机作业实验报告
一．实验任务用MA TLAB 语言编写连续函数最佳平方逼近的算法程序（函数式M 文件）。
并用此程序进行数值试验，写出实验报告。
二．实验方法最佳平方逼近方法采用基于正交多项式的最佳平方逼近，选择Lengendre 多项式做基。
计算组合系数时，函数的积分采用变步长复化梯形求积法。
三．程序功能和使用说明1.采用基于正交多项式的最佳平方逼近，选择Lengendre 多项式做基利用递推关系0112()1,()()(21)()(1)()/2,3,.....n n n P x P x xP x n xP x n P x n n --===---⎡⎤⎣⎦=可构造出用户需要的任意次数的最佳平方逼近多项式。
2. 用M 文件建立数学函数，实现程序通过修改建立数学函数的M 文件以适用不同的被逼近函数。
3.已经考虑一般的情况]1,1[],[)(+-≠∈b a x f ，程序有变量代换的功能。
4.计算组合系数时，函数的积分采用变步长复化梯形求积法5.可根据需要，求出二次、三次、。
最佳平方逼近函数）x s (。
6.最后作出逼近函数）x s (和被逼近函数)(x f 的曲线图可进行比较，分别用绘图函数plot 和fplot 绘图。
7.在matlab 的命令窗口，输入[c,sx]=leastp(@func1,a,b,n)，func1是被逼近函数，b 和a 分别是逼近函数的上、下区间，n 为最佳平方逼近的次数，可为任意次数。
四．程序代码（含注释）1. 最佳平方逼近主函数function [c,sx]=leastp(func,a,b,n)%LEASTP.m:least-square fitting with legendre polynomials%func 指被逼近函数，调用需要用句柄%a,b 分别指被逼近函数的区间上下限%n 指最佳平方逼近的次数syms t;syms x;%以Lengendre 多项式为基，构造任意次数的最佳平方逼近多项式p(2)=t;p(1)=1;if n>1for j=3:1:(n+1)p(j)=((2*j-3)*t*p(j-1)-(j-2)*p(j-2))/(j-1);endend%变量代换，区间调整为[-1,1]f=feval(func,(b-a)/2*t+(b+a)/2);%计算组合系数，其中调用变步长复化梯形求积函数trapzfor j=1:1:(n+1)c(j)=(2*j-1)/2*trapz(f*p(j),-1,1);end%将组合系数与对应的最佳平方多项式相乘然后求和，得到最佳逼近函数sx=0;for j=1:1:(n+1)sx=sx+c(j)*p(j);end%将变量替换还原sx=subs(sx,(2*x-a-b)/(b-a));%使用fplot绘制原函数图像f1=feval(func,x);f1=inline(f1);[x,y]=fplot(f1,[a,b]);plot(x,y,'r-','linewidth',1.5);hold on;%使用plot绘制最佳平方逼近函数图像g=linspace(a,b,(b-a)*300);fsx=subs(sx,g);plot(g,fsx,'b-','linewidth',1.5);str=strcat(num2str(n),'次最佳平方逼近');legend('原函数',str);end2. 计算组合系数，变步长复化梯形求积法function To1=trapz(func,a,b)%半分区间复化梯形公式计算定积分%func指需要求积分的原函数%a,b分别指积分上下区间%初值h=b-a;To=(subs(func,a)+subs(func,b))*(b-a)/2;e=1;while e>10^-6%迭代终止条件，前后两次积分值差小于10^-6 H=0;x=a+h/2;while x<bH=H+subs(func,x);%计算出所有二分新出现的值的和x=x+h;endTo1=0.5*(To+h*H);%计算出新的积分值e=abs(To1-To);h=h/2;%继续半分区间，进行迭代计算To=To1;endend3. 以.m文件定义被逼近函数function y=func1(x)y=x*cos(x);end五．实验结果1. 一次最佳平方逼近c =-1.1702 -2.4235sx=1.253290 - 1.211752*x2. 二次最佳平方逼近c =-1.1702 -2.4235 -0.4265sx=-0.159939*x^2 - 0.571997*x + 0.8267873. 三次最佳平方逼近c =-1.1702 -2.4235 -0.4265 1.2216sx=0.381759*x^3 - 2.450495*x^2 + 3.092892*x - 0.3948434. 四次最佳平方逼近c =-1.1702 -2.4235 -0.4265 1.2216 0.3123sx =0.085392*x^4 - 0.301375*x^3 - 0.693864*x^2 + 1.531443*x - 0.082553六．分析与讨论从次数从1到4的最佳平方逼近图像对比可以发现，次数越高，图像拟合效果越好。
清华大学数值分析报告某实验报告材料
数值分析实验报告一、 实验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 ki k ik k i na a ≤≤=，当k i k ≠时，对B 中第k 行与第k i 行交换，然后按照和顺序消去法相同的步骤进行。
数值分析第二次大作业
《数值分析》计算实习报告第二题院系：机械工程及自动化学院学号：姓名：2017年11月一、题目要求试求矩阵A =[a ij ]10×10的全部特征值，并对其中的每一个实特征值求相应的特征向量，已知a ij ={sin (0.5i +0.2j ) i ≠j1.52cos (i +1.2j ) i =j(i,j =1,2, (10)说明：1.用带双步位移的QR 方法求矩阵特征值，要求迭代的精度水平为ε=10−12。
2.打印以下内容： （1）全部源程序；（2）矩阵A 经过拟上三角化后所得的矩阵A (n−1)； （3）对矩阵A (n−1)实行QR 方法迭代结束后所得的矩阵；（4）矩阵A 的全部特征值λi =(R i ,I i ) (i =1,2,⋯,10)，其中R i =Re(λi )，I i = Im(λi ) 。
若λi 是实数，则令I i =0； （5）A 的相应于实特征值的特征向量。
3.采用e 型数输出实型数，并且至少显示12位有效数字。
二、算法设计思路和方案1. 将矩阵A 拟上三角化得到矩阵A (n−1)为了减少计算量，一般先利用Householder 矩阵对矩阵A 作相似变换，把A 化为拟上三角矩阵A (n−1)，然后用QR 方法计算A (n−1)的全部特征值，而A (n−1)的特征值就是A 的特征值。
具体算法如下：记(1)A A =，()r A 的第r 列至第n 列的元素为(r)(1,2,,;,1,,)ij a i n j r r n ==+。
对于1,2,,2r n =-执行（1）若()(2,3,,)r ira i r r n =++全为零，则令(1)()r r A A +=，转(5)；否则转（2）。
（2）计算r d =()()1,sgn r r r r r c a d +=- (若()1,0r r r a +=，则取r r c d =)2()1,r r r r r r h c c a +=-（3）令()()()()1,2,0,,0,,,,T r r r n r r rr r rnru ac aaR ++=-∈。
数值分析大作业
第二次计算实验：SVD及其应用梁杰存2014310739航博1431.方法求矩阵A奇异值分解一个途径是求解A T A的特征值，但因为舍入误差容易丢掉小奇异值。
因此通常先将矩阵上双对角化，即构造正交阵Q和W,使得Q T AW=B(upper−bidiagnal)。
这一过程可以通过逐次Household变换或逐次Given’s变换完成，还有一种基于待定系数法思想的Lanczos算法。
由于Linpack中SVD算法需要输入上双对角矩阵，本文采用Lanczos 算法实现上双对角化。
1.1.隐式零移位QR法（implicit zero-shift QR）与传统移位QR迭代算法不同，隐式零移位QR算法不进行移位，并且第一步构造右乘Given’s变换矩阵GR(1,2)将上双对角矩阵B的(1,2)位置上的元素12b消零，而不是传统方法中引入一个非零元素21b。
但这一步可能会使原来为零的b12变为非零。
第二步左乘Given’s阵GL(1,2)使得12b为0，但可能会使为零b13变为非零。
与上述步骤类似，将b13变为0后可能会使b23非零。
如下图所示，重复上述步骤最终将恢复为上双对角矩阵，即完成一步隐式零移位QR迭代。
反复迭代，矩阵B将趋近与对角阵阵，对角元即特征值。
图1隐式QR迭代1.2.分而治之(Divide-and-conquer)分而治之算法将上双对角阵B分成有两个互相独立对角块矩阵与另一矩阵之和，即：B=B100B2+b m vvT=Q1Σ1Q1T00Q2Σ1Q2T+b m vv T =Q100Q2(Σ100Σ1+b m uuT)Q1T00Q2T所以矩阵B的特征值与矩阵D+ρu u T的特征值相同，其中D=Σ100Σ1为对角阵，又：det D+ρu u T−λI=det((D−λI)(I+ρD−λ−1u u T))由于D−λI非奇异，则det I+ρD−λ−1u u T=1+ρu T D−λ−1u=1+ρu i2d i−λ=0ni=1在每个d i与d i+1之间分布着一个特征值，可用牛顿法快速找到该特征值。
数值分析上机作业2
数值实验数值实验综述：线性代数方程组的解法是一切科学计算的基础与核心问题。
求解方法大致可分为直接法和迭代法两大类。
直接法——指在没有舍入误差的情况下经过有限次运算可求得方程组的精确解的方法，因此也称为精确法。
当系数矩阵是方的、稠密的、无任何特殊结构的中小规模线性方程组时，Gauss消去法是目前最基本和常用的方法。
如若系数矩阵具有某种特殊形式，则为了尽可能地减少计算量与存储量，需采用其他专门的方法来求解。
Gauss消去等同于矩阵的三角分解，但它存在潜在的不稳定性，故需要选主元素。
对正定对称矩阵，采用平方根方法无需选主元。
方程组的性态与方程组的条件数有关，对于病态的方程组必须采用特殊的方法进行求解。
实验一一、实验名称：矩阵的LU分解.二、实验目的：用不选主元的LU分解和列主元LU分解求解线性方程组Ax=b, 并比较这两种方法.三、实验内容与要求（1）用所熟悉的计算机语言将不选主元和列主元LU分解编成通用的子程序，然后用编写的程序求解下面的84阶方程组将计算结果与方程组的精确解进行比较，并就此谈谈你对Gauss消去法的看法。
（2）写出追赶法求解三对角方程组的过程，并编写程序求该实验中的方程组（1）①不选主元高斯消去法求解方程组function [L,U]=gauss1(A,b)n=length(A);for k=1:n-1A(k+1:n,k)=A(k+1:n,k)/A(k,k);A(k+1:n,k+1:n)=A(k+1:n,k+1:n)-A(k+1:n,k)*A(k,k+1:n);endL=tril(A(:,1:n),-1)+eye(n);U=triu(A);主程序为：Clear；clc；a=ones([84,1])*6;b=ones([83,1]);c=ones([83,1])*8;A=diag(a)+diag(b,1)+diag(c,-1); d=ones([82,1])*15;b=[7;d;14];[L U]=gauss1(A,b);n=length(A);y=ones(n,1);y=L\b;x=ones(n,1);x=U\y解为：x=11.000000000000001.000000000000001.000000000000000.9999999999999981.000000000000000.9999999999999931.000000000000010.9999999999999721.000000000000060.9999999999998861.000000000000230.9999999999995451.000000000000910.9999999999981811.000000000003640.9999999999927251.000000000014550.9999999999708981.000000000058200.9999999998835921.000000000232820.9999999995343671.000000000931270.9999999981374691.000000003725060.9999999925498741.000000014900250.9999999701994971.000000059601010.9999998807979861.000000238404030.9999995231919461.000000953616110.9999980927677831.000003814464440.9999923710711301.000015257857740.9999694842845201.000061031430960.9998779371380811.000244125723840.9995117485523231.000976502895350.9980469942092911.003906011581410.9921879768371871.015624046325570.9687519073490881.062496185300920.8750076294018071.249984741181830.5000305176945351.99993896437811 -0.999877927824961 4.99975585192486 -6.99951168894947 16.9990233182979 -30.9980463981918 64.9960918427676 -126.992179871071 256.984344484284 -510.968627937136 1024.93701174855 -2046.87304699420 4096.74218797682 -8190.46875190732 16383.8750076293 -32764.5000305175 65531.0001220701 -131055.000488281 262097.001953124 -524127.007812498 1048001.03125000 -2094975.124999994185857.49999999-8355328.9999999716645128.9999999-33028126.999999965007744.9999998-125821439.000000234866688.999999-402628606.999999536838144.999998可见，这是一个病态方程，从56个跟开始发散。
清华大学 李津老师 数值分析第二次实验作业
就不再赘述了。 二、实际计算 生成十个不同的(最好属不同类型或有不同性质的)的 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
清华大学第五版【数值分析】习题答案
第一章 绪论(1)1．设0x >,x 的相对误差为δ，求ln x 的误差。
解：近似值*x 的相对误差为*****r e x x e x x δ-=== 而ln x 的误差为()1ln *ln *ln **e x x x e x =-≈ 进而有(ln *)x εδ≈2．设x 的相对误差为2%，求nx 的相对误差。
解：设()n f x x =，则函数的条件数为'()||()p xf x C f x = 又1'()n f x nx-=, 1||n p x nx C n n-⋅∴== 又((*))(*)r p r x n C x εε≈⋅且(*)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题所给的数。
解：*41*32*13*34*151()1021()1021()1021()1021()102x x x x x εεεεε-----=⨯=⨯=⨯=⨯=⨯***124***1244333(1)()()()()1111010102221.0510x x x x x x εεεε----++=++=⨯+⨯+⨯=⨯ ***123*********123231132143(2)()()()()1111.10210.031100.031385.610 1.1021385.6102220.215x x x x x x x x x x x x εεεε---=++=⨯⨯⨯+⨯⨯⨯+⨯⨯⨯≈**24****24422*4335(3)(/)()()110.0311056.430102256.43056.43010x x x x x x xεεε---+≈⨯⨯+⨯⨯=⨯=5计算球体积要使相对误差限为1，问度量半径R 时允许的相对误差限是多少？ 解：球体体积为343V R π=则何种函数的条件数为23'4343p R V R R C V R ππ===(*)(*)3(*)r p r r V C R R εεε∴≈=又(*)1r V ε=故度量半径R 时允许的相对误差限为1(*)10.333r R ε=⨯≈6．设028Y =，按递推公式1n n Y Y -=（n=1,2,…）计算到100Y 27.982≈（5位有效数字），试问计算100Y 将有多大误差？解：1n n Y Y -=-10099Y Y ∴=-9998Y Y =9897Y Y =-……10Y Y =-依次代入后，有1000100Y Y =-即1000Y Y =27.982, 100027.982Y Y ∴=-*310001()()(27.982)102Y Y εεε-∴=+=⨯100Y ∴的误差限为31102-⨯。
北航数值分析全部三次大作业
北航数值分析全部三次大作业第一次大作业是关于解线性方程组的数值方法。
我们被要求实现各种常用的线性方程组求解算法，例如高斯消元法、LU分解法和迭代法等。
我首先学习了这些算法的原理和实现方法，并借助Python编程语言编写了这些算法的代码。
在实验中，我们使用了不同规模和条件的线性方程组进行测试，并比较了不同算法的性能和精度。
通过这个作业，我深入了解了线性方程组求解的原理和方法，提高了我的编程和数值计算能力。
第二次大作业是关于数值积分的方法。
数值积分是数值分析中的重要内容，它可以用于计算曲线的长度、函数的面积以及求解微分方程等问题。
在这个作业中，我们需要实现不同的数值积分算法，例如矩形法、梯形法和辛普森法等。
我学习了这些算法的原理和实现方法，并使用Python编写了它们的代码。
在实验中，我们计算了不同函数的积分值，并对比了不同算法的精度和效率。
通过这个作业，我深入了解了数值积分的原理和方法，提高了我的编程和数学建模能力。
第三次大作业是关于常微分方程的数值解法。
常微分方程是数值分析中的核心内容之一，它可以用于描述众多物理、化学和生物现象。
在这个作业中，我们需要实现不同的常微分方程求解算法，例如欧拉法、龙格-库塔法和Adams法等。
我学习了这些算法的原理和实现方法，并使用Python编写了它们的代码。
在实验中，我们解决了一些具体的常微分方程问题，并比较了不同算法的精度和效率。
通过这个作业，我深入了解了常微分方程的原理和方法，提高了我的编程和问题求解能力。
总的来说，北航数值分析课程的三次大作业非常有挑战性，但也非常有意义。
通过这些作业，我在数值计算和编程方面得到了很大的提升，也更加深入地了解了数值分析的理论和方法。
虽然这些作业需要大量的时间和精力，但我相信这些努力将会对我未来的学习和工作产生积极的影响。
清华大学高等数值分析实验设计及答案
高等数值分析实验一工物研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，由于上实验结果可以看出实际叠代次数都比上限值要小较多。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
与 与 与 Arnoldi与 与 与 与 与 与 与 与 m与 与 与 与 与 900
与 与 与 Arnoldi与 与
800
与 与 与 GMRES与 与
700
600
与与与与与
500
400
300
200
100 0
100 200 300 400 500 600 700 800 900 1000 m
注：上图标题有误，应为两种基本算法迭代步长与负特征值个数 m 的关系曲线
4. 结论 从上述结果中分析可得如下结论： 1) 对于这两种算法来说，重启动算法的总迭代步长要大于基本的算法，这个很容易理解，这
是由于在重启动时，会丢掉一些之前算过的信息，从另一个点我们认为更好（其实不一定 好）的初值x0重新开始计算。 2) 这两种重启动算法，在两次重启动的交接点处都会产生不连续。GMRES在两次重启动之间 不满足残差单调不增，然而在单次重启的内部计算过程中，则满足单调不增的原则。 Arnodli过程由于不具有残差二范数最小的性质，所以在以上两个位置处都可能会发生跳 动。 3) 一般来说，重启次数随着m的增大不断减小，当m非常大的时候，就过度到了基本的不启 动算法。 4) 随着m的增大，迭代总次数并不是单调下降的。重新启动的Arnoldi方法在m=40时迭代次数 最少，重新启动的GMRES方法m=30或40时迭代次数最少。 5) 计算代价（用运算时间估算）与m不成线性关系，存在一个最优值，当m大到一定程度， 虽然重启次数很小，但是每步代价很大，随着m的进一步增加，代价逐渐增大。 6) 重启的算法如果收敛的话，代价可能会远远小于不重启的算法。但是最合适的m的选取因 问题不同而不同，不好把握。 7) 对于同样的矩阵A，GMRES方法的收敛速度一般比相应的Arnoldi方法快。这是由于其残差 具有最优性，而且在迭代步数相同的情况下，GMRES方法的相对残差比相应的Arnoldi方法 的相对残差要小。
与 与 与 与 与 与 与 与 与 与 与 与 m与 与 与 与 与 与 与
与 与 GMRES 与 与 Aronldi
1600
1500
与与与与与
1400
1300
1200
1100
1000
900 0
50
100
150
m
从表中数据可以看出，结果相差并不大 2. 重启次数与 m 之间的关系
与 与 与 与 与 与 与 与 与 与 与 m与 与 与 与 与 与 与 250
与 与 与 Arnoldi与 与 与 与 与 GMRES与 与
100
10-2
10-4
10-6
10-8 0
102 100
100
200
300
400
500
600
700
与与与与
与 与 与 与 与 与 与 ||rk||与 与 与 与 (与 与 与 与 与 与 与 20与 )
与 与 与 Arnoldi与 与 与 与 与 GMRES与 与
第二题。对于 1 中的矩阵，将特征值进行平移，使得实部有正有负，和 1 的结果进行比较，方法的收敛速度会如 何？基本的 Arnoldi 算法有无峰点？若有，基本的 GMRES 算法相应地会怎样？ 答：
1. 使用如下方法对第一题中的 A 矩阵进行平移
A' AI
则得到的矩阵的特征值为 i' i 。
此外，从收敛曲线上看，由于特征值均处在右半平面，收敛曲线平滑，收敛速度（收敛因子）比较均匀。
3、重启动的 GMRES 和 Arnoldi 算法 对上述 A、b、x0 使用重启动的 Arnoldi 和 GMRES 算法。变化 m 经过多次实验，得到如下结果： 1. 总迭代步长与 m 之间的关系
1900 1800 1700
2 2
A
2 2O
OO
n / 2 n / 2 n / 2 n / 2 nn
其中，A 由 n/2 个块形成，每个对角块具有如下形式，对应一对特征向量i ii
A
i i
i i
、
这里，取 n=1000，得到矩阵 A。经过验证，A 的特征值分布均在右半平面，如下图所示
||rk||/||b||
||rk||/||b||
||rk||/||b||
与 与 与 与 与 与 与 ||rk||与 与 与 与 (与 与 与 与 与 与 与 1与 ) 102
与 与 与 Arnoldi与 与 与 与 与 GMRES与 与
100
10-2
10-4
10-6
10-8 0
104 102
100
200
与 与 GMRES 与 与 Aronldi
200
150
与与与与
100
50
0
0
50
3. 计算时间与 m 之间的关系
与 与 与 Arnoldi与 与 与 与 与 与 与 与 与 m与 与 与 与 与 5.5
100
150
m
与 与 与 GMRES与 与 与 与 与 与 与 与 与 m与 与 与 与 与 10
第二次实验作业 清华大学 高等数值分析 贾仲孝老师作业 第一题：构造例子特征值全部在右半平面时，观察基本的 Arnoldi 方法和 GMRES 方法的数值性态，和相应重新 启动算法的收敛性。 答： 1、计算初始条件 1) 矩阵 A 的生成 根据实 Schur 分解，构造矩阵如下形式
1 1
1 1
100
与 与 与 GMRES与 与
10-1
10-1
10-2
10-2
||rk||/||b||
||rk||/10-4
10-4
10-5
10-5
10-6
10-6
10-7 0
101 100
20 40 60 80 100 120 140 160 180 200 与与与与
与 与 与 与 与 与 与 ||rk||与 与 与 与 (与 与 与 与 与 与 与 950与 )
与与与与
与 与 与 与 与 与 与 ||rk||与 与 与 与 (与 与 与 与 与 与 与 600与 ) 101
与 与 与 Arnoldi与 与
100
与 与 与 GMRES与 与
与 与 与 与 与 与 与 ||rk||与 与 与 与 (与 与 与 与 与 与 与 900与 ) 101
与 与 与 Arnoldi与 与
104 102
100
200
300
400
500
600
700
800
与与与与
与 与 与 与 与 与 与 ||rk||与 与 与 与 (与 与 与 与 与 与 与 400与 )
与 与 与 Arnoldi与 与 与 与 与 GMRES与 与
100
10-2
10-4
10-6
10-8 0
50
100
150
200
250
300
400
500
600
700
与与与与
与 与 与 与 与 与 与 ||rk||与 与 与 与 (与 与 与 与 与 与 与 5与 )
与 与 与 Arnoldi与 与 与 与 与 GMRES与 与
100
10-2
10-4
10-6 0
102 100
100
200
300
400
500
600
700
与与与与
与 与 与 与 与 与 与 ||rk||与 与 与 与 (与 与 与 与 与 与 与 50与 )
10-1
10-1
10-2
10-2
||rk||/||b||
||rk||/||b||
10-3
10-3
10-4
10-4
10-5
10-5
10-6
10-6
10-7 0
20
40
60
80
100 120 140 160 180
与与与与
10-7 0
20
40
60
80
100 120 140 160 180
与与与与
当负特征值个数逐渐增多时，迭代次数与负特征值个数 m 的关系曲线如下图
与 与 Re(x)
500 400 300 200 100
0 -100 -200 -300 -400 -500
0
与 与 与 与 A与 与 与 与 与 与 与 与
与与与
50 100 150 200 250 300 350 400 450 500 与 与 Im(x)
2) b 的初值为 b=(1,1,1,1..1)T 3) 迭代初值为 x0=0 4) 停机准则为
10-4
10-5
10-6
10-7 0
100
200
300
400
500
600
与与与与
结果讨论：从图中可以看出，基本的 Arnoldi 方法经过 554 步收敛，基本的 GMRES 方法经过 535 步收敛。这 是由于 GMRES 具有残差最优性，会略快于 Arnoldi 方法，但是，由于两种方法的基本原理近似，GMRES 方法不 会实质性的提速。
对第一题的 A，分别取、、、、、、、、、、，得到 的负特征值分别为、、、、、、、、、、个。分别求解上述矩阵，得到的结 果如下
||rk||/||b||
10-2
10-4
10-6
10-8 0
104 102
100
200
300
400
500
600
700
与与与与
与 与 与 与 与 与 与 ||rk||与 与 与 与 (与 与 与 与 与 与 与 100与 )