(整理)高等数值分析作业-第二次实验

合集下载

数值分析第二次作业及答案

数值分析第二次作业及答案1. 用矩阵的直接三角分解法（LU 分解）解方程组1234102050101312431701037x x x x ⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦。

212223243132333441424344212223243132333441424344110201020101011124310103110201101,112110101l uu u l l u u l l l u luu u l l u u l l l u ⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⇒=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎢⎥⎣⎦⎣⎦⎣⎦解：设111122223333444410201012125115102053101310136212117216420101724y y x x y y x x y y x x y y x x ⎡⎤⎢⎥⎢⎥=⎢⎥⎢⎥⎥⎣⎦==⎡⎤⎧⎡⎤⎧⎡⎤⎡⎤⎡⎤⎡⎤⎢⎥⎪⎢⎥⎪⎢⎥⎢⎥⎢⎥⎢⎥==⎪⎪⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥=⇒=⇒⎨⎨⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥==⎪⎪⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎪⎪==⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦⎣⎦⎣⎦⎩⎣⎦⎩由2.矩阵第一行乘以一数，成为211A λλ⎡⎤=⎢⎥⎣⎦,证明当23λ=±时，()cond A ∞有最小值。

112131213 21121223A AA A λλλλλλλλ--∞∞⎧⎛⎫≥- ⎪⎪⎛⎫⎪=⇒==⇒=+ ⎪⎨⎪ ⎪⎝⎭⎪<- ⎪⎪⎩⎝⎭解：123632() min ()7.22343cond A AA cond A λλλλλ-∞∞∞∞⎧+≥⎪⎪∴====⎨⎪+<⎪⎩故当时， 3. 设有方程组AX b =，其中121011221,,302223A b ⎡⎤⎢⎥-⎡⎤⎢⎥⎢⎥⎢⎥==⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎢⎥-⎢⎥⎣⎦已知它有解12130X ⎡⎤⎢⎥⎢⎥⎢⎥=-⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦。

北航数值分析大作业二(纯原创,高分版)

(R_4 ,I_4 )=( 1.590313458807e+000, 0.000000000000e+000)
(R_5 ,I_5 )=(-1.493147080915e+000, 0.000000000000e+000)
(R_6 ,I_6 )=(-9.891143464723e-001, 1.084758631502e-001)
-0.8945216982
-0.0993313649
-1.0998317589
0.9132565113
-0.6407977009
0.1946733679
-2.3478783624
2.3720579216
1.8279985523
-1.2630152661
0.6790694668
-0.4672150886
6.220134985374e-001
-1.119962139645e-001
-2.521344456568e+000
-1.306189420531e+000
-3.809101150714e+000
8.132800093357e+000
-1.230295627285e+000
-6.753086301215e-001
而其本质就是
1.令以及最大迭代步数L；
2.若m≤0，则结束计算，已求出A的全部特征值，判断或或m≤2是否成立，成立则转3，否则转4；
3.若，则得一个特征值，m=m-1，降阶；若，则计算矩阵：
的特征值得矩阵A的两个特征值，m=m-2，降阶，转2.；
4.若k≤L，成立则令
k=k+1，转2，否则结束计算，为计算出矩阵A的全部特征值；

数值分析第二次大作业——编写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的最佳平方逼近图像对比可以发现，次数越高，图像拟合效果越好。

数值分析第二次大作业

《数值分析》计算实习报告第二题院系：机械工程及自动化学院学号：姓名：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 ++=-∈。

清华大学李津老师数值分析第二次实验作业

就不再赘述了。二、实际计算生成十个不同的(最好属不同类型或有不同性质的)的 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

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

高等数值分析第二次实验作业注：矩阵阶数均为1000阶T1. 构造例子特征值全部在右半平面时，观察基本的Arnoldi方法和GMRES方法的数值性态，和相应重新启动算法的收敛性。

Answer：➢关于特征值均在右半平面的矩阵A：首先构造对角矩阵D，其中D矩阵的对角是由以下形式的2×2的子矩阵组成其对角元：S=[a−bb a]，其中a，b为实数此时，S矩阵的特征值分别为a+bi和a-bi，这样D=diag(S1,S2,S3……S n)矩阵的特征值均分布在右半平面。

