技術標籤：數值分析-matlab程式matlab演算法

題目如圖所示：

要求程式中不存係數矩陣A,分別對不同的階數取w=1.1, 1.2, …,1.9進行迭代；
首先編寫根據矩陣的階數生成係數矩陣的函式，然後再SOR迭代程式中呼叫此函式；觀察係數矩陣的規律，編寫函式SOR_A如下：
程式碼塊1：根據階數生成係數矩陣：

function [A,b]=SOR_A(n)
A1=(-4)*eye(n);
for i=1:n
    for j=1:n
        if i==j+1
            L(i,j)=1;
        else L(i,j)=0;
        end%
 
生成下三角矩陣L
        if j==i+1
            U(i,j)=1;
        else U(i,j)=0;
        end%生成上三角矩陣U
    end
end
A=A1+L+U;%生成係數矩陣A
b1(1)=-3;
for i=2:n-1
    b1(i)=-2;
end
b1(n)=-3;
b=b1;%生成結果向量b

SOR迭代法和Jacobi迭代和Gauss-Seidel迭代法構造思路一致，只是迭代矩陣的構造不同，外加了一個鬆弛因子w，當w=1時，此迭代法和高斯塞德爾迭代法一致，SOR迭代法程式碼塊如下：
程式碼塊2：SOR迭代法：

%函式SOR_A根據輸入的方程階數自動生成矩陣A
clc;
 
clear;
n=input('請輸入係數矩陣A的階數:n=');
%n=15;
[A,b1]=SOR_A(n);%上機實習題5中的矩陣A；b為行向量，計算時需要轉置
%A=input('請輸入係數矩陣A:A=');%把上行註釋掉後自由度更大，自己定義矩陣A；
%b=input('請輸入結果向量b:b=');%把上行註釋掉後自由度更大，自己定義矩陣b；
x1=input('請輸入初始向量x1:x1=');%x1為行向量，計算時需要轉置
%x1=[3 3 3 3 3 3 3 3 3 3 3 3 3 3 3];
%eps=input('請輸入停止精度要求:eps=');
eps=1e-6;
w=input
 
('請輸入鬆弛因子:w=');
for i=1:n
    for j=1:n
        if i==j
            D(i,j)=A(i,j);
        else D(i,j)=0;
        end%生成矩陣D
        if i>j
            L(i,j)=-A(i,j);
        else L(i,j)=0;
        end%生成矩陣L
        if i<j
            U(i,j)=-A(i,j);
        else U(i,j)=0;
        end%生成矩陣U
    end
end
D1=inv(D-w*L);%矩陣(D-w*L)的逆矩陣記作D1
D2=((1-w)*D+w*U);
Lw=D1*D2;
f=w*D1*b1';%f為行向量
k=2;
xk=Lw*x1'+f;%xk為列向量，儲存時需要轉置
x=[x1;xk'];
bk=A*xk;
b=[b1;bk'];
while max(abs(x(k,:)-x(k-1,:)))>eps%
    if k>500||max(abs(eig(Lw)))>1%如果迭代矩陣B的譜半徑大於1，則迭代不收斂
        error('迭代失敗!可能是迭代精度過高或迭代矩陣譜半徑大於1的緣故,迭代矩陣譜半徑為:%0.8f\n',max(abs(eig(Lw))))
    else
        k=k+1;%第一次為k=4
        xk=Lw*x(k-1,:)'+f;%x(k-1,:)為X矩陣的第k-1行的行元素
        x=[x;xk'];%第一次後X矩陣為4行3列；
        bk=A*xk;
        b=[b;bk'];
    end
end
fprintf('經過%d次迭代，方程組的根的近似解為:\n',k-1)
disp(vpa(x(k,:),8))
fprintf('SOR迭代法迭代矩陣的譜半徑為：%0.8f\n',max(abs(eig(Lw))))

執行並輸入相應變數，結果如下圖：

分別輸入不同的階數n和不同的鬆弛因子w，得到不同的結果，並對比迭代次數，觀察迭代因子w與迭代矩陣譜半徑和迭代速度的關係