线性方程组的迭代解法

提示：代码仅供参考

文章目录

• 前言
• 一、Jacobi迭代
• 二、Gauss-Seidel迭代
• 三、松弛法

---

前言

  算法原理可参考金一庆的《数值计算》。

---

提示：以下是本篇文章正文内容，下面案例可供参考

一、Jacobi迭代

tic
t1=clock;
%A=input('请输入系数矩阵A：');
A=[1.5 0.25 -0.3;0.15 1.28 -0.4;-0.54 0.46 1.17]
b=[4 5 6]
x=size(A)
%b=input('请输入常向量b：');
b=b'
if all(diag(A))==1 || x(1)==x(2)
    disp('该矩阵可以使用Jacobi迭代')
else
    disp('该矩阵不能使用Jacobi迭代')
end
%若矩阵可以使用Jacobi迭代法则进行以下迭代计算
d=diag(A).*eye(size(A))
d1=inv(d)
%input('请输入初始向量')
x0=[2.7 4.9 4.4]'
for i=1:20
    x=(eye(size(A))-d1*A)*x0+d1*b
    if normest(x-x0)<eps
        disp('Jacobi迭代收敛，且收敛值为:')
        break
        x0=x
    end
end
x,normest(x-x0)
toc
disp(['toc计算最后一次循环运行时间',num2str(toc)])
disp(['etime程序总运行时间：',num2str(etime(clock,t1))])

二、Gauss-Seidel迭代

tic
t1=clock;
%A=input('请输入系数矩阵A：');
A=[1.5 0.25 -0.3;0.15 1.28 -0.4;-0.54 0.46 1.17]
b=[4 5 6]
x=size(A)
L=tril(A,-1)
U=triu(A,1)
d=diag(A).*eye(size(A))
%b=input('请输入常向量b：');
b=b'
if det(d+L)~=0 || x(1)==x(2)
    disp('该矩阵可以使用Gauss-Seidel迭代')
else
    disp('该矩阵不能使用Gauss-Seidel迭代')
end
%若矩阵可以使用Gauss-Seidel迭代,则进行以下迭代计算
K=inv(d+L)
%input('请输入初始向量')
x0=[2.7 4.9 4.4]'
for i=1:20
    x=-K*U*x0+K*b
    if normest(x-x0)<eps
        disp('Gauss-Seidel迭代收敛，且收敛值为:')
        break
    end
    x0=x
end
x,normest(x-x0)
toc
disp(['toc计算最后一次循环运行时间',num2str(toc)])
disp(['etime程序总运行时间：',num2str(etime(clock,t1))])

三、松弛法 

tic
t1=clock;
%A=input('请输入系数矩阵A：');
A=[4 -1 0;-1 4 -1;0 -1 4]
b=[1 4 -3]
x=size(A)
%b=input('请输入常向量b：');
b=b'
L=tril(A,-1)
U=triu(A,1)
d=diag(A).*eye(size(A))
w=1.1%1.03 1  1.1
H=inv(d+w*L)*((1-w)*d-w*U)
ro=max(eig(H'*H))
if ro<1
    disp('可使用松弛法，且松弛法收敛')
else
    disp('不建议使用松弛法')
end
%判断松弛法是否适用
x0=[0 0 0]'
for i=1:20
    x=H*x0+inv(d+w*L)*w*b
    
    if normest(x-x0)<eps
        disp('松弛法迭代收敛，且收敛值为:')
        break
    end
    x0=x
end
x
toc
disp(['toc计算最后一次循环运行时间',num2str(toc)])
disp(['etime程序总运行时间：',num2str(etime(clock,t1))])

---

总结

数值计算实验大概是我本学期学起来最轻松的一门课了。