另矩阵A=Q T AQ，则A矩阵的特征值也均在右半平面，生成矩阵的过程代码如下所示：N=500 %生成的矩阵为2N阶A=zeros(2*N);%delta控制特征值分布的密集程度，即控制条件数大小，delta越大条件数越大delta=0.1;%构造D矩阵for j=1:NA(2*j-1,2*j-1)=N+j*delta;A(2*j-1,2*j)=-N-j*delta;A(2*j,2*j-1)=N+j*delta;A(2*j,2*j)=N+j*delta;endU = orth(rand(2*N,2*N));A = U'*A*U;➢首先进行观察基本的Arnoldi方法和GMRES方法的数值性态，考虑N=500，即矩阵为1000阶，取X0=zeros(N,1)，b=ones(N,1)，收敛准则为e=10-6;Arnoldi方法函数如下：function [xm,error,num]=Arnoldi(A,x0,b,e) N=size(A,1);r=b-A*x0;belta=norm(r);v=r/belta;V=v;j=0;while norm(r)>e & j<Nj=j+1;num(j)=j;w=A*V(1:N,j);for i=1:j h(i,j)=V(1:N,i)'*w; w=w-h(i,j)*V(1:N,i);endh(j+1,j)=norm(w);if h(j+1,j) ~=0v=w/h(j+1,j);V=[V v];ende1=zeros(j,1);e1(1)=1;be1=belta*e1;try[L,U,S]=lu(h(1:j,1:j)); be1=S*be1;lym=LX(L,be1);ym=UX(U,lym);xm=x0+V(1:N,1:j)*ym; r=b-A*xm;enderror(j)=log10(norm(r)); endendGMRES方法的函数如下：function [xm,error,num] = GMRES(A,x0,b,e) %ARNOLDI Summary of this function goes here % Detailed explanation goes hereN=size(A,1);r=b-A*x0;belta=norm(r);v=r/belta;V=v;j=0;er=1000;while er > e & j<Nj=j+1;num(j)=j;w=A*V(1:N,j);for i=1:jh(i,j)=V(1:N,i)'*w;w=w-h(i,j)*V(1:N,i);endh(j+1,j)=norm(w);if h(j+1,j) ~=0v=w/h(j+1,j);V=[V v];endtry[Q,R]=qr(h(1:j+1,1:j)); gk=Q'*belta*eye(j+1,1); ym=minresYK(R,gk);xm=x0+V(1:N,1:j)*ym;er=gk(j+1);enderror(j)=log10(er);endend其中UX，LX和minresYK是自己编写的求解（上、下）三角矩阵的函数，代码如下：function yk=LX(TK,b)n=size(TK,1);yk=zeros(n,1);yk(1)=b(1)/TK(1,1);for i=2:nyk(i)=b(i);for j=1:i-1yk(i)=yk(i)-TK(i,j)*yk(j);endyk(i)=yk(i)/TK(i,i); endendfunction yk=UX(TK,b)n=size(TK,2);yk=zeros(n,1);yk(n)=b(n)/TK(n,n);for i=1:n-1k=n-i;yk(k)=b(k);for j=k+1:nyk(k)=yk(k)-TK(k,j)*yk(j);endyk(k)=yk(k)/TK(k,k);endendfunction yk=minresYK(TK,b)n=size(TK,2);yk=zeros(n,1);yk(n)=b(n)/TK(n,n);for i=1:n-1k=n-i;yk(k)=b(k);for j=k+1:nyk(k)=yk(k)-TK(k,j)*yk(j);endyk(k)=yk(k)/TK(k,k);endend两种方法的计算结果如下所示：Delta=1 即条件数较大 Delta=0.1 即条件数较小可以看到delta=0.1时两种方法收敛都很快，delta=1时，两种方法略有差别，GMRES 方法收敛速度略快一些，下面就去delta=1的情况进行第一题的讨论。

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

清华大学高等数值分析课程作业第二次实验作业第一题：构造例子特征值全部在右半平面时，观察基本的Arnoldi 方法和GMRES 方法的数值性态，和相应重新启动算法的收敛性。

答：1、计算初始条件1) 矩阵A 的生成根据实Schur 分解，构造矩阵如下形式11112222/2/2/2/2n nA n n n n ⨯-⎛⎫ ⎪ ⎪ ⎪- ⎪= ⎪ ⎪ ⎪- ⎪⎪⎝⎭其中，A 由n/2个块形成，每个对角块具有如下形式，对应一对特征向量i i i αβ+ii i i A αββα-⎛⎫= ⎪⎝⎭、这里，取n=1000，得到矩阵A 。

经过验证，A 的特征值分布均在右半平面，如下图所示50100150200250300350400450500-500-400-300-200-1000100200300400500复平面中A 的特征值分布情况实部 Im(x)虚部 R e (x )特征值2) b 的初值为 b=(1,1,1,1..1)T 3) 迭代初值为 x 0=0 4) 停机准则为 ε=10-62、基本的Arnoldi 和GMRES 方法代入前面提到的初始值A 、b 、x0，得到的收敛结果如下10020030040050060010-710-610-510-410-310-210-110两种基本算法的||r k ||收敛曲线 (阶数n=1000)迭代次数||r k ||/||b ||基本的Arnoldi 算法基本的GMRES 算法结果讨论：从图中可以看出，基本的Arnoldi 方法经过554步收敛，基本的GMRES 方法经过535步收敛。

