病态问题之（改进）Jacobi迭代法（附matlab源码）

病态方程的解法有三类：
第一类，病态观测方程的直接解法，如阶段奇异值分解法，
第二类，病态法方程的修正，如有偏估计：（广义）岭估计，遗传算法；
第三类，迭代算法，如改进的雅克比迭代算法【孔祥强，2013】，谱修正迭代算法【王新洲，2001】

本文主要是针对（改进）雅克比迭代算法，解算病态方程；
附：Jacobi迭代计算公式和Modified Jacobi计算公式

非病态方程解：
系数阵A的条件数:cond(A)=2.7244
方法1：传统Jacobi迭代法

%%Jacobi迭代法为整体代换法，n为系数矩阵A的阶数，b为右端常向量，x0为迭代初始值，

%jacobi[A,b,x0,n,maxN];
%A：系数矩阵
%b：常数项
%x0:迭代初值
%n:系数矩阵的阶数；
%maxn:最大迭代次数
%reference：【数值分析,book,P36】

%2018/10/27    executive！

A=[20.0 2.0 3.0;
    1.0 8.0 1.0;
    2.0 -3.0 15.0];
b=[24.0 12.0 30.0]';
% b=[24.0 12.0 30.0];
n=size(A,1);
%初始值
for  i=1:n
    x0(i)=100000;
end
maxN=100;%最大迭代值

%调用jacobi函数
% x=jacobi(A,b,x0,n,maxN);
%调用函数不管用，直接写吧
fid=fopen('Jacobi_Result.txt','w');
for i=1:n
    fprintf(fid,'%f ',x0(i));
    %fprintf(fid,'\n');
end
fprintf(fid,'\n');

for k=1:maxN
    %Jacobi迭代一次
    for i=1:n
        sum=0.0;
        for j=1:n
            if(j~=i)
                sum=sum+A(i,j)*x0(j);
                x(i)=-1.0/A(i,i)*sum+b(i)/A(i,i);
            end
        end

    end
     for i=1:n
        x0(i)=x(i);
        fprintf(fid,'%.10f ',x0(i));
%         fprintf(fid,'\n');
     end
     fprintf(fid,'\n');
end
    
    fclose(fid);

迭代结果（迭代19次收敛）：

方法2：改进的Jacobi迭代法

**%针对雅克比迭代法，无法有效解算病态方程，**
%本文则对jacobi方法进行改进
%reference:【孔祥强，2013】【MyNote,P63】

%A：系数矩阵
%b：常数项
%x0:迭代初值
%n:系数矩阵的阶数；
%maxn:最大迭代次数
%2018/10/27   

A=[20.0 2.0 3.0;
    1.0 8.0 1.0;
    2.0 -3.0 15.0];
b=[24.0 12.0 30.0]';
% b=[24.0 12.0 30.0];
n=size(A,1);
%初始值
for  i=1:n
    x0(i)=100000;
end
maxN=100;%最大迭代值

%求F矩阵对角元素f(i)
FF(:)=[];
for i=1:n
    if(A(i,i)>0)
        for j=1:n
            if(i~=j)
             temp=abs(A(i,:));
             temp(i)=0;
            end
        end
        f(i)=max(temp);
    else
        f(i)=-max(temp);
    end
    FF(i,i)=f(i);
        
end

%w如何取值呢？？慢慢试，一般在（0,2）之间
w=0.5;

fid=fopen('Modified_Jacobi_Result.txt','w');
for i=1:n
    fprintf(fid,'%f ',x0(i));
    %fprintf(fid,'\n');
end
fprintf(fid,'\n');

for k=1:maxN
    %Jacobi迭代一次
    for i=1:n
        sum1=0.0;
        sum2=0.0;
        for j=1:n
            if(j~=i)
                 sum1=sum1+A(i,j)*x0(j);
                 sum2=w*f(i)*x0(i);
                 x(i)=(-1.0*sum1+sum2+b(i))/(A(i,i)+w*f(i));
            end
        end

    end
     for i=1:n
        x0(i)=x(i);%把x(k)赋给x(0)
        fprintf(fid,'%.10f ',x0(i));
%         fprintf(fid,'\n');
     end
     fprintf(fid,'\n');
end
    
fclose(fid);
    

迭代结果（迭代22次收敛）：

发现上面两个算法无法解算病态方程；问题还没解决？？？
病态方程解：

方法1：传统Jacobi迭代