迭代法

這裡一共提供了四種迭代法：
+ 雅可比迭代法
+ 高斯賽德迭代法
+ 超鬆弛迭代法（SOR）
+ 共軛迭代法

隨機生成方程組

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

演算法設計

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

程式碼實現：

%n代表的是維度
b = randn([n(k) 1]);
Q = diag(rand([1 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;`

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

  x = B*x0 + f;

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

程式碼實現：

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

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

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

    count = 1;

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

    tic;
    while norm(x-x0)> eps
        x0 = x;
        %y(count) = norm(x-A\b);
 

        x = B*x0 + f;
        count = count + 1;
        if count > 2000
            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;

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

     x = B*x0 + f;

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

程式碼實現：

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

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

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

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

    tic;
    while norm(x-x0) > eps

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

        if count > 2000
            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) \ ((1-w) * D + w*U);
    f = w * inv(D - w*L) * b;

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

    x = B*x0 + f;

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

程式碼實現：

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

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

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

    count = 1;

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

    tic;
    while norm(x-x0) > eps

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

        if count > 2000
            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(1000,1);
    eps = 1.0e-6;
    x0 = zeros(n,1);
    r0 = b - A * x0;
    p0 = r0;

    count = 1;

    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 + 1;

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

    disp(count);
end

相對誤差計算

演算法設計

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

程式碼實現

具體的程式碼實現，已經穿插在四種迭代法中