解线性方程组的迭代法

1.实验目的

通过上机实验，理解利用计算机迭代求解线性方程组的整个过程，加深对所学计算方法的理论及算法特点的理解。

2.实验题目

2.1

利用算法2.1（Jacobi迭代法），编制MATLAB程序，求线性方程组
(1)

\[\begin{pmatrix} 14&4&4&4\\ 4&14&4&4\\ 4&4&14&4\\ 4&4&4&14\\ \end{pmatrix} \begin{pmatrix} x_1\\ x_2\\ x_3\\ x_4\\ \end{pmatrix} = \begin{pmatrix} -4\\ 16\\ 36\\ 56\\ \end{pmatrix} \]

(2)

\[\begin{pmatrix} 10.9&1.2&2.1&0.9\\ 1.2&11.2&1.5&2.5\\ 2.1&1.5&9.8&1.3\\ 0.9&2.5&1.3&12.3\\ \end{pmatrix} \begin{pmatrix} x_1\\ x_2\\ x_3\\ x_4\\ \end{pmatrix} = \begin{pmatrix} -7.0\\ 5.3\\ 10.3\\ 24.6\\ \end{pmatrix} \]

的近似解，取初值 \(x = (0,0,0,0)^T\).

2.2

利用算法2.2（Gauss-Seidel迭代法），编制MATLAB程序，求线性方程组
(1)

\[\begin{pmatrix} 6&-2&-1&-1\\ -2&12&-1&-1\\ -1&-1&6&2\\ -1&-1&-1&12\\ \end{pmatrix} \begin{pmatrix} x_1\\ x_2\\ x_3\\ x_4\\ \end{pmatrix} = \begin{pmatrix} -16\\ 6\\ 8\\ 54\\ \end{pmatrix} \]

(2)

\[\begin{pmatrix} 0.78&-0.02&-0.12&-0.14\\ -0.02&0.86&-0.04&-0.06\\ -0.12&-0.04&-0.72&-0.08\\ -0.14&-0.06&-0.08&0.74\\ \end{pmatrix} \begin{pmatrix} x_1\\ x_2\\ x_3\\ x_4\\ \end{pmatrix} = \begin{pmatrix} 0.76\\ 0.08\\ 1.12\\ 0.68\\ \end{pmatrix} \]

的近似解，取初值 \(x = (0,0,0,0)^T\).

2.3

利用算法2.3（SOR法），编制MATLAB程序，求线性方程组
(1)

\[\begin{pmatrix} -4&1&1&1\\ 1&-4&1&1\\ 1&1&-4&1\\ 1&1&1&-4\\ \end{pmatrix} \begin{pmatrix} x_1\\ x_2\\ x_3\\ x_4\\ \end{pmatrix} = \begin{pmatrix} 1\\ 1\\ 1\\ 1\\ \end{pmatrix} \]

(2)

\[\begin{pmatrix} 1&0&-0.25&-0.25\\ 0&1&-0.25&-0.25\\ -0.25&-0.25&1&0\\ -0.25&-0.25&0&1\\ \end{pmatrix} \begin{pmatrix} x_1\\ x_2\\ x_3\\ x_4\\ \end{pmatrix} = \begin{pmatrix} 0.5\\ 0.5\\ 0.5\\ 0.5\\ \end{pmatrix} \]

的近似解，取初值 \(x = (0,0,0,0)^T\)，松弛因子分别为 \(\omega = 1.3 , \omega = 1.1\).

3.程序文本

Jacobi迭代

function x = jac(A,b,x0,ep,N)
n = length(b);
if nargin < 5,N = 500;end
if nargin < 4,ep = 1e-6;end
if nargin < 3,x0 = zeros(n,1);end
x = zeros(n,1);k = 0;
while k < N
    for i = 1:n
        x(i) = (b(i) - A(i,[1:i-1,i+1:n]) * x0([1:i-1,i+1:n])) / A(i,i);
    end
    if norm(x-x0,inf) < ep,break;end
    x0 = x;k = k+1;