这是由于GMRES 具有残差最优性，会略快于Arnoldi 方法，但是，由于两种方法的基本原理近似，GMRES 方法不会实质性的提速。

此外，从收敛曲线上看，由于特征值均处在右半平面，收敛曲线平滑，收敛速度（收敛因子）比较均匀。

3、重启动的GMRES 和Arnoldi 算法对上述A 、b 、x0使用重启动的Arnoldi 和GMRES 算法。

西南交大数值分析第二次大作业(可以运行)

数值分析第二次大作业1（1）用Lagrange插值法程序：function f=Lang(x,y,x0)syms t;f=0;n=length(x);for(i=1:n)l=y(i);for(j=1:i-1)l=l*(t-x(j))/(x(i)-x(j));end;for(j=i+1:n)l=l*(t-x(j))/(x(i)-x(j));endf=f+l;simplify(f);if(i==n)if(nargin==3)f=subs(f,'t',x0);elsef=collect(f);f=vpa(f,6);endendendx=[1,2,3,-4,5];y=[2,48,272,1182,2262]; t=[-1];disp('插值结果')yt=Lang(x,y,t)disp('插值多项式')yt=Lang(x,y)ezplot(yt,[-1,5]);运行结果：插值结果：Yt= 12.0000插值多项式：yt =4.0*t^4 - 2.0*t^3 + t^2 - 3.0*t + 2.0（2）构造arctan x在[1,6]基于等距节点的10次插值多项式程序：function f=New(x,y,x0)syms t;if(length(x)==length(y))n=length(x);c(1:n)=0.0;elsedisp('xºÍyÎ¬Êý²»µÈ£¡');return;endf=y(1);y1=0;xx=linspace(x(1),x(n),(x(2)-x(1)));for(i=1:n-1)for(j=1:n-i)y1(j)=y(j+1)-y(j);endc(i)=y1(1);l=t;for(k=1:i-1)l=l*(t-k);end;f=f+c(i)*l/factorial(i);simplify(f);y=y1;if(i==n-1)if(nargin==3)f=subs(f,'t',(x0-x(1))/(x(2)-x(1)));elsef=collect(f);f=vpa(f,6);endendend>>x=[1,1.5,2,2.5,3,3.5,4,4.5,5,5.5,6];y=[atan(1),atan(1.5),atan(2),atan(2.5),atan(3),atan(3.5),atan(4),atan(4.5),atan(5),atan (5.5),atan(6)];disp('插值多项式')yt=New(x,y)ezplot(yt,[1,6]);hold onezplot('atan(t)',[1,6])grid on运行结果：插值多项式yt =1.34684*10^(-10)*t^10 - 6.61748*10^(-9)*t^9 + 1.28344*10^(-7)*t^8 - 0.00000104758*t^7 - 0.00000243837*t^6 + 0.000149168*t^5 - 0.00176296*t^4 + 0.0125826*t^3 - 0.0640379*t^2 + 0.250468*t + 0.7853982（1）用MATLAB自带spline函数用于进行三次样条插值程序：>>x=[-5,-4,-3,-2,-1,0,1,2,3,4,5];y=[0.03846,0.05882,0.10000,0.20000,0.50000,1.00000,0.50000,0.20000,0.10000,0.0 5882,0.03846];xi=linspace(-5,5)yi=spline(x,y,xi);plot(x,y,'rp',xi,yi);hold on;syms xfx=1/(1+x^2);ezplot(fx);grid on运行结果：由图可知，三次样条插值多项式图像与原函数图像基本一致。

河北工业大学数值分析实验二

实验二拉格朗日插值法（2学时）一、目的与要求：熟悉拉格朗日插值多项式和牛顿插值多项式，注意其不同特点；二、实验内容：通过拉格朗日插值和牛顿插值多项式的两个实例的计算，了解两种求解方法，分析他们的优缺点。

三、程序与实例算法 1．输入x i ,y i (i=0,1,2,⋯,n),令L(x n )=0; 2．对=0,1,2,⋯,n 计算ll i (x)=∏≠=--nij i jijx xx x 0L n ← L n +l i (x)y i 程序与实例例1 已知函数表用三次拉格朗日多项式求x=0.5635的函数近似值。

牛顿插值多项式算法1.输入n,xi ,yi(i=0,1,2⋯,n);2.对k=1,2,3⋯,n, i=1,2, ⋯,k计算各阶差商f(x0,x1⋯,xk);3.计算函数值Nn (x)=f(x)+f[x, x1](x- x)+⋯+f[x, x1,⋯,xn](x- x)(x- x1)⋯(x-x1-n)程序与实例例2已知函数表用牛顿插值多项式求Nn (0.596)和Nn(0.895)。

