gmres源程序matlab

GM （1,1）模型中的MATLAB 程序一、GM （1,1）模型的建立：（1）、一次累加生成序列的MA TLAB 命令：>> X0=[142 340 200 500 900 800 490 980 463 1100];>> X1(1)=X0(1)X1 =142>> for k=2:10X1(k)=X1(k-1)+X0(k)endX1 =142 482 682 1182 20822882 3372 4352 4815 5915（2）、由一次累加生成序列紧邻均值生成的)1(Z 的MA TLAB 命令：>> X0=[142 340 200 500 900 800 490 980 463 1100];>> X1(1)=X0(1);>> for k=2:10X1(k)=X1(k-1)+X0(k)Z(k)=(1/2)*(X1(k)+X1(k-1))EndX1 =142 482 682 1182 20822882 3372 4352 4815 5915Z =1.0e+003 *0 0.3120 0.5820 0.9320 1.6320 2.4820 3.12703.86204.58355.3650（3）、GM （1，1）的灰微分方程模型为：b k aZ k X =+)()()1()0(。

设∧α为待估计参数向量，⎥⎦⎤⎢⎣⎡=∧b a α。

利用最小二乘法得到Y B B B ')'(1-∧=α，MA TLAB 程序如下：>> B=[-Z(2:10)',ones(9,1)];>> Y=(X0(2:10))';>> alfa=inv(B'*B)*B'*Yalfa =-0.1062371.6018（4）、GM （1，1）的灰微分方程模型 b k aZ k X =+)()()1()0(的时间相应序列为：ab e a b X k X ak +⋅-=+-∧))1(()1()0()1( 由.6018.371,1062.0=-=b a 令.)1(,)0(u X v a b u -== 计算得到1.3499-=u ， 1.3641=v 。

function [x，flag,relres,iter,resvec] = gmres(A,b,restart,tol,maxit,M1，M2,x,varargin)%GMRES Generalized Minimum Residual Method。

％X = GMRES(A,B) attempts to solve the system of linear equations A*X = B% for X。

The N-by—N coefficient matrix A must be square and the right％hand side column vector B must have length N。

This uses the unrestarted% method with MIN(N,10) total iterations.%％X = GMRES（AFUN,B）accepts a function handle AFUN instead of the matrix% A。

AFUN(X）accepts a vector input X and returns the matrix-vector% product A*X. In all of the following syntaxes, you can replace A by% AFUN.％％X = GMRES(A，B,RESTART）restarts the method every RESTART iterations.％If RESTART is N or [］then GMRES uses the unrestarted method as above.％% X = GMRES(A，B,RESTART,TOL）specifies the tolerance of the method。

If% TOL is ［] then GMRES uses the default, 1e-6.％% X = GMRES（A，B,RESTART，TOL，MAXIT）specifies the maximum number of outer % iterations。

2.1set 与get 函数 (1)2.2callback函数 (2)2.3元胞数组 (4)2.4结构数组 (6)2.5矩阵操作 (9)2.6字符串操作 (13)2.7判断函数使用大全 (16)2.11打开外部程序 (21)2.11程序运行时间 (22)2.14动画 (23)2.12动画 (24)2.23显示多行内容 (26)2.24 uitable 使用 (26)2.27鼠标操作 (27)2.28键盘操作 (27)2.32粘贴板 (28)2.1set 与get 函数set(edit_handle,'String','my value!'); %String为Edit控件的属性%%%2.1-1%创建figure对象hfig=figure(1);%创建坐标轴对象，指定其父对象为figure 1haxes1=axes('parent',hfig);prop.Color='b';prop.FontSize=12;set(haxes1,prop);%%%2.1-2hfig=figure(1);%查询其Units属性值get(hfig,'units')%其Units属性值为pixels（像素）% ans=% pixels%%%2.1-3%figure的Pointer属性标识了鼠标指针的形状set(gcf,'pointer');% 返回值为：[ crosshair | fullcrosshair | {arrow} | ibeam | watch | topl | topr | botl | botr | left | top | right | bottom | circle | cross | fleur | custom | hand ]%%%2.1-4%首先取得标识电脑屏幕大小的度量单位get(0,'units')% ans =% pixels%取得屏幕的尺寸get(0,'screensize')% ans =% 1 1 1280 8002.2callback函数%定义M文件的主函数名称为DefineCallback，不带输入和输出参数function DefineCallbackhFig= figure('units','normalize',...'position',[0.4 0.4 0.3 0.2]);%在窗口中创建按钮控件，并定义其Callback属性uicontrol('parent',hFig,...'style','pushbutton',...'String','Execute Callback',...'units','normalize',...'position',[0.4 0.4 0.3 0.2],...'callback',['figure;',...'x = 0:pi/20:2*pi;',...'y = sin(x);',...'plot(x,y);']);%定义M文件的主函数名称为DefineCallback，不带输入和输出参数function DefineCallback%创建界面窗口hFig= figure('units','normalize',...'position',[0.4 0.4 0.3 0.2]);%在窗口中创建按钮控件hpush=uicontrol('parent',hFig,...'style','pushbutton',...'String','Execute Callback',...'units','normalize',...'position',[0.4 0.4 0.3 0.2]);%设置按钮的Callback属性set(hpush,'callback',@mycallback);%定义回调函数为子函数function mycallback(hobj,event)figure;x = 0:pi/20:2*pi;y = sin(x);plot(x,y);2.3元胞数组a={'hello' [1 2 3;4 5 6];1 {'1''2'}}a ='hello' [2x3 double][ 1] {1x2 cell }%示例2：将元胞数组a中的元胞逐一赋值>> a{1,1}='hello';a{1,2}=[1 2 3;4 5 6];a{2,1}=1;a{2,2}={'1' '2'};>> aa ='hello' [2x3 double][ 1] {1x2 cell }%示例3：使用cell函数来创建元胞数组%生成2x3的元素为空数组的元胞数组>> a=cell(2,3)a =[] [] [][] [] []%示例4：判断数组A是否为元胞数组%定义一个元胞数组A>> A={1 2 3};%判断A是否为元胞数组，如果为元胞数组，则函数>> tf = iscell(A)tf =1%示例5：显示元胞数组C中的内容>> clear>> C={'Smith' [1 2;3 4] [12]};%直接显示元胞数组C中的内容>> celldisp(C)C{1} =SmithC{2} =1 23 4C{3} =12%显示元胞数组C中的内容，数组的名称用cellcontent代替>> celldisp(C,'cellcontent')cellcontent{1} =Smithcellcontent{2} =1 23 4cellcontent{3} =12%示例6：将字符数组转换为元胞数组>> S = ['abc '; 'defg'; 'hi m'];>> cellstr(S)ans ='abc'%原先abc后面的空格被清除'defg''hi m'%i和m之间的空格仍然保留%示例7：显示元胞数组S中的内容（包括空格和字符）>> S = {'abc ', 'defg','hi m'};>> cellplot(S)%示例8：将数字数组A按行或按列转换为元胞数组%A是4x3的数组>> A=[1 2 3;4 5 6;7 8 9;10 11 12];%把A的每一列转换为一个元胞，得到的C是1×3的元胞数组>> C=num2cell(A,1)C =[4x1 double] [4x1 double] [4x1 double] %把A的每一行转换为一个元胞，得到的C是4×1的元胞数组>> C=num2cell(A,2)C =[1x3 double][1x3 double][1x3 double][1x3 double]2.4结构数组%示例1：使用直接法来创建结构数组>> A(1).name = 'Pat';A(1).number = 176554;A(2).name = 'Tony';A(2).number = 901325;>> AA =1x2 struct array with fields:namenumber%示例2：利用struct函数来创建结构数组>> A(1)=struct('name','Pat','number',176554);A(2)=struct('name','Tony','number',901325);>> AA =1x2 struct array with fields:namenumber%示例3：使用deal函数来得到结构体中各结构域的值%定义结构数组A>> = 'Pat'; A.number = 176554;A(2).name = 'Tony';A(2).number = 901325;%得到结构数组中所有name结构域的数据>>[name1,name2] = deal(A(:).name)name1 =Patname2 =Tony%示例4：使用getfield函数来取得结构体中结构域的值%定义mystr结构数组>> mystr(1,1).name = 'alice';mystr(1,1).ID = 0;mystr(2,1).name = 'gertrude';mystr(2,1).ID = 1;%取得mystr(2,1)的结构域name的值>> f = getfield(mystr, {2,1}, 'name')f =gertrude%示例5：删除结构数组中的指定结构域%定义结构数组s>> s.field1=[1 2 3];s.field2='string';s.field3={1 2 3;4 5 6}; %删除结构域field1>> s=rmfield(s,'field1')s =field2: 'string'field3: {2x3 cell}%删除结构域'field2','field3'>> s=rmfield(s,{'field2','field3'})s =field1: [1 2 3]%示例6：%定义结构数组s>> s.field1=[1 2 3];s.field2='string';s.field3={1 2 3;4 5 6}; >> ss =field1: [1 2 3]field2: 'string'field3: {2x3 cell}%将结构数组转换为元胞数组>> c=struct2cell(s)c =[1x3 double]'string'{2x3 cell }%示例7：>>c = {'birch', 'betula', 65; 'maple', 'acer', 50}c ='birch''betula' [65]'maple''acer' [50]>>fields = {'name', 'genus', 'height'}; %fields包含struct中的结构域名>>s = cell2struct(c, fields, 2); %dim=2表示把c中的各行转换为struct数组s =2x1 struct array with fields:namegenusheight>> s(1)ans =name: 'birch'genus: 'betula'height: 65>> s(2)ans =name: 'maple'genus: 'acer'height: 50>> fields = {'field1', 'field2'};>> s = cell2struct(c, fields, 1); %dim=1表示把c中的各列转换为struct数组>> s(1)ans =field1: 'birch'field2: 'maple'>> s(2)ans =field1: 'betula'field2: 'acer'>> s(3)ans =field1:65field2:502.5矩阵操作%示例1： find函数的使用方法。

1、实验目的（1）学习MATLAB的各种二维绘图；（2）学习MATLAB的三维绘图；（3）M ATLAB的绘图修饰（多种、图形注释、绘图颜色、色图矩阵）。

2、实验内容（1）基本二维绘图1）向量绘图x=0:2*pi/100:2*pi;y1=sin(2*x);y2=cos(2*x);plot(x,y1)plot(x,y2)plot(x,y1,x,y2)plot(x,y1);hold on;plot(x,y2);hold off;plot(x',[y1',y2'])设定颜色与线型plot(x,y1,'c:',x,y2,'wo')多窗口绘图figure(1);plot(x,y1);figure(2);plot(x,y2);子图绘图subplot(221);plot(x,y1);subplot(222);plot(x,y2);subplot(223);plot(x,y1,x,y1+y2);subplot(224);plot(x,y2,x,y1-y2);复变函数绘图w=0.01:0.01:10;G=1./(1+2*w*i);subplot(121);plot(G);subplot(122);plot(real(G),imag(G));插值绘图x=0:2*pi/8:2*pi;y=sin(x); plot(x,y,'o');hold on;xi=0:2*pi/100:2*pi;yi=spline(x,y,xi);plot(xi,yi,'m');反白绘图与绘图背景颜色设定whitebgwhitebg('b')whitebg('k')2）函数绘图fplot('sin',[0 4*pi])f='sin(x)';fplot(f,[0 4*pi])fplot('sin(1/x)',[0.01 0.1],1e-3)fplot('[tan(x),sin(x),cos(x)]',[-2*pi,2*pi,-2*pi,2*pi])3)符号函数快捷绘图syms xf=exp(-0.5*x)*sin(x);ezplot(f,[0,10])（2）多种二维绘图1）半对数绘图（频率特性绘图）w=logspace(-1,1);g=20*log10(1./(1+2*w*i));p=angle(1./(1+2*w*i))*180/pi;subplot(211);semilogx(w,g);grid;subplot(212);semilogx(w,p);grid;2）极坐标绘图t=0:2*pi/180:2*pi;mo=cos(2*t);polar(t,mo);t=0:2*pi/8:2*pi; y=sin(t);bar(t,y);t=0:2*pi/8:2*pi; y=sin(t);stem(t,y);t=0:2*pi/8:2*pi; y=sin(t);stairs(t,y);t=-pi:pi/200:pi;comet(t,tan(sin(t))-sin(tan(t)));（3）图形注释y1=dsolve('D2u+2*Du+10*u=0','Du(0)=1,u(0)=0','x');y2=dsolve('D2u+2*Du+10*u=1','Du(0)=1,u(0)=0','x');ezplot(y1,[0,5]);hold on;ezplot(y2,[0,5]);axis([0,5,-0.1,0.2]);title('二阶系统时间响应');axis([0,5,-0.1,0.2]);title('二阶系统时间响应');xlabel('时间t');xlabel('时间t');title('二阶系统时间响应');xlabel('时间t');ylabel('响应幅值y');gtext('零输入响应');gtext('零状态响应');grid;hold off;（4）三维绘图1）三维线图t=0:pi/50:10*pi;plot3(sin(t),cos(t),t);comet3(sin(t),cos(t),t);Z2=[1 1;1 -1];Z4=[Z2 Z2;Z2 -Z2]; Z8=[Z4 Z4;Z4 -Z4]; mesh(Z8);x=-4:0.5:4;y=x; [X,Y]=meshgrid(x,y); Z=X.^2-Y.^2;mesh(X,Y,Z);4)圆锥面网线图t1=0:0.1:0.9;t2=1:0.1:2;r=[t1 -t2+2];[X,Y,Z]=cylinder(r,40);mesh(X,Y,Z);5）视角修饰t1=0:0.1:0.9;t2=1:0.1:2;r=[t1 -t2+2];[X,Y,Z]=cylinder(r,30);mesh(X,Y,Z);subplot(2,2,1);mesh(X,Y,Z);view(0,0);subplot(2,2,2);mesh(X,Y,Z);view(-20,20);subplot(2,2,3);mesh(X,Y,Z);view(-30,30);subplot(2,2,4);mesh(X,Y,Z);view(-40,40);6）MATLAB暖色（hot）函数色图peaks(20);z = 3*(1-x).^2.*exp(-(x.^2) - (y+1).^2) ...- 10*(x/5 - x.^3 - y.^5).*exp(-x.^2-y.^2) ...- 1/3*exp(-(x+1).^2 - y.^2)axis('off')colormap(hot);colorbar('horiz');7）光照修饰surfl(peaks(20));colormap(hot);shading interp;8）表面修饰subplot(131);peaks(20);shading flat;z = 3*(1-x).^2.*exp(-(x.^2) - (y+1).^2) ...- 10*(x/5 - x.^3 - y.^5).*exp(-x.^2-y.^2) ...- 1/3*exp(-(x+1).^2 - y.^2)subplot(132);peaks(20);shading interp;z = 3*(1-x).^2.*exp(-(x.^2) - (y+1).^2) ... - 10*(x/5 - x.^3 - y.^5).*exp(-x.^2-y.^2) ... - 1/3*exp(-(x+1).^2 - y.^2)subplot(133);peaks(20);shading faceted;z = 3*(1-x).^2.*exp(-(x.^2) - (y+1).^2) ... - 10*(x/5 - x.^3 - y.^5).*exp(-x.^2-y.^2) ... - 1/3*exp(-(x+1).^2 - y.^2)9）透视与消隐P=peaks(30);subplot(121);mesh(P);hidden off;subplot(122);mesh(P);hidden on;10）裁剪修饰P=peaks(30);P(20:23,9:15)=NaN*ones(4,7);subplot(121);meshz(P);P=peaks(30);subplot(121);mesh(P);hidden off;subplot(122);mesh(P);hidden on;P=peaks(30);P(20:23,9:15)=NaN*ones(4,7);subplot(121);meshz(P);subplot(122);meshc(P);11）水线修饰waterfall(peaks(20))colormap([1 0 0])12)等高线修饰contour(peaks(20),6)contour3(peaks(20),10)clabel( contour(peaks(20),4))clabel( contour3(peaks(20),3))。

以下是matlab源程序=========================================== %%%%%%%%%%%%%%人造地震波程序%%%%%%%%%%%%%% 本程序利用反应谱转换为功率谱方法模拟人造地震波% 变量说明：% s0: 谱强度因子dt: 地震波采样时间间隔% kao: 地震波总时间长度% Ts: 平稳地震动的持续时间% t1,t2,C: 控制主震段首末时间和衰减快慢的参数% keceg,omigag: 地表覆盖层的阻尼比和卓越频率% omigah: 反映基岩特性的谱参数%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%xg(1:1500,1)=zeros(1500,1);dt=0.02;kao=21.9;t1=1.5;t2=15;c=0.18;Domiga=2*pi*(25-1/6)/100;Tg=0.25;kece=0.05;alfa=0.9;for i=1:100omiga=(2*pi/6)+Domiga*(i-0.5);T=2*pi/omiga;if T<0.1r=5.5*alfa*T+0.45*alfa;elseif (T>=0.1&T<=Tg)r=alfa;elser=alfa*(Tg/T)^0.9;endgama=omiga/pi;M1=(t2-0.3*t1)/kao;M2=(exp(-2*kece*omiga*t2)-exp(-2*kece*omiga*t1))/(2*kece*omiga*kao);M3=(1-exp(-c*(kao+0.5*t1-t2)))/(c*kao);M=M1+M2+M3;Gomiga=4*kece*r^2/(pi*omiga*M*(sqrt(2*log(gama*kao))+0.5772/sqrt(2*log(gama*kao)))^2); ak=2*sqrt(Gomiga*Domiga);qrand=rand(1)*2*pi;for j=1:1500xg(j)=xg(j)+ak*sin(omiga*j*dt+qrand);endendfor i=1:1500if i*dt<t1fai=(i*dt/t1)^2;elseif (i*dt>=t1&i*dt<t2)fai=1;elsefai=exp(-c*(i*dt-t2));endxg(i)=fai*xg(i);endsubplot(2,1,1);plot(xg);%xg=xg./(max(abs(xg)));save ren xg -ascii;fans=fft(xg);fans=fans.*conj(fans);subplot(2,1,2);plot(fans);。

正弦波的源程序：(一)，用到的函数1,f2t函数function x=f2t(X)global dt df t f T N%x=f2t(X)%x为时域的取样值矢量%X为x的傅氏变换%X与x长度相同并为2的整幂%本函数需要一个全局变量dt(时域取样间隔) X=[X(N/2+1:N),X(1:N/2)];x=ifft(X)/dt;end2,t2f函数。

function X=t2f(x)global dt df N t f T%X=t2f(x)%x为时域的取样值矢量%X为x的傅氏变换%X与x长度相同,并为2的整幂。

%本函数需要一个全局变量dt(时域取样间隔) H=fft(x);X=[H(N/2+1:N),H(1:N/2)]*dt;end(二),主程序。

1，%（1）绘出正弦信号波形及频谱global dt df t f Nclose allk=input('取样点数＝2^k, k取10左右');if isempty(k), k=10; endf0=input('f0=取1(kz)左右');if isempty(f0), f0=1; endN=2^k;dt=0.01; %msdf=1/(N*dt); %KHzT=N*dt; %截短时间Bs=N*df/2; %系统带宽f=[-Bs+df/2:df:Bs]; %频域横坐标t=[-T/2+dt/2:dt:T/2]; %时域横坐标s=sin(2*pi*f0*t); %输入的正弦信号S=t2f(s); %S是s的傅氏变换a=f2t(S); %a是S的傅氏反变换a=real(a);as=abs(S);subplot(2,1,1) %输出的频谱plot(f,as,'b');gridaxis([-2*f0,+2*f0,min(as),max(as)])xlabel('f (KHz)')ylabel('|S(f)| (V/KHz)') %figure(2)subplot(2,1,2)plot(t,a,'black') %输出信号波形画图gridaxis([-2/f0,+2/f0,-1.5,1.5])xlabel('t(ms)')ylabel('a(t)(V)')gtext('频谱图')最佳基带系统的源程序：(一)，用到的函数f2t函数和t2f函数。

第二章习题参考答案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 =是一种向量范数。

gmresGeneralized minimum residual method (with restarts)expand all in page Syntaxx = gmres(A,b)gmres(A,b,restart)gmres(A,b,restart,tol)gmres(A,b,restart,tol,maxit)gmres(A,b,restart,tol,maxit,M)gmres(A,b,restart,tol,maxit,M1,M2)gmres(A,b,restart,tol,maxit,M1,M2,x0)[x,flag] = gmres(A,b,...)[x,flag,relres] = gmres(A,b,...)[x,flag,relres,iter] = gmres(A,b,...)[x,flag,relres,iter,resvec] = gmres(A,b,...)Descriptionx = gmres(A,b) attempts to solve the system of linear equations A*x = b for x. The n-by-n coefficient matrix A must be square and should be large and sparse. The column vector b must have length n. A can be a function handle, afun, such that afun(x) returns A*x. For this syntax, gmres does not restart; the maximum number of iterations is min(n,10).Parameterizing Functions explains how to provide additional parameters to the function afun, as well as the preconditioner function mfun described below, if necessary.If gmres converges, a message to that effect is displayed. If gmres fails to converge after the maximum number of iterations or halts for any reason, a warning message is printed displaying the relative residual norm(b-A*x)/norm(b) and the iteration number at which the method stopped or failed. gmres(A,b,restart) restarts the method every restart inner iterations. The maximum number of outer iterations is min(n/restart,10). The maximum number of total iterationsis restart*min(n/restart,10). If restart is n or [], then gmres does not restart and the maximum number of total iterations is min(n,10).gmres(A,b,restart,tol) specifies the tolerance of the method. If tol is [], then gmres uses the default, 1e-6.gmres(A,b,restart,tol,maxit) specifies the maximum number of outer iterations, i.e., the total number of iterations does not exceed restart*maxit. If maxit is [] then gmres uses the default, min(n/restart,10). If restart is n or [], then the maximum number of total iterationsis maxit (instead of restart*maxit).gmres(A,b,restart,tol,maxit,M) and gmres(A,b,restart,tol,maxit,M1,M2) use preconditioner M or M = M1*M2and effectively solve the system inv(M)*A*x = inv(M)*b for x.If M is [] then gmres applies no preconditioner. M can be a function handle mfun suchthat mfun(x) returns M\x.gmres(A,b,restart,tol,maxit,M1,M2,x0) specifies the first initial guess. If x0 is [],then gmres uses the default, an all-zero vector.[x,flag] = gmres(A,b,...) also returns a convergence flag:flag = 0gmres converged to the desired tolerance tol within maxit outer iterations.flag = 1gmres iterated maxit times but did not converge.flag = 2Preconditioner M was ill-conditioned.flag = 3gmres stagnated. (Two consecutive iterates were the same.)Whenever flag is not 0, the solution x returned is that with minimal norm residual computed over all the iterations. No messages are displayed if the flag output is specified.[x,flag,relres] = gmres(A,b,...) also returns the relative residual norm(b-A*x)/norm(b).If flag is 0, relres <= tol. The third output, relres, is the relative residual of the preconditioned system.[x,flag,relres,iter] = gmres(A,b,...) also returns both the outer and inner iterationnumbers at which x was computed, where 0 <= iter(1) <= maxit and 0 <= iter(2) <= restart.[x,flag,relres,iter,resvec] = gmres(A,b,...) also returns a vector of the residual norms at each inner iteration. These are the residual norms for the preconditioned system.ExamplesUsing gmres with a Matrix InputA = gallery('wilk',21);b = sum(A,2);tol = 1e-12;maxit = 15;M1 = diag([10:-1:1 1 1:10]);x = gmres(A,b,10,tol,maxit,M1);displays the following message:gmres(10) converged at outer iteration 2 (inner iteration 9) toa solution with relative residual 3.3e-013Using gmres with a Function HandleThis example replaces the matrix A in the previous example with a handle to a matrix-vector productfunction afun, and the preconditioner M1 with a handle to a backsolve function mfun. The example iscontained in a function run_gmres that∙Calls gmres with the function handle @afun as its first argument.∙Contains afun and mfun as nested functions, so that all variables in run_gmres are available to afun and mfun.The following shows the code for run_gmres:function x1 = run_gmresn = 21;b = afun(ones(n,1));tol = 1e-12; maxit = 15;x1 = gmres(@afun,b,10,tol,maxit,@mfun);function y = afun(x)y = [0; x(1:n-1)] + ...[((n-1)/2:-1:0)'; (1:(n-1)/2)'].*x + ...[x(2:n); 0];endfunction y = mfun(r)y = r ./ [((n-1)/2:-1:1)'; 1; (1:(n-1)/2)'];endendWhen you enterx1 = run_gmres;MATLAB® software displays the messagegmres(10) converged at outer iteration 2 (inner iteration 10) to a solution with relative residual 1.1e-013.Using a Preconditioner without RestartThis example demonstrates the use of a preconditioner without restarting gmres.Load west0479, a real 479-by-479 nonsymmetric sparse matrix.load west0479;A = west0479;Set the tolerance and maximum number of iterations.tol = 1e-12; maxit = 20;Define b so that the true solution is a vector of all ones.b = full(sum(A,2));[x0,fl0,rr0,it0,rv0] = gmres(A,b,[],tol,maxit);fl0 is 1 because gmres does not converge to the requested tolerance 1e-12 within the requested 20 iterations. The best approximate solution that gmres returns is the last one (as indicated by it0(2) = 20). MATLAB stores the residual history in rv0.Plot the behavior of gmres.semilogy(0:maxit,rv0/norm(b),'-o');xlabel('Iteration number');ylabel('Relative residual');The plot shows that the solution converges slowly. A preconditioner may improve the outcome.Use ilu to form the preconditioner, since A is nonsymmetric.[L,U] = ilu(A,struct('type','ilutp','droptol',1e-5));Error using iluThere is a pivot equal to zero. Consider decreasingthe drop tolerance or consider using the 'udiag' option.Note MATLAB cannot construct the incomplete LU as it would result in a singular factor, which is useless as a preconditioner.As indicated by the error message, try again with a reduced drop tolerance.[L,U] = ilu(A,struct('type','ilutp','droptol',1e-6));[x1,fl1,rr1,it1,rv1] = gmres(A,b,[],tol,maxit,L,U);fl1 is 0 because gmres drives the relative residual to 9.5436e-14 (the value of rr1). The relative residual is less than the prescribed tolerance of 1e-12 at the sixth iteration (the value of it1(2)) when preconditioned by the incomplete LU factorization with a drop tolerance of 1e-6. The output, rv1(1), is norm(M\b), where M = L*U. The output, rv1(7), is norm(U\(L\(b-A*x1))).Follow the progress of gmres by plotting the relative residuals at each iteration starting from the initial estimate (iterate number 0).semilogy(0:it1(2),rv1/norm(b),'-o');xlabel('Iteration number');ylabel('Relative residual');Using a Preconditioner with RestartThis example demonstrates the use of a preconditioner with restarted gmres.Load west0479, a real 479-by-479 nonsymmetric sparse matrix.load west0479;A = west0479;Define b so that the true solution is a vector of all ones.b = full(sum(A,2));Construct an incomplete LU preconditioner as in the previous example.[L,U] = ilu(A,struct('type','ilutp','droptol',1e-6));The benefit to using restarted gmres is to limit the amount of memory required to execute the method. Without restart, gmres requires maxit vectors of storage to keep the basis of the Krylov subspace. Also,gmres must orthogonalize against all of the previous vectors at each step. Restarting limits the amount of workspace used and the amount of work done per outer iteration. Note that even though preconditioned gmres converged in six iterations above, the algorithm allowed for as many as twenty basis vectors and therefore, allocated all of that space up front.Execute gmres(3), gmres(4), and gmres(5)tol = 1e-12; maxit = 20;re3 = 3;[x3,fl3,rr3,it3,rv3] = gmres(A,b,re3,tol,maxit,L,U);re4 = 4;[x4,fl4,rr4,it4,rv4] = gmres(A,b,re4,tol,maxit,L,U);re5 = 5;[x5,fl5,rr5,it5,rv5] = gmres(A,b,re5,tol,maxit,L,U);fl3, fl4, and fl5 are all 0 because in each case restarted gmres drives the relative residual to less than the prescribed tolerance of 1e-12.The following plots show the convergence histories of eachrestarted gmres method. gmres(3) converges at outer iteration 5, inner iteration 3 (it3 = [5, 3]) which would be the same as outer iteration 6, inner iteration 0, hence the marking of 6 on the final tick mark.figuresemilogy(1:1/3:6,rv3/norm(b),'-o');set(gca,'XTick',[1:1/3:6],...'XTickLabel',...['1';' ';' ';'2';' ';' ';'3';' ';' ';'4';' ';' ';'5';' ';' ';'6';])title('gmres(3)')xlabel('Iteration number');ylabel('Relative residual');figuresemilogy(1:1/4:3,rv4/norm(b),'-o');set(gca,'XTick',[1:1/4:3],...'XTickLabel',['1';' ';' ';' ';'2';' ';' ';' ';'3'])title('gmres(4)')xlabel('Iteration number');ylabel('Relative residual');figuresemilogy(1:1/5:2.8,rv5/norm(b),'-o');set(gca,'XTick',[1:1/5:2.8],...'XTickLabel',['1';' ';' ';' ';' ';'2';' ';' ';' ';' ']) title('gmres(5)')xlabel('Iteration number');ylabel('Relative residual');ReferencesBarrett, R., M. Berry, T. F. Chan, et al., Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods, SIAM, Philadelphia, 1994.Saad, Youcef and Martin H. Schultz, "GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems," SIAM J. Sci. Stat. Comput., July 1986, Vol. 7, No. 3, pp. 856-869.。

GMRES在matlab中的描述参考文献：MATLAB函数库查询辞典徐东艳孟晓刚编著[函数描述]X=gmres(A,b)视图求解线性方程组A*x=b的解x。

nXn的稀疏矩阵A必须是方程且应是大型稀疏矩阵。

列向量b的长度必须为n。

参数A可以是一个函数afun以使得afun(x)返回A*x，对于这一语法格式，gmres并不重新启动，迭代的最大次数为min(n,10)。

如果gmres收敛，则显示这一结果的信息。

如果gmres在最大迭代步后没有收敛或者由于某种原因而停止，则打印警告信息，显示相对残差范数norm(b-A*x)/norm(b)并且显示该方法停止或者失效时的迭代步数。

Gmres(A,b,restart)在每一次重新开始内部迭代时才重新开始方法。

外部迭代的最大数目为min(n/restart,10),总迭代的最大步数为restart*min(n/restart,10).如果restart为n或者[]，则gmres并不重新开始，则总迭代的最大数目为min(n，10)Gmres(A,b,restart,tol)指定方法的误差。

如果tol为[]，则函数gmres将使用默认值1e-6.Gmres(A,b,restart,tol,maxit)指定外部迭代的最大次数，也就是说，迭代的总次数不超过restart*matrix。

如果参数maxit为[]，则矩阵gmres使用默认值min(n/restart,10);如果参数restart为n或者[],则总迭代的最大次数为maxit（而非restart*matxit）。

Gmres(A,b,restart,tol,maxit,M)和Gmres(A,b,restart,tol,maxit,M1,M2)使用预处理矩阵M或者M=M1*M2有效地求解方程组inv(M)*A*x=inv(M)*b.如果M为[],则函数gmres不适应预处理矩阵，M可以是一个函数，它返回M\xGmres(A,b,restart,tol,maxit,M1,M2,x0)指定初始的猜测值。

GMRES在电力系统潮流计算中的应用研究作者：高成来源：《科学与技术》2018年第23期摘要：本文阐述一种基于matlab的GMRES实现方法，主要针对于大规模跨区域电力系统分析中潮流计算的修正方程，在求解这个高维数超稀疏的线性方程组中避免了采用直接法进行求解，满足大型电力系统潮流计算的需要。

该方法采用基于Krylov子空间的Arnoldi过程生成单位正交化基和系数矩陣，然后采用基于Given变换的回代法求解最小二乘问题求解向量，结合重启动技术，极大地提高计算的效率和收敛性。

该方法具有算法稳定、简单及通用性强等特点，能够为后续电力工作者应用于潮流计算和电力系统分析中进一步扩展应用。

算例结果验证文中实现的正确性和有效性。

关键词：Matlab GMRES方法潮流计算修正方程大型稀疏线性方程组引言所谓电力系统潮流，是指系统中所有运行参数的总体，包括各个母线电压的大小和相位、各个发电机和负荷的功率及电流，以及各个变压器和线路等元件所通过的功率、电流和其损耗。

潮流问题是电力系统分析的基础和核心。

随着我国经济的发展，电力工业发展很快，而且向着大电网、智能化的方向发展，根据我国制定“西电东送、南北互供、全国联网”的战略构想，已基本上实现除西藏、新疆、海南、台湾外全国联网。

Krylov子空间法是二十世纪十大算法之一，具有存储量少，计算量小且易于并行等优点，非常适合并行求解高维稀疏线性方程组。

根据不同的极小化方法Krylov子空间法主要有CG （Conjugate Gradient）法、MINRES（Minimal Residual）法、GMRES 法（Generalised Minimal Residual）、BiCG 法（Bi-Conjugate Gradient）、CGS法（Conjugate Gradient Squared）。

其中，GMRES适合于求解大型非对称线性方程组，结合预条件处理技术的GMRES法具有良好的收敛特性和较高的数值稳定性，同时采用重启动技术可以有效的降低对内存占有的需求。