end
if k == N ,warning('111');end
disp(['k=',num2str(k)])

Gauss-Seidel迭代

function x = gau(a,b,x0,ep,N)
n = length(b);
if nargin < 5,N = 500;end
if nargin < 4,ep = 1e-6;end
if nargin < 3,x0 = zeros(n,1);end
x = zeros(n,1);k = 0;
while k < N
    for i = 1:n
        if i == 1
            x(1) = (b(1) - a(1,2:n) * x0(2:n)) / a(1,1);
        else 
            if i == n
                x(n) = (b(n) - a(n,1:n-1) * x(1:n-1)) / a(n,n);
            else
                x(i) = (b(i) - a(i,1:i-1) * x(1:i-1) - a(i,i+1:n) * x0(i+1:n)) / a(i,i);
            end
        end
    end
    if norm(x-x0,inf) < ep,break;end
    x0 = x;k = k+1;
end
if k == N,warning('111');end
disp(['k=',num2str(k)])

SOR法

function x = sor(a,b,ome,x0,ep,N)
n = length(b);
if nargin < 6,N = 500;end
if nargin < 5,ep = 1e-6;end
if nargin < 4,x0 = zeros(n,1);end
if nargin < 3,ome = 1.5;end
x = zeros(n,1);k = 0;
while k < N
    for i = 1:n
        if i == 1
            x1(1) = (b(1) - a(1,2:n) * x0(2:n)) / a(1,1);
        else
            if i == n
                x1(n) = (b(n) - a(n,1:n-1) * x(1:n-1)) / a(n,n);
            else
                x1(i) = (b(i) - a(i,1:i-1) * x(1:i-1) - a(i,i+1:n) * x0(i+1:n)) / a(i,i);
            end
        end
        x(i) = (1 - ome) * x0(i) + ome * x1(i);
    end
    if norm(x-x0,inf) < ep,break;end
    x0 = x;k = k+1;
end
if k == N,warning('111');end
disp(['k=',num2str(k)])

4.运行结果与分析

>> a = [14 4 4 4;4 14 4 4;4 4 14 4;4 4 4 14];
>> b = [-4 16 36 56]';
>> jac(a,b)
k=94

ans =

   -2.0000
    0.0000
    2.0000
    4.0000

>> a = [10.9 1.2 2.1 0.9;1.2 11.2 1.5 2.5;2.1 1.5 9.8 1.3;0.9 2.5 1.3 12.3];
>> b = [-7.0 5.3 10.3 24.6]';
>> jac(a,b)
k=17

ans =

   -0.9986
    0.0071
    1.0032
    1.9656

>> a = [6 -2 -1 -1;-2 12 -1 -1;-1 -1 6 -2;-1 -1 -1 12];
>> b = [-16 6 8 54]';
>> gau(a,b)
k=11

ans =

   -1.0782
    0.9553
    2.8897
    4.7306

>> a = [0.78 -0.02 -0.12 -0.14;-0.02 0.86 -0.04 -0.06;-0.12 -0.04 0.72 -0.08;-0.14 -0.06 -0.08 0.74];
>> b = [0.76 0.08 1.12 0.68]';
>> gau(a,b)
k=8

ans =

    1.5503
    0.3232
    1.9934
    1.4539

>> a = [-4 1 1 1;1 -4 1 1;1 1 -4 1;1 1 1 -4];
>> b = [1 1 1 1]';
>> sor(a,b,1.3)
k=13

ans =

   -1.0000
   -1.0000
   -1.0000
   -1.0000

>> a = [1 0 -0.25 -0.25;0 1 -0.25 -0.25;-0.25 -0.25 1 0;-0.25 -0.25 0 1];
>> b = [0.5 0.5 0.5 0.5]';
>> sor(a,b,1.1)
k=6

ans =

    1.0000
    1.0000
    1.0000
    1.0000