上机实验作以下两题：1．按下列数据作二次插值，并求x1=-2，x2=0，x3=2.75时的函数近似值#include<stdio.h>void main(void){float xi[3],x,l,l1,l11,l2,l22,l3,l33,l4,l44,l5,l55;int i;for(i=0;i<3;i++)scanf("%f",&xi[i]);for(i=0;i<3;i++){x=xi[i];l1=(x+1.0)*(x-1.0)*(x-2.0)*(x-3.0)*(1.0);l11=(-3.0+1.0)*(-3.0-1.0)*(-3.0-2.0)*(-3.0-3.0);l1=l1/l11;l2=(x+3.0)*(x-1.0)*(x-2.0)*(x-3.0)*(1.5);l22=(-1.0+3.0)*(-1.0-1.0)*(-1.0-2.0)*(-1.0-3.0);l2=l2/l22;l3=(x+3.0)*(x+1.0)*(x-2.0)*(x-3.0)*(2.0);l33=(1.0+3.0)*(1.0+1.0)*(1.0-2.0)*(1.0-3.0);l3=l3/l33;l4=(x+3.0)*(x+1.0)*(x-1.0)*(x-3.0)*(2.0);l44=(2.0+3.0)*(2.0+1.0)*(2.0-1.0)*(2.0-3.0);l4=l4/l44;l5=(x+3.0)*(x+1.0)*(x-1.0)*(x-2.0)*(1.0);l55=(3.0+3.0)*(3.0+1.0)*(3.0-1.0)*(3.0-2.0);l5=l5/l55;printf("%f\n",l1+l2+l3+l4+l5);}2．}按下列数据作五次插值，并求x1=0.46，x2=0.55，x3=0.60时的函数近似值#include<stdio.h>void main(void){float xi[3],x,l,l1,l11,l2,l22,l3,l33,l4,l44,l5,l55,l6,l66;int i;for(i=0;i<3;i++)scanf("%f",&xi[i]);for(i=0;i<3;i++){x=xi[i];l1=(x-0.42)*(x-0.50)*(x-0.58)*(x-0.66)*(x-0.72)*(1.044030);l11=(0.30-0.42)*(0.30-0.50)*(0.3-0.58)*(0.3-0.66)*(0.3-0.72);l1=l1/l11;l2=(x-0.3)*(x-0.5)*(x-0.58)*(x-0.66)*(x-0.72)*(1.084620);l22=(0.42-0.3)*(0.42-0.5)*(0.42-0.58)*(0.42-0.66)*(0.42-0.72);l2=l2/l22;l3=(x-0.3)*(x-0.42)*(x-0.58)*(x-0.66)*(x-0.72)*(1.118030);l33=(0.5-0.3)*(0.5-0.42)*(0.5-0.58)*(0.5-0.66)*(0.5-0.72);l3=l3/l33;l4=(x-0.3)*(x-0.42)*(x-0.5)*(x-0.66)*(x-0.72)*(1.156030);l44=(0.58-0.3)*(0.58-0.42)*(0.58-0.5)*(0.58-0.66)*(0.58-0.72);l4=l4/l44;l5=(x-0.3)*(x-0.42)*(x-0.5)*(x-0.58)*(x-0.72)*(1.198170);l55=(0.66-0.3)*(0.66-0.42)*(0.66-0.5)*(0.66-0.58)*(0.66-0.72);l5=l5/l55;l6=(x-0.3)*(x-0.42)*(x-0.5)*(x-0.58)*(x-0.66)*(1.23223);l66=(0.72-0.3)*(0.72-0.42)*(0.72-0.5)*(0.72-0.58)*(0.72-0.66);l6=l6/l66;l=l1+l2+l3+l4+l5+l6;printf("%f\n",l);}}实验报告1.实验目的2.实验内容及源程序3.实验结果及分析4.实验编码调试过程中出现的问题总结5.心得体会。

数值分析-第二次作业答案

2020级数值分析第二次作业(非线性方程求根)参考答案和评分标准班级学号姓名一(20分)用二分法求方程3()10f x x x =--=在区间[1.0,1.5]内的一个实根，且要求有3位有效数字。

试完成：(1)估计需要二分的次数；(8分)解：容易知道方程在[1.0,1.5]有且仅有一个实根。

记此实根为*x ，根据二分法误差估计公式有12()1*22k k k b a x x ++--≤=要使得近似解有3位有效数字，只需要有22111022k -+≤⨯从而可得6k ≥，即满足精度要求的二分次数为6次。

(2)将计算过程中数据填入表1.(中间过程填写到小数点后面3位)(12分，每个k x 得2分，其它空不计分)表1题1计算过程kk a k b kx 0 1.0 1.5 1.251 1.25 1.5 1.3752 1.25 1.375 1.3163 1.313 1.375 1.3494 1.313 1.344 1.3285 1.313 1.328 1.32061.3201.3281.324二.(10分)为了计算方程()3sin 2120f x x x =--=的根，某同学将()0f x =改写为14sin 23x x =+，并建立迭代公式114sin 23k k x x +=+。

