关于MATLAB误差不变的问题

Question

我在验证其他人的方法时，出现了这样的问题，人家文章中的结果是收敛的，但是我无论怎么改都不收敛，所以想麻烦大家帮我一起看看，那篇论文的基本思路是用HSS解AX+XB=C，即转化为解下面的两个方程：（图1）其中H(A)=1/2(A+A')，S(A)=1/2(A-A')，H(B)与S(B)同理。然后我写的程序如下：%用HSS解Sylester方程,其中内迭代用MATLAB自带的lyap函数function [ite_num re] = HSS_2(A,B,alpha,beta,tol,max_num)n=size(A,2);m=size(B,2);I=eye(n);E=eye(m);F=eye(m,n);X0 = zeros(m,n);H_A=(1/2)*(A+A');S_A=(1/2)*(A-A');H_B=(1/2)*(B+B');S_B=(1/2)*(B-B');A1=alpha*I+H_A;A2=alpha*I+S_A;A3=alpha*I-S_A;A4=alpha*I-H_A;B1=beta*E+H_B;B2=beta*E+S_B;B3=beta*E-S_B;B4=beta*E-H_B;R_k = norm(F-A*X0-X0*B,'fro');eps_k=1;R_0=R_k;re(1) = eps_k;ite_num = 1;while ite_num < max_num && eps_k > tol C1=A3*X0+X0*B3+F; X1=lyap(A1,B1,C1);%计算第二个式子 C2=A4*X1+X1*B4+F; X0=lyap(A2,B2,C2); R_k = norm(F-A*X0-X0*B,'fro'); eps_k=R_k/R_0; ite_num = ite_num + 1; re(ite_num) = eps_k;endite_num = ite_num - 1;主程序是：clear allclcclose alltol = 1e-6; max_num = 200;A=SPD2(16,0.1);B=A;n=size(A,2);m=size(B,2);for k = 1:20 a(k)=0.1*k; b(k)=a(k); [ite(k) re] = HSS_2(A,B,a(k),b(k),tol,max_num);end %Select the alpha corresponding to the minimum iteration number [~,num] = min(ite);alpha = a(num);beta = b(num);[ite_num re] = HSS_2(A,B,alpha,beta,tol,max_num);disp([' '])disp([ 'Output result'])disp([' '])disp([' The optimum acceleration parameter a: ' num2str(alpha) ])disp([' '])disp([' The optimum acceleration parameter b: ' num2str(beta) ])disp([' '])disp([' The fastest number of iteration : ' num2str(ite_num) ])% Plot the k-step iteration residualssemilogy(re,'k.-','LineWidth',0.8)title([ 'n = ',num2str(n),', \alpha = ',num2str(alpha)]);xlabel('Iterative numbers ');ylabel('Residual');legend({'HSS'},1)% Setting of rectangular coordinate system parameterset(get(gca,'title'), 'FontSize', 10,'FontName','Times New Roman');set(get(gca,'XLabel'),'FontSize',10,'FontName','Times New Roman');set(get(gca,'YLabel'),'FontSize',10,'FontName','Times New Roman');xlim([1 ite_num])shg然后里面涉及到的例子是：% A = M + 2rN + (100/((n+1)^2)) I% M = tridiag(-1,2,-1)% N = tridiag(0.5,0,-0.5)function A = SPD2(n,r)A = diag((-1+r)*ones(n-1,1),-1) + diag((-1-r)*ones(n-1,1),1) + (2+100/((n+1)^2))*eye(n);我得到的相对误差的图像是 （图2）我自己研究了好久多不知道我错在哪，但肯定是某一点错了，所以想麻烦各位帮我看看了，我曾经以为是F的问题，然后我试着换了F，但最后得到的R_k居然是一样的，这就更让我纳闷了，所以。。。辛苦大家了，谢谢！！！

Answer 1

1设计到误差问题，就是精度问题，控制算法本身就是为了提高精度而提出的，老师一定会问你的精度能达到多少，和其他算法比起来属于一个什么样的标准
2涉及到校正肯定会问你具体是什么校正，串联还是并联，滞后还是超前
3如果是控制问题的话在matlab仿真最关键的一点就是仿真图了，形象具体，老师会问你，控制性能如何（包括超调和各种响应时间快慢的问题）

Answer 2

help lyap看下lyap函数的用法，

lyap(A, B, C)求解的是AX+XB+C=0，而不是AX+XB=C。所以你HSS_2函数中循环的两行应该改成：

...
X1 = lyap(A1, B1, -C1);

...
X0 = lyap(A2, B2, -C2);