实验十二.线性方程迭代法求解实现

实验内容	1. 编写MATLAB函数实现高斯消去法,并求解以下方程组 2. 分别编写MATLAB函数实现用迭代法(雅可比、高斯-塞德尔法、SOR方法)求解线性方程组,并通过求解以下方程组比较三种方法。

实 验 步 骤、过 程

1.

A=[1,1,1;0,4,-1;2,-2,1];
b=[6,5,1];% 系数矩阵
y=inv(A).*b;% n维向量
n=length(b);% matlab的计算结果
x=zeros(n,1);
%-------------消去-----------
for k=1:n-1
%     if A(k,k)==0;
%         error('Error');
%     end
    for i=k+1:n
%         A(i,k)=A(i,k)/A(k,k); 
         Aik=A(i,k)/A(k,k)
         for j=k:n
              A(i,j)=A(i,j)-Aik*A(k,j);
         end
         A
         b(i)=b(i)-Aik*b(k)
    end
end
%-------------回代-----------
x(n)=b(n)/A(n,n)
for k=n-1:-1:1
    S=b(k);
    for j=k+1:n
        S=S-A(k,j)*x(j);
    end
    x(k)=S/A(k,k)
end
x         %程序的计算结果

2.

(1)雅可比

function[x,number] = su(A,b,x0,er)
%Jacobi 迭代法矩阵形式
%x迭代向量列，x0迭代初值，er误差，number迭代次数
D = diag(diag(A));
D = inv(D);
U = triu(A,1);
L = tril(A,-1);
B = -D.*(L+U);
f = D.*b;
number = 0;
x = B.*x0 + f;
number = number + 1;
while norm(x-x0) > er
    x0 = x
    x = B.*x0 + f;
    number = number + 1;
end
命令行窗口输入
A=[8,-3,2;4,11,-1;6,3,12];
b=[20,33,36]';
er = 1e-5;
[x,n] = su(A,b,[0,0,0]',er)

(2)高斯-塞德尔法

function[x,number] = su(A,b,x0,er)
%高斯——赛德尔迭代法矩阵形式
%x迭代向量列，x0迭代初值，er误差，number迭代次数
D = diag(diag(A));
U = triu(A,1);
L = tril(A,-1);
D = inv(D+L);
B = -D.*U;
f = D.*b;
number = 0;
x = B.*x0 + f;
number = number + 1;
while norm(x-x0)>er
    x0 = x;
    x = B.*x0 + f;
    number = number + 1;
end

命令行窗口输入

A=[8,-3,2;4,11,-1;6,3,12];
b=[20,33,36];
er = 1e-5;
[x,n] = su(A,b,[0,0,0]',er)

(3)SOR方法

function[x,number] = su(A,b,x0,er,w)
%SOR方法
%x迭代向量列，x0迭代初值，er误差，number迭代次数，w松弛因子
[m,n] = size(A);
x = ones(m,1);
number = 0;
x = su2(A,b,x0,er,w);
number = number + 1;
while norm(x-x0) > er
    x0 = x;
    x = su2(A,b,x0,er,w);
    number = number + 1;
end
function x = su2(A,b,x0,er,w)
[m,n] = size(A);
x = ones(m,1);
for k = 1:n
    if k == 1
        ss = 0;
        for j = k:n
            ss = ss + A(k,j).*x0(j);
        end
        x(k) = x0(k) + w.*(b(k)-ss)./A(k,k);
    else
        ss = 0;
        for j = 1:(k-1)
            ss = ss + A(k,j).*x(j);
        end
        for j = k:n
            ss = ss + A(k,j).*x0(j);
        end
        x(k) = x0(k) + w.*(b(k) - ss)./A(k,k);
    end
end

命令行窗口输入

A=[8,-3,2;4,11,-1;6,3,12];
b=[20,33,36];
er = 1e-5;
[x,n] = su(A,b,[0,0,0]',er,1.12)

答案供参考,如有错误请自行改正，本人只是随意分享之前的作业，代码均来源于网络+个人稍微修改，认为有用可参考哈～