请问此迭代公式在R 上是否全局收敛的吗？说明理由。

证明：(1)对任意的x R ∈，有11113()4sin 2,333x x R ϕ⎡⎤=+∈⊆⎢⎣⎦；(2)对任意的x R ∈，有22'()cos 2133x x ϕ=≤<;从而可知，迭代格式在R 上全局收敛。

三.(20分)设有方程3()10f x x x =--=，试回答下列问题：(1)确定方程3()10f x x x =--=实根的数目；(4分)解：由2'()31f x x =-可知函数3()1f x x x =--的单调递增区间是,33⎛⎛⎫-∞-⋃+∞ ⎪ ⎪⎝⎭⎝⎭，单调递减区间是,33⎛- ⎝⎭。

数值分析第二次大作业SOR最优松弛因子选取方法及SOR迭代法的改进

《数值分析》第二次大作业题目：SOR最优松弛因子选取方法及SOR迭代法的改进内容：1.SOR最优松弛因子选取方法2.SOR迭代法的改进（SSOR迭代法）3.SSOR迭代法的Matlab程序4.举例比较jacobi，Gauss-Seidel，SOR及SSOR 迭代法的收敛速度姓名：合肥工业大学学号：2011班级：信息与计算科学11-1班参考资料：1.《确定SOR最优松弛因子的一个实用算法》李春光等《计算力学学报》2.《数值分析与实验》，薛毅,北京工业大学出版社.3.《数值分析中的迭代法解线性方程组》，马云,科学出版社4.《非线性互补问题的改进超松弛迭代算法》，段班祥等,江西师范大学出版社5.《迭代法解线性方程组的收敛性比较》，郑亚敏,江西科学出版社.一、SOR最优松弛因子选取方法SOR迭代法迭代公式：x(k+1)i=(1-ω)xi+(k) bi-∑aijxjaii⎝j=1ω⎛ i-1(k+1)-j=i+1∑axijn(k)j⎫⎪ (i=1,2,..n.), ⎪⎭1.二分比较法将松弛因子1/2，ω的区间（1,2）进行二分，每个小区间的长度为ω去中间值3/2，按照SOR 迭代法迭代公式，求出跌代次数k,如果k不超过指定的发散常数，则可确定ω的值；否则将（1,2）四等分，每个区间长度为1/4,ω取各分点值，继续迭代，一般地，将1区间（1,2）二分M次，每次二分步长为，ω一次取取各分点值，2M按照SOR迭代法迭代公式，求出跌代次数k,如果k不超过指定的发散常数，则可确定的ω的值，这样总能找到一个不超过指定发散常数ω值。

2.逐步搜索法将1+ω的取值区间（1,2）进行M等分，ω分别取ω的值。

12M-1,1+,...,1+，通过迭代公式依次对同意精度要求求出迭代MMM次数k的值，并从中选出最优松弛因子3.黄金分割法依据黄金分割比的思想，通过计算机主动选取最优松弛因子的近似值，步骤如下a.对（1，2）区间进行第一次0.618的分割，区间边界a1=1,b1=2，在区间(a1,b1)分割出黄金点p1=a1+0.618(b1-a1)，进行SOR迭代法的迭代，求出迭代次数k的值，如果没有超过规定的发散常数，迭代结束，否则做步骤b。

数值分析第二次作业答案answer2

