疊代法——Matlab中實作

疊代法

這裡一共提供了四種疊代法：

+ 雅可比疊代法

+ 高斯賽德疊代法

+ 超松弛疊代法（SOR）

+ 共轭疊代法

随機生成方程組

此處随機生成特征值服從獨立同分布的[0,1]間的均勻分布的A矩陣，跟服從獨立同分布的正态分布的b向量

算法設計

• A矩陣的實作，首先需要用rand得到一組特征值，将該組特征值通過diag函數生成對角陣Q，之後通過orth函數生成對稱矩陣U，再通過UQU’即可得到對稱正定矩陣
• b向量的實作同第一題，用randn函數即可

代碼實作：

%n代表的是次元
b = randn([n(k) ]);
Q = diag(rand([ n(k)]));
U = orth(rand([n(k) n(k)]));
           

Jacobi 疊代法函數

算法設計

• 傳入參數為矩陣A，向量b，以及次元n
• 傳出參數為結果與時間（or相對誤差數組），具體輸出視操作而定
• 通過diag函數，tril函數，triu函數得到A矩陣的對角陣D，上三角矩陣U跟下三角矩陣L
• 通過以下公式得到疊代所需的B跟f

B = D(L+U);
  f = D\b;`
           

• 通過以下公式進行疊代：

• 疊代終止的判斷，求出前後兩次疊代的向量差的範數，直至小于1^-6
• 如果疊代次數超過2000，則視為發散，彈出錯誤

代碼實作：

function [x,y] = Jacobi(n, A ,b)

    %y = zeros(1000,1);
    eps = ;
    D = diag(diag(A));
    L = -tril(A, -);
    U = -triu(A, );

    B = D\(L+U);
    f = D\b;

    count = ;

    x0 = zeros(n,);
    x = B*x0 + f;

    tic;
    while norm(x-x0)> eps
        x0 = x;
        %y(count) = norm(x-A\b);
        x = B*x0 + f;
        count = count + ;
        if count > 
            disp('error:該矩陣不收斂');
            return;
        end
    end
    toc;
    y = toc;
    disp(count);
end
           

Gauss-Seidel 疊代法函數

算法設計

• 傳入參數為矩陣A，向量b，以及次元n
• 傳出參數為結果與時間（or相對誤差數組），具體輸出視操作而定
• 通過diag函數，tril函數，triu函數得到A矩陣的對角陣D，上三角矩陣U跟下三角矩陣L
• 通過以下公式得到疊代所需的B跟f

B = (D - L)\U;
    f = (D - L)\b;
           

• 通過以下公式進行疊代：

• 疊代終止的判斷，求出前後兩次疊代的向量差的範數，直至小于1^-6
• 如果疊代次數超過2000，則視為發散，彈出錯誤

代碼實作：

function [x, y]= Gauss_Seidel(n, A ,b)

    %y = zeros(1000,1);
    eps = ;
    D = diag(diag(A));
    L = -tril(A, -);
    U = -triu(A, );

    B = (D - L)\U;
    f = (D - L)\b;

    count = ;
    x0 = zeros(n,);
    x = B*x0 + f;

    tic;
    while norm(x-x0) > eps

        x0 = x;
       % y(count) = norm(x-A\b);
        x = B*x0 + f;
        count = count + ;

        if count > 
            disp('error:該矩陣不收斂');
            return;
        end
    end
    toc;
    y = toc;
    disp(count);
end
           

逐次超松弛疊代法函數

算法設計

• 傳入參數為矩陣A，向量b，以及次元n
• 傳出參數為結果與時間（or相對誤差數組），具體輸出視操作而定
• 通過diag函數，tril函數，triu函數得到A矩陣的對角陣D，上三角矩陣U跟下三角矩陣L
• 通過以下公式得到疊代所需的B跟f

B = (D - w*L) \ ((-w) * D + w*U);
    f = w * inv(D - w*L) * b;
           

• 通過以下公式進行疊代：

• 疊代終止的判斷，求出前後兩次疊代的向量差的範數，直至小于1^-6
• 如果疊代次數超過2000，則視為發散，彈出錯誤

代碼實作：

function [x, y] = SOR(n, A ,b, w)

    y = zeros(,);
    eps = ;
    D = diag(diag(A));
    L = -tril(A, -);
    U = -triu(A, );

    B = (D - w*L) \ ((-w) * D + w*U);
    f = w * inv(D - w*L) * b;

    count = ;

    x0 = zeros(n,);
    x = B*x0 + f;

    tic;
    while norm(x-x0) > eps

        x0 = x;
        y(count) = norm(x-A\b);
        x = B*x0 + f;
        count = count + ;

        if count > 
            disp('error:該矩陣不收斂');
            return;
        end
    end
    toc;
    y = toc;
    disp(count);
end
           

共轭梯度法函數

算法設計

• 傳入參數為矩陣A，向量b，以及次元n
• 傳出參數為結果與時間（or相對誤差數組），具體輸出視操作而定
• 通過以下公式得到疊代所需的r0跟p0

r0 = b - A * x0;
    p0 = r0;
           

• 通過以下公式進行疊代：

a = r0'r0/(p0'A*p0);
      x = x0 + a*p0;
      r = r0 - a*A*p0;
      B = r'*r/(r0'*r0);
      p = r + B*p0;
           

• 疊代終止的判斷，求出前後兩次疊代的向量差的範數，直至小于1^-6
• 如果疊代次數超過2000，則視為發散，彈出錯誤

代碼實作：

function [x, y] = CG(n, A , b)

    y = zeros(,);
    eps = ;
    x0 = zeros(n,);
    r0 = b - A * x0;
    p0 = r0;

    count = ;

    tic;
    while norm(p0) > eps
        a = r0'*r0/(p0'*A*p0);
        x = x0 + a*p0;
        r = r0 - a*A*p0;
        B = r'*r/(r0'*r0);
        p = r + B*p0;

        p0 = p;
        r0 = r;
        x0 = x;
       % y(count) = norm(x-A\b);
        count = count + ;

        if count > 
            disp('error:該矩陣不收斂');
            return;
        end
    end
    toc;
    y = toc;

    disp(count);
end
           

相對誤差計算

算法設計

• 通過A\b得到近似解
• 通過求出近似解與得到的解的差向量的範數作為相對誤差

代碼實作

具體的代碼實作，已經穿插在四種疊代法中