S4 = 0.11157238253891，S8 = 0.11157181325263。同学们根据自己理解计算 S4 ，S8 都可。复合梯形公式和复合 Simpson 公式的代码已重复多次，同学们自己整理。 3. 用 Simpson 公式计算积分误差为 |R(f )| = | − η ∈ (0, 1)。 4. 推导下列三种矩形求积公式： ∫b f (x)dx ∫a b f (x)dx ∫a b a f (x)dx = (b − a)f (a) + = (b − = (b −
14.7 53.63 从而 a = −7.855048，b = 22.25376。 2. 已知实验数据如下：。 xi 19 25 31 38
44
yi 19.0 32.3 49.0 73.3 97.8 用最小二乘法求形如 y = a + bx2 的经验公式。答案：两个待定常数，只能两个 φ。 φ0 ，φ1 也必须形如 y = a + bx2 。可设 φ0 = 1，φ1 = x2 。法方程为： ( 5 5327 )( a b ) = ( 271.4 369321.5 )
第三章函数逼近 1. 观测物体的直线运动，得出以下数据：时间 t(s) 0 0.9 1.9 3.0 3.9 5.0 距离 s(m) 0 求运动方程。 ( 10 φ0 = 1，φ1 = t。法方程为： 6 14.7 )( a b ) = ( 280 1078 )
6
1. 用 LU 分解及列主元高斯消去法解线性方程组 8 10 −7 0 1 x1 −3 2.099999 6 2 x 5.900001 2 = 5 5 − 1 5 − 1 x 3 x4 1 2 1 0 2 输出 Ax = b 中系数 A = LU 分解的矩阵 L 及 U ，解向量 x 及 det A；列主元法的行交换次序，解向量 x 及 det A；比较两种方法所得的结果。代码： A=[10,-7,0,1;-3,2.099999,6,2;5,-1,5,-1;2,1,0,2]; b=[8,5.900001,5,1]'; x=A\b;x(1) 结果：1.7764e-016 LU分解代码： A=[10,-7,0,1;-3,2.099999,6,2;5,-1,5,-1;2,1,0,2]; b=[8,5.900001,5,1]'; [m,n] = size(A); if m~=n, error('A matrix needs to be square'); end for i=1:n-1 pivot = A(i,i); if abs(pivot)<50*eps, error('zero pivot encountered'); end for k = i+1:n A(k,i) = A(k,i)/pivot; A(k,i+1:n) = A(k,i+1:n) - A(k,i)*A(i,i+1:n); end end 7

高等数值分析第二次实验作业(贾忠孝)

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 分解，构造如下形式的矩阵：

清华大学贾仲孝高等数值分析第二次作业

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. 证明：反证法。

高等数值分析第二次实验

一.1.基本的Arnolidi方法数值性态1.1.构造矩阵Aa1=linspace(60,300,500)b1=linspace(-300,300,500)U=zeros(1000)U=sparse(U)for i=1:500U(2*i-1,2*i-1)=a1(i)U(2*i,2*i)=a1(i)U(2*i-1,2*i)=b1(i)U(2*i,2*i-1)=-b1(i)endA=Q*U*Q' %Q为已经构造好的1000阶正交阵A的特征值均在坐标轴右半平面1.2. Arnoldi方法程序:b=randn(1000,1)V=[]M=40tol =10^-6x0=zeros(n,1)r=[]r0=b-A3*x0normr0=norm(r0)r(1)=normr0V(:,1)=r0/normr0m=1while m<Mw=A3*V(:,m)for i=1:mH(i,m)=V(:,i)'*ww=w-H(i,m)*V(:,i)endif norm(w)<10*epse1=[1;zeros(m-1,1)]ym=H(1:m,1:m)\(mormr0*e1')breakelseH(i+1,m)=norm(w)V(:,m+1)=w/H(i+1,m) ende1=[1;zeros(m-1,1)]ym=H(1:m,1:m)\(normr0*e1)emT=[zeros(m-1,1);1]r(m+1)=H(m+1,m)*abs(emT'*ym)m=m+1endr=log(r/norm(b))1.3收敛曲线以步数为横坐标，以相对残量的10为底的对数为纵坐标，得到下图的收敛曲线：可以看出，此时收敛曲线单调下降，有限步即能达到收敛要求。

2.基本的GMRES方法数值性态2.1构造A的方法同1.12.2用MATLAB自带的gmres函数进行求解，具体程序如下：b=randn(1000,1)tol=10^-18maxit=1restart=60r=[][x,flag,relres,iter,revec]=gmres(A3,b,restart,tol,maxit)for i=1:iter(2)r(1,i)=log(revec(i)/norm(b))endplot(1:iter(2) ,r,'-r')xlabel('step')ylabel('lg(morm(rk))')2.3以步数为横坐标，以相对残量的10为底的对数为纵坐标，得到下图的收敛曲线：收敛曲线单调下降，有限步即能达到收敛精度。

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

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)。

Answer：关于特征值均在右半平面的矩阵A：首先构造对角矩阵D，其中D矩阵的对角是由以下形式的2×2的子矩阵组成其对角元：，其中，为实数此时，S矩阵的特征值分别为a+bi和a-bi，这样D=diag(S1,S2,S3……S n)矩阵的特征值均分布在右半平面。

另矩阵A=Q T AQ，则A矩阵的特征值也均在右半平面，生成矩阵的过程代码如下所示：N=500 %生成的矩阵为2N阶A=zeros(2*N);%delta控制特征值分布的密集程度，即控制条件数大小，delta越大条件数越大delta=0.1;%构造D矩阵for j=1:NA(2*j-1,2*j-1)=N+j*delta;A(2*j-1,2*j)=-N-j*delta;A(2*j,2*j-1)=N+j*delta;A(2*j,2*j)=N+j*delta;endU = orth(rand(2*N,2*N));A = U'*A*U;首先进行观察基本的Arnoldi方法和GMRES方法的数值性态，考虑N=500，即矩阵为1000阶，取X0=zeros(N,1)，b=ones(N,1)，收敛准则为e=10-6;Arnoldi方法函数如下：function [xm,error,num]=Arnoldi(A,x0,b,e) N=size(A,1);r=b-A*x0;belta=norm(r);v=r/belta;V=v;j=0;while norm(r)>e & j<N j=j+1;num(j)=j;w=A*V(1:N,j);for i=1:jh(i,j)=V(1:N,i)'*w; w=w-h(i,j)*V(1:N,i);endh(j+1,j)=norm(w);if h(j+1,j) ~=0v=w/h(j+1,j);V=[V v];ende1=zeros(j,1);e1(1)=1;be1=belta*e1;try[L,U,S]=lu(h(1:j,1:j)); be1=S*be1;lym=LX(L,be1);ym=UX(U,lym);xm=x0+V(1:N,1:j)*ym; r=b-A*xm;enderror(j)=log10(norm(r)); endendGMRES方法的函数如下：function [xm,error,num] = GMRES(A,x0,b,e) %ARNOLDI Summary of this function goes here % Detailed explanation goes hereN=size(A,1);r=b-A*x0;belta=norm(r);v=r/belta;V=v;j=0;er=1000;while er > e & j<Nj=j+1;num(j)=j;w=A*V(1:N,j);for i=1:jh(i,j)=V(1:N,i)'*w;w=w-h(i,j)*V(1:N,i);endh(j+1,j)=norm(w);if h(j+1,j) ~=0v=w/h(j+1,j);V=[V v];endtry[Q,R]=qr(h(1:j+1,1:j)); gk=Q'*belta*eye(j+1,1); ym=minresYK(R,gk);xm=x0+V(1:N,1:j)*ym;er=gk(j+1);enderror(j)=log10(er);endend其中UX，LX和minresYK是自己编写的求解（上、下）三角矩阵的函数，代码如下：function yk=LX(TK,b)n=size(TK,1);yk=zeros(n,1);yk(1)=b(1)/TK(1,1);for i=2:nyk(i)=b(i);for j=1:i-1yk(i)=yk(i)-TK(i,j)*yk(j) ;endyk(i)=yk(i)/TK(i,i);endend function yk=UX(TK,b)n=size(TK,2);yk=zeros(n,1);yk(n)=b(n)/TK(n,n);for i=1:n-1k=n-i;yk(k)=b(k);for j=k+1:nyk(k)=yk(k)-TK(k,j)*yk(j);endyk(k)=yk(k)/TK(k,k);endendfunction yk=minresYK(TK,b)n=size(TK,2);yk=zeros(n,1);yk(n)=b(n)/TK(n,n);for i=1:n-1k=n-i;yk(k)=b(k);for j=k+1:nyk(k)=yk(k)-TK(k,j)*yk(j);endyk(k)=yk(k)/TK(k,k);endend两种方法的计算结果如下所示：Delta=1 即条件数较大 Delta=0.1 即条件数较小可以看到delta=0.1时两种方法收敛都很快，delta=1时，两种方法略有差别，GMRES 方法收敛速度略快一些，下面就去delta=1的情况进行第一题的讨论。

delta=1时，计算结果如下：Arnoldi 方法良好，Arnoldi 方法会有平台段，但是GMRES 方法的则平稳地收敛，这也是最后迭代次数GMRES 较Arnoldi 方法少一次的原因。

对于重启的Arnoldi 方法和GMRES 方法代码如下所示：A 、首先对之前的标准Arnoldi 方法和GMRES 方法进行修改，代码如下（m 为间隔m 步重启）： function [xm,error]=Arnoldi_m(A,x0,b,e,m) N=size(A,1); r=b-A*x0; belta=norm(r); v=r/belta; V=v; j=0; while abs(norm(r))>e & j<m j=j+1; w=A*V(1:N,j); for i=1:j h(i,j)=V(1:N,i)'*w; w=w-h(i,j)*V(1:N,i); end h(j+1,j)=norm(w); if h(j+1,j) ~=0 v=w/h(j+1,j); V=[V v]; end e1=zeros(j,1); e1(1)=1; be1=belta*e1; try [L,U,S]=lu(h(1:j,1:j)); be1=S*be1; lym=LX(L,be1); ym=UX(U,lym); xm=x0+V(1:N,1:j)*ym; r=b-A*xm; end end error=abs(norm(r)); end function [xm,error]=GMRES_m(A,x0,b,e,m) N=size(A,1); r=b-A*x0; belta=norm(r); v=r/belta; V=v; j=0; er=1000; while abs(er) > e & j<m j=j+1;迭代次数l g （|r m |）迭代次数l g （|r m |）w=A*V(1:N,j); for i=1:jh(i,j)=V(1:N,i)'*w; w=w-h(i,j)*V(1:N,i); endh(j+1,j)=norm(w); if h(j+1,j) ~=0 v=w/h(j+1,j); V=[V v]; endtry[Q,R]=qr(h(1:j+1,1:j)); gk=Q'*belta*eye(j+1,1); ym=minresYK(R,gk); xm=x0+V(1:N,1:j)*ym; er=gk(j+1); end error=abs(er); end endB 、调用上面更改后的代码，编写重启的算法代码如下：function[xm,err,num]=Arnoldim(A,x0,b,e,m) N=size(A,1); jmax=N/m; j=0; xm=x0; error=1000;while j<jmax & error>e j=j+1;[xm,error]=Arnoldi_m(A,xm,b,e,m); num(j)=j;err(j)=log10(error); end jlog10(error) endfunction [xm,err,num]=GMRESm(A,x0,b,e,m) N=size(A,1); jmax=N/m; j=0; xm=x0; error=1000;while j<jmax & error>e j=j+1;[xm,error]=GMRES_m(A,xm,b,e,m); num(j)=j;err(j)=log10(error); end jlog10(error) end计算结果如下所示（计算时取m=5，即每隔5步进行算法的重启）：delta=1 m=5（每隔5步重启） delta=1 m=10（每隔10步重启）重启次数l g （|r m |）重启次数l g （|r m |）delta=10 m=5（每隔5步重启） delta=10 m=10（每隔10步重启）当delta=100 m=5（每隔5步重启），结果如下：从结果可以看出：1. 重启算法的总迭代次数（重启次数×m ）要多于基本的算法；2. 重启的Arnoldi 方法收敛速度要比GMRES 方法快；3. 对于这个问题，重启算法也可以稳定地收敛，及每次重启后的x0都能有所改善；4. 当条件数增大时，明显重启次数会增加；5. 总的来说，对题中的矩阵A ，四种方法均能稳定地收敛。

T2. 对于1中的矩阵，将特征值进行平移，使得实部有正有负，和1的结果进行比较，方法的收敛速度会如何？基本的Arnoldi 算法有无峰点？若有，基本的GMRES 算法相应地会怎样？ Answer ：对A 的特征值进行平移，由于只要变换A 矩阵的实部，这里将上面的矩阵生成做如下改造：A=zeros(2*N); delta=0.1;s=1.5 %s 用于控制实部为负的特征值的个数，s 越大，实部为负的特征值个数越多。

for j=1:NA(2*j-1,2*j-1)=-s*delta+j*delta; A(2*j-1,2*j)=s*delta-j*delta; A(2*j,2*j-1)=-s*delta+j*delta; A(2*j,2*j)=-s*delta+j*delta; end重启次数l g （|r m |）重启次数l g （|r m |）重启次数l g （|r m |）U = orth(rand(2*N,2*N)); A = U'*A*U;基本算法的结果如下所示：2个实部为负的特征值10个实部为负的特征值100个实部为负的特征值200个实部为负的特征值500个实部为负的特征值900个实部为负的特征值迭代次数l g （|r m |）迭代次数l g （|r m |）迭代次数l g （|r m |）迭代次数l g （|r m |）迭代次数l g （|r m |）迭代次数l g （|r m |）迭代次数l g （|r m |）迭代次数l g （|r m |）990个实部为负的特征值998个实部为负的特征值2个实部为负的特征值，每10步重启 2个实部为负的特征值，每100步重启迭代次数l g （|r m |）重启次数l g （|r m |）重启次数l g （|r m |）10个实部为负的特征值，每10步重启10个实部为负的特征值，每100步重启100个实部为负的特征值，每10步重启 100个实部为负的特征值，每100步重启500个实部为负的特征值，每10步重启 500个实部为负的特征值，每100步重启998个实部为负的特征值，每10步重启 500个实部为负的特征值，每100步重启重启次数l g （|r m |）重启次数l g （|r m |）重启次数l g （|r m |）重启次数l g （|r m |）重启次数l g （|r m |）重启次数l g （|r m |）重启次数l g （|r m |）重启次数l g （|r m |）结论：1. 对于上面的矩阵，基本的Arnoldi 方法与GMRES 方法均能收敛，通知GMRES 方法收敛速度相对较快，而重启的算法则均不收敛； 2. 可以看到单特征值均匀地分布在左右两边平面时（即有一半的特征值实部为负，另一半实部为正），问题的病态程度最高；3. 问题越病态，重启的Arnoldi 算法得到的结果越不好，而重启的GMRES 对于结果基本就没有改进，这两种方法均不收敛；4. 基本的Arnoldi 算法有峰点，这是由于实部为负的特征值的存在，使得H 矩阵发生近似奇异导致发生近似中断而引起的，而GMRES 算法则相对稳定，当然这两种算法都